9
$\begingroup$

Motivation: I'm writing a state estimator in MATLAB (the unscented Kalman filter), which calls for the update of the (upper-triangular) square-root of a covariance matrix $S$ at every iteration (that is, for a covariance matrix $P$, it is true that $P=SS^{T}$). In order for me to perform the requisite calculations, I need to do a Rank-1 Cholesky Update and Downdate using MATLAB's cholupdate function.

Problem: Unfortunately, during the course of the iterations, this matrix $S$ can sometimes lose positive definiteness. The Cholesky downdate fails on non-PD matrices.

My question is: are there any simple and reliable ways in MATLAB to make $S$ positive-definite?

(or more generally, is there a good way of making any given covariance $X$ matrix positive-definite?)

Notes:

  • $S$ is full rank
  • I've tried the eigendecomposition approach (which did not work). This basically involved finding $S = VDV^{T}$, setting all negative elements of $V,D = 1 \times 10^{-8}$, and reconstructing a new $S' = V' D' V'^{T}$ where $V',D'$ are matrices with only positive elements.
  • I am aware of the Higham approach (which is implemented in R as nearpd), but it seems to only project to the nearest PSD matrix. I need a PD matrix for the Cholesky update.
Improve this question
$\endgroup$
5
  • $\begingroup$ I think maybe you wanted $S = PP^\top$, where $S$ is the covariance, and $P$ is the Cholesky factor. $\endgroup$ – shabbychef Jan 19 '11 at 1:39
  • $\begingroup$ Actually I want $S$ (the square root, or in this case, the Cholesky factor) to be positive definite. I have clarified the question though; thanks! $\endgroup$ – Gilead Jan 19 '11 at 1:50
  • $\begingroup$ I am having the same problem. I did try sqrtm(x) function but it did work just for few iterations. Did you find the solution? $\endgroup$ – user36937 Jan 6 '14 at 20:26
  • $\begingroup$ I tried various methods, but ended up using the standard (non-rigorous) trick of perturbing the diagonals, i.e. $S + kI$, where $k$ is a constant sufficiently large to make $S$ positive-definite. There are perhaps better approaches, but this one was quick and easy. $\endgroup$ – Gilead Jan 7 '14 at 3:50
  • $\begingroup$ A bit late, but pivoting is a helpful strategy for numerically unstable. $\endgroup$ – probabilityislogic Jul 10 '14 at 13:04
4
$\begingroup$

Here is code I've used in the past (using the SVD approach). I know you said you've tried it already, but it has always worked for me so I thought I'd post it to see if it was helpful. 

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
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for your effort -- unfortunately, it did not work. (I was doing something very similar in my 3-line program: [V,D] = eig(A); D(D <= 1e-10) = 1e-6; Apd = V*A*V';). This approach is similar to the one by Rebonato and Jackel, and it seems to fail for pathological cases like mine. $\endgroup$ – Gilead Jan 19 '11 at 20:35
  • $\begingroup$ Thats too bad. I'd be interested in an example matrix you've found that causes this (and other methods you've tried) to fail if you have time to post one. This is such an aggravating problem to keep running into, I hope you find a solution. $\endgroup$ – Nick Jan 19 '11 at 21:02
2
$\begingroup$

in Matlab:

help cholupdate

I get 

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.
Improve this answer
$\endgroup$
3
  • $\begingroup$ I am using cholupdate but my question is about making R (in this case) positive definite. I have a case where my R is non-pd, and cholupdate(R,X,'-') (a downdate) fails. $\endgroup$ – Gilead Jan 19 '11 at 1:51
  • 1
    $\begingroup$ all online algorithms of this form (update & downdate) suffer from precision issues like this. I had similar issues in 1d resulting in negative estimates of variance. My suggestion would be to keep a circular buffer of the last k vectors observed, and when cholupdate fails, recompute the covariance based on that circular buffer and eat the cost. If you have the memory, and can stand the occasional time hit when this happens, you will not find a better method in terms of accuracy and ease of implementation. $\endgroup$ – shabbychef Jan 19 '11 at 5:57
  • $\begingroup$ Thanks, that's something to think about. Unfortunately, my covariance matrix goes through so many transformations that it is not clear at which point I have to a recomputation from the circular buffer. Nevertheless, if all else fails, I should be able to use the last known PD covariance matrix -- with the hope it doesn't produce a bias in my estimates. $\endgroup$ – Gilead Jan 19 '11 at 20:38
2
$\begingroup$

One alternative way to compute the Cholesky factorisation is by fixing the diagonal elements of S to 1, and then introducing a diagonal matrix D, with positive elements.

This avoids the need to take square roots when doing the computations, which can cause problems when dealing with "small" numbers (i.e. numbers small enough so that the rounding which occurs due to floating point operations matters). The wikipedia page has what this adjusted algorithm looks like.

So instead of $P=SS^T$ you get $P=RDR^T$ with $S=RD^{\frac{1}{2}}$

Hope this helps!

Improve this answer
$\endgroup$
1
1
$\begingroup$

Effectively the Cholesky factorization can fail when your matrix is not "really" positif definite. Two cases appears, or you have a negative eingen value, or your smallest eingen value is positive, but close to zero. The second case must theorically give a solution, but numerically difficult. If have just intuitive add a small constant to the diagonal of my matrix for solving the problem. But this way is not rigourous because it modify slightly the solution. If you must to compute a really hight accuracy solution, try some research on modified Cholesky factorization.

Improve this answer
$\endgroup$
0
$\begingroup$

If you try to estimate with P not positive definite you are asking for problems and chalenging algorithms, you should avoid this situation. If your problem is numeric : P is positive definite but the numerical eigenvalue are too small - try a new scalling for your states If you problem is indeed non positive definite - try different set of state variables. I hope the advis eis not too late Regards, Zeev

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.