[Feature request] Add Matrix Multiplication

[ajos6183]

Firstly, thanks for your wonderful effort.

My question is if it possible to include matrix multiplication in CSparse.NET.
I currently have a loop using CSparse.NET to solve equations and Mathnet.Numerics to do multiplication. So i am having to convert back and forth between the the different formats of MathNet Vector and CSparse.Net double().

I don't know if there will be any performance improvements. This is my code.

Dim KeGlobalFree As SparseMatrix = MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.Build.Sparse(FreeDOFs.Count, FreeDOFs.Count)
Dim KgGlobalFree As SparseMatrix = MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.Build.Sparse(FreeDOFs.Count, FreeDOFs.Count)

........Missing Code

Dim EigenVectorEstimate = MathNet.Numerics.LinearAlgebra.Double.SparseVector.Build.Dense(KgGlobalFree.RowCount, 1)
Dim Errors As Double = 10 ' = MathNet.Numerics.LinearAlgebra.Double.SparseVector.Build.Dense(KgGlobalFree.RowCount, 1)
Dim EigenValueEstimate As Double = 1.0
Dim y2 = CSparse.Double.Vector.Create(KeGlobalFree.ColumnCount, 0)
While Errors >= 10 ^ -8
   Dim y1 = -KgGlobalFree * EigenVectorEstimate
   Dim storage = DirectCast(KeGlobalFree.Storage, MathNet.Numerics.LinearAlgebra.Storage.SparseCompressedRowMatrixStorage(Of Double)) ' Get CSR storage.
   Dim A = New CSparse.Double.SparseMatrix(KeGlobalFree.ColumnCount, KeGlobalFree.ColumnCount) With {.ColumnPointers = storage.RowPointers, .RowIndices = storage.ColumnIndices, .Values = storage.Values} ' Create CSparse matrix and Assign storage arrays.
   Dim SC = CSparse.Double.Factorization.SparseCholesky.Create(A, CSparse.ColumnOrdering.MinimumDegreeAtPlusA) ' Apply column ordering to A to reduce fill-in.
   SC.Solve(y1.ToArray, y2)
   Dim y2vec = MathNet.Numerics.LinearAlgebra.Double.SparseVector.Build.DenseOfArray(y2)
   Dim norm = y2vec * -KgGlobalFree * y2vec
   EigenVectorEstimate = y2vec.Divide(Math.Sqrt(norm))
   Dim NewEigenValueEstimate = EigenVectorEstimate * KeGlobalFree * EigenVectorEstimate
   Errors = Math.Abs((NewEigenValueEstimate - EigenValueEstimate) / EigenValueEstimate)
   EigenValueEstimate = NewEigenValueEstimate
End While

[epsi1on]

Are you multiplying a sparse matrix with a sparse vector, or sparse matrix with a dense vector?
I think dense vector multiplication is already in the library.

[ajos6183]

I am multiplying sparse matrices with dense vectors.
How do i use it? Do i just use the multiplication operator like A * B? where A is double(,) and B is double()

Thanks.

[wo80]

1. Actually, you shouldn't use operator overloads at all. Using MathNet, you are tempted to write code like -KgGlobalFree * EigenVectorEstimate, but this hides some important internals from the user. For example, the unary minus operator for sparse matrices will create a copy of the data. Doing this inside a loop is a bad idea, if you care about performance (lots of unnecessary memory allocations). The same argument holds for all operator overloads (matrix-vector etc.).

2. Since KeGlobalFree isn't modified inside the loop, you should create the factorization outside the loop.

3. Why don't you use CSparse all the way? The following C# code wasn't tested. It's just to get an idea:

using CSparse.Double;
using CSparse.Double.Factorization;
using CSparse.Storage;

