1
$\begingroup$

From what I know, for any square real matrix A, a matrix generated with the following should be a positive semidefinite (PSD) matrix:

Q = A @ A.T

I have this matrix A, which is sparse and not symmetric. However, regardless of the properties of A, I think the matrix Q should be PSD.

However, upon using np.linalg.eigvals, I get the following:

np.sort(np.linalg.eigvals(Q))
>>>array([-1.54781185e+01+0.j, -7.27494242e-04+0.j,  2.09363431e-04+0.j, ...,
        3.55351888e+15+0.j,  5.82221014e+17+0.j,  1.78954577e+18+0.j])

I think the complex eigenvalues result from the numerical instability of the operation. Using scipy.linalg.eigh, which takes advantage of the fact that the matrix is symmetric, gives,

np.sort(eigh(Q, eigvals_only=True))
>>>array([-3.10854357e+01, -6.60108485e+00, -7.34059692e-01, ...,
        3.55351888e+15,  5.82221014e+17,  1.78954577e+18])

which again, contains negative eigenvalues.

My goal is to perform Cholesky decomposition on the matrix Q, however, I keep getting this error message saying that the matrix Q is not positive definite, which can be again confirmed with the negative eigenvalues shown above.

Does anyone know why the matrix is not PSD? Thank you.

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ This is the result of floating point imprecision: the value -1.54...e+01 is effectively zero compared to the largest eigenvalue of 1.78...e+18. There are many ways to fix the problem, but which one(s) to recommend depend on how you compute this matrix, what it represents, and what you want to do with the Cholesky decomposition. Maybe you could edit in that information? $\endgroup$
    – whuber ♦
    Commented Mar 22, 2022 at 18:47
  • 2
    $\begingroup$ For useful posts, see stats.stackexchange.com/questions/269323 and stats.stackexchange.com/questions/58527. Do either (or both) solve your problem? $\endgroup$
    – whuber ♦
    Commented Mar 22, 2022 at 18:58
  • $\begingroup$ Putting it simply: the matrix has a ridiculous condition number so it can be effectively factorised. This is most likely a combination of the sparseness and the columns having different scales. First of all, ensure that all the non-zero values are of similar scales. Second why the Cholesky? If we are trying to solve a linear system then there are better solvers we can use. Third, if we absolutely need that Cholesky we can look into different algorithms, maybe $LDL^T$? Finally, if we totally need this, look into CHOLMOD that is for sparse matrices. $\endgroup$
    – usεr11852
    Commented Mar 22, 2022 at 19:00
  • 2
    $\begingroup$ I see a fundamental problem that is easily fixed: don't do the calculation with $Q$! Perform the matrix decomposition with $A$ directly. This will double the precision of the result. cc @usεr11852 $\endgroup$
    – whuber ♦
    Commented Mar 22, 2022 at 19:05
  • 1
    $\begingroup$ @whuber: Agreed. That is my second point. That said, sometimes Cholesky is "very helpful"; for example, we might calculate once and then use it to solve our main system in a way a "standard" QR would fail - case in point: lme4::lmer at the end of the day is (computationally) a CHOLMOD. You get that decomposition done and the rest almost glides through). $\endgroup$
    – usεr11852
    Commented Mar 22, 2022 at 19:13

0

Reset to default

Your Answer

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