Instabilité numérique du calcul de la matrice de covariance inverse

8

J'ai 65 échantillons de données en 21 dimensions (collées ici ) et j'en construis la matrice de covariance. Lorsque calculé en C ++, je récupère la matrice de covariance collée ici . Et lorsque calculé dans matlab à partir des données (comme indiqué ci-dessous), je reçois la matrice de covariance collée ici

Code Matlab pour calculer la cov à partir des données:

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

Comme vous pouvez le voir, la différence dans les matrices de covariance est minuscule (~ e-07), ce qui est probablement dû à des problèmes numériques dans le compilateur utilisant l'arithmétique à virgule flottante.

Cependant, lorsque je calcule la matrice de covariance pseudo-inverse à partir de la matrice de covariance produite par matlab et celle produite par mon code C ++, j'obtiens des résultats très différents. Je les calcule de la même manière c'est à dire:

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);

La différence est si énorme que lorsque je calcule la distance mahalanobis d'un échantillon (collé ici ) à la distribution des 65 échantillons par:

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

en utilisant les différentes matrices de covariance inverse ( ) j'obtiens des résultats très différents c'est-à-dire: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

Est-il normal que les petites différences (e-7) de matrice de covariance aient un tel effet sur le calcul de la matrice pseudo-inverse? Et si oui, que puis-je faire pour atténuer cet effet?

À défaut, y a-t-il d'autres mesures de distance que je peux utiliser qui n'impliquent pas la covariance inverse? J'utilise la distance de Mahalanobis comme nous le savons pour n échantillons, elle suit une distribution bêta, que j'utilise pour les tests d'hypothèse

Merci d'avance

