[ANSWERED] Performing Cholesky Decomposition of Complex Hermitian Matrices...

Question

I wish to create a function to calculate Cholesky Decomposition of a Complex Hermitian Matrix. I am using the Cholesky-Banachiewicz Algorithm from here : https://en.wikipedia.org/wiki/Cholesky_decomposition#Example

Here is my function : https://code.sololearn.com/cqEyMBcBBX8t/#cpp

But now if the matrix to be decomposed is :

5 2+3i 5i
2-3i 7 1+7i
-5i 1-7i 12

The program is returning :

2.2361 0 0
0.89443-1.3416i 2.6458 0
-2.2361i 0.37796-2.6458i 3.4641

When the expected answer is :

2.2361 0 0
0.89443-1.3416i 2.098 0
-2.2361i -0.953-2.384i 0.64

Note that all the functions used inside the Cholesky function work as expected and return correct answers. Then, why do I get a wrong answer? Kindly help.

Answer 1

197# - (L1.CTranspose().mat.at(j).at(k))
197# + (L1.mat.at(j).at(k).Conjugate());

204# - CTranspose().mat.at(j).at(k));
204# + mat.at(j).at(k).Conjugate());

Just to add a few tips
1) your pair is using copy constructor -> your Matrices are copied into the pair...
Better use make_pair(std::move(L1),std::move(L2)) so your matrices are going to be moved into the pair
2) There is no need to use .at(...) here (no need for boundary checks), use operator []
3) use size_t instead of int

Best Regards!

Answer 2

@Ozren Zaric
Thank You very much! I assigned the transposed elements, and thats why the answer was different. Silly me.

Answer 3

While I won't pretend to understand the topic itself and the language isn't my favored flavor, the rules of debugging all boil down to the same thing. If you can't see the problem easily, step by step is the only way. If you can step through and check values as you go you can find where something is not being calculated or passed properly. If stepping isn't an option, loading it up with writes is always an option for tracking values through loops.

One thing I have learned in 18 years is that expected values can sometimes be incorrect as well. I debugged a financial calculation for 2 days only to discover that the math others had used to give me the example to model was where the flaw was and my code was right.

Sorry I can't provide more direct help. Hopefully someone can. But if not, brute force debugging is hard to beat.

Answer 4

Post a minimal complete example. Your code is incomplete, this way nobody can really help you!
BTW you should really use blas/ lapack routines! I think they even have choleskey!

Answer 5

@Ozren Zaric

I have posted a minimal working code in the same link. Kindly try it now and let me know if you find the bug. Thank you for your time.
I don't know anything about blas/lapack routines, and so didn't use them.

Answer 6

@Adam

Thank You for your reply. I have tried debugging the code, by visualizing the matrix's state after each operation, but still couldn't find why the answer is printed incorrectly.

It may be that the problem is with the formula itself, or I misinterpreted it somehow while making the code. I am trying to find the bug myself, but two heads are sometimes better than one, and maybe I can't see the problem in the aspect other professional programmers can, and that may be the reason the error still persists.