Hacer que la raíz cuadrada de la matriz de covarianza sea positiva-definida (Matlab)

Motivación : estoy escribiendo un estimador de estado en MATLAB (el filtro de Kalman sin perfume), que requiere la actualización de la raíz cuadrada (triangular superior) de una matriz de covarianza en cada iteración (es decir, para una matriz de covarianza , es cierto que ). Para poder realizar los cálculos necesarios, necesito hacer una actualización y una actualización de Cholesky de rango 1 usando la función de MATLAB .$S$P = S S $P$$P=SS^{T}$cholupdate

Problema : Desafortunadamente, durante el curso de las iteraciones, esta matriz veces puede perder definición positiva. La actualización de Cholesky falla en matrices que no son PD.$S$

Mi pregunta es : ¿hay alguna manera simple y confiable en MATLAB para hacer que definitivo?$S$

( o más en general, ¿hay una buena manera de hacer que una matriz de covarianza positiva-definida?$X$ )

Notas :

• $S$ es rango completo
• He intentado el enfoque de descomposición propia (que no funcionó). Básicamente, esto implicaba encontrar , establecer todos los elementos negativos de , y reconstruir un nuevo donde son matrices con solo 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'$
• Soy consciente del enfoque de Higham (que se implementa en R como nearpd), pero parece que solo se proyecta a la matriz PSD más cercana. Necesito una matriz PD para la actualización de Cholesky.

Aquí está el código que he usado en el pasado (usando el enfoque SVD). Sé que dijiste que ya lo has probado, pero siempre me ha funcionado, así que pensé en publicarlo para ver si 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

en Matlab:

help cholupdate

yo obtengo

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.

Una forma alternativa de calcular la factorización de Cholesky es fijando los elementos diagonales de S a 1, y luego introduciendo una matriz diagonal D, con elementos positivos.

Esto evita la necesidad de tomar raíces cuadradas al hacer los cálculos, lo que puede causar problemas al tratar con números "pequeños" (es decir, números lo suficientemente pequeños como para que el redondeo que se produce debido a operaciones de punto flotante sea importante). La página de wikipedia tiene el aspecto de este algoritmo ajustado.

Entonces, en lugar de , obtienes con P = R D R T S = R D $P=SS^T$$P=RDR^T$$S=RD^{\frac{1}{2}}$

¡Espero que esto ayude!

Efectivamente, la factorización de Cholesky puede fallar cuando su matriz no es "realmente" positiva. Aparecen dos casos, o tiene un valor de eingen negativo, o su valor de eingen más pequeño es positivo, pero cercano a cero. El segundo caso debe dar una solución teórica, pero numéricamente difícil. Si tiene solo intuitiva, agregue una pequeña constante a la diagonal de mi matriz para resolver el problema. Pero de esta manera no es riguroso porque modifica ligeramente la solución. Si debe calcular una solución de alta precisión, intente investigar sobre la factorización de Cholesky modificada.

Si intenta estimar con P no positivo definido, está solicitando problemas y algoritmos desafiantes, debe evitar esta situación. Si su problema es numérico: P es positivo definido pero el valor propio numérico es demasiado pequeño - intente una nueva escala para sus estados Si su problema es realmente no positivo definido - intente diferentes conjuntos de variables de estado. Espero que el consejo no sea demasiado tarde Saludos, Zeev