Ile regulacji należy dodać, aby zapewnić stabilność SVD?

Korzystałem z SVD Intela MKL ( dgesvdprzez SciPy) i zauważyłem, że wyniki są znacząco różne, kiedy zmieniam precyzję pomiędzy float32i float64kiedy moja matryca jest źle uwarunkowana / nie w pełni ustawiona. Czy istnieje przewodnik dotyczący minimalnej ilości regularyzacji, którą powinienem dodać, aby wyniki były niewrażliwe na float32-> float64zmianę?

W szczególności robienie $A=UDV^{T}$, Widzę to $L_\infty$ norma $V^{T}X$przesuwa się o około 1, gdy zmieniam precyzję pomiędzy float32i float64.$L_2$ norma $A$ jest $10^5$ i ma około 200 zerowych wartości własnych z 784 ogółem.

Włączanie SVD $\lambda I + A$ z $\lambda=10^{-3}$ sprawiło, że różnica zniknęła.

Chociaż pytanie ma świetną odpowiedź, oto ogólna zasada dla małych pojedynczych wartości z fabułą.

Jeśli pojedyncza wartość jest niezerowa, ale bardzo mała, należy zdefiniować jej odwrotność jako zero, ponieważ jej pozorna wartość jest prawdopodobnie artefaktem błędu zaokrąglenia, a nie znaczącą liczbą. Prawdopodobna odpowiedź na pytanie „jak mały jest mały?” polega na edycji w ten sposób wszystkich pojedynczych wartości, których stosunek do największej jest mniejszy niż$N$ razy dokładność maszyny $\epsilon$ .

$\qquad$- Przepisy numeryczne str. 795

Dodano: poniższe kilka linii oblicza tę zasadę.

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

---

Macierz Hilberta wydaje się być szeroko stosowana jako przypadek testowy dla błędu zaokrąglenia:

Tutaj bity niskiego rzędu w mantysach macierzy Hilberta są zerowane A.astype(np.float__).astype(np.float64), a następnie np.linalg.svdsą wprowadzane float64. (Wyniki dla svdwszystkich float32są prawie takie same.)

Po prostu obcięcie do float32może być nawet przydatne do odrzucania danych wielowymiarowych, np. Do klasyfikacji pociągu / testu.

Prawdziwe przypadki testowe byłyby mile widziane.

Dekompozycja wartości osobliwych dla macierzy symetrycznej $A=A^{T}$ jest taki sam jak jego kanoniczny skład eigend (tj. z ortonormalną matrycą wektorów własnych), podczas gdy to samo dotyczy macierzy niesymetrycznej $M=U \Sigma V^T$ jest tylko kanonicznym rozkładem wartości własnej macierzy symetrycznej

$$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}$$ Zatem bez utraty ogólności rozważmy ściśle powiązane pytanie: jeśli dwie macierze symetryczne są w przybliżeniu takie same, to czy powinniśmy oczekiwać, że ich kanoniczne składy eigend będą w przybliżeniu takie same?

Odpowiedź brzmi zaskakująco nie. Pozwolić$\epsilon>0$ bądź mały i rozważ dwie matryce

$$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}$$ oba mają wartości własne $\Lambda_{\epsilon}=\mathrm{diag}(1+\epsilon,1-\epsilon)$, ale których wektorami własnymi są $$V=\frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1\\ 1 & -1 \end{bmatrix},\qquad U=\begin{bmatrix}1 & 0\\ 0 & 1 \end{bmatrix}.$$ Podczas gdy matryce $A_{\epsilon} \approx B_{\epsilon}$ są w przybliżeniu takie same, ich matryce wektorów własnych $V$ i $U$są bardzo różne. Rzeczywiście, ponieważ kompozycje eigend są wyjątkowe$\epsilon>0$, naprawdę nie ma wyboru $U,V$ takie, że $U\approx V$

Teraz, stosując tę ​​wiedzę z powrotem do SVD ze skończoną precyzją, napiszmy $M_{0}=U_{0}\Sigma_{0}V_{0}^{T}$jako float64 precyzyjna matryca , oraz$M_{\epsilon}=U_{\epsilon}\Sigma_{\epsilon}V_{\epsilon}^{T}$ jak ta sama matryca w float32precyzji. Jeśli założymy, że same SVD są dokładne, to liczby osobliwe$\Sigma_{0},\Sigma_{\epsilon}$ musi różnić się nie więcej niż małym stałym współczynnikiem wynoszącym $\epsilon\approx10^{-7}$, ale pojedyncze wektory $U_{0},U_{\epsilon}$ i $V_{0},V_{\epsilon}$ może różnić się dowolnie dużą ilością. Dlatego, jak pokazano, nie ma sposobu, aby SVD była „stabilna” w sensie pojedynczych wektorów.