EDIT: Ajout de code C ++ pour calculer la matrice de covariance ci-dessous: Le vector<vector<double> >représente la collection de lignes du fichier collé.

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
la source
Inverser les matrices ..... C'est une chose dangereuse! Habituellement, il est préférable de trouver des alternatives à cela (par exemple pseudoinverse)
Ander Biguri
1
@Aly: les matrices que vous cherchez à inverser ne sont pas des matrices de covariances "valides" car elles ne sont pas définies positives; numériquement, ils ont même des valeurs propres négatives (mais proches de zéro). J'ajouterais probablement juste une très petite constante le long de la diagonale; c'est vraiment une forme de correction de Tikhonov ( ). N'utilisez pas non plus de flottants, utilisez des doubles pour stocker votre matrice de covariance. (Et en plus vous utilisez déjà OpenCV, vous pourriez aussi bien utiliser Eigen ou Armadillo ..)Χ+λI
usεr11852
1
@Aly: Wikipedia, vraiment. (c'est le lemme: régularisation de Tikhonov). La méthode que Whuber a mentionnée en utilisant le SVD vous donnera une matrice définie non négative si vous définissez de petites valeurs propres à zéro; vous devrez toujours ajouter une petite constante à toutes vos valeurs propres pour les rendre définitives positives. Pratiquement les deux méthodes font de même. Je me suis contenté de ne pas utiliser le SVD mais d'affecter directement les valeurs propres des échantillons en ajoutant à chacun d'eux. Je n'ai trouvé aucune référence, les deux méthodes sont assez intuitives je crois. λ
usεr11852
1
@ user11852 S'il vous plaît pouvez-vous faire une réponse à vos commentaires, j'expérimente toujours, mais si prometteur acceptera. De plus, si d'autres font leurs suggestions, je voterai car elles ont été très utiles / utiles à ma compréhension du problème
Aly
1
J'ai commenté dans votre autre thread que le fait d'avoir des variables totalisant 1 , comme le fait votre ensemble de données, encourage l'instabilité et contient une variable redondante. Veuillez supprimer une colonne. Vous n'avez même pas besoin du pinv: la matrice de covariance n'est plus singulière.
Cam.Davidson.Pilon

Réponses:

7

Les matrices que vous cherchez à inverser ne sont pas des matrices de covariances "valides" car elles ne sont pas définies positives; numériquement, ils ont même des valeurs propres négatives (mais proches de zéro). Cela est probablement dû aux zéros de la machine, par exemple la dernière valeur propre de votre matrice "matlab_covariance" est -0,000000016313723. Pour corriger au positif défini, vous pouvez faire deux choses:

  1. Ajoutez simplement une très petite constante le long de la diagonale; une forme de correction de Tikhonov vraiment ( avec ).Χ+λIλ0
  2. (Sur la base de ce qui est proposé) Utilisez SVD, définissez les valeurs propres "problématiques" sur une petite valeur fixe (non nulle), reconstruisez votre matrice de covariance, puis inversez-la. De toute évidence, si vous définissez certaines de ces valeurs propres à zéro, vous vous retrouverez avec une matrice non négative (ou semi-positive), qui ne sera toujours pas inversible.

Une matrice non négative n'a pas d'inverse mais elle a un pseudo inverse (toutes les matrices avec des entrées réelles ou complexes ont un pseudo-inverse), néanmoins le pseudo-inverse de Moore – Penrose est plus coûteux en calcul qu'un vrai inverse et si l'inverse existe il est égal au pseudo-inverse. Alors optez pour l'inverse :)

Les deux méthodes essaient pratiquement de gérer les valeurs propres évaluées à zéro (ou en dessous de zéro). La première méthode est un peu ondulée mais probablement beaucoup plus rapide à mettre en œuvre. Pour quelque chose d'un peu plus stable, vous voudrez peut-être calculer le SVD, puis définir le égal à l'absolu de la plus petite valeur propre (donc vous obtenez non négatif) plus quelque chose de très petit (donc vous obtenez positif). Faites juste attention à ne pas imposer de positivité à une matrice qui est évidemment négative (ou déjà positive). Les deux méthodes modifieront le nombre de conditionnement de votre matrice.λ

En termes statistiques, ce que vous faites en ajoutant que sur la diagonale de votre matrice de covariance, vous ajoutez du bruit à vos mesures. (Parce que la diagonale de la matrice de covariance est la variance de chaque point et en ajoutant quelque chose à ces valeurs, vous dites simplement que "la variance aux points pour lesquels j'ai des lectures est en fait un peu plus grande que ce que je pensais à l'origine".)λ

Un test rapide pour déterminer le caractère positif d'une matrice est l'existence (ou non) de sa décomposition Cholesky.

Également comme note de calcul:

  1. N'utilisez pas de flottants, utilisez des doubles pour stocker votre matrice de covariance.
  2. Utilisez des bibliothèques d'algèbre linéaire numérique en C ++ (comme Eigen ou Armadillo) pour obtenir des inverses de matrices, produits matriciels, etc. C'est plus rapide, plus sûr et plus concis.

EDIT: Étant donné que vous avez une décomposition Cholesky de votre matrice telle que (vous devez le faire pour vérifier que vous avez une matrice Pos.Def.), Vous devriez pouvoir résoudre immédiatement le système . Vous venez de résoudre Ly = b pour y par substitution directe, puis L ^ Tx = y pour x par substitution inverse. (Dans eigen, utilisez simplement la méthode .solve (x) de votre objet Cholesky) Merci à bnaul et Zen d'avoir souligné que je me suis tellement concentré sur l'obtention du be Pos.Def. que j'ai oublié pourquoi nous nous en soucions en premier lieu :)KLLTKx=bK

usεr11852
la source
3
+1. En utilisant Mathematica et en l'appliquant aux données (plutôt qu'à la matrice de covariance affichée, qui peut avoir été présentée avec trop peu de précision), je ne trouve aucune valeur propre négative. C'est comme ça: quand une matrice de covariance est calculée exactement, elle est garantie semi-définie positive, donc toute valeur propre négative doit être attribuée à l'imprécision dans les calculs. Toute procédure inverse généralisée décente doit «reconnaître» ces minuscules valeurs négatives comme des zéros et les traiter en conséquence.
whuber
Merci les gars pour l'effort, comme indiqué, j'ai voté et je vais les essayer et les commenter ou les accepter en conséquence.
Aly
Désolé, je suis un peu confus, comment la résolution du Cholesky donne-t-elle une utilisation de la distance de Mahalanobis?
Aly
Vérifiez le lien dans le message original de bnaul. Mais n'utilisez pas mais Cholesky (c'est ce qu'ils entendent par LDL *). LU
usεr11852
2

Les réponses et les commentaires publiés mettent tous en évidence les dangers de l'inversion de matrices presque singulières. Cependant, pour autant que je sache, personne n'a mentionné que le calcul de la distance de Mahalanobis ne nécessite pas en fait d'inverser la covariance de l'échantillon. Voir cette question StackOverflow pour une description de la façon de procéder à l'aide de la décomposition .LU

Le principe est le même que pour résoudre un système linéaire: lorsque l'on essaie de résoudre pour tel que , il existe des méthodes beaucoup plus efficaces et numériquement stables que de prendre .xAx=bx=A1b

Edit: cela va sans dire, mais cette méthode produit la valeur exacte de la distance, alors que l'ajout de à et l'inversion ne donnent qu'une approximation.λIS

bnaul
la source
1
Tu as raison, @bnaul. Cependant, sans une sorte de régularisation, la LUdécomposition ne fonctionnera pas non plus. J'ajouterai un commentaire à ce sujet dans ma réponse.
Zen
@bnaul: pourquoi faire le LU quand vous le faites au Cholesky pour vérifier le caractère positif? En supposant que vous avez une matrice de covariance valide résolvant pour y par substitution directe, puis pour x par substitution inverse sera plus rapide. Bon point cependant, je me concentre définitivement sur l'obtention d'une définition positive, car j'ai oublié pourquoi je m'en souciais à l'origine! : DK=LLTLy=bLTx=y
usεr11852
0

(Des années plus tard) un tout petit exemple: avec rang déficient en A, valeurs propres de seront de 0 à la précision de la machine - et environ la moitié de ces "zéros" peuvent être :Ar<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
la source