Quanta regularização adicionar para tornar o SVD estável?

Eu tenho usado o SVD da Intel MKL ( dgesvdatravés do SciPy) e notei que os resultados são significativamente diferentes quando eu mudo a precisão entre float32e float64quando minha matriz está mal condicionada / não está na classificação completa. Existe um guia sobre a quantidade mínima de regularização que devo adicionar para tornar os resultados insensíveis a float32-> float64alteração?

Em particular, fazendo , eu ver que L ∞ norma de V T X move-se em cerca de 1 quando mudar entre precisão e . A norma L 2 de A é 10 5 e tem cerca de 200 autovalores zero de um total de 784.$A=UDV^{T}$$L_\infty$$V^{T}X$float32float64$L_2$$A$$10^5$

Fazer SVD em com λ = 10 - 3 fez a diferença desaparecer.$\lambda I + A$$\lambda=10^{-3}$

Embora a pergunta tenha uma ótima resposta, aqui está uma regra prática para pequenos valores singulares, com um gráfico.

Se um valor singular for diferente de zero, mas muito pequeno, você deve definir seu recíproco como zero, pois seu valor aparente provavelmente é um artefato de erro de arredondamento, não um número significativo. Uma resposta plausível à pergunta "quão pequeno é pequeno?" é editar desta forma todos os valores singulares cuja razão para o maior é menor que vezes a precisão da máquina ϵ .$N$$\epsilon$

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

Adicionado: as duas linhas a seguir calculam essa regra de ouro.

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

---

A matriz Hilbert parece ser amplamente utilizada como um caso de teste para erro de arredondamento:

Aqui, os bits de baixa ordem nas mantissas da matriz Hilbert são zerados A.astype(np.float__).astype(np.float64)e depois np.linalg.svdexecutados float64. (Os resultados com svdtodos float32são praticamente os mesmos.)

Simplesmente truncar float32pode até ser útil para excluir dados de alta dimensão, por exemplo, para classificação de trem / teste.

Casos de teste reais seriam bem-vindos.

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