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

10

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=UDVTLVTXfloat32float64L2A105

Fazer SVD em com λ = 10 - 3 fez a diferença desaparecer.λI+Aλ=103

Yaroslav Bulatov
fonte
O que é o tamanho de um N × N matriz A para que o exemplo (se ele mesmo uma matriz quadrada)? 200 autovalores zero ou valores singulares? Uma norma de Frobenius | | Um | | F para um exemplo representativo também seria útil. NN×NA||A||F
Anton Menshov
Neste caso 784 x 784 matriz, mas eu estou mais interessado em técnica geral para encontrar um bom valor de lambda
Yaroslav Bulatov
Então, a diferença em apenas nas últimas colunas corresponde aos valores singulares zero? V
Nick Alger #
2
Se houver vários valores singulares iguais, o svd não é exclusivo. No seu exemplo, acho que o problema vem dos múltiplos valores singulares de zero e que uma precisão diferente leva a uma escolha diferente da base para o respectivo espaço singular. Eu não sei por que isso muda quando você regularizar ...
Dirk
11
... o que é ? X
Federico Poloni

Respostas:

1

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ϵ

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

insira a descrição da imagem aqui

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.

denis
fonte
btw, scipy parece acrescentar um fator de 1e3 para float32 e 1E6 para float64, curioso de onde estes vieram
Yaroslav Bulatov
@Yaroslav Bulatov, numpye scipy.linalg.svdchame LAPACK gesdd , consulte o parâmetro JOBRem dgejsv: "Especifica o RANGE para os valores singulares. Emite a licença para definir zero valores singulares positivos pequenos se eles estiverem fora ..." ( scipy.sparse.linalg.svdsquebra ARPACK e tem um parâmetro tol, Tolerance para valores singulares.)
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
fonte
11
Essa é uma ótima referência. Eu não sei, eu aprendi neste exemplo em particular, há muitos anos na aula de matemática :-)
Richard Zhang