How Does Tikhonov Regularization Handle Underflows in IEEE Double Precision?

[vibe3]

The least-squares solution of A x = b using Tikhonov regularization with a matrix \mu^2 I has the solution:

<br /> x = \sum_i \left( \frac{\sigma_i^2}{\sigma_i^2 + \mu^2} \right) \left( \frac{u_i^T b}{\sigma_i} \right) v_i<br />

where A = U S V^T is the SVD of A and u_i,v_i are the columns of U,V.

For ill-conditioned matrices, the singular values \sigma_i could be tiny leading to problems in computing the quantity \left( \frac{\sigma_i^2}{\sigma_i^2 + \mu^2} \right) since \sigma_i^2 could underflow.

Does anyone know how to compute this solution safely and efficiently in IEEE double precision?

 

[AlephZero]

If $\sigma_i^2 \ll \mu^2$, then the corresponding term of the sum is approximately
$$\frac{\sigma_i}{\mu^2} (u_i^T b) v_i$$

But in practice you are unlikely to have a problem, because the ratio of $\sigma_\min / \sigma_\max$ will be limited by the numerical precision of the matrix $A$, and unless $A$ has some pathological properties $\sigma_\min / \sigma_\max$ is unlikely to be less than about $10^{-16}$.

If all the SVs and $\mu$ are very small (e.g. less than $10^{-100}$) so underflows are likely to affect everything, the best fix would be to rescale the original problem to make them of order 1.