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:

$$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$$

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.