9
$\begingroup$

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!

$\endgroup$
0

3 Answers 3

18
$\begingroup$

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.

$\endgroup$
0
24
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Thanks a lot! That was awesome. But I was looking for methods that do not use functions like hypot(a, b) as I have some limits at the moment. But It'll be really useful to me in the future. Thanks again. $\endgroup$ Oct 31, 2023 at 14:00
  • 2
    $\begingroup$ @IlkayBurak This answer is brilliant and does not really force you to use hypot in any way. That's just his suggestion to avoid overflows. If you want, you can actually copy the implementation of hypot from an (appropriately-licensed) open source language into your own code since the implementation is quite simple. Alternatively, you can also use a cases based solution and approximations as in @WolfgangBangerth's good answer. $\endgroup$ Oct 31, 2023 at 14:42
  • 3
    $\begingroup$ Unless I'm missing something obvious, there's no need to assume positive $a$: $\sqrt{a^2+b^2}-|a| = \frac{(\sqrt{a^2+b^2}-|a|)(\sqrt{a^2+b^2}+|a| )}{\sqrt{a^2+b^2}+|a|}=\frac{b^2}{\sqrt{a^2+b^2}+|a|}$ $\endgroup$
    – Luca Citi
    Nov 1, 2023 at 1:13
  • 1
    $\begingroup$ @LucaCiti I agree; that assumption is just there to save a few absolute values and have simpler-looking formulas. $\endgroup$ Nov 1, 2023 at 7:22
  • $\begingroup$ Sure, I understand. However my derivation is equally simple, just requiring multiplying and dividing by the same quantity and then noticing that you obtain a sum times a difference of two terms, which is a difference of squares. Then it's easy from there. $\endgroup$
    – Luca Citi
    Nov 1, 2023 at 23:36
5
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.