Lightweight fast matrix class in C# (Strassen algorithm, LU decomposition)

This is lightweight (19 kB) matirx class written in C#, that includes basic matrix operations.

All this is written in C#. Parse() and ToString() methods included. Matrix class can throw exceptions (MException).

Implemented matrix operations

  • exponentiation by integer
  • LU decomposition
  • determinant
  • inversion
  • solve system of linear equations
  • transpose, …

Overridden infix operators

  • array-like access (a[2,4])
  • addition (a + b)
  • subtraction (a – b)
  • unary minus (-a)
  • multiplication (a * b)
  • constant multiplication (5.4 * a)

Strassen algorithm (fast multiplication of large matrices)

It contains matrix multiplication with Strassen algorithm, which is memory efficient and cahce-oblivious, it works with still the same array in each level of recursion.

It is about 3 times faster in multiplication of two 500×500 matrices than the standard algorithm (tested). Used only when multiplying matrices larger than 32.

Documentation

Show documentationHide documentation

Example

// solving system of linear equations ... A * x = b
Matrix A = Matrix.Parse("1 2\r\n-3 4.5"); // \r\n is newline in textbox
Matrix b = Matrix.Parse("-1\r\n-3);
Matrix x = A.SolveWith(b);

// LU decomposition (P is permutation matrix)
Matrix PLU = A.GetP() * A.L * A.U;   // -> PLU == A
double d = A.Det();
Matrix I = A.Invert();
Matrix Powered = Matrix.Power(A, 10);
Matrix S = A.Duplicate();
Matrix Z = Matrix.ZeroMatrix(4,3);
Matrix I = Matrix.IdentityMatrix(6,7);
Matrix R = Matrix.RandomMatrix(4, 4, 80); // rows, cols, dispersion

Download

CS file

EXE example (.NET – Windows)

Old comments (closed because of spam)

25 Comments

  1. Griggs said:

    Amazing class, right what I needed! Thanks!

    May 2nd, 2011
  2. Nick said:

    Mate, that is really great. Thank you for the code. Ideal!

    May 18th, 2011
  3. Stephan said:

    Great! I was looking for some performance boosts in my .NET application and found the strassen-algorithm in your class. This gave my application a good performance improvement.

    Maybe I will implement the coppersmith–winograd algorithm too.

    Thanks for your work.

    Best regards!

    August 14th, 2011
  4. Alexis said:

    Ivan, thanks for you work! Your lightweight class with Strassen algorithm helps me very mutch.

    Best regards!

    March 31st, 2012
  5. Andy said:

    Hi Ivan, thanks so much! This is excellent. I’m writing a little app to help my son with his linear algebra homework and I found your matrix library – it saved me a lot of work.

    I have one bug fix and two new features to share with you (since you were kind enough to share your library with others):

    Bug Fix:
    The transpose method only worked on square matrices, I changed the last line from
    t[i,j] = m[j,i];
    to
    t[j, i] = m[i, j];

    This seemed to fix it.

    Feature:
    I added some code to dummy-proof the string passed to the parse method.

    // Try to dummy proof the string before we pass it to Matrix.Parse().
    // The Matrix.Parse() method expects rows to end with newline (except last
    // row) and expects columns (row elements) to be delimited with a single
    // space and no trailing spaces at end of row. public static string NormalizeMatrixString(string matStr)
    {
    // Remove any double spaces
    while (matStr.IndexOf(” “) != -1)
    {
    matStr = matStr.Replace(” “, ” “);
    }

    // Remove any spaces before or after newlines
    matStr = matStr.Replace(” \r\n”, “\r\n”);
    matStr = matStr.Replace(“\r\n “, “\r\n”);

    // If the data ends in a newline, remove the trailing newline.
    // Make it easier by first replacing \r\n’s with |’s then
    // restore the |’s with \r\n’s
    matStr = matStr.Replace(“\r\n”, “|”);
    while (matStr.LastIndexOf(“|”) == (matStr.Length – 1))
    {
    matStr = matStr.Substring(0, matStr.Length – 1);
    }
    matStr = matStr.Replace(“|”, “\r\n”);

    return matStr;
    }

    Then, in the Parse Method, pass the input string to the NormalizeMatrixString method before parsing it.

    New Feature 2:
    I found an algorithm for reduced row echelon form at
    http://rosettacode.org/wiki/Reduced_row_echelon_form, I added this method to your library with very little tweaking (I just had to rename a few variables).

    Thanks again for sharing your work!

    Andy

    May 25th, 2012
  6. Bernard Mensah said:

    Theres an error in your transpose method which causes out of bounds exceptions

    June 14th, 2012
  7. Ivan Kuckir said:

    Thank you guys!

    I corrected Transpose method, now it works with any rectangular matrices.

    I also added matrix string normalization. Thanks Andy!

    June 14th, 2012
  8. Blue said:

    The method invert does not work fine for matrices 12×12. It gives different result than for example inverting in open calc. However, using your method of getting determinant made my code a lot faster.

    July 29th, 2012
  9. Blue said:

    OK, I have made some more tests and your method does work :). Unfortunately, all online calculators I spotted able to invert matrices (and also your class) happen to give different results when determinant is very small (open calc too). So I guess it would be a brilliant idea to add a method which let us operate on a matrix with a very small determinant.

    August 1st, 2012
  10. doubt said:

    where can i find Main method??

    September 27th, 2012
  11. Thomas said:

    Hi Ivan,

    your class looks like it could be very useful. Could you be convinced to include a license snippet at the top?

    Thank you!

    Thomas

    December 8th, 2012
  12. hrobeers said:

    Hi Ivan,

    Can you please include a license? So we know what our rights are?

    I suppose that you intend to publish this under public domain? If so, please mention.

    Without a license, nobody is actually allowed to use this code.

    Keep up the good work!

    January 23rd, 2013
  13. VCCGeek said:

    Ivan, you’ve done some great work here. This is actually the only generalised C# implementation of a matrix that I can find anywhere. All of the others only work on 3×3 or 4×4 matrices, which is practically useless to me.

    One thing to note is that by modifying your code to use jagged arrays, I was able to speed up matrix multiplication by about 30% (21 seconds versus 30 seconds on a 2012 MBP) on a 30×30 matrix. I’m still trying to get some speedups, though, because I’m dealing with matrices on the order of 2000×2000 (Matlab does that in 1.3 seconds, but that’s not a production language). I think the problem may be in the number of copies the code makes during Strassen’s algorithm. I get an out-of-memory error if I get much larger. If we could find a way to do it in-place – or at least with fewer copies – we might get some impressive speed. I’m convinced that we could match Matlab’s performance, or at least come close.

    Again, great work. It’s saved me a great deal of trouble, and it’s well-coded. Much appreciated!

    March 18th, 2013
  14. The_Keeper said:

    Great class, thank you very much!

    I believe there’s a mistake in the comment beside RandomMatrix method:
    public static Matrix RandomMatrix(int iRows, int iCols, int dispersion) // Function generates the zero matrix

    March 20th, 2013
  15. Ivan Kuckir said:

    Thomas, hrobeers : I attached the MIT licence.

    VCCGeek : thank you :) I am planning to rewrtie it into 1D array, think it will be even faster than Jagged arrays.

    The_Keeper : thanks for reporting it. There really was a “bug” in comment :)

    March 20th, 2013
  16. VCCGeek said:

    According to my experimentation, it’s not (much) faster with a 1-D, but it’s a bit easier to implement. I think the reason 2-D arrays are so slow is because they have the overhead of several extra function calls. Badly implemented, IMHO. (Meaning the 2d arrays, of course, not your Matrix class.)

    Here’s how I rewrote it to a jagged array. I’ll change it to a 1-D, since that’s what you want to use. Hope it saves you some time. I hereby give this under the MIT licence as stated in the above source, yada yada. :)

    Accessor:

    public double this[int iRow, int iCol] // Access this matrix as a 2D array
    {
    get { return mat[iRow*rows + iCol]; }
    set { mat[iRow*rows + iCol] = value; }
    }

    And change your constructor:
    public Matrix(int iRows, int iCols) // Matrix Class constructor
    {
    rows = iRows;
    cols = iCols;
    mat = new double[rows*cols];
    }

    Column accessor methods:
    public Matrix GetCol(int k)
    {
    Matrix m = new Matrix(rows, 1);
    for (int i = 0; i < rows; i++) m[i, 0] = mat[i*cols + k];
    return m;
    }
    public void SetCol(Matrix v, int k)
    {
    for (int i = 0; i < rows; i++) mat[i*cols + k] = v[i, 0];
    }

    Duplicator:
    public Matrix Duplicate() // Function returns the copy of this matrix
    {
    Matrix matrix = new Matrix(rows, cols);
    for (int i = 0; i < rows; i++)
    for (int j = 0; j < cols; j++)
    matrix[i, j] = mat[i*cols + j];
    return matrix;
    }

    ToString() override:
    public override string ToString() // Function returns matrix as a string
    {
    string s = "";
    for (int i = 0; i < rows; i++)
    {
    for (int j = 0; j 1000, so that would partially explain my results. I think it would be almost impossible to implement this without making any copies, but I do think we could reduce it. That would really produce a performance boost, particularly on the huge matrices I’m throwing around. Still, it’s great work. Keep it up!

    March 21st, 2013
  17. VCCGeek said:

    Yuck. It got rid of my formatting. Does this posting system have a code block tag?

    March 21st, 2013
  18. VCCGeek said:

    On a second look, it also got rid of a bunch of my code and part of my post. Here’s the ToString() method again (going to try an html tag here):

    public override string ToString() // Function returns matrix as a string
    {
    string s = “”;
    for (int i = 0; i < rows; i++)
    {
    for (int j = 0; j < cols; j++) s += String.Format("{0,5:0.00}", mat[i][j]) + " ";
    s += "\r\n";
    }
    return s;
    }

    March 21st, 2013
  19. Alexander Popov said:

    Maybe you should add some type of exception handling in public double this[int iRow, int iCol] when iRow or iCol is outside the matrix range.

    July 8th, 2013
  20. Nepi Nepuansky said:

    I’ve tested your code and i tried 1000 multiplications with random generated matrix and their size. After that i compared results in matlab and everything is perfect!!!

    August 27th, 2013
  21. Kadir Güngör said:

    Thanks bro, you saved me.
    Great job, great class, it’s matrix revolution!

    October 16th, 2013
  22. fvdl said:

    Hi,

    Great class, saved me a lot of time.

    October 26th, 2013
  23. Recommendation for C# Matrix Library | Ask Programming & Technology said:

    [...] Lightweight fast matrix class in C# (Strassen algorithm, LU decomposition) [...]

    November 2nd, 2013
  24. marc jacobs et louis vuitton said:

    louis vuitton totally mm

    January 12th, 2014
  25. Darrell Pittman said:

    Thanks for posting this class! I had been looking for a C# implementation of the Strassen algorithm. I will likely use it in special cases.

    Compared to the old three-nested-loop naïve algorithm, Strassen is much faster, but only for square matrices whose dimensions are powers of two. In cases of non-square matrices, or matrices whose dimensions are not powers of two, it is much slower (sometimes 2-3x slower) than the naïve algorithm.

    Perhaps in your cutover to Strassen multiplication, in addition to the size being over 32, you should also check for square matrices whose dimensions are powers of two.

    I put the Matrix class I downloaded here into a new VS2010 project and made it a class library, then I wrote a console app to call that class’s matrix multiplication algorithm numerous times with various sizes, along with doing the multiplication the naïve way on the exact same data, printing the execution times as it goes.

    The following is the output of the test timings. Note that it was run on a Intel Core vPro, 2x3Ghz, running Windows 7 64-bit, with 12 GB RAM. The code was compiled in release mode for x64 with optimizations enabled:

    m1[ 32, 32] * m2[ 32, 32] : 1 ms for naive, 6 ms for Strassen
    m1[ 64, 64] * m2[ 64, 64] : 13 ms for naive, 11 ms for Strassen
    m1[ 128, 128] * m2[ 128, 128] : 104 ms for naive, 97 ms for Strassen
    m1[ 256, 256] * m2[ 256, 256] : 875 ms for naive, 498 ms for Strassen
    m1[ 512, 512] * m2[ 512, 512] : 6863 ms for naive, 3758 ms for Strassen
    m1[1024, 1024] * m2[1024, 1024] : 60495 ms for naive, 23855 ms for Strassen
    m1[ 128, 128] * m2[ 128, 1024] : 817 ms for naive, 23165 ms for Strassen
    m1[ 256, 256] * m2[ 256, 128] : 408 ms for naive, 480 ms for Strassen
    m1[ 256, 256] * m2[ 256, 512] : 1635 ms for naive, 3633 ms for Strassen
    m1[ 75, 75] * m2[ 75, 75] : 19 ms for naive, 87 ms for Strassen
    m1[ 150, 150] * m2[ 150, 150] : 151 ms for naive, 468 ms for Strassen
    m1[ 300, 300] * m2[ 300, 300] : 1217 ms for naive, 3637 ms for Strassen
    m1[ 600, 600] * m2[ 600, 600] : 11204 ms for naive, 23240 ms for Strassen
    m1[1100, 1100] * m2[1100, 1100] : 75918 ms for naive, 169865 ms for Strassen
    m1[ 36, 34] * m2[ 34, 34] : 2 ms for naive, 9 ms for Strassen
    m1[ 69, 65] * m2[ 65, 67] : 13 ms for naive, 91 ms for Strassen
    m1[ 129, 131] * m2[ 131, 129] : 97 ms for naive, 467 ms for Strassen
    m1[ 261, 257] * m2[ 257, 258] : 783 ms for naive, 3649 ms for Strassen
    m1[ 514, 516] * m2[ 516, 514] : 6999 ms for naive, 23473 ms for Strassen
    m1[1029, 1029] * m2[1029, 1025] : 61644 ms for naive, 162917 ms for Strassen
    m1[ 79, 76] * m2[ 76, 77] : 21 ms for naive, 87 ms for Strassen
    m1[ 152, 154] * m2[ 154, 155] : 168 ms for naive, 484 ms for Strassen
    m1[ 305, 305] * m2[ 305, 304] : 1373 ms for naive, 3646 ms for Strassen
    m1[ 603, 605] * m2[ 605, 604] : 11771 ms for naive, 23337 ms for Strassen
    m1[1101, 1102] * m2[1102, 1105] : 76688 ms for naive, 164019 ms for Strassen

    In all cases, the outputs of the Strassen and naïve algorithms matched.

    You can see that for square matrices whose sizes are powers of two, Strassen is much faster than the naïve algorithm. For other cases, it is significantly slower.

    The code that produced the timings is given below. As I mentioned, the Matrix class downloaded here was used verbatim, except for enclosing it in a namespace declaration.

    Code listing follows. I hope it formats okay…

    using System;
    using System.Collections.Generic;
    using System.Diagnostics;
    using System.Drawing;
    using TestStrassen;

    namespace StrassenApp
    {
    class Program
    {
    static Random random = new Random();

    const double Epsilon = 1e-6;

    static void Main(string[] args)
    {
    // Sizes to run.
    // Each pair gives the size if m1 and m2.
    // Negative values indicate a that small random size variation
    // should be applied so that the matrices are not square.
    var sizes = new List
    {
    // Square, powers of two
    new Point(32, 32), new Point(64, 64), new Point(128, 128),
    new Point(256, 256), new Point(512, 512), new Point(1024, 1024),
    // Non-square, powers of two
    new Point(128, 1024), new Point(256, 128), new Point(256, 512),
    // Square, non-powers of two
    new Point(75, 75), new Point(150, 150), new Point(300, 300),
    new Point(600, 600), new Point(1100, 1100),
    // Non-square, random variations
    new Point(-32, -32), new Point(-64, -64), new Point(-128, -128),
    new Point(-256, -256), new Point(-512, -512), new Point(-1024, -1024),
    new Point(-75, -75), new Point(-150, -150), new Point(-300, -300),
    new Point(-600, -600), new Point(-1100, -1100),
    };

    // Run an experiment for each size
    foreach (var size in sizes)
    {
    // Sizes of the matrices
    var m1Size = Math.Abs(size.X);
    var m2Size = Math.Abs(size.Y);

    // Whether or not to use square matrices
    var m1Randomize = size.X < 0;
    var m2Randomize = size.Y < 0;

    // Matrix sizes
    var m1Rows = m1Size + (m1Randomize ? 1 + random.Next(5) : 0);
    var m1Cols = m1Size + (m1Randomize ? 1 + random.Next(5) : 0);
    var m2Rows = m1Cols;
    var m2Cols = m2Size + (m2Randomize ? 1 + random.Next(5) : 0);

    // Allocate random matrices of the appropriate sizes
    var m1 = Matrix.RandomMatrix(m1Rows, m1Cols, 1);
    var m2 = Matrix.RandomMatrix(m2Rows, m2Cols, 1);

    // Start a timer on the naive matrix multiplication algorithm
    var stopwatchNaive = Stopwatch.StartNew();

    // Multiply matrices m1 and m2 the old-fashioned way.
    // Their product will be used later as the "expected" value.
    var expected = NaiveMatrixMultiplication(m1, m2);

    // Stop the timer for the naive algorithm, get the milliseconds elapsed
    stopwatchNaive.Stop();
    var naive = stopwatchNaive.ElapsedMilliseconds;

    // Start a different timer for the Strassen algorithm
    var stopwatchStrassen = Stopwatch.StartNew();

    // Do the Strassen multiplication
    var actual = m1 * m2;

    // Stop the Strassen timer, get the milliseconds elapsed
    stopwatchStrassen.Stop();
    var strassenMsec = stopwatchStrassen.ElapsedMilliseconds;

    // Write a message giving the timings to the trace log
    Console.WriteLine(
    string.Format("m1[{0,4}, {1,4}] * m2[{2,4}, {3,4}] : {4,6} ms for naive, {5,6} ms for Strassen",
    m1Rows, m1Cols, m2Rows, m2Cols, naive, strassenMsec));

    Debug.Assert(expected.rows == actual.rows);
    Debug.Assert(expected.cols == actual.cols);

    // Compare the two products element by element.
    var equal = true;
    for (var row = 0; equal && row < expected.rows; row++)
    {
    for (var col = 0; equal && col Epsilon)
    {
    Console.WriteLine(“Products not equal!”);
    equal = false;
    }
    }
    }
    }
    }

    private static Matrix NaiveMatrixMultiplication(Matrix a, Matrix b)
    {
    var product = new Matrix(a.rows, b.cols);

    // Matrix multiplication the old-fashioned way
    for (var row = 0; row < product.rows; row++)
    {
    for (var col = 0; col < product.cols; col++)
    {
    for (var inner = 0; inner < a.cols; inner++)
    {
    product[row, col] += a[row, inner] * b[inner, col];
    }
    }
    }

    return product;
    }
    }
    }

    March 9th, 2014