Wie viel Regularisierung muss hinzugefügt werden, um SVD stabil zu machen?

10

Ich habe die SVD von Intel MKL ( dgesvdüber SciPy) verwendet und festgestellt, dass die Ergebnisse erheblich voneinander abweichen, wenn ich die Genauigkeit zwischen float32und float64wenn meine Matrix schlecht konditioniert ist / nicht den vollen Rang hat. Gibt es einen Leitfaden zum Mindestmaß an Regularisierung, das ich hinzufügen sollte, um die Ergebnisse unempfindlich gegenüber float32-> float64Änderungen zu machen?

Insbesondere wenn ich mache , sehe ich, dass sich die L -Norm von V T X um ungefähr 1 bewegt, wenn ich die Genauigkeit zwischen und ändere . Die L 2 -Norm von A ist 10 5 und hat ungefähr 200 Null-Eigenwerte von insgesamt 784.A=UDVTLVTXfloat32float64L2A105

Durch SVD auf mit λ = 10 - 3 verschwand der Unterschied.λI+Aλ=103

Jaroslaw Bulatow
quelle
Was ist die Größe einer N × N- Matrix A für dieses Beispiel (ist es sogar eine quadratische Matrix)? 200 Null-Eigenwerte oder Singularwerte? Eine Frobenius-Norm | | A | | F für ein repräsentatives Beispiel wäre ebenfalls hilfreich. NN×NA||A||F
Anton Menshov
In diesem Fall 784 x 784 Matrix, aber ich bin mehr an der allgemeinen Technik interessiert, um einen guten Wert von Lambda zu finden
Jaroslaw Bulatow
Entspricht die Differenz in nur in den letzten Spalten den Null-Singularwerten? V
Nick Alger
2
Wenn es mehrere gleiche Singularwerte gibt, ist die SVD nicht eindeutig. In Ihrem Beispiel denke ich, dass das Problem von den mehreren Null-Singularwerten herrührt und dass eine andere Genauigkeit zu einer anderen Wahl der Basis für den jeweiligen Singularraum führt. Ich weiß nicht, warum sich das ändert, wenn Sie regulieren ...
Dirk
1
X

Antworten:

1

Obwohl die Frage eine gute Antwort hat, ist hier eine Faustregel für kleine Singularwerte mit einer Darstellung.

Nϵ

- Numerische Rezepte p. 795

Hinzugefügt: Die folgenden Zeilen berechnen diese Faustregel.

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


Die Hilbert-Matrix scheint als Testfall für Rundungsfehler weit verbreitet zu sein:

Geben Sie hier die Bildbeschreibung ein

Hier Bits niedriger Ordnung in die Mantissen der Hilbert - Matrix werden auf Null gesetzt, A.astype(np.float__).astype(np.float64)und dann np.linalg.svdzulaufen float64. (Die Ergebnisse sind bei svdallen float32ungefähr gleich.)

Das einfache Abschneiden auf float32kann sogar nützlich sein, um hochdimensionale Daten zu entrauschen, z. B. für die Zug- / Testklassifizierung.

Echte Testfälle wären willkommen.

denis
quelle
Übrigens scheint scipy einen Faktor von 1e3 für float32 und 1e6 für float64 hinzuzufügen, neugierig, woher diese kamen
Yaroslav Bulatov
@Yaroslav Bulatov, numpyund scipy.linalg.svdrufen Sie LAPACK gesdd auf , siehe Parameter JOBRin dgejsv: "Gibt den BEREICH für die Singularwerte an. Gibt die Lizenz aus, kleine positive Singularwerte auf Null zu setzen, wenn sie außerhalb liegen ..." ( scipy.sparse.linalg.svdsumschließt ARPACK und hat den Parameter tolToleranz für singuläre Werte.)
denis
13

A=ATM=UΣVT

H=[0MMT0]=[U00V][0ΣΣ0][U00V]T

ϵ>0

Aϵ=[1ϵϵ1]=VΛϵVT,Bϵ=[1+ϵ001ϵ]=UΛϵUT
Λϵ=diag(1+ϵ,1ϵ)
V=12[1111],U=[1001].
AϵBϵVUϵ>0U,VUV

M0=U0Σ0V0Tfloat64Mϵ=UϵΣϵVϵTfloat32Σ0,Σϵϵ107U0,UϵV0,Vϵ

Richard Zhang
quelle
1
Das ist eine großartige Referenz. Ich weiß nicht, ich habe dieses spezielle Beispiel vor vielen Jahren im Matheunterricht gelernt :-)
Richard Zhang