Numerical accuracy of expression involving norm squared

Question

I am computing the following quantity: $$ \text{lhs} := ||a+b||^2 = ||a||^2 + 2a^\top b + ||b||^2 =: \text{rhs} $$ for $a=c-d$, where $a,b,c,d$ are $n$-vectors. Is there a rule of thumb for when I should have my program compute $\text{lhs}$ vs $\text{rhs}$ in floating point arithmetic? It seems that if $a,b,c,d$ are all of the same scale, then either way is fine (based on this question). On the other hand, what if some of $a,b,c,d$ are of very different scales? I will add some thoughts below.

---

1. if $c+b$ is of the same scale as $d$, then $\text{lhs}$ will be better since the final step in $\text{rhs}$ involves adding something where $2a^\top b$ could be of opposite sign as both $||a||^2$ and $||b^2||$ but of the same scale, resulting in less precision.

Answer 1

The LHS is stable, the RHS is not.

The RHS is liable to numerical cancellation, because you subtract two computed quantities that have a smaller (relative) difference. The difference is visible even in dimension 1:

julia> a::Float64 = 10^12 + 1
1000000000001

julia> b::Float64 = -10^12
-1000000000000

julia> (a+b)^2, a^2 + 2*a*b + b^2
(1.0, -1.34217728e8)

You are computing values of magnitude about $10^{24}$, so you incur in an error of magnitude about eps(1e24).

Isn't there a subtraction in the LHS, too? Yes, but it applies to the exact inputs $a$ and $b$ (or their entries if they are vectors). If they are given exactly, then the algorithm will return the exact result perturbed by a relative error of the size of the machine precision eps, since it contains no subtractions; you can write down a detailed error analysis to confirm it. The returned result will have an expression of the form $\|a+b\|^2(1+\delta)$, where $|\delta|\leq n\mathsf{u} + O(\mathsf{u}^2)$, with $n$ the length of the vectors and $\mathsf{u}$ the machine precision eps(1).

If you pass to the LHS not the exact inputs $a$ and $b$, but perturbed version $\tilde{a}$ and $\tilde{b}$, then the computed result can still change by a large quantity. But this is not the fault of the algorithm, since even the best algorithm in the world to compute $\|a-b\|^2$ will have the same issue. In other words, evaluating the formula in the LHS gives a stable algorithm for a problem ill-conditioned for certain $a$ and $b$.

In general, a good principle to get stable algorithms is the "no inaccurate cancellation" principle: try to avoid all subtractions (as in "summing numbers with different signs"), unless they involve input values. In your case, the subtractions in the LHS involve input values, those in the RHS do not.