# Making square-root of covariance matrix positive-definite (Matlab)

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.
-
I think maybe you wanted $S = PP^\top$, where $S$ is the covariance, and $P$ is the Cholesky factor. –  shabbychef Jan 19 '11 at 1:39
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! –  Gilead Jan 19 '11 at 1:50

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.


-
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. –  Gilead Jan 19 '11 at 1:51
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. –  shabbychef Jan 19 '11 at 5:57
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. –  Gilead Jan 19 '11 at 20:38

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

-
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. –  Gilead Jan 19 '11 at 20:35
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. –  Nick Jan 19 '11 at 21:02

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!

-
Thanks, that's a technique I could potentially use. It seems like it's been implemented here: infohost.nmt.edu/~borchers/ldlt.html –  Gilead Jan 19 '11 at 20:40