Quadratwurzel der Kovarianzmatrix positiv-definit machen (Matlab)

9

Motivation : Ich schreibe einen Zustandsschätzer in MATLAB (dem nicht parfümierten Kalman-Filter), der die Aktualisierung der (oberen Dreiecks-) Quadratwurzel einer Kovarianzmatrix bei jeder Iteration ( dh für eine Kovarianzmatrix P ) fordert ist es wahr, dass P = S S T ). Damit ich die erforderlichen Berechnungen durchführen kann, muss ich mithilfe der MATLAB- Funktion ein Cholesky-Update und -Downdate vom Rang 1 durchführen .SPP=SSTcholupdate

Problem : Leider kann diese Matrix im Verlauf der Iterationen manchmal an positiver Bestimmtheit verlieren. Das Cholesky-Downdate schlägt bei Nicht-PD-Matrizen fehl.S

Meine Frage ist : Gibt es in MATLAB einfache und zuverlässige Möglichkeiten, positiv-definitiv zu machen ?S

( oder allgemeiner, gibt es eine gute Möglichkeit, eine bestimmte Kovarianz- Matrix positiv-definitiv zu machen?X )


Anmerkungen :

  • ist voller RangS
  • Ich habe den Eigendekompositionsansatz ausprobiert (der nicht funktioniert hat). Dies beinhaltete im Wesentlichen das Finden von , das Setzen aller negativen Elemente von V , D = 1 × 10 - 8 und das Rekonstruieren eines neuen S ' = V ' D ' V ' T, wobei V ' , D ' nur Matrizen mit sind positive Elemente.S=VDVTV,D=1×108S=VDVTV,D
  • Ich kenne den Higham-Ansatz (der in R as implementiert ist nearpd), aber er scheint nur auf die nächste PSD-Matrix zu projizieren. Ich benötige eine PD-Matrix für das Cholesky-Update.
Gilead
quelle
S=PPSP
S
Ich habe das gleiche Problem. Ich habe die Funktion sqrtm (x) ausprobiert, aber sie hat nur für wenige Iterationen funktioniert. Haben Sie die Lösung gefunden?
S+kIkS
Ein bisschen spät, aber das Schwenken ist eine hilfreiche Strategie für numerisch instabile.
Wahrscheinlichkeitslogik

Antworten:

4

Hier ist Code, den ich in der Vergangenheit verwendet habe (unter Verwendung des SVD-Ansatzes). Ich weiß, dass du gesagt hast, du hast es bereits versucht, aber es hat immer bei mir funktioniert, also dachte ich, ich würde es posten, um zu sehen, ob es hilfreich ist.

function [sigma] = validateCovMatrix(sig)

% [sigma] = validateCovMatrix(sig)
%
% -- INPUT --
% sig:      sample covariance matrix
%
% -- OUTPUT --
% sigma:    positive-definite covariance matrix
%

EPS = 10^-6;
ZERO = 10^-10;

sigma = sig;
[r err] = cholcov(sigma, 0);

if (err ~= 0)
    % the covariance matrix is not positive definite!
    [v d] = eig(sigma);

    % set any of the eigenvalues that are <= 0 to some small positive value
    for n = 1:size(d,1)
        if (d(n, n) <= ZERO)
            d(n, n) = EPS;
        end
    end
    % recompose the covariance matrix, now it should be positive definite.
    sigma = v*d*v';

    [r err] = cholcov(sigma, 0);
    if (err ~= 0)
        disp('ERROR!');
    end
end
Nick
quelle
Vielen Dank für Ihre Mühe - leider hat es nicht funktioniert. (Ich habe in meinem 3-Zeilen-Programm etwas sehr Ähnliches gemacht :) [V,D] = eig(A); D(D <= 1e-10) = 1e-6; Apd = V*A*V';. Dieser Ansatz ähnelt dem von Rebonato und Jackel und scheint in pathologischen Fällen wie meinem fehlzuschlagen.
Gilead
Das ist sehr schade. Ich würde mich für eine Beispielmatrix interessieren, die dazu führt, dass diese (und andere Methoden, die Sie ausprobiert haben) fehlschlägt, wenn Sie Zeit haben, eine zu veröffentlichen. Ich hoffe, Sie finden eine Lösung für dieses Problem.
Nick
2

in Matlab:

help cholupdate

Ich bekomme

