Numerische Instabilität der Berechnung der inversen Kovarianzmatrix

8

Ich habe 65 Stichproben von 21-dimensionalen Daten ( hier eingefügt ) und konstruiere daraus die Kovarianzmatrix. Bei der Berechnung in C ++ wird hier die Kovarianzmatrix eingefügt . Und wenn ich in Matlab aus den Daten berechnet werde (wie unten gezeigt), wird die Kovarianzmatrix hier eingefügt

Matlab-Code zur Berechnung von cov aus Daten:

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

Wie Sie sehen können, sind die Unterschiede in den Kovarianzmatrizen winzig (~ e-07), was wahrscheinlich auf numerische Probleme im Compiler mit Gleitkomma-Arithmetik zurückzuführen ist.

Wenn ich jedoch die pseudo-inverse Kovarianzmatrix aus der von matlab und der von meinem C ++ - Code erzeugten Kovarianzmatrix berechne, erhalte ich sehr unterschiedliche Ergebnisse. Ich berechne sie auf die gleiche Weise, dh:

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

Der Unterschied ist so groß, dass bei der Berechnung des Mahalanobis-Abstands von einer Probe ( hier eingefügt ) zur Verteilung der 65 Proben durch:

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

Mit den verschiedenen inversen Kovarianzmatrizen ( 1 ) erhalte ich sehr unterschiedliche Ergebnisse, dh:

 (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

Ist es normal, dass die kleinen (e-7) Unterschiede in der Kovarianzmatrix einen solchen Effekt auf die Berechnung der pseudoinversen Matrix haben? Und wenn ja, was kann ich tun, um diesen Effekt abzuschwächen?

Wenn dies nicht gelingt, kann ich andere Entfernungsmetriken verwenden, die nicht die inverse Kovarianz beinhalten? Ich verwende den Mahalanobis-Abstand, da wir wissen, dass er für n Proben einer Beta-Verteilung folgt, die ich für Hypothesentests verwende

Vielen Dank im Voraus

BEARBEITEN: Hinzufügen von C ++ - Code zur Berechnung der Kovarianzmatrix unten: Die vector<vector<double> >repräsentiert die Sammlung von Zeilen aus der eingefügten Datei.

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
quelle
Matrizen invertieren ..... Das ist eine gefährliche Sache! Normalerweise ist es vorzuziehen, Alternativen dazu zu finden (z. B. Pseudoinverse)
Ander Biguri
1
@Aly: Die Matrizen, die Sie invertieren möchten, sind keine "gültigen" Kovarianzmatrizen, da sie nicht eindeutig positiv sind. numerisch haben sie sogar einige Eigenwerte, die negativ sind (aber nahe bei Null liegen). Ich würde wahrscheinlich nur eine sehr kleine Konstante entlang der Diagonale hinzufügen; es ist wirklich eine Form der Tichonow-Korrektur ( ). Verwenden Sie auch keine Floats, sondern doppelte, um Ihre Kovarianzmatrix zu speichern. (Und abgesehen davon, dass Sie bereits OpenCV verwenden, können Sie auch Eigen oder Armadillo verwenden.)Χ+λI
usεr11852
1
@Aly: Wikipedia, wirklich. (Es ist das Lemma: Tichonow-Regularisierung). Die Methode, die mit der SVD erwähnt wurde, liefert eine nicht negative definitive Matrix, wenn Sie kleine Eigenwerte auf Null setzen. Sie müssen immer noch eine kleine Konstante zu all Ihren Eigenwerten hinzufügen, um sie positiv zu definieren. Praktisch beide Methoden machen das gleiche. Ich habe nur versucht, die SVD nicht zu verwenden, sondern die Eigenwerte der Stichproben direkt zu beeinflussen, indem ich allen hinzufüge . Ich habe keine Referenzen gefunden, beide Methoden sind meiner Meinung nach sehr intuitiv. λ
usεr11852
1
@ user11852 Bitte können Sie Ihre Kommentare eine Antwort geben, ich experimentiere noch, aber wenn vielversprechend zu akzeptieren. Auch wenn andere ihre Vorschläge beantworten, werde ich abstimmen, da sie für mein Verständnis des Problems sehr hilfreich / nützlich waren
Aly
1
Ich habe in Ihrem anderen Thread kommentiert , dass Variablen , die wie Ihr Datensatz 1 ergeben , die Instabilität fördern und eine redundante Variable enthalten. Bitte versuchen Sie, eine Spalte zu löschen. Sie brauchen nicht einmal den Pinv: Die Kovarianzmatrix ist nicht mehr singulär.
Cam.Davidson.Pilon

Antworten:

7

Die Matrizen, die Sie invertieren möchten, sind keine "gültigen" Kovarianzmatrizen, da sie nicht eindeutig positiv sind. numerisch haben sie sogar einige Eigenwerte, die negativ sind (aber nahe bei Null liegen). Dies ist höchstwahrscheinlich auf Maschinennullen zurückzuführen, beispielsweise ist der letzte Eigenwert Ihrer Matrix "matlab_covariance" -0.000000016313723. Um positiv positiv zu korrigieren, können Sie zwei Dinge tun:

  1. Fügen Sie einfach eine sehr kleine Konstante entlang der Diagonale hinzu. eine Form der Tichonow-Korrektur wirklich ( mit ).Χ+λIλ0
  2. (Basierend auf dem, was whuber vorgeschlagen hat) Verwenden Sie SVD, setzen Sie die "problematischen" Eigenwerte auf einen festen kleinen Wert (nicht Null), rekonstruieren Sie Ihre Kovarianzmatrix und invertieren Sie diese dann. Wenn Sie einige dieser Eigenwerte auf Null setzen, erhalten Sie eindeutig eine nicht negative (oder halbpositive) Matrix, die immer noch nicht invertierbar ist.

Eine nicht negative Matrix hat keine Inverse, aber eine Pseudo-Inverse (alle Matrizen mit reellen oder komplexen Einträgen haben eine Pseudo-Inverse), dennoch ist die Moore-Penrose-Pseudo-Inverse rechenintensiver als eine echte Inverse und wenn Das Inverse existiert, es ist gleich dem Pseudo-Inversen. Also mach einfach das Gegenteil :)

