Instabilidade numérica do cálculo da matriz de covariância inversa

8

Eu tenho 65 amostras de dados 21-dimensionais (colados aqui ) e estou construindo a matriz de covariância a partir dele. Quando computado em C ++, recebo a matriz de covariância colada aqui . E quando computado no Matlab a partir dos dados (como mostrado abaixo), recebo a matriz de covariância colada aqui

Código Matlab para computação cov a partir de dados:

data = csvread('path/to/data');
matlab_cov = cov(data);

Como você pode ver, as diferenças nas matrizes de covariância são mínimas (~ e-07), o que provavelmente se deve a problemas numéricos no compilador usando aritmética de ponto flutuante.

No entanto, quando computo a matriz de covariância pseudo-inversa da matriz de covariância produzida pelo matlab e a produzida pelo meu código C ++, obtenho resultados amplamente diferentes. Estou computando-os da mesma maneira, ou seja:

data = csvread('path/to/data');
matlab_cov = cov(data);
my_cov = csvread('path/to/cov_file');
matlab_inv = pinv(matlab_cov);
my_inv = pinv(my_cov);

A diferença é tão grande que, quando estou computando a distância dos mahalanobis de uma amostra (colada aqui ) até a distribuição das 65 amostras por:

(65/642)×((samplemean)×1×(samplemean))

usando as diferentes matrizes de covariância inversa ( ), obtenho resultados amplamente diferentes, ou seja:1

 (65/(64^2))*((sample-sample_mean)*my_inv*(sample-sample_mean)')
ans =

   1.0167e+05

(65/(64^2))*((sample-sample_mean)*matlab_inv*(sample-sample_mean)')
ans =

  109.9612

É normal que as pequenas diferenças (e-7) na matriz de covariância tenham esse efeito no cálculo da matriz pseudo-inversa? E se sim, o que posso fazer para mitigar esse efeito?

Caso contrário, existem outras métricas de distância que eu possa usar que não envolvam a covariância inversa? Eu uso a distância de Mahalanobis como sabemos para n amostras segue uma distribuição beta, que eu uso para testes de hipóteses

Muito obrigado antecipadamente

EDIT: Adicionando código C ++ para calcular a matriz de covariância abaixo: O vector<vector<double> >representa a coleção de linhas do arquivo colado.

Mat covariance_matrix = Mat(21, 21, CV_32FC1, cv::Scalar(0));
    for(int j = 0; j < 21; j++){
        for(int k = 0; k < 21; k++){
            for(std::vector<vector<double> >::iterator it = data.begin(); it!= data.end(); it++){
                covariance_matrix.at<float>(j,k) += (it->at(j) - mean.at(j)) * (it->at(k) - mean[k]);
            }
            covariance_matrix.at<float>(j,k) /= 64; 
        }
    }
Aly
fonte
Inverter matrizes ..... Isso é perigoso! Geralmente é preferível encontrar alternativas para isso (por exemplo, pseudoinverso)
Ander Biguri
1
@ Aly: as matrizes que você deseja inverter não são matrizes de covariâncias "válidas" porque não são definidas positivamente; numericamente, eles ainda têm alguns autovalores negativos (mas próximos de zero). Eu provavelmente adicionaria uma constante muito pequena ao longo da diagonal; é realmente uma forma de correção de Tikhonov ( ). Além disso, não use flutuadores, use duplos para armazenar sua matriz de covariância. (E, além disso você já usa OpenCV, assim como você pode usar Eigen ou Armadillo ..)Χ+λI
usεr11852
1
@ Aly: Wikipedia, realmente. (é o lema: regularização de Tikhonov). O método mencionado pelo whuber usando o SVD fornecerá uma matriz definida não negativa se você definir valores próprios pequenos para zero; você ainda precisará adicionar uma pequena constante a todos os seus autovalores para torná-los positivos definidos. Praticamente os dois métodos fazem o mesmo. Acabei de não usar o SVD, mas afetou diretamente os autovalores das amostras adicionando a todos eles. Não encontrei nenhuma referência, acredito que ambos os métodos sejam bastante intuitivos. λ
usεr11852
1
@ user11852 Você pode fazer dos seus comentários uma resposta, ainda estou experimentando, mas se a promessa aceitar. Além disso, se outros fazem suas sugestões respostas que vai até voto como têm sido muito útil / útil para o meu entendimento do problema
Aly
1
Comentei em seu outro segmento que ter variáveis ​​que somam 1 , como o conjunto de dados, incentiva a instabilidade e contém uma variável redundante. Tente soltar uma coluna. Você nem precisa do pinv: a matriz de covariância não é mais singular.
Cam.Davidson.Pilon

Respostas:

7

As matrizes que você deseja inverter não são matrizes de covariâncias "válidas" porque não são definidas positivamente; numericamente, eles ainda têm alguns autovalores negativos (mas próximos de zero). Provavelmente devido a zeros da máquina, por exemplo, o último valor próprio da matriz "matlab_covariance" é -0,000000016313723. Para corrigir a definição positiva, você pode fazer duas coisas:

  1. Basta adicionar uma constante muito pequena ao longo da diagonal; realmente uma forma de correção de Tikhonov ( com ).Χ+λIλ0
  2. (Com base no que o proponente propôs) Use SVD, defina os autovalores "problemáticos" para um valor pequeno fixo (não zero), reconstrua sua matriz de covariância e depois inverta isso. Claramente, se você definir alguns desses autovalores para zero, você terminará com uma matriz não negativa (ou semi positiva), que ainda não será invertível.

Uma matriz não negativa não tem um inverso, mas possui um pseudo inverso (todas as matrizes com entradas reais ou complexas têm um pseudo inverso); no entanto, o pseudo inverso de Moore – Penrose é mais computacionalmente caro do que um inverso verdadeiro e, se o inverso existe, é igual ao pseudo-inverso. Então, basta ir para o inverso :)

Ambos os métodos praticamente tentam manipular os valores próprios que são avaliados como zero (ou abaixo de zero). O primeiro método é um pouco ondulado, mas provavelmente muito mais rápido de implementar. Para algo um pouco mais estável, convém calcular o SVD e definir o igual ao absoluto do menor valor próprio (para que você não seja negativo) mais algo muito pequeno (para que seja positivo). Apenas tome cuidado para não impor positividade a uma matriz que é obviamente negativa (ou já positiva). Ambos os métodos alterarão o número de condicionamento da sua matriz.λ

Em termos estatísticos, o que você faz adicionando na diagonal de sua matriz de covariância adiciona ruído às suas medições. (Como a diagonal da matriz de covariância é a variação de cada ponto e, ao adicionar algo a esses valores, você apenas diz "a variação nos pontos em que tenho leituras é realmente um pouco maior do que eu pensava originalmente".)λ

Um teste rápido para a definição positiva de uma matriz é a existência (ou não) da decomposição de Cholesky.

Também como uma nota computacional:

  1. Não use carros alegóricos, use duplos para armazenar sua matriz de covariância.
  2. Use bibliotecas de álgebra linear numérica em C ++ (como Eigen ou Armadillo) para obter inversões de matrizes, produtos matriciais, etc. É mais rápido, seguro e conciso.

EDIT: Dado que você tem uma decomposição de Cholesky da sua matriz tal forma que (você precisa fazer isso para verificar se está tendo uma matriz Pos.Def.), Você poderá resolver imediatamente o sistema . Você apenas resolve Ly = b para y por substituição direta e, em seguida, L ^ Tx = y para x por substituição inversa. (Em si, basta usar o método .solve (x) do seu objeto Cholesky) Obrigado a bnaul e Zen por apontar que eu me concentrei tanto em obter o ser Pos.Def. que eu esqueci por que nos importamos com isso em primeiro lugar :)KLLTKx=bK

