How much regularization to add to make SVD stable?

Question

I've been using Intel MKL's SVD (dgesvd through SciPy) and noticed that results are are significantly different when I change precision between float32 and float64 when my matrix is badly conditioned/not full rank. Is there a guide on the minimum amount of regularization I should add to make results insensitive to float32->float64 change?

In particular, doing $A=UDV^{T}$, I see that $L_\infty$ norm of $V^{T}X$ moves by about 1 when I change precision between float32 and float64. $L_2$ norm of $A$ is $10^5$ and it has about 200 zero eigenvalues out of 784 total.

Doing SVD on $\lambda I + A$ with $\lambda=10^{-3}$ made the difference vanish.

Answer 1

The singular value decomposition for a symmetric matrix $A=A^{T}$ is one and the same as its canonical eigendecomposition (i.e. with an orthonormal matrix-of-eigenvectors), while the same thing for a nonsymmetric matrix $M=U \Sigma V^T$ is just the canonical eigenvalue decomposition for the symmetric matrix $$ H=\begin{bmatrix}0 & M\\ M^{T} & 0 \end{bmatrix}=\begin{bmatrix}U & 0\\ 0 & V \end{bmatrix}\begin{bmatrix}0 & \Sigma\\ \Sigma & 0 \end{bmatrix}\begin{bmatrix}U & 0\\ 0 & V \end{bmatrix}^{T} $$ Hence, without loss of generality, let us consider a closely related question: If two symmetric matrices are approximately the same, then should we expect their canonical eigendecompositions to also be approximately the same?

The answer is a surprising no. Let $\epsilon>0$ be small, and consider the two matrices $$ A_{\epsilon}=\begin{bmatrix}1 & \epsilon\\ \epsilon & 1 \end{bmatrix}=V\Lambda_{\epsilon}V^{T},\qquad B_{\epsilon}=\begin{bmatrix}1+\epsilon & 0\\ 0 & 1-\epsilon \end{bmatrix}=U\Lambda_{\epsilon}U^{T} $$ both of which have eigenvalues $\Lambda_{\epsilon}=\mathrm{diag}(1+\epsilon,1-\epsilon)$, but whose eigenvectors are $$ V=\frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1\\ 1 & -1 \end{bmatrix},\qquad U=\begin{bmatrix}1 & 0\\ 0 & 1 \end{bmatrix}. $$ While the matrices $A_{\epsilon} \approx B_{\epsilon}$ are approximately the same, their matrix-of-eigenvectors $V$ and $U$ are very different. Indeed, since the eigendecompositions are unique for $\epsilon>0$, there really exists no choice of $U,V$ such that $U\approx V$

Now, applying this insight back to the SVD under finite precision, let us write $M_{0}=U_{0}\Sigma_{0}V_{0}^{T}$ as your matrix in float64 precision, and $M_{\epsilon}=U_{\epsilon}\Sigma_{\epsilon}V_{\epsilon}^{T}$ as the same matrix in float32 precision. If we assume that the SVDs themselves are exact, then the singular values $\Sigma_{0},\Sigma_{\epsilon}$ must differ by no more than a small constant factor of $\epsilon\approx10^{-7}$, but the singular vectors $U_{0},U_{\epsilon}$ and $V_{0},V_{\epsilon}$ can differ by an arbitrarily large quantity. Hence, as shown, there is no way to make the SVD "stable" in the sense of the singular vectors.

Answer 2

Although the question has a great answer, here's a rule of thumb for small singular values, with a plot.

If a singular value is nonzero but very small, you should define its reciprocal to be zero, since its apparent value is probably an artifact of roundoff error, not a meaningful number. A plausible answer to the question "how small is small ?" is to edit in this fashion all singular values whose ratio to the largest is less than $N$ times the machine precision $\epsilon$ .

$\qquad$ — Numerical Recipes p. 795

Added: the following couple of lines calculate this rule-of-thumb.

#!/usr/bin/env python2

from __future__ import division
import numpy as np
from scipy.sparse.linalg import svds  # sparse, dense or LinOp

#...............................................................................
def howsmall( A, singmax=None ):
    """ singular values < N float_eps sing_max  may be iffy, questionable
        "How small is small ?"
        [Numerical Recipes p. 795](http://apps.nrbook.com/empanel/index.html?pg=795)
    """
        # print "%d singular values are small, iffy" % (sing < howsmall(A)).sum()
        # small |eigenvalues| too ?
    if singmax is None:
        singmax = svds( A, 1, return_singular_vectors=False )[0]  # v0=random

    return max( A.shape ) * np.finfo( A.dtype ).eps * singmax

---

The Hilbert matrix seems to be widely used as a test case for roundoff error:

Here low-order bits in the mantissas of the Hilbert matrix are zeroed, A.astype(np.float__).astype(np.float64), then np.linalg.svd is run in float64. (Results with svd all float32 are about the same.)

Simply truncating to float32 might even be useful for denoising high-dimensional data, e.g. for train / test classification.

Real test cases would be welcome.