wo80 / CSparse.NET

A concise library for solving sparse linear systems with direct methods.
GNU Lesser General Public License v2.1
58 stars 25 forks source link

Rectangular QR decomposition with MathNet.Numerics #33

Closed kenryot closed 3 years ago

kenryot commented 3 years ago

Hi,

Thank you very much for publishing this library. Really great to have an open source sparse solver in .NET. I've got just one question. I encountered an error with feeding a rectangular matrix defined in MathNet.Numerics to CSparse for QR decomposition. I followed the tutorial from this page (https://github.com/wo80/CSparse.NET/wiki/Math.NET-Numerics-and-CSparse) but didn't manage to figure out the cause. It seems to work fine for square matrices.

Have anyone made it work? If it is supported, could you elaborate the procedure?

Best, Kenryo

wo80 commented 3 years ago

Maybe you forgot to switch the roles of rows and columns (CSC vs CSR storage). Here's the SparseQR.Create(...) method. If this doesn't help, I'll post the full code as a gist.

public static SparseQR Create(SparseMatrix matrix, CSparse.ColumnOrdering ordering)
{
    int rows = matrix.RowCount;
    int cols = matrix.ColumnCount;

    // Get CSR storage.
    var storage = (SparseCompressedRowMatrixStorage<double>)matrix.Storage;

    // Create CSparse matrix.
    var A = new CSparseMatrix(cols, rows);

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

    return new SparseQR(CSparseQR.Create(A, ordering), cols, rows);
}
kenryot commented 3 years ago

It seems I'm getting an error for not being a square matrix in AMD class.. Is it possible ?

private static SymbolicColumnStorage ConstructMatrix(SymbolicColumnStorage A, ColumnOrdering order) { SymbolicColumnStorage result = null;

        // Compute A'
        var AT = A.Transpose();

        int m = A.RowCount;
        int n = A.ColumnCount;

        if (order == ColumnOrdering.MinimumDegreeAtPlusA)
        {
            if (n != m)
            {
                throw new ArgumentException(Resources.MatrixSquare, "A");
            }

            // Return A+A'
            result = A.Add(AT);
        }
wo80 commented 3 years ago

If the matrix isn't square, you cannot compute A^t + A. Use ColumnOrdering.MinimumDegreeAtA instead.

kenryot commented 3 years ago

Right! Now the code runs. Thanks a lot.

Kenryo,