CODEWITHSUNDEEP
divisionnumerical-methods
vibe
vibe

Reputation: 385

Division in double precision

I have two double variables:

a > 0
b >= 0

which could be tiny numbers. 'a' represents singular values of a matrix and 'b' represents the Tikhonov regularization constant. As part of the Tikhonov least squares solution, it is necessary to compute the quantity:

c = a*a / (a*a + b)

However if a is really small (ie small singular values of the matrix), a*a may not be representable in double precision. How can I compute this quotient c in a numerically stable way for the given ranges of a,b?

Upvotes: 1

Views: 131

Answers (1)

TypeIA
TypeIA

Reputation: 17248

The best I can come up with is:

c = 1 / (1 + b / a / a)

To derive this equivalency, note that 1/c is (a^2 + b)/c and then decompose the fraction. This form might be more numerically stable since it doesn't require a^2 to be calculated at any point. It'll still lose precision if both b and a are very small. If that case must be handled too, you might look at a Taylor series expansion (may or may not work for this case).

Upvotes: 1

Related Questions