2
$\begingroup$

I'm following Peter Jackel's book "Monte Carlo Methods in Finance", where an algorithm for transforming malformed correlation matrices into acceptable correlation matrices (positive semi-definite) to perform a Cholesky decomposition later on.

The basic approach is transforming the negative eigenvalues into 0 and then transforming the matrix to adjust to these new values. In http://www.mathworks.com/matlabcentral/fileexchange/26475-salvaging-a-linear-correlation-matrix/content/SpectralDecomp.m there's an implementation for matlab of this method.

The problem in this case is in floating point arithmetics. Suppose I have this matrix:

1  1  1
1  1  1
1  1  1

The eigenvalues for this matrix are:

-9.11271925974066E-17
1.11022302462516E-16
3

This matrix is not positive semi-definite, because of the first eigenvalue. However, when I try to adjust it to a new matrix that is positive semi-definite, I end up with the very same matrix as above! And the reason is quite clear: the first 2 eigenvalues are too close to zero already - if they're substituted by zero, nothing really changes at all. In fact, by inputting this matrix in a system with less precision (ie http://www.mathstools.com/section/main/system_equations_solver#.UYrN5rV-jK0), they do appear to be zero. In the final matrix multiplication, I end up with a sum that's not mathematically correct, but a consequence of floating point arithmetic:

3.57e-33 + 1 = 1

I'm using VBA for these calculations, with the main operations made mainly by the functions in the excellent Alglib library, but Excel only stores 15 precision digits for floating point numbers. Maybe I could try the hypershpere decomposition approach instead of the spectral decomposition, but I think I'll probably bump into the same issue. Do R, Matlab or any other mathematical software deal with this issue?

To reproduce the problem, try to multiply the following matrix by its own transpose:

0    7.45058059692383E-09    1
0   -7.45058059692383E-09    1   
0            0               1

The result shouldn't be a matrix of 1's, but that's what's returned by the calculations.

Is there any way to circumvent this limitation? Or any other method to obtain a positive semi-definite matrix out of an arbitrary given correlation matrix? I need a Cholesky decomposition of this matrix to generate random samples for a copula, perhaps there's another method for this that does not involve the Cholesky decomposition and the consequent constraint of the matrix being positive semi-definite?

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ use the svd decomposition instead of the spectral decomposition. The singular values of a PSD matrix are guaranteed to always be positive. For your matrix (which is PSD) these are: 3.000000e+00 1.338312e-16 8.967126e-49. $\endgroup$
    – user603
    Commented May 9, 2013 at 6:30
  • 1
    $\begingroup$ I don't get it...this process tries to convert a non-semi positive matrix into one, to apply a Cholesky decomposition later on for use in a copula. How can the SVD decomposition fit into this? $\endgroup$
    – mr_mouse
    Commented May 22, 2013 at 15:39
  • 4
    $\begingroup$ Because this exposition wanders and never really seems to get to a single clear question, it is hard to see what you are asking. Most of it reads like an exploration of the (well-known) consequences of floating-point imprecision. And what kind of simulation are you running where imprecision on the order of $10^{-16}$ is going to make any difference? Would it be fair to sum up your question by characterizing it as "How can I obtain a Cholesky-like decomposition for a degenerate correlation matrix using (only) double-precision arithmetic"? $\endgroup$
    – whuber ♦
    Commented May 22, 2013 at 16:35

1 Answer 1

Reset to default
2
$\begingroup$

I'm also trying to solve a similar problem: transform a matrix that is not positive definite into a positive definite one for then applying a Cholesky Decomposition.

The best way to solve the problem would be change the equations to use a Cholesky decomposition $LDL^T$ that doesn't take the square root of the matrix (see 1 and 2). This way, you might be able to work directly with the non-definite matrix.

Anyway, here's a simpler method to "convert" the matrix to a positive definite one and it seems to be better at reconstructing the original matrix than the one you pointed out.

Results of a test with a matrix that is already positive definite, so the matrix should be preserved:

Before: 

   1.00114  0.0015438  0.0055613  0.0083755  0.0114398  0.0271191
 0.0015438    1.00209 0.00751505  0.0113179  0.0154588  0.0366463
 0.0055613 0.00751505    1.02707  0.0407709  0.0556878   0.132013
 0.0083755  0.0113179  0.0407709     1.0614  0.0838676   0.198815
 0.0114398  0.0154588  0.0556878  0.0838676    1.11455   0.271555
 0.0271191  0.0366463   0.132013   0.198815   0.271555    1.64375

After, using the method you mentioned: 

         1 0.00154131 0.00548439 0.00812499  0.0108298  0.0211402
0.00154131          1 0.00740762  0.0109742  0.0146276  0.0285536
0.00548439 0.00740762          1   0.039049  0.0520486   0.101601
0.00812499  0.0109742   0.039049          1  0.0771088   0.150519
 0.0108298  0.0146276  0.0520486  0.0771088          1   0.200628
 0.0211402  0.0285536   0.101601   0.150519   0.200628          1

After using the simpler method: 

   1.00114  0.0015438  0.0055613  0.0083755  0.0114398  0.0271191
 0.0015438    1.00209 0.00751505  0.0113179  0.0154588  0.0366463
 0.0055613 0.00751505    1.02707  0.0407709  0.0556878   0.132013
 0.0083755  0.0113179  0.0407709     1.0614  0.0838676   0.198815
 0.0114398  0.0154588  0.0556878  0.0838676    1.11455   0.271555
 0.0271191  0.0366463   0.132013   0.198815   0.271555    1.64375

Test done using C++ Eigen library with double precision.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Because this answer lacks an explanation, here it is: the function used here merely zeros out any small negative eigenvalues. Sometimes that works, but there are other ways (such as applying a tiny amount of regularization). I have found that automatically but selectively changing results due to floating point rounding error can introduce additional systematic error that may propagate through subsequent calculations, frequently (and ironically) creating worse results. $\endgroup$
    – whuber ♦
    Commented Mar 22, 2022 at 19:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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