The results of linear algebra functions might not be consistent

Algorithms
1

I’m trying to implement a model employing some linear algebra functions. I find the results is confusing. I generated a matrix (Phat) and use QR decomposition, SVD, and Cholesky decomposition on it. Here is the code and results.

print("1-1 P-hat: ", Phat)
print("    QR - Q: ", qr_Q(Phat))
print("    QR - R: ", qr_R(Phat))
print("    determinant: ", determinant(Phat))
print("    eigonvalue: ", eigenvalues_sym(Phat))
print("    decompose: ", cholesky_decompose(Phat))

1-1 P-hat: [[1e+06,0,0,0,0,0,0,0,0],[0,416667,0,0,-83333.3,-83.3333,-83.3333,-83333.3,-83.3333],[0,0,1e+06,0,0,0,0,0,0],[0,0,0,1e+06,0,0,0,0,0],[0,-83333.3,0,0,416667,-83.3333,-83.3333,-83333.3,-83.3333],[0,-83.3333,0,0,-83.3333,0.416667,-0.0833333,-83.3333,-0.0833333],[0,-83.3333,0,0,-83.3333,-0.0833333,0.416667,-83.3333,-0.0833333],[0,-83333.3,0,0,-83333.3,-83.3333,-83.3333,416667,-83.3333],[0,-83.3333,0,0,-83.3333,-0.0833333,-0.0833333,-83.3333,0.416667]]
QR - Q: [[1,-0,0,0,-0,-0,-0,-0,-0],[0,0.96225,0,0,0.136083,-0.2357,-0.0007071,0.000408247,-0.00057735],[0,-0,1,0,-0,-0,-0,-0,-0],[0,-0,0,1,-0,-0,-0,-0,-0],[0,-0.19245,0,0,0.952579,-0.2357,-0.0007071,0.000408247,-0.00057735],[0,-0.00019245,0,0,-0.000272165,0.00329981,-0.7071,0.408247,-0.57735],[0,-0.00019245,0,0,-0.000272165,-0.000942803,0.707107,0.408247,-0.57735],[0,-0.19245,0,0,-0.272165,-0.942803,-0.0028284,0.00163299,-0.00057735],[0,-0.00019245,0,0,-0.000272165,-0.000942803,-2.8284e-06,-0.816496,-0.57735]]
QR - R: [[1e+06,0,0,0,0,0,0,0,0],[-0,433013,-0,-0,-144338,-48.1126,-48.1126,-144338,-48.1126],[0,0,1e+06,0,0,0,0,0,0],[0,0,0,1e+06,0,0,0,0,0],[-0,-0,-0,-0,408248,-68.0414,-68.0414,-204124,-68.0414],[-0,-0,-0,-0,-0,117.852,117.85,-353551,117.85],[-0,-0,-0,-0,-0,-0,0.707104,-1060.65,0.353549],[-0,-0,-0,-0,-0,-0,-0,612.372,-0.612372],[-0,-0,-0,-0,-0,-0,-0,-0,1.92624e-14]]
determinant: -5.20417e+18
eigonvalue: [5.22073e-12,0.5,0.5,250000,500000,500000,1e+06,1e+06,1e+06]
decompose:
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: cholesky_decompose: Matrix m is not positive definite

The results are hard to read, so I put them in Excel and format them.

屏幕快照 2020-05-13 下午8.54.54
屏幕快照 2020-05-13 下午8.54.541322×1228 459 KB

The Cholesky decomposition causes error and says the matrix is not positive definite. The determinant is negative, but the eigenvalues are all positive. The multiplication of the diagonal elements of the matrix R of QR (1.74E+20) is different from the multiplication of the eigenvalues (8.16E+22).

The results are confusing for me. Am I using anything wrong? I use PyStan 2.19.1.1.

2 Likes
2

This happens all the time when a matrix is closed to singular. Basically rounding error in the cholesky algorithm leads it to fail. One solution is to add a small constant fo the diagonal (GP people call this a jitter). Or you can use a more expensive but more stable algorithm (Like a QR decomposition or a full eigendecomposition)

