I estimated the sample covariance matrix C of a sample and get a symmetric matrix. With C, I would like to create n-variate normal distributed r.n. but therefore I need the Cholesky decomposition of C. What should I do if C is not positive definite?
Solution Method A:
- If C is not symmetric, then symmetrize it. D <– 0.5(C+CT)
- Add a multiple of the Identity matrix to the symmetrized C sufficient to make it positive definite with whatever margin, m, is desired, i.e., such that smallest eigenvalue of new matrix has minimum eigenvalue = m. Specifically, D <– D+(m−min(eigenvalue(D)))I, where I is the identity matrix. D contains the desired positive definite covariance matrix.
In MATLAB, the code would be
D = 0.5 * (C + C'); D = D + (m - min(eig(CD)) * eye(size(D));
Solution Method B:
Formulate and solve a Convex SDP (Semidefinite Program) to find the nearest matrix D to C according to the frobenius norm of their difference, such that D is positive definite, having specified minimum eigenvalue m.
Using CVX under MATLAB, the code would be:
n = size(C,1); cvx_begin variable D(n,n) minimize(norm(D-C,'fro')) D -m *eye(n) == semidefinite(n) cvx_end
Comparison of Solution Methods: Apart from symmetrizing the initial matrix, solution method A adjusts (increases) only the diagonal elements by some common amount, and leaves the off-diagonal elements unchanged. Solution method B finds the nearest (to the original matrix) positive definite matrix having the specified minimum eigenvalue, in the sense of minimum frobenius norm of the difference of the positive definite matrix D and the original matrix C, which is based on the sums of squared differences of all elements of D – C, to include the off-diagonal elements. So by adjusting off-diagonal elements, it may reduce the amount by which diagonal elements need to be increased, and diagoanl elements are not necessarily all increased by the same amount.