Rendre positif-défini la racine carrée de la matrice de covariance (Matlab)

9

Motivation : j'écris un estimateur d'état dans MATLAB (le filtre de Kalman non parfumé), qui appelle à la mise à jour de la racine carrée (triangulaire supérieure) d'une matrice de covariance à chaque itération (c'est-à-dire pour une matrice de covariance , il est vrai que ). Pour que je puisse effectuer les calculs requis, je dois faire une mise à jour et une mise à jour Cholesky de rang 1 à l'aide de la fonction MATLAB .SPP=SSTcholupdate

Problème : Malheureusement, au cours des itérations, cette matrice peut parfois perdre en précision positive. La mise à jour de Cholesky échoue sur les matrices non PD.S

Ma question est : existe-t-il des moyens simples et fiables dans MATLAB pour rendre positif-défini?S

( ou plus généralement, existe-t-il un bon moyen de rendre positive une matrice de covariance donnée ?X )


Remarques :

  • S est le rang complet
  • J'ai essayé l'approche eigendecomposition (qui n'a pas fonctionné). Cela impliquait essentiellement de trouver , de définir tous les éléments négatifs de , et de reconstruire un nouveau où sont des matrices avec uniquement des éléments positifs.S=VDVTV,=1×dix-8S=VVTV,
  • Je connais l'approche Higham (qui est implémentée dans R as nearpd), mais elle ne semble se projeter que vers la matrice PSD la plus proche. J'ai besoin d'une matrice PD pour la mise à jour Cholesky.
Gilead
la source
Je pense que tu voulais peut-être S=PP, où S est la covariance, et Pest le facteur Cholesky.
shabbychef
En fait je veux S(la racine carrée, ou dans ce cas, le facteur de Cholesky) pour être définie positive. J'ai cependant clarifié la question; Merci!
Gilead
Je rencontre le même problème. J'ai essayé la fonction sqrtm (x) mais cela n'a fonctionné que pour quelques itérations. Avez-vous trouvé la solution?
J'ai essayé différentes méthodes, mais j'ai fini par utiliser l'astuce standard (non rigoureuse) de perturber les diagonales, à savoir S+kje, où k est une constante suffisamment grande pour faire Sdéfinie positive. Il existe peut-être de meilleures approches, mais celle-ci a été rapide et facile.
Gilead
Un peu tard, mais le pivotement est une stratégie utile pour l'instabilité numérique.
probabilitéislogic

Réponses:

4

Voici le code que j'ai utilisé dans le passé (en utilisant l'approche SVD). Je sais que vous avez dit que vous l'aviez déjà essayé, mais cela a toujours fonctionné pour moi, donc j'ai pensé le poster pour voir si c'était utile.

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
pseudo
la source
Merci pour vos efforts - malheureusement, cela n'a pas fonctionné. (Je faisais quelque chose de très similaire dans mon programme à 3 lignes:) [V,D] = eig(A); D(D <= 1e-10) = 1e-6; Apd = V*A*V';. Cette approche est similaire à celle de Rebonato et Jackel, et elle semble échouer pour des cas pathologiques comme le mien.
Gilead
C'est dommage. Je serais intéressé par un exemple de matrice que vous avez trouvé qui fait échouer cela (et d'autres méthodes que vous avez essayées) si vous avez le temps d'en poster une. C'est un problème tellement aggravant à continuer à courir, j'espère que vous trouverez une solution.
Nick
2

dans Matlab:

help cholupdate

Je reçois

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
la source
J'utilise cholupdatemais ma question est de rendre R(dans ce cas) définitif positif. J'ai un cas où mon Rest non-pd et cholupdate(R,X,'-')(une mise à jour) échoue.
Gilead
1
tous les algorithmes en ligne de ce formulaire (mise à jour et mise à jour) souffrent de problèmes de précision comme celui-ci. J'ai eu des problèmes similaires en 1d, ce qui a entraîné des estimations négatives de la variance. Ma suggestion serait de garder un tampon circulaire des k derniers vecteurs observés et, en cas d' cholupdateéchec, de recalculer la covariance sur la base de ce tampon circulaire et d'en manger le coût. Si vous avez la mémoire et que vous pouvez supporter le temps occasionnel lorsque cela se produit, vous ne trouverez pas de meilleure méthode en termes de précision et de facilité de mise en œuvre.
shabbychef
Merci, c'est quelque chose à penser. Malheureusement, ma matrice de covariance passe par tellement de transformations qu'il n'est pas clair à quel point je dois effectuer un recalcul à partir du tampon circulaire. Néanmoins, si tout le reste échoue, je devrais pouvoir utiliser la dernière matrice de covariance PD connue - dans l'espoir qu'elle ne produise pas de biais dans mes estimations.
Gilead
2

Une autre façon de calculer la factorisation de Cholesky consiste à fixer les éléments diagonaux de S à 1, puis à introduire une matrice diagonale D, avec des éléments positifs.

Cela évite d'avoir à prendre des racines carrées lors des calculs, ce qui peut poser des problèmes lors du traitement de "petits" nombres (c'est-à-dire des nombres suffisamment petits pour que l'arrondi qui se produit en raison des opérations en virgule flottante soit important). La page wikipedia a à quoi ressemble cet algorithme ajusté.

Donc au lieu de P=SST vous obtenez P=RRT avec S=R12

J'espère que cela t'aides!

probabilitéislogique
la source
1
Merci, c'est une technique que je pourrais potentiellement utiliser. Il semble que cela ait été mis en œuvre ici: infohost.nmt.edu/~borchers/ldlt.html
Gilead
1

En effet, la factorisation de Cholesky peut échouer lorsque votre matrice n'est pas "vraiment" positive définie. Deux cas apparaissent, ou vous avez une valeur eingen négative, ou votre plus petite valeur eingen est positive, mais proche de zéro. Le second cas doit théoriquement donner une solution, mais numériquement difficile. Si j'ai juste intuitif, ajoutez une petite constante à la diagonale de ma matrice pour résoudre le problème. Mais cette voie n'est pas rigoureuse car elle modifie légèrement la solution. Si vous devez calculer une solution de très haute précision, essayez des recherches sur la factorisation de Cholesky modifiée.

ctNGUYEN
la source
0

Si vous essayez d'estimer avec P non positif défini que vous demandez des problèmes et des algorithmes de défi, vous devriez éviter cette situation. Si votre problème est numérique: P est défini positif mais les valeurs propres numériques sont trop petites - essayez un nouveau scalling pour vos états Si votre problème est en effet non positif défini - essayez différents ensembles de variables d'état. J'espère que le conseil n'est pas trop tard Cordialement, Zeev

Zeev Berman
la source