More stable algorithm to calculate `sqrt(a^2 + b^2) - abs(a)` in MatLab

Question

Suppose we want to calculate $\sqrt{a^2+b^2}-|a|$ in MatLab. Using sqrt(a^2 + b^2) - abs(a) will have some problems:

• If a or b are too large, it may cause overflow
• if a and b are too small, it may cause underflow
• if a>>b, a^2 + b^2 may be equal to a^2 in machine

What's a more stable algorithm to calculate this without leading to overflow, underflow or omitting b^2?

First thing that came to my mind was to use if clauses that if one error is going to happen, then stop. But I think there can be a way to calculate the answer anyway, because we don't need a^2 or b^2, we need sqrt(a^2 + b^2) - abs(a). But I couldn't find a way to do the calculation when for example b is too large. Second thing that came to my mind was to change the formula. I only could reach sqrt((a + b)^2 - 2*a*b) - abs(a) and couldn't expand the formula anymore. This would help in 3rd problem case, but won't help overflow or underflow I think. Any help is so much appreciated!

Answer 1

A manipulation that may help is the following. Assume for simplicity $a>0$. We have the identity $$b^2 = (a^2+b^2)-a^2 = (\sqrt{a^2+b^2}-a)(\sqrt{a^2+b^2}+a),$$ hence $$ \sqrt{a^2+b^2}-a = \frac{b^2}{\sqrt{a^2+b^2}+a}, $$ and the cancellation is gone.

Note that $\sqrt{a^2+b^2}$ can be computed without the overflow associated to squaring with hypot(a,b). If one rearranges the computation as (b / (hypot(a,b)+a)) * b I think that most of your corner cases are solved. The only danger is that hypot(a,b)+a may overflow, but that will happen only if $a,b$ are of the order of $10^{300}$.

EDIT: you may be interested in checking out methods for computing accurate solutions to scalar quadratic equations. This is a closely related problem, as one encounters similar issues when computing the numerator of the classical formula $b \pm \sqrt{b^2-4ac}$; and indeed your problem can be reformulated as finding the positive solution of the quadratic equation $x^2 + 2ax - b^2=0$. For instance, there is this note by one of the founding fathers of floating-point computing, and this SO answer with useful algorithms.

Answer 2

You can break down the domain of your function into three distinct cases:

• $|a|\gg |b|$: In this case, $\sqrt{a^2+b^2} \approx |a|$ and a naive application of the formula will likely result in poor accuracy. But you can transform $$ \sqrt{a^2+b^2} - |a| = \sqrt{a^2(1+b^2/a^2)} - |a| = |a| \left(\sqrt{1+b^2/a^2}-1\right) $$ and at this point you use a Taylor expansion. To first order, because $\sqrt{1+x}\approx 1+\frac{1}{2}x$ for small $x$, you get $$ \sqrt{a^2+b^2} - |a| \approx|a| \left(\frac 12 b^2/a^2\right) = \frac{b^2}{2|a|}. $$ If you want better accuracy, you can always higher order Taylor expansions.

• $|a|\ll|b|$: This also results in poor accuracy because the effect of $a$ is lost. But again we can do the following: $$ \sqrt{a^2+b^2} - |a| = \sqrt{b^2(1+a^2/b^2)} - |a| = |b|\left(\sqrt{1+a^2/b^2} - \frac{|a|}{|b|}\right). $$ By the same Taylor expansion, you can approximate $$ \sqrt{a^2+b^2} - |a| \approx |b|\left(1+\frac{a^2}{2b^2} - \frac{|a|}{|b|}\right). $$ This can again be stably evaluated in terms of the small parameter $a/b$. As before, you can use a higher order Taylor approximation.

• $|a|\approx |b|$: In this case, the two quantities are approximately the same and no truncation error should happen when evaluating the original formula.

Answer 3

This is basically a variant of Federico Poloni's excellent answer, based on the principle that one might want to try proven building blocks first. I do not know if Matlab provides this functionality, but various computational platforms offer a function sqrt1pm1(), which computes $\sqrt{x+1}-1$ accurately while avoiding subtractive cancellation. This is an excellent building block for situations like this one. Where a function sqrt1pm1() is not provided, we could easily build it ourselves:

/* compute sqrt (x+1)-1 */
double sqrt1pm1 (double a) { 
    return a / (1 + sqrt (a + 1)); 
}

With this building block, the computation from the question turns into:

/* compute sqrt(a**2 + b**2) - |a| */
double func (double a, double b) 
{
    double ratio = b / a;
    return fabs (a) * sqrt1pm1 (ratio * ratio);
}

While this addresses accuracy issues arising from subtractive cancellation, it does not address all potential issues of potential overflow and underflow in intermediate computation. In particular it lacks the range-mitigating effect provided by the use of hypot(). Obviously this approach also requires that $a \ne 0$, but $a=0$ is easily handled separately.