¿Cuánta regularización agregar para hacer SVD estable?

He estado usando SVD de Intel MKL (a dgesvdtravés de SciPy) y noté que los resultados son significativamente diferentes cuando cambio la precisión entre float32y float64cuando mi matriz está mal condicionada / no tiene rango completo. ¿Existe una guía sobre la cantidad mínima de regularización que debo agregar para que los resultados sean insensibles a float32-> float64cambios?

En particular, haciendo $A=UDV^{T}$, Veo que $L_\infty$ norma de $V^{T}X$se mueve aproximadamente 1 cuando cambio la precisión entre float32y float64.$L_2$ norma de $A$ es $10^5$ y tiene alrededor de 200 valores propios cero de un total de 784.

Haciendo SVD en $\lambda I + A$ con $\lambda=10^{-3}$ hizo que la diferencia se desvaneciera.

Aunque la pregunta tiene una gran respuesta, aquí hay una regla general para pequeños valores singulares, con un gráfico.

Si un valor singular no es cero pero es muy pequeño, debe definir su recíproco como cero, ya que su valor aparente es probablemente un artefacto de error de redondeo, no un número significativo. Una respuesta plausible a la pregunta "¿qué tan pequeño es pequeño?" es editar de esta manera todos los valores singulares cuya relación con el mayor es menor que$N$ veces la precisión de la máquina $\epsilon$ .

$\qquad$- Recetas numéricas p. 795

Agregado: las siguientes dos líneas calculan esta regla general.

#!/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

---

La matriz de Hilbert parece ser ampliamente utilizada como un caso de prueba para el error de redondeo:

Aquí los bits de orden inferior en las mantisas de la matriz de Hilbert se ponen a cero A.astype(np.float__).astype(np.float64)y luego np.linalg.svdse introducen float64. (Los resultados con svdtodos float32son casi iguales).

Simplemente truncar float32podría incluso ser útil para eliminar datos de alta dimensión, por ejemplo, para la clasificación de trenes / pruebas.

Casos de prueba reales serían bienvenidos.

La descomposición del valor singular para una matriz simétrica. $A=A^{T}$ es uno y el mismo que su descomposición propia canónica (es decir, con una matriz de vectores propios ortonormales), mientras que lo mismo para una matriz no simétrica $M=U \Sigma V^T$ es solo la descomposición del valor propio canónico para la matriz simétrica

$$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}$$ Por lo tanto, sin pérdida de generalidad, consideremos una pregunta estrechamente relacionada: si dos matrices simétricas son aproximadamente iguales, ¿deberíamos esperar que sus descomposiciones propias canónicas también sean aproximadamente las mismas?

La respuesta es un sorprendente no. Dejar$\epsilon>0$ ser pequeño y considerar las dos 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}$$ ambos tienen valores propios $\Lambda_{\epsilon}=\mathrm{diag}(1+\epsilon,1-\epsilon)$, pero cuyos vectores propios son $$V=\frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1\\ 1 & -1 \end{bmatrix},\qquad U=\begin{bmatrix}1 & 0\\ 0 & 1 \end{bmatrix}.$$ Mientras las matrices $A_{\epsilon} \approx B_{\epsilon}$ son aproximadamente iguales, sus vectores de matriz de vectores $V$ y $U$son muy diferentes. De hecho, dado que las descomposiciones propias son únicas para$\epsilon>0$, realmente no hay elección de $U,V$ tal que $U\approx V$

Ahora, aplicando esta información al SVD con precisión finita, escribamos $M_{0}=U_{0}\Sigma_{0}V_{0}^{T}$como su matriz en float64 precisión, y$M_{\epsilon}=U_{\epsilon}\Sigma_{\epsilon}V_{\epsilon}^{T}$ como la misma matriz en float32precisión. Si suponemos que los SVD son exactos, entonces los valores singulares$\Sigma_{0},\Sigma_{\epsilon}$ debe diferir en no más que un pequeño factor constante de $\epsilon\approx10^{-7}$, pero los vectores singulares $U_{0},U_{\epsilon}$ y $V_{0},V_{\epsilon}$ puede diferir en una cantidad arbitrariamente grande. Por lo tanto, como se muestra, no hay forma de hacer que la SVD sea "estable" en el sentido de los vectores singulares.