Matrix must be positive definite.

[AlbertLiDesign]

Hello, when I attempted to solve a least sqaure question using Cholesky decomposition, I got an error, 'Matrix must be positive definite'. I have ensured that the matrix A is column full rank. Besides, I have attempted to output this matrix and I found that it can be calculated by numpy. I don't know where the problem is. This problem is urgent for me. Please help me. Thank you.

C# code: var X = A.TransposeThisAndMultiply(A).Cholesky();

Matrix A can be downloaded here https://drive.google.com/file/d/1yz5N4iDocHfS9YCstEpMvkBpyNrwQu1f/view?usp=sharing

Edited

I find the reason why I fail to decompose using Cholesky is using sparse matrix. But my matrix could be very large. It means if I use a dense matrix, it will cause serious time and memory consumption. Are there any good solutions here?

[wo80]

You might want to give CSparse.NET a try (sparse Cholesky example).

[AlbertLiDesign]

@wo80 Hi, I've tired CSparse.NET and convert the sparse matrix from Math.NET to CSparse. But I'm not sure whether my implement is correct because I don't know how to get the results in CSparse.NET. Could you please help me to check the code?

        public static void CSparseCholesky(Matrix<double> sparseMatrix, Matrix<double> input, Matrix<double> result)
        {
            int n = sparseMatrix.ColumnCount;
            int m = sparseMatrix.RowCount;
            // Apply column ordering to A to reduce fill-in.
            var order = CSparse.ColumnOrdering.MinimumDegreeAtPlusA;
            // Get CSR storage.
            var storage = (SparseCompressedRowMatrixStorage<double>)sparseMatrix.Storage;

            // Create CSparse matrix.
            var A = new CSparse.Double.SparseMatrix(m, n);

            // Assign storage arrays.
            A.ColumnPointers = storage.RowPointers;
            A.RowIndices = storage.ColumnIndices;
            A.Values = storage.Values;

            var AT = A.Transpose();
            var ATA = AT.Multiply(A);

            var b = input.Storage as DenseColumnMajorMatrixStorage<double>;
            var x = result.Storage as DenseColumnMajorMatrixStorage<double>;

            var c = new CSparse.Double.DenseMatrix(m,n);
            AT.Multiply(b.Data, c.Values);
            var cholesky = CSparse.Double.Factorization.SparseCholesky.Create(ATA, order);
            cholesky.Solve(c.Values, x.Data);
        }

[wo80]

1) CSparse doesn't support solving for rhs matrix, it expects a simple vector. You will have to extract the columns one by one. 2) CSparse uses compressed sparse column (CSC) storage instead of MathNet CSR, so you cannot simply interchange the storage, unless your matrix is symmetric.

static void CSparseCholesky(Matrix<double> A, Matrix<double> B, Matrix<double> X)
{
    int rows = A.RowCount;
    int columns = A.ColumnCount;

    if (rows < columns)
    {
        throw new Exception("underdetermined system");
    }

    if (B.ColumnCount != X.ColumnCount)
    {
        throw new Exception("X and B must have the same numnber of columns.");
    }

    var At = A.Transpose();
    var AtA = At * A;
    var C = At * B;

    // Get CSR storage of A'A (same as CSC storage, since matrix is symmetric).
    var storage = (SparseCompressedRowMatrixStorage<double>)AtA.Storage;

    // Create CSparse matrix.
    var csc = new CSparse.Double.SparseMatrix(columns, columns);

    // Assign storage arrays.
    csc.ColumnPointers = storage.RowPointers;
    csc.RowIndices = storage.ColumnIndices;
    csc.Values = storage.Values;

    // Apply column ordering to A to reduce fill-in.
    var order = CSparse.ColumnOrdering.MinimumDegreeAtPlusA;

    var cholesky = CSparse.Double.Factorization.SparseCholesky.Create(csc, order);

    var tmp = new DenseVector(columns);
    var x = new double[columns];

    for (int i = 0; i < B.ColumnCount; i++)
    {
        // Get column i of A'b.
        C.Column(i, tmp);

        cholesky.Solve(tmp.AsArray(), x);

        X.SetColumn(i, x);
    }
}

[AlbertLiDesign]

@wo80 Thank you sooooooooo much. It works now. I notice that you also developed parallel solving method. Looks amazing! I like your work!!