CHOLUPDATE Rank 1 update to Cholesky factorization.
    If R = CHOL(A) is the original Cholesky factorization of A, then
    R1 = CHOLUPDATE(R,X) returns the upper triangular Cholesky factor of A + X*X',
    where X is a column vector of appropriate length.  CHOLUPDATE uses only the
    diagonal and upper triangle of R.  The lower triangle of R is ignored.

    R1 = CHOLUPDATE(R,X,'+') is the same as R1 = CHOLUPDATE(R,X).

    R1 = CHOLUPDATE(R,X,'-') returns the Cholesky factor of A - X*X'.  An error
    message reports when R is not a valid Cholesky factor or when the downdated
    matrix is not positive definite and so does not have a Cholesky factorization.

    [R1,p] = CHOLUPDATE(R,X,'-') will not return an error message.  If p is 0
    then R1 is the Cholesky factor of A - X*X'.  If p is greater than 0, then
    R1 is the Cholesky factor of the original A.  If p is 1 then CHOLUPDATE failed
    because the downdated matrix is not positive definite.  If p is 2, CHOLUPDATE
    failed because the upper triangle of R was not a valid Cholesky factor.

    CHOLUPDATE works only for full matrices.

    See also chol.
shabbychef
quelle
Ich benutze, cholupdateaber meine Frage ist, R(in diesem Fall) positiv definitiv zu machen. Ich habe einen Fall, in dem mein Rnicht pd ist und cholupdate(R,X,'-')(ein Downdate) fehlschlägt.
Gilead
1
Alle Online-Algorithmen dieser Form (Update & Downdate) leiden unter solchen Präzisionsproblemen. Ich hatte ähnliche Probleme in 1d, was zu negativen Varianzschätzungen führte. Mein Vorschlag wäre, einen zirkulären Puffer der letzten k beobachteten Vektoren cholupdatebeizubehalten und, wenn dies fehlschlägt, die Kovarianz basierend auf diesem zirkulären Puffer neu zu berechnen und die Kosten zu essen. Wenn Sie über den Arbeitsspeicher verfügen und in diesem Fall gelegentlich Zeit sparen können, werden Sie keine bessere Methode in Bezug auf Genauigkeit und einfache Implementierung finden.
Shabbychef
Danke, das ist etwas zum Nachdenken. Leider durchläuft meine Kovarianzmatrix so viele Transformationen, dass nicht klar ist, an welchem ​​Punkt ich eine Neuberechnung aus dem Kreispuffer durchführen muss. Wenn alles andere fehlschlägt, sollte ich dennoch in der Lage sein, die letzte bekannte PD-Kovarianzmatrix zu verwenden - mit der Hoffnung, dass meine Schätzungen keine Verzerrung hervorrufen.
Gilead
2

Eine alternative Methode zur Berechnung der Cholesky-Faktorisierung besteht darin, die diagonalen Elemente von S auf 1 zu fixieren und dann eine diagonale Matrix D mit positiven Elementen einzuführen.

Dies vermeidet die Notwendigkeit, bei den Berechnungen Quadratwurzeln zu ziehen, was zu Problemen beim Umgang mit "kleinen" Zahlen führen kann (dh Zahlen, die klein genug sind, damit die Rundung aufgrund von Gleitkommaoperationen von Bedeutung ist). Auf der Wikipedia-Seite sehen Sie, wie dieser angepasste Algorithmus aussieht.

P=SSTP=RDRTS=RD12

Hoffe das hilft!

Wahrscheinlichkeitslogik
quelle
1
Danke, das ist eine Technik, die ich möglicherweise anwenden könnte. Es scheint, als wäre es hier implementiert worden: infohost.nmt.edu/~borchers/ldlt.html
Gilead
1

Tatsächlich kann die Cholesky-Faktorisierung fehlschlagen, wenn Ihre Matrix nicht "wirklich" positiv ist. Es treten zwei Fälle auf, oder Sie haben einen negativen Eingen-Wert oder Ihr kleinster Eingen-Wert ist positiv, aber nahe Null. Der zweite Fall muss theoretisch eine Lösung geben, ist aber numerisch schwierig. Wenn Sie nur intuitiv sind, fügen Sie der Diagonale meiner Matrix eine kleine Konstante hinzu, um das Problem zu lösen. Dieser Weg ist jedoch nicht streng, da er die Lösung geringfügig verändert. Wenn Sie eine wirklich hochgenaue Lösung berechnen müssen, versuchen Sie es mit Untersuchungen zur modifizierten Cholesky-Faktorisierung.

ctNGUYEN
quelle
0

Wenn Sie versuchen, mit P nicht positiv definitiv zu schätzen, fragen Sie nach Problemen und Chalenging-Algorithmen, sollten Sie diese Situation vermeiden. Wenn Ihr Problem numerisch ist: P ist positiv definit, aber der numerische Eigenwert ist zu klein - versuchen Sie eine neue Abwertung für Ihre Zustände. Wenn Ihr Problem tatsächlich nicht positiv definit ist - versuchen Sie es mit einem anderen Satz von Zustandsvariablen. Ich hoffe, der Rat kommt nicht zu spät. Grüße, Zeev

Zeev Berman
quelle