3 Likes
3

Thanks for answering me. I got another case.

屏幕快照 2020-05-13 下午9.53.16
屏幕快照 2020-05-13 下午9.53.161048×1032 357 KB

For the matrix P[1], the Cholesky decomposition also causes error and it is consistent to the results of the determinant and eigenvalues. But the diagonal elements multiplication of QR-R is opposite to the multiplication of eigenvalues. They have the same quantity but one is positive and the other is negative.

Does it mean that I should only use the qr-R function？

1-1 P[1]: [[1e+06,0,0,0,0,0,0],[0,416667,0,0,-83345.4,-12.893,-83333.6],[0,0,1e+06,0,0,0,0],[0,0,0,1e+06,0,0,0],[0,-83345.4,0,0,-2.00465e+06,-12.9034,-139838],[0,-12.893,0,0,-12.9034,0.00997369,-12.8932],[0,-83333.6,0,0,-139838,-12.8932,415348]]
QR - Q: [[1,-0,0,0,0,-0,0],[0,0.962245,0,0,-0.203347,-0.180926,0.00015472],[0,-0,1,0,0,-0,0],[0,-0,0,1,0,-0,0],[0,-0.192477,0,0,-0.97836,0.0759231,-6.02439e-05],[0,-2.97748e-05,0,0,-1.51316e-06,0.000698502,1],[0,-0.19245,0,0,-0.0382326,-0.980562,0.000679137]]
QR - R: [[1e+06,0,0,0,0,0,0],[-0,433015,-0,-0,332562,-7.44128,-133205],[0,0,1e+06,0,0,0,0],[0,0,0,1e+06,0,0,0],[0,0,0,0,1.98356e+06,15.7389,137878],[-0,-0,-0,-0,-0,13.9956,-402814],[0,0,0,0,0,0,264.716]]
determinant: -3.18215e+33
eigonvalue: [-2.01589e+06,0.00919742,343255,500000,1e+06,1e+06,1e+06]
decompose:
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: cholesky_decompose: Matrix m is not positive definite

4

The diagonal of the R matrix is always positive. In this case the determinant of Q is -1

1 Like
5

Thanks a lot.

6

Just to be clear: similar issues as you noted (minor inconsistencies) are to be somewhat expected due to floating point inaccuracies and cannot be resolved completely. But sometimes they may signal that there is some problem with Stan’s implementation, so it is good to let us know and thanks for that. A good way to see if this might be an issue with Stan (and not just unavoidable numerical imprecisions) is to see if other software (e.g. R, numpy) have similar inconsistencies with the same matrices. My bet would be that in this case, other software would have similar issues, but if you find an implementation that doesn’t, let us know - it might be worthwhile for us to check if we can learn something from this implementation to improve Stan.

Best of luck with your modelling :-)

1 Like
7

Thank you guys for developing Stan for us.

BTW. Is there any function for checking whether a matrix is positive definite or not? I only find related posts discussing C++ implementation, not for Stan users. If we can check a matrix first, we can walk around by add small amount when it’s the problem of floating point inaccuracies or stop the program when it’s a real non-positive definite.

My matrix is defined in the model block, so it cannot be declared as cov_matrix and take advantage of the auto checking.

8

Add one more case.

I did the Cholesky decomposition by hand. Here is a matrix (Phat) that Stan can do decomposition successfully. The right triangular matrix is the handmade decomposition. They are the same. I do this comparison to make sure I did it right.

屏幕快照 2020-05-14 下午6.23.28
屏幕快照 2020-05-14 下午6.23.281960×560 342 KB

This is the bad case. The matrix (P[3]) cannot be decomposed by Stan and the results of the qr-R function, the determinant function, and the eigenvalue function are consistent although different.
I think it is in the situation of floating point inaccuracies. But the handmade Cholesky decomposition, the right side triangular matrix, seems ok.

屏幕快照 2020-05-14 下午6.24.05
屏幕快照 2020-05-14 下午6.24.051960×1068 1.06 MB

1-1 Phat: [[1.00001e+06,986595,-2.35392,-0.0695205,-6.20572e-12,-5.54236e-14,-0.909642],[986595,1.72653e+06,-87516.3,-2584.7,100245,1550.16,361955],[-2.35392,-87516.3,613221,-11429.7,242552,-4.54747e-13,-149469],[-0.0695205,-2584.7,-11429.7,999746,-4949.32,-2138.21,970175],[-6.20572e-12,100245,242552,-4949.32,256374,522.352,59759.1],[-5.54236e-14,1550.16,-3.97904e-13,-2138.21,522.352,30.9068,-1001.52],[-0.909642,361955,-149469,970175,59759.1,-1001.52,1.77257e+06]]
decompose: [[1000,0,0,0,0,0,0],[986.59,867.854,0,0,0,0,0],[-0.00235391,-100.839,776.565,0,0,0,0],[-6.95201e-05,-2.97819,-15.105,999.754,0,0,0],[-6.20569e-15,115.509,327.339,0.339224,368.62,0,0],[-5.54233e-17,1.7862,0.231944,-2.12991,0.653323,4.76437,0],[-0.000909638,417.07,-138.317,969.566,153.36,52.573,783.041]]

1-1 P[3]: [[1.00001e+06,986593,-6.86094,-0.361009,-4.85556,-0.00959368,-2.7535],[986593,1.68123e+06,-187386,-9043.72,-7348.08,1337.58,321097],[-6.86094,-187386,393025,-25670.8,5326.86,-468.712,-239554],[-0.361009,-9043.72,-25670.8,998825,-20291.7,-2168.52,964348],[-4.85556,-7348.08,5326.86,-20291.7,803.547,17.393,-37291.7],[-0.00959368,1337.58,-468.712,-2168.52,17.393,29.9091,-1193.28],[-2.7535,321097,-239554,964348,-37291.7,-1193.28,1.73572e+06]]
QR - Q: [[0.711866,-0.565053,-0.309915,0.146129,-0.23782,0.00144911,4.20798e-16],[0.702315,0.572737,0.314135,-0.148113,0.241051,-0.00147016,-4.79519e-06],[-4.88402e-06,-0.299169,0.873058,0.298058,-0.24378,0.000201814,1.17229e-16],[-2.56988e-07,-0.0144385,-0.0926745,0.771569,0.629185,0.00310373,-9.38503e-16],[-3.45648e-06,-0.0117217,0.00378779,-0.0106546,0.0146979,-0.272303,-0.961961],[-6.82935e-09,0.00213561,0.000825099,-0.00233018,-0.00171833,0.96218,-0.272388],[-1.96011e-06,0.512674,-0.185574,0.52195,-0.655603,-0.00678922,-0.020854]]
QR - R: [[1.40477e+06,1.88308e+06,-131610,-6353.75,-5164.08,939.399,225507],[-0,626326,-347406,482708,-24634.3,325.716,1.13194e+06],[-0,-0,331125,-296855,11147.8,433.468,-519892],[-0,-0,-0,1.26791e+06,-32454.1,-2634.07,1.53146e+06],[-0,-0,-0,-0,8624.42,-145.19,-395934],[-0,-0,-0,-0,-0,23.3517,-305.012],[-0,-0,-0,-0,-0,-0,5.67482e-10]]
determinant Q: -1
determinant R: 4.22171e+19
determinant: -4.80413e+19
eigonvalue: [-4.68585e-11,24.3854,192406,247021,530233,2.19332e+06,2.64664e+06]
decompose:
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: cholesky_decompose: Matrix m is not positive definite

9

The most efficient way to check is to try a cholesky decomposition. If it fails the matrix isn’t postitive definite (in floating point). That’s what all the stan functions do

10

Tht matrix has two negative Eigen values. It isn’t positive definite.

It would probably also be useful to scale your matrix so if doesn’t have such large entries.

11

Thanks a lot. I might misunderstand the relationship between the Eigen values and the Cholesky decomposition. I thought as the sqrt is available when computing the diagonal values of the lower triangular matrix in the Cholesky decomposition, the Eigen values might be all positive.