void Example()
{
    int N = FreeDOFs.Count;

    // Use coordinate storage to assemble the matrices.
    var cKeGlobalFree = new CoordinateStorage<double>(N, N, 2 * N);
    var cKgGlobalFree = new CoordinateStorage<double>(N, N, 2 * N);
    
    // Convert to sparse matrix.
    var KeGlobalFree = CSparse.Converter.ToCompressedColumnStorage(cKeGlobalFree);
    var KgGlobalFree = CSparse.Converter.ToCompressedColumnStorage(cKgGlobalFree);

    double Errors = 10;
    double EigenValueEstimate = 1;

    var EigenVectorEstimate = Vector.Create(N, 1);

    var y1 = Vector.Create(N, 0);
    var y2 = Vector.Create(N, 0);
    
    var SC = SparseCholesky.Create(KeGlobalFree, CSparse.ColumnOrdering.MinimumDegreeAtPlusA);

    while (Errors >= 10e-8)
    {
        // y1 = -KgGlobalFree * EigenVectorEstimate
        KgGlobalFree.Multiply(-1, EigenVectorEstimate, 0, y1);

        // EigenVectorEstimate = KeGlobalFree \ y1
        SC.Solve(y1, EigenVectorEstimate);

        // y2 = -KgGlobalFree * EigenVectorEstimate
        KgGlobalFree.Multiply(-1, EigenVectorEstimate, 0, y2);

        // norm = EigenVectorEstimate * -KgGlobalFree * EigenVectorEstimate
        double norm = Vector.DotProduct(EigenVectorEstimate, y2);

        // EigenVectorEstimate = EigenVectorEstimate / sqrt(norm)
        Vector.Scale(Math.Sqrt(norm), EigenVectorEstimate);

        // y2 = KeGlobalFree * EigenVectorEstimate
        KeGlobalFree.Multiply(1, EigenVectorEstimate, 0, y2);

        // NewEigenValueEstimate = EigenVectorEstimate * KeGlobalFree * EigenVectorEstimate
        double NewEigenValueEstimate = Vector.DotProduct(y2, EigenVectorEstimate);

        Errors = Math.Abs((NewEigenValueEstimate - EigenValueEstimate) / EigenValueEstimate);

        EigenValueEstimate = NewEigenValueEstimate;
    }
}

[ajos6183]

I didn't know that a minus sign creates a copy of the data.
Thanks for your help! I will implement your suggestions.

[ajos6183]

wo80, I would like to thank you again! Your algorithm is a lot faster, about 70 times faster for a 700x700 sparse matrix! I still haven't utilized your idea of using CoordinateStorage because when i assign cKeGlobalFree it is always in an additive manner (cKeGlobalFree (i,j)=cKeGlobalFree (i,j)+x) and i couldn't figure out how to do this with CoordinateStorage.

For the benefit of someone who might be looking at this later, the following line of code should be changed:
From: Vector.Scale(Math.Sqrt(norm), EigenVectorEstimate);
To: Vector.Scale(1.0/Math.Sqrt(norm), EigenVectorEstimate);

I realize i am going off topic here, but do you know of any fast C# EigenProblem Solvers?
The EigenValue solver in Math.Net is really slow as it finds all eigenvalues. In reality, someone might be interested in the first couple EigenValues. The code above uses the inverse iteration method and is for solving the buckling problem. There are faster, more complicated methods out there such as the subspace iteration method. Any suggestions would be appreciated.

[wo80]

The CoordinateStorage will automatically add up all values you pass in, so

// Create a 2x2 matrix.
var C = new CoordinateStorage<double>(2, 2, 4);

// Add some values.
C.At(0, 0, 1.0);
C.At(1, 1, 1.0);
C.At(1, 1, 0.5); // Again position (1 ,1)

// Convert to sparse matrix.
var A = CSparse.Converter.ToCompressedColumnStorage(C);

will produce the matrix A with entries

[ 1.0  0.0 ]
[ 0.0  1.5 ]

Unfortunately, I don't know any sparse Eigenproblem solvers for .NET. I've used ARPACK / ARPACK++, calling the native code from C# using P/Invoke.

[ajos6183]

Thanks!
As you might have already noticed, i am not a programmer by training. Do you have any example using ARPACK in C# ?

I looked at http://sweb.cityu.edu.hk/kincau/blog/files/0e9a2b646554191e9ddc831b697cf04d-1.html and tried it in visual studio but i am left with 3 errors!

[wo80]

Even if the C# code would compile, you'd still need to build the native ARPACK DLL, which is rather involved.

I've used a lot of native solvers from C# and I'm planning to release all the code here on Github (including ARPACK), but at the moment, I don't have much time to work on it.