Combien de régularisation ajouter pour rendre SVD stable?

J'ai utilisé le SVD d'Intel MKL ( dgesvdvia SciPy) et j'ai remarqué que les résultats sont considérablement différents lorsque je change de précision entre float32et float64lorsque ma matrice est mal conditionnée / n'est pas pleine. Existe-t-il un guide sur le montant minimum de régularisation à ajouter pour rendre les résultats insensibles aux float32-> float64changements?

En particulier, en faisant , je vois que la norme L ∞ de V T X se déplace d'environ 1 lorsque je change de précision entre et . La norme L 2 de A est 10 5 et elle a environ 200 valeurs propres nulles sur 784 au total.$A=UDV^{T}$$L_\infty$$V^{T}X$float32float64$L_2$$A$$10^5$

Faire SVD sur avec λ = 10 - 3 a fait disparaître la différence.$\lambda I + A$$\lambda=10^{-3}$

Bien que la question ait une excellente réponse, voici une règle de base pour les petites valeurs singulières, avec un tracé.

Si une valeur singulière est non nulle mais très petite, vous devez définir sa réciproque à zéro, car sa valeur apparente est probablement un artefact d'erreur d'arrondi, pas un nombre significatif. Une réponse plausible à la question "comment petit est petit?" consiste à éditer de cette façon toutes les valeurs singulières dont le rapport à la plus grande est inférieur à fois la précision de la machine ϵ .$N$$\epsilon$

$\qquad$- Recettes numériques p. 795

Ajouté: les deux lignes suivantes calculent cette règle empirique.

#!/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 matrice de Hilbert semble être largement utilisée comme cas de test pour l'erreur d'arrondi:

Ici, les bits de poids faible dans les mantisses de la matrice de Hilbert sont mis à zéro A.astype(np.float__).astype(np.float64), puis np.linalg.svdexécutés float64. (Les résultats avec svdtous float32sont à peu près les mêmes.)

Une simple troncature à float32pourrait même être utile pour débruiter des données de grande dimension, par exemple pour la classification des trains / tests.

De vrais cas de test seraient les bienvenus.

$A=A^{T}$$M=U \Sigma V^T$

$$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}$$

$\epsilon>0$

$$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}$$ $\Lambda_{\epsilon}=\mathrm{diag}(1+\epsilon,1-\epsilon)$ $$V=\frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1\\ 1 & -1 \end{bmatrix},\qquad U=\begin{bmatrix}1 & 0\\ 0 & 1 \end{bmatrix}.$$ $A_{\epsilon} \approx B_{\epsilon}$$V$$U$$\epsilon>0$$U,V$$U\approx V$

$M_{0}=U_{0}\Sigma_{0}V_{0}^{T}$float64$M_{\epsilon}=U_{\epsilon}\Sigma_{\epsilon}V_{\epsilon}^{T}$float32$\Sigma_{0},\Sigma_{\epsilon}$$\epsilon\approx10^{-7}$$U_{0},U_{\epsilon}$$V_{0},V_{\epsilon}$