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

10

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=UDVT, Veo que L norma de VTXse mueve aproximadamente 1 cuando cambio la precisión entre float32y float64.L2 norma de A es 105 y tiene alrededor de 200 valores propios cero de un total de 784.

Haciendo SVD en λI+A con λ=103 hizo que la diferencia se desvaneciera.

Yaroslav Bulatov
fuente
Cual es el tamaño N de una N×N matriz Apara ese ejemplo (¿es incluso una matriz cuadrada)? 200 valores propios cero o valores singulares? Una norma de Frobenius||A||Fpara un ejemplo representativo también sería útil.
Anton Menshov
En este caso matriz de 784 x 784, pero estoy más interesado en la técnica general para encontrar un buen valor de lambda
Yaroslav Bulatov
Entonces, es la diferencia en Vsolo en las últimas columnas correspondientes a los valores singulares cero?
Nick Alger
2
Si hay varios valores singulares iguales, el svd no es único. En su ejemplo, supongo que el problema proviene de los múltiples valores singulares cero y que una precisión diferente conduce a una elección diferente de la base para el espacio singular respectivo. No sé por qué eso cambia cuando se regulariza ...
Dirk
1
...que es X?
Federico Poloni

Respuestas:

1

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 queN veces la precisión de la máquina ϵ .

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

ingrese la descripción de la imagen aquí

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.

denis
fuente
por cierto, scipy parece agregar un factor de 1e3 para float32 y 1e6 para float64, curioso de dónde provienen
Yaroslav Bulatov
@Yaroslav Bulatov, numpyy scipy.linalg.svdllame a LAPACK gesdd , vea el parámetro JOBRen dgejsv: "Especifica el RANGO para los valores singulares. Emite la licencia para poner a cero pequeños valores singulares positivos si están fuera ..." ( scipy.sparse.linalg.svdsenvuelve ARPACK y tiene un parámetro tol, Tolerancia para valores singulares.)
denis
13

La descomposición del valor singular para una matriz simétrica. A=AT 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ΣVT es solo la descomposición del valor propio canónico para la matriz simétrica

H=[0MMT0]=[U00V][0ΣΣ0][U00V]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ϵ>0 ser pequeño y considerar las dos matrices

Aϵ=[1ϵϵ1]=VΛϵVT,Bϵ=[1+ϵ001ϵ]=UΛϵUT
ambos tienen valores propios Λϵ=diag(1+ϵ,1ϵ), pero cuyos vectores propios son
V=12[1111],U=[1001].
Mientras las matrices AϵBϵ son aproximadamente iguales, sus vectores de matriz de vectores V y Uson muy diferentes. De hecho, dado que las descomposiciones propias son únicas paraϵ>0, realmente no hay elección de U,V tal que UV

Ahora, aplicando esta información al SVD con precisión finita, escribamos M0=U0Σ0V0Tcomo su matriz en float64 precisión, yMϵ=UϵΣϵVϵT como la misma matriz en float32precisión. Si suponemos que los SVD son exactos, entonces los valores singularesΣ0,Σϵ debe diferir en no más que un pequeño factor constante de ϵ107, pero los vectores singulares U0,Uϵ y V0,Vϵ 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.

Richard Zhang
fuente
1
Esa es una gran referencia. No sé, aprendí este ejemplo en particular hace muchos años en la clase de matemáticas :-)
Richard Zhang