usεr11852
fonte
3
+1. Usando Mathematica e aplicando-à dados (em vez da matriz de covariância postado, que podem ter sido apresentados com muito pouca precisão) I encontrar valores próprios negativos. É assim que deve ser: quando uma matriz de covariância é calculada exatamente, é garantida uma semi-definida positiva; portanto, quaisquer autovalores negativos devem ser atribuídos à imprecisão nos cálculos. Qualquer procedimento inverso generalizado decente deve "reconhecer" esses minúsculos valores negativos como zeros e tratá-los de acordo.
whuber
Obrigado pessoal pelo esforço, como declarado. Eu votei e vou testá-los e comentar ou aceitar de acordo.
Aly
Desculpe, estou um pouco confuso, como resolver o Cholesky usa a distância de Mahalanobis?
Aly
Verifique o link na postagem original do bnaul. Mas não use mas Cholesky (é isso que eles querem dizer com LDL *). LU
usεr11852
2

As respostas e comentários postados são bons pontos sobre os perigos da inversão de matrizes quase singulares. No entanto, até onde eu sei, ninguém mencionou que calcular a distância de Mahalanobis não requer, na verdade, inverter a covariância da amostra. Consulte esta pergunta StackOverflow para obter uma descrição de como fazer isso usando a decomposição da .LU

O princípio é o mesmo que resolver um sistema linear: ao tentar resolver modo que , existem métodos muito mais eficientes e numericamente estáveis ​​do que usar .xAx=bx=A1b

Edit: provavelmente escusado será dizer, mas esse método produz o valor exato da distância, enquanto a adição de a e a inversão produz apenas uma aproximação.λIS

bnaul
fonte
1
Você está certo, @bnaul. No entanto, sem algum tipo de regularização, a LUdecomposição também não funcionará. Vou adicionar um comentário sobre isso na minha resposta.
Zen
@ Bnaul: por que a LU quando você faz com o Cholesky para verificar a certeza positiva? Supondo que você tenha uma matriz de covariância válida, resolvendo para y por substituição direta e, em seguida, para x por substituição traseira será mais rápida. Bom ponto, porém, definitivamente me concentro em obter uma certeza positiva que esqueci por que me importava com isso originalmente! : DK=LLTLy=bLTx=y
usεr11852 28/02
0

(Anos depois), um pequeno exemplo: com deficiente em classificação, autovalores de serão de 0 a precisão da máquina - e cerca de metade desses "zeros" pode ser :r < n , n - r A T A < 0Ar<n, nrATA<0

#!/usr/bin/env python2
""" many eigenvalues of A'A are tiny but < 0 """
# e.g. A 1 x 10: [-1.4e-15 -6.3e-17 -4e-17 -2.7e-19 -8.8e-21  1e-18 1.5e-17 5.3e-17 1.4e-15  7.7]

from __future__ import division
import numpy as np
from numpy.linalg import eigvalsh  # -> lapack_lite
# from scipy.linalg import eigvalsh  # ~ same
from scipy import __version__

np.set_printoptions( threshold=20, edgeitems=10, linewidth=140,
        formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
print "versions: numpy %s  scipy %s \n" % (
        np.__version__, __version__  )

np.random.seed( 3 )

rank = 1
n = 10
A = np.random.normal( size=(rank, n) )
print "A: \n", A
AA = A.T.dot(A)
evals = eigvalsh( AA )
print "eigenvalues of A'A:", evals
denis
fonte