# 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 documentation –

Hide documentation
### Properties

int rows | Number of rows in matrix |

int cols | Number of columns in matrix |

Matrix L | Matrix L for LU decomposition |

Matrix U | Matrix U for LU decomposition |

### Methods of an Instance

Matrix(int r, int c) : void | Class constructor
r – number of rows, c – number of columns |

this[int r, int c] : double | You can access matrix values as 2D array
r – number of row, c – number of column |

GetCol(int k) : Matrix | Returns a k-th column of matrix
k – number of column that you want |

SetCol(Matrix m, int k) : void | Sets a k-th column of matrix
m – one-column matrix, k – index of target column |

MakeLU() : void | Makes a LU decomposition |

SolveWith(Matrix v) : Matrix | Computes a set of linear equations: A*x = v, returns “x”
v – vector of right results of equations |

Invert() : Matrix | Returns inverted matrix |

Det() : double | Returns a determinant |

GetP() : Matrix | Returns a permutation matrix P of PLU decomposition |

Duplicate() : Matrix | Returns a copy of the matrix |

ToString() : string | Converts the matrix to a “nice” string |

### Static methods

ZeroMatrix(int r, int c) : Matrix | Returns a matrix of zeros
r – number of rows, c – number of columns |

IdentityMatrix(int r, int c) : Matrix | Returns an identity matrix (ones on diagonal)
r – number of rows, c – number of columns |

Parse(string s) : Matrix | Parses a string into a matrix
s – source, spaces between numbers, \r\n as a new line |

Transpose(Matrix m) : Matrix | Transposes the matrix
m – source matrix |

Power(Matrix m, int pow) : Matrix | Powers the matrix to an integer
m – source matrix, pow – exponent (can be negative) |

### Operators

-Matrix A | Negation of matrix (multiplies the values by -1) |

Matrix A + Matrix B | Matrix addition |

Matrix A – Matrix B | Matrix subtraction |

Matrix A * Matrix B | Matrix multiplication |

double n * Matrix B | Multiplies the values by “n” |

## 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

Amazing class, right what I needed! Thanks!

May 2nd, 2011Mate, that is really great. Thank you for the code. Ideal!

May 18th, 2011Great! 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, 2011Ivan, thanks for you work! Your lightweight class with Strassen algorithm helps me very mutch.

Best regards!

March 31st, 2012Hi 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, 2012Theres an error in your transpose method which causes out of bounds exceptions

June 14th, 2012Thank you guys!

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

I also added matrix string normalization. Thanks Andy!

June 14th, 2012The 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, 2012OK, 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, 2012where can i find Main method??

September 27th, 2012Hi 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, 2012Hi 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, 2013Ivan, 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, 2013Great class, thank you very much!

I believe there’s a mistake in the comment beside RandomMatrix method:

March 20th, 2013public static Matrix RandomMatrix(int iRows, int iCols, int dispersion) // Function generates the zero matrix

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, 2013According 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:

March 21st, 2013public 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!

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

March 21st, 2013On 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

March 21st, 2013{

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;

}

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, 2013I’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, 2013Thanks bro, you saved me.

October 16th, 2013Great job, great class, it’s matrix revolution!

Hi,

Great class, saved me a lot of time.

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

November 2nd, 2013louis vuitton totally mm

January 12th, 2014Thanks 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}

}

}