Tornando a raiz quadrada da matriz de covariância positiva-definida (Matlab)

Motivação : estou escrevendo um estimador de estado no MATLAB (o filtro Kalman sem cheiro), que exige a atualização da raiz quadrada (triangular superior) de uma matriz de covariância a cada iteração (ou seja, para uma matriz de covariância , é verdade que ). Para que eu possa executar os cálculos necessários, preciso fazer uma Atualização e Downdate Rank-1 de Cholesky usando a função do MATLAB .P P = S S $S$$P$$P=SS^{T}$cholupdate

Problema : Infelizmente, durante o curso das iterações, essa matriz às vezes pode perder uma definição positiva. O downdate de Cholesky falha em matrizes não PD.$S$

Minha pergunta é : existem formas simples e confiáveis ​​no MATLAB de tornar positivo-definido?$S$

( ou, mais geralmente, existe uma boa maneira de tornar positiva qualquer matriz covariância positiva$X$ ) ?

Notas :

• $S$ é classificação completa
• Eu tentei a abordagem eigendecomposition (que não funcionou). Isso basicamente envolveu encontrar , definir todos os elementos negativos de e reconstruir um novo que são matrizes apenas com elementos positivos. V , D = 1 × 10 - 8 S ′ = V ′ D ′ V ′ T V ′ , D $S = VDV^{T}$$V,D = 1 \times 10^{-8}$$S' = V' D' V'^{T}$$V',D'$
• Estou ciente da abordagem de Higham (que é implementada em R as nearpd), mas parece projetar apenas a matriz PSD mais próxima. Eu preciso de uma matriz PD para a atualização Cholesky.

Aqui está o código que eu usei no passado (usando a abordagem SVD). Sei que você disse que já tentou, mas sempre funcionou para mim, então pensei em publicá-lo para ver se era útil.

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

em Matlab:

help cholupdate

eu recebo

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.

Uma maneira alternativa de calcular a fatoração de Cholesky é fixando os elementos diagonais de S em 1 e, em seguida, introduzindo uma matriz diagonal D, com elementos positivos.

Isso evita a necessidade de criar raízes quadradas ao realizar os cálculos, o que pode causar problemas ao lidar com números "pequenos" (ou seja, números pequenos o suficiente para que o arredondamento que ocorre devido a operações de ponto flutuante seja importante). A página da wikipedia tem a aparência desse algoritmo ajustado.

Então, em vez de você obtém com P = R D R T S = R D $P=SS^T$$P=RDR^T$$S=RD^{\frac{1}{2}}$

Espero que isto ajude!

Efetivamente, a fatoração de Cholesky pode falhar quando sua matriz não é "realmente" positivamente definida. Aparecem dois casos, ou você tem um valor de eingen negativo ou seu menor valor de eingen é positivo, mas próximo de zero. O segundo caso deve teoricamente dar uma solução, mas numericamente difícil. Se tiver apenas intuitivo, adicione uma pequena constante à diagonal da minha matriz para resolver o problema. Mas esse caminho não é rigoroso porque modifica um pouco a solução. Se você precisar calcular uma solução de alta precisão, tente algumas pesquisas sobre a fatoração de Cholesky modificada.

Se você tentar estimar com P não positivo positivo, está solicitando problemas e algoritmos de desafio, evite esta situação. Se o seu problema for numérico: P é definitivo positivo, mas o autovalor numérico é muito pequeno - tente uma nova ameaça para os seus estados. Espero que o conselho não seja tarde demais Atenciosamente, Zeev