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 involves allocating memory for the transposed matrix, and as you will see in a moment, there's an easy way to avoid this.

I will now present a convenient approach, which closely follows the design of Math.NET's dense solvers (only the LU case will be covered here). Two classes will be introduced: SparseLU, which wraps the CSparse implementation, and SparseMatrixExtensions, which will add extension methods to Math.NET's SparseMatrix class.

SparseLU.cs

// 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);
        }
    }
}

SparseMatrixExtensions.cs

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);
        }

        /// <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></returns>
        public static LU SparseLU(this SparseMatrix matrix, ColumnOrdering ordering,
            double tol = 1.0)
        {
            return LU.Create(matrix, ordering, tol);
        }
    }
}