Code covered by the BSD License  

Highlights from

5.0 | 4 ratings Rate this file 109 Downloads (last 30 days) File Size: 23.2 KB File ID: #42885 Version: 1.1



John D'Errico (view profile)


30 Jul 2013 (Updated )

Finding the nearest positive definite matrix

| Watch this File

File Information

This tool saves your covariance matrices, turning them into something that really does have the property you will need. That is, when you are trying to use a covariance matrix in a tool like mvnrnd, it makes no sense if your matrix is not positive definite. So mvnrnd will fail in that case.

But sometimes, it appears that users end up with matrices that are NOT symmetric and positive definite (commonly abbreviated as SPD) and they still wish to use them to generate random numbers, often in a tool like mvnrnd. A solution is to find the NEAREST matrix (minimizing the Frobenius norm of the difference) that has the desired property of being SPD.

I see the question come up every once in a while, so I looked in the file exchange to see what is up there. All I found was nearest_posdef. While this usually almost works, it could be better. It actually failed completely on most of my test cases, and it was not as fast as I would like, using an optimization. In fact, in the comments to nearest_posdef, a logical alternative was posed. That alternative too has its failures, so I wrote nearestSPD, partly based on what I found in the work of Nick Higham.

nearestSPD works on any matrix, and it is reasonably fast. As a test, randn generates a matrix that is not symmetric nor is it at all positive definite in general.
U = randn(100);

nearestSPD will be able to convert U into something that is indeed SPD, and for a 100 by 100 matrix, do it quickly enough.

tic,Uj = nearestSPD(U);toc
Elapsed time is 0.008964 seconds.

The ultimate test of course, is to use chol. If chol returns a second argument that is zero, then MATLAB (and mvnrnd) will be happy!
[R,p] = chol(Uj);
p =


Nearest Positive Semi Definite Covariance Matrix inspired this file.

Required Products MATLAB
MATLAB release MATLAB 8.1 (R2013a)
Other requirements Nothing fancy in here, so much older MATLAB releases will works easily.
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (10)
05 Sep 2015 John D'Errico

John D'Errico (view profile)

Please send me an example case that has this problem.

Comment only
05 Sep 2015 Yi He

Yi He (view profile)

Thanks John. But when I run some 125*125 covariance matrices, the progress stands at ' mineig = min(eig(Ahat));' for pretty long time (actually almost over 10 hours). What could I do with this issue?
Thanks again.

Comment only
24 Nov 2014 Aeimit Lakdawala  
09 Nov 2014 Meng

Meng (view profile)

Thank you. Very useful!

08 Nov 2014 John D'Errico

John D'Errico (view profile)

I actually nudge the matrix just a bit at the very end if necessary. The code should ensure that chol applied to the result will always yield a valid factorization, and that is essentially the test in MATLAB to be truly SPD. So while the Higham algorithm will ensure positive semi-definite, if chol should always work, then it will indeed be positive definite.

Comment only
07 Nov 2014 Ioannis P

Could you please explain if this code is giving a positive definite or a semi-positive definite matrix? You have written the following:

"From Higham: "The nearest symmetric positive semidefinite matrix in the
Frobenius norm to an arbitrary real matrix A is shown to be (B + H)/2,
where H is the symmetric polar factor of B=(A + A')/2."
arguments: (input)
A - square matrix, which will be converted to the nearest Symmetric
Positive Definite Matrix."

Could you please clarify this? Thanks!

Comment only
29 Oct 2014 Ahmed  
13 Nov 2013 Eric

Eric (view profile)

Kudos to you, John, mostly for calling attention to Higham's paper. Trying to use the other files you mentioned was driving me crazy, because of their high probability of failure.

31 Jul 2013 John D'Errico

John D'Errico (view profile)

Sorry about that. Frobenius norm is minimized.

Comment only
31 Jul 2013 Petr Pošík

Hi John. I miss in the description how the "nearness" of the 2 matrices, U and Uj, is measured. Could you comment on that? Thanks, Petr

Comment only
31 Jul 2013 1.1

Documentation fix

Contact us