Accuracy of Cholesky factorization

[diluculo]

To get covariance matrix in polynomial fitting, I tried to do Cholesky factorization. but I found something wrong I can't understand. Please see the following code.

           // independent variable x
           var x = Enumerable.Range(0, 21).Select(Convert.ToDouble).ToArray();
           // polynomial order
           int order = 5; 
           // number of points
           int Ns = x.Length; // 21
           // degree of freedom
           int DF = Ns - (order + 1); // 21 - 6 = 15

           // design Matrix X
           Matrix<double> X = Matrix<double>.Build.Dense(Ns, order + 1, (i, j) => Math.Pow(x[i], j));
           // Transposed matrix of X
           Matrix<double> Xt = X.Transpose();
           // XtX = X'X is a symmetric square matrix
           Matrix<double> XtX = Xt * X;

           // Let's do Cholesky factorization
           var chol = XtX.Cholesky();
           var cholL = chol.Factor;
           var cholNorm = (XtX - cholL * cholL.ConjugateTranspose()).L2Norm(); // 0.00390625
           // A is not LL'. Something wrong !!!

           // Let's see how DoCholesky method works? 
           // (1) 1st column
           var l11 = Math.Sqrt(XtX[0, 0]);
           var l21 = XtX[1, 0] / l11;
           var l31 = XtX[2, 0] / l11; // and so on 
           // (1a) doCholeskyStep
           var l22 = XtX[1, 1] - l21 * l21;
           var l32 = XtX[2, 1] - l21 * l31;
           var l33 = XtX[2, 2] - l31 * l31; // and so on
           // (2) 2nd column
           l22 = Math.Sqrt(l22); 
           l32 = l32 / l22; // and so on
           // (2a) doCholeskyStep
           l33 = l33 - l32 * l32; // and so on

           // so far all calculated values are exactly same to my reference values.

           // (3) 3rd column
           l33 = Math.Sqrt(l33); // 149.77538738613421,
           // but l33 should be     149.775387386134.
           //                                       ^
           // This small errors propagates through next steps.
           var l66 = cholL[5, 5]; // 21011.779008808317
           // it should be           21011.779009604055
           //                                   ^  

           // I think the error is related with the fact the followings give different values!!!
           var l33a = Math.Sqrt(XtX[2, 2] - l31 * l31 - l32 * l32); // 149.77538738613421
           var l33b = Math.Sqrt(XtX[2, 2] - (l31 * l31 + l32 * l32)); // 149.775387386134

[diluculo]

I think we can get || A - LL' || = 0, if we can rearrange the calculation steps like l33 = Math.Sqrt(A[2, 2] - (l31 * l31 + l32 * l32)).

[ggraham-412]

I implemented a re-arrangement of the Cholesky Factor routine. It is in the issue552 branch of my fork at https://github.com/ggraham-412/mathnet-numerics.git. I'll generate a pull request today, The example code you provided runs against that fork with the following results:

    static void Main(string[] args)
    {
        // independent variable x
        var x = Enumerable.Range(0, 21).Select(Convert.ToDouble).ToArray();
        // polynomial order
        int order = 5;
        // number of points
        int Ns = x.Length; // 21
                           // degree of freedom
        int DF = Ns - (order + 1); // 21 - 6 = 15

        // design Matrix X
        Matrix<double> X = Matrix<double>.Build.Dense(Ns, order + 1, (i, j) => Math.Pow(x[i], j));
        // Transposed matrix of X
        Matrix<double> Xt = X.Transpose();
        // XtX = X'X is a symmetric square matrix
        Matrix<double> XtX = Xt * X;

        // Let's do Cholesky factorization
        var chol = XtX.Cholesky();
        var cholL = chol.Factor;
        var cholNorm = (XtX - cholL * cholL.ConjugateTranspose()).L2Norm();   // 0.0 so A is now LL'
        var boolTest = (XtX.Equals(cholL * cholL.ConjugateTranspose()));      // true, A is LL' elementwise as well

    }