Beide Methoden versuchen praktisch, mit den Eigenwerten umzugehen, die sich zu Null (oder unter Null) ergeben. Die erste Methode ist etwas handgewellt, aber wahrscheinlich viel schneller zu implementieren. Für etwas etwas stabileres möchten Sie vielleicht die SVD berechnen und dann das gleich dem Absoluten des kleinsten Eigenwerts (damit Sie nicht negativ werden) plus etwas sehr Kleines (damit Sie positiv werden) setzen. Achten Sie nur darauf, dass eine Matrix, die offensichtlich negativ (oder bereits positiv) ist, nicht positiv beeinflusst wird. Beide Methoden ändern die Konditionierungsnummer Ihrer Matrix.λ

In statistischer Hinsicht fügen Sie Ihren Messungen Rauschen hinzu, indem Sie dieses über die Diagonale Ihrer Kovarianzmatrix addieren. (Da die Diagonale der Kovarianzmatrix die Varianz jedes Punktes ist und Sie diesen Werten etwas hinzufügen, sagen Sie einfach "die Varianz an den Punkten, für die ich Messwerte habe, ist tatsächlich etwas größer als ursprünglich angenommen".)λ

Ein schneller Test für die positive Bestimmtheit einer Matrix ist das Vorhandensein (oder Nichtvorhandensein) der Cholesky-Zersetzung derselben.

Auch als rechnerische Anmerkung:

  1. Verwenden Sie keine Floats, sondern Double, um Ihre Kovarianzmatrix zu speichern.
  2. Verwenden Sie numerische lineare Algebra-Bibliotheken in C ++ (wie Eigen oder Armadillo), um Inversen von Matrizen, Matrixprodukten usw. zu erhalten. Es ist schneller, sicherer und präziser.

EDIT: Wenn Sie eine Cholesky-Zerlegung Ihrer Matrix so dass (Sie müssen dies tun, um zu überprüfen, ob Sie eine Pos.Def. Matrix haben), sollten Sie in der Lage sein, das System sofort zu lösen . Sie lösen einfach Ly = b für y durch Vorwärtssubstitution und dann L ^ Tx = y für x durch Rückwärtssubstitution. (Verwenden Sie im Eigen einfach die .solve (x) -Methode Ihres Cholesky-Objekts.) Vielen Dank an bnaul und Zen für den Hinweis, dass ich mich so sehr darauf konzentriert habe, das be Pos.Def zu erhalten. dass ich vergessen habe, warum uns das überhaupt interessiert hat :)KLLTKx=bK

usεr11852
quelle
3
+1. Wenn ich Mathematica verwende und es auf die Daten anwende (anstelle der veröffentlichten Kovarianzmatrix, die möglicherweise mit zu geringer Genauigkeit dargestellt wurde), finde ich keine negativen Eigenwerte. Das ist so, wie es sein sollte: Wenn eine Kovarianzmatrix genau berechnet wird, ist sie garantiert positiv semidefinit, daher müssen alle negativen Eigenwerte der Ungenauigkeit in den Berechnungen zugeordnet werden. Jedes anständige verallgemeinerte inverse Verfahren sollte diese winzigen negativen Werte als Nullen "erkennen" und sie entsprechend behandeln.
whuber
Vielen Dank für die Mühe, wie gesagt, ich habe abgestimmt und werde diese ausprobieren und entweder kommentieren oder entsprechend akzeptieren.
Aly
Entschuldigung, ich bin ein bisschen verwirrt. Wie nutzt das Lösen des Cholesky die Mahalanobis-Distanz?
Aly
Überprüfen Sie den Link im Originalbeitrag von bnaul. Verwenden Sie jedoch nicht sondern Cholesky (das ist es, was sie mit LDL * meinen). LU
usεr11852
2

Die veröffentlichten Antworten und Kommentare machen alle gute Hinweise auf die Gefahren der Invertierung nahezu singulärer Matrizen. Soweit ich das beurteilen kann, hat jedoch niemand erwähnt, dass für die Berechnung der Mahalanobis-Entfernung keine Umkehrung der Stichproben-Kovarianz erforderlich ist. In dieser StackOverflow-Frage finden Sie eine Beschreibung der Verwendung der Zerlegung.LU

Das Prinzip ist dasselbe wie das Lösen eines linearen Systems: Wenn versucht wird, nach zu lösen, so dass , gibt es viel effizientere und numerisch stabilere Methoden als .A x = b x = A - 1 bxAx=bx=A1b

Bearbeiten: Es ist wahrscheinlich selbstverständlich, aber diese Methode liefert den genauen Abstandswert, während das Addieren von zu und das Invertieren nur eine Annäherung ergibt.S.λIS

bnaul
quelle
1
Du hast recht, @bnaul. Ohne irgendeine Art von Regularisierung LUfunktioniert die Zerlegung jedoch auch nicht. Ich werde in meiner Antwort einen Kommentar dazu hinzufügen.
Zen
@bnaul: Warum die LU, wenn Sie mit dem Cholesky tun, um die positive Bestimmtheit zu überprüfen? Angenommen, Sie haben eine gültige Kovarianzmatrix , die für y durch Vorwärtssubstitution löst , und dann ist für x durch Rückwärtssubstitution schneller. Guter Punkt, ich konzentriere mich definitiv darauf, eine positive Bestimmtheit zu bekommen, die ich vergessen habe, warum ich mich ursprünglich darum gekümmert habe! : D L y = b L T x = yK=LLTLy=bLTx=y
usεr11852
0

(Jahre später) ein winziges Beispiel: Mit Rangmangel werden Eigenwerte von innerhalb der Maschinengenauigkeit 0 sein - und ungefähr die Hälfte dieser "Nullen" kann :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
quelle