What's the problem with my Cholesky decomposition C++ code?

[BRN]

Hi at all,
I have to calculate the Cholesky decomposition of a symmetric matrix and this is the C ++ code I wrote:

boost::numeric::ublas::matrix<double> Math::cholesky(const boost::numeric::ublas::matrix<double> &MatrixA)
{
    int dim = MatrixA.size1();

    boost::numeric::ublas::matrix<double> MatrixL (dim,dim);

    double sum;
    for (int i = 0; i < dim; ++i)
    {
        for (int j = 0; j <= i; ++j)
        {
            sum = 0;
            if (j == i)
            {
                for (int k = 0; k < j - 1; ++k)
                {
                    sum += (MatrixL(j,k) * MatrixL(j,k));
                    MatrixL(j,j) = sqrt(MatrixA(j,j) - sum);
                }
            }
            else
                {
                    for (int k = 0; k < j; ++k)
                    {
                        sum += (MatrixL(i,k) * MatrixL(j,k));
                        MatrixL(i,j) = (MatrixA(i,j) - sum) / MatrixL(j,j);
                    }
                }
        }
    }

    return MatrixL;
}

But this don't work!
My output is wrong:
[CODE title="MatrixA"]25 15 -5
15 18 0
-5 0 11[/CODE]

[CODE title="MatrixL"]4.6396e-310 0 0
0 3.2387e-319 2.122e-314
0 0 3.31662 [/CODE]

what's wrong with my code?
thank you!

 

[jedishrfu]

Try stepping through it line by line printing key values. Key values might be your indices too so you know what cell of your matrix is being referenced. It should become obvious once you see the output.

The ++i and ++j seems odd as most programmers use it as i++ and j++ ie increment after the fact.

The other thing if the if block inside your second for loop and the i==j condition could be refactored to be outside of the second for loop right? which eliminates the if block too ie the else will be executed inside the loop and the i==j code will be executed after the second loop has completed. This also means the j<=i loop condition can be changed to i<j and the i==j code can be run following the second loop.

I noticed that you are changing the values of the L matrix when computing the A matrix won’t that affect your answer adversely.

 

[Filip Larsen]

The ++i and ++j seems odd as most programmers use it as i++ and j++ ie increment after the fact.

Just a small comment on that.

Pre-increment is to my knowledge still recommended over post-increment in "pure" increment statements of general iterators to avoid construction of a temporary object associated with post-increment, and for consistency this idiom then carries over to primitive data types as well.

I am not aware if modern compilers in general can optimize this temporary object away, but judging from the Notes section of https://en.cppreference.com/w/cpp/language/operator_incdec I would think not. As it is, I am still using this idiom for my iterators.

 

[jim mcnamara]

@jedishrfu 's point is the interchange of post-increment and pre-increment operators changes the logic of the program. On a casual look we lose "expected output" of "1" at line #8 as follows:

Owner@Owner-PC ~
$ cat -n t.cpp
     1  #include <iostream>
     2
     3  int main()
     4  {
     5      int i = 0;
     6
     7      std::cout << "i = " << i++ << '\n';
     8      std::cout << "i = " << ++i << '\n';
     9  }
Owner@Owner-PC ~
$ ./t
i = 0
i = 2

@Filip Larsen I am sure you are aware of this, but your wording, IMO, is a bit misleading. Newbies (to C++) may encounter problems with your suggestion. Those two operators are not 100% interchangeable.

 

[jim mcnamara]

@BRN -
If you write C++ you will have to learn a little about how to debug code.

If you are on linux, you probably have gdb (it is free to download), Windows 10 has WSL as well. Introduction to each debugging environment:



This uses GDB in Windows 10 on C++ code written using WSL linux

 

[jedishrfu]

There were some programming puzzle books where many of the problems used post and preincrement in the context of referencing array elements and it was quite challenging to see how things worked out.

In some cases, at the time, it was compiler dependent and so I would follow the common convention of post increment and to never use it in the context of any expression to avoid confusion.

Unless of course I wanted to insert some tricky hard to visually debug code for some future programmer.

 

[jedishrfu]

Here's the C Puzzle book for those who want to plumb the depths from C to shining C!

https://www.amazon.com/dp/0201604612/?tag=pfamazon01-20

 

[Filip Larsen]

I am sure you are aware of this, but your wording, IMO, is a bit misleading.

Well, in the context of this thread I would think that calling it a "pure increment statement" (as opposed to expression) should be hint enough that the effective result of "++i" and "i++" is the same, even if the compiler (for general iterator classes) may produce slightly different code. If a person new to C++ reads my comment and then, without any additional understanding as to what difference of pre- and post-increment is, applies a pre-increment in an expression where post-increment semantics is required, then I would claim he lacks the proper understanding of basic C++ in the first place.

Another way to deliver my point could also be to note that if one put value in selecting the operation that most closely follow the meaning of a pure increment statement (namely the semantics that after this statement the iterator is now pointing to next element), then the pre-increment operator is the right choice whereas the semantics of post-increment introduce additional elements that is irrelevant. In fact, for some container iterators the post-increment is directly implemented using pre-increment (and a temporary copy of the original iterator value).

Sorry for hammering out such minor detail in the context of this thread; now back to the regular show :)

 

[jedishrfu]

Yes, let's limit the discussion to the code and agree to disagree.

I would NOT recommend prefix operators to be used along with postfix operators as it can cause unwanted complication for future programmers working with the code. Using them by themselves usually is no problem, it's when we try to mix them in with expressions or method/function calls that their side-effects become known.

 

[BRN]

I use "++i" because my IDE (Clion) is added by default. Also in my programming course it was said to prefer "++ i"...

@jim mcnamara now I have installed gdb and will try to learn how to use it. Thanks for your advice!

@jedishrfu and @Filip Larsen thanks for the book and for the explanations!

I rewrote the code by removing the "if/else" statement and using a conditional statement in this way:

// Cholesky Decomposition
boost::numeric::ublas::matrix<double> Math::cholesky(const boost::numeric::ublas::matrix<double> &MatrixA)
{
    int dim = MatrixA.size1();

    boost::numeric::ublas::matrix<double> MatrixL(boost::numeric::ublas::zero_matrix<double>(dim,dim));

    double sum;
    for (int i = 0; i < dim ; ++i)
    {
        for (int k = 0; k < i + 1; ++k)
        {
            sum = 0;
            for (int j = 0; j < k; ++j)
            {
                sum += MatrixL(i,j) * MatrixL(k,j);
            }

            MatrixL(i,k) = (i == k) ? sqrt(MatrixA(i,i) - sum) :
                                                        ((1./MatrixL(k,k)) * (MatrixA(i,k) - sum));
        }
    }
    return MatrixL;
}

also preloading the matrix with zeros.
Now, it works fine!

Thank you all!