Math.NET Numerics and CSparse

Math.NET Numerics (version 3.x) doesn't provide any direct solvers for sparse matrices. On this page you will learn how to solve large systems in Math.NET Numerics using CSparse.

One difficulty you will face, is that Math.NET and CSparse use different storage schemes. While Math.NET uses Compressed Sparse Row storage (CSR), CSparse uses Compressed Sparse Column storage (CSC). Both schemes are closely related: If a sparse matrix is stored as CSR, its transpose will be the original matrix represented in CSC.

So the approach which immediately comes to mind, is to transpose the Math.NET matrix and feed the storage data to the CSparse solver routines. This is a valid approach, but it involves allocating memory for the transposed matrix, which you want to avoid when dealing with large matrices.

With CSparse.NET, you can then proceed as follows:

• SparseCholesky: for symmetric matrices both storage schemes are equivalent, so just call the Solve method without transposing the storage.
• SparseLU: use the SolveTransposed method without transposing the storage (see the example code below).
• SparseQR: transpose the storage and then call the Solve method.

I will now present a convenient approach for the LU case, which closely follows the design of Math.NET's dense solvers:

// Extend existing namespace.
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
{
    using CSparse;
    using CSparse.Double;
    using MathNet.Numerics.LinearAlgebra.Storage;
    using MathNet.Numerics.Properties;
    using System;

    // Create an alias for CSparse's SparseLU class.
    using CSparseLU = CSparse.Double.Factorization.SparseLU;

    public class SparseLU
    {
        int n;
        CSparseLU lu;

        private SparseLU(CSparseLU lu, int n)
        {
            this.n = n;
            this.lu = lu;
        }

        /// <summary>
        /// Compute the sparse LU factorization for given matrix.
        /// </summary>
        /// <param name="matrix">The matrix to factorize.</param>
        /// <param name="ordering">The column ordering method to use.</param>
        /// <param name="tol">Partial pivoting tolerance (form 0.0 to 1.0).</param>
        /// <returns>Sparse LU factorization.</returns>
        public static SparseLU Create(SparseMatrix matrix, ColumnOrdering ordering,
            double tol = 1.0)
        {
            int n = matrix.RowCount;

            // Check for proper dimensions.
            if (n != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.MatrixMustBeSparse);
            }

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

            // Create CSparse matrix.
            var A = new CompressedColumnStorage(n, n);

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

            return new SparseLU(new CSparseLU(A, ordering, tol), n);
        }

        /// <summary>
        /// Solves a system of linear equations, <c>Ax = b</c>, with A LU factorized.
        /// </summary>
        /// <param name="input">The right hand side vector, <c>b</c>.</param>
        /// <param name="result">The left hand side vector, <c>x</c>.</param>
        public void Solve(Vector<double> input, Vector<double> result)
        {
            // Check for proper arguments.
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            // Check for proper dimensions.
            if (input.Count != result.Count)
            {
                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
            }

            if (input.Count != n)
            {
                throw new ArgumentException("Dimensions don't match", "input");
            }

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

            if (b == null || x == null)
            {
                throw new NotSupportedException("Expected dense vector storage.");
            }

            // Copy the contents of input to result.
            Buffer.BlockCopy(b.Data, 0, x.Data, 0, b.Data.Length * Constants.SizeOfDouble);

            lu.SolveTranspose(x.Data);
        }
    }
}

You can use extension methods to extend the original SparseMatrix class:

namespace MathNet.Numerics.LinearAlgebra.Double
{
    using CSparse;

    using LU = MathNet.Numerics.LinearAlgebra.Double.Factorization.SparseLU;

    public static class SparseMatrixExtensions
    {
        /// <summary>
        /// Compute the sparse LU factorization for given matrix.
        /// </summary>
        /// <param name="matrix">The matrix to factorize.</param>
        /// <param name="tol">Partial pivoting tolerance (form 0.0 to 1.0).</param>
        /// <returns></returns>
        public static LU SparseLU(this SparseMatrix matrix, double tol = 1.0)
        {
            return LU.Create(matrix, ColumnOrdering.MinimumDegreeAtPlusA, tol);
        }
    }
}