11,710,472 members (90,360 online)

# Efficient Matrix Programming in C#

, 11 Dec 2002 341.2K 9.1K 146
 Rate this:
Fast matrix expressions evaluation, based on dynamic code generation and partial evaluation

## Introduction

This article shows a technique to write clear and efficient matrix math code in C# language. The matrix operations are expressed using operator redefinition, but the code is generated dynamically, following the principle of partial evaluation. This gives a solution efficient and elegant as the one obtained by the Expression Templates in C++. This article is the prosecution of my studies in dynamic code generation and operator overloading with C# (here).

## The problem

An object oriented language as C#, lets you define very elegant matrix classes, that encapsulate element access, and with the mechanism of operator overloading you can write easily:

```Matrix m = new Matrix(10,5);
Matrix n = new Matrix(10,5);
Matrix r;
double d = 2;

r = m + n;
r = n*2;
```

The classic implementation of these classes is to create a temporary `Matrix` object for each sub expression, and the problem is much more evident in the following code where we generate two temporary objects.

```Matrix a,b,c,r; // with m and n initialized

r = a+3*(b+c);
```

This approach introduces a big overhead in object based matrix code, and usually the math programmers use arrays and indices losing all the object orientation. This kind of code can be called C#Tran (from C#/Fortran) like the JavaTran coined by Veldhuizen.

```double [] a, b, c, r;

for(int i = 0; i < a.length; i++)
r[i] = a[i]+3*(b[i]+c[i]);
```

## Solutions

A possible solution to this problem is the defered evaluation: when writing b+c we generate an object `MatrixAdd` that incapuslates the operation. It has overloaded operators too and we can write 3*(b+c) generating a `MatrixAddMul` object that evaluates the operation efficiently. This time the overhead is relative to all possible combinations of operations and the fact that if our specific expression is not covered by a deferred operation, we revert to the inefficient temporary based code.

My proposal is to use a generative programming approach: generate specialized code for specialized operations. I use the dynamic code capabilities of the .NET framework to generate on the fly the for-based code directly in IL. A compiled expression can be evaluated one time or stored for reuse into a delegate of type `Evaluator`.

These are examples of usage:

```Vector a = new Vector(3);
Vector b = new Vector(3);
Vector r = new Vector(3);

(a+b*3).AssignTo(r)();
```

The IL code generated is equivalent to the following. You should notice that the n size is known at the generation moment, so the for loop has a constant value: this is the partial evaluation, because the code is specialized to specific values and the JIT can generate optimized code.

```int n = a.rows();
for(int i = 0; i < <<n>>; i++)
{
r[i] = a[i]+b[i]*3;
}
```

This is another example with `Matrix` objects:

```Matrix a, b, c, r;

// evaluate directly
(a + Expr.Cos(b*c*3)).AssignTo(r)();

// store and reuse the compiled expression ...
Evaluator ev = (a + Expr.Cos(b*c*3)).AssignTo(r);

// and then evaluate
ev();
```

This is the equivalent code:

```int n = a.rows();
int m = a.cols();

for(int i = 0; i < <<n>>; i++)
for (int j = 0; j < <<m>>; j++)
{
double d = 0;
for(int k = 0; k < <<q>>; k++)
d += b[i,k]*q[k,j];
r[i,j] = a[i,j]+cos(d*3);
}
```

## The implementation

Internally each expression is stored as a tree of objects derived from the class `Expr`, and assigned to an object derived from `LeftExpr`. We can evaluate the expression slowly, or generate dynamically the IL code, that evaluates the expression. We consider operations with vectors and two dimensional matrices but the system can be extended to more dimensions. Before generating the code we check the sizes of the matrices involved in the operations, throwing the exception `SizeMismatchException` when they are not correct.

An expression of class `Expr` has the following virtual methods:

```// evaluates the expression against the two indices
float Eval(int i, int j)
// compile the expression using the context
void  Compile(ILGenerator g, CompilerContext cc)
// a property to find the size of the expression
OpSize Size { get }
```

The `LeftExpr` is a specialized `Expr` that handles assignments:

```// assign a value v to the Expression
void Assign(int i,  int j, float v);
// compile the assignment (two phase compilation)
void CompileAssign(ILGenerator g, CompilerContext cc, bool post);
```

The hierarchy of the classes that I've defined is shown in the following image:

The class generated dynamically is more complex than the one described in the previous article, this time we have to store all the array objects, the parameter values, and the indices `i`, `j`. To analyze the generated code, set the `SafeMode` variable in the `Compiler` class before generating the first Expression, and then, at the end, generate the assembly with `Compiler.Save()`.

## SubMatrix Access

One of the features obtained with this system is the fast access to submatrices. You can evaluate expression relative to parts of a matrix, like rows, columns or the main diagonal in case of a square matrix. The access to the different parts of the matrix is efficient because of the dynamic compilation. I use also an easy to read notation (I think), based upon the `Tag` Class approach.

The code is now much clearer than the standard array manipulation. The following code shows some submatrix accesses and their MATLAB equivalent (remember that MATLAB uses 1-based indexing).

```Matrix m;

// a row, like MATLAB m(3,:)
m[2][Matrix.All]

// a range of rows like m(3:4, :)
m[Matrix.Range(2,3)][Matrix.All]

// a column, like m(:,5)
m[Matrix.All][4]

// columns, like m(:,3:6)
m[Matrix.Range(2,5)[Matrix.All]

// a submatrix m(2:3,4:5)
m[Matrix.Range(1,2)][Matrix.Range(3,4)]

// a single element m(3,4)
m[2][3]

// the third element taking the matrix as a row vector
// there is no MATLAB equivalent because MATLAB
//uses column ordered matrices
m[2]

// some elements taking the matrix as a row vector
m[Matrix.Range(2,3)]

// the diagonal
m.Diagonal()
```

These operations are based on `MatrixRowRef` structures that stores the first array access and generates a `SubMatrix` object. It's a `LeftExpr` derived class so can be used in the assignments.

```Vector v = new Vector(3);
Matrix m = new Matrix(3,5);
(v + 22).AssignTo(m[Matrix.All][2])();
```

## Speed Considerations

The code generated by this system is as efficient as the `for` based code giving high speed math evaluation, but there is a little overhead during the generation phase, so this approach is effective when the expression is evaluated many times. Specifically the first expression dynamically evaluated is a little slower than the others because of the creation of the dynamic Assembly.

## Conclusion

The dynamic code generation of .NET, combined with the C# custom operators, let us build elegant and efficient math expressions, that let use C# as a Math Language. Citing Jim Blinn, "we have to put as much intelligence in the class libraries as possible, to make things as easy as possible for users of those classes". The final user of this code should just know that it's fast and easy to use.

### References

A list of licenses authors might use can be found here

## Share

 Software Developer (Senior) Scuola Superiore S.Anna Italy
Assistant Professor in Applied Mechanics working in Virtual Reality, Robotics and having fun with Programming

## You may also be interested in...

 First PrevNext
 Generating symbols using matrix in C# Namithasujir14-Feb-10 19:23 Namithasujir 14-Feb-10 19:23
 Does anyone use this? prictor1-Oct-09 15:10 prictor 1-Oct-09 15:10
 Re: Does anyone use this? aspdotnetdev17-Oct-09 16:39 aspdotnetdev 17-Oct-09 16:39
 Re: Does anyone use this? prictor9-Nov-09 12:04 prictor 9-Nov-09 12:04
 Bug: natrual Log function prictor1-Oct-09 15:07 prictor 1-Oct-09 15:07
 Bug - MatrixMulOp DajMiSciagnacPlik20-Dec-08 1:10 DajMiSciagnacPlik 20-Dec-08 1:10
 Matrix Multiplicaion Locking Raul Rangel24-Nov-08 19:35 Raul Rangel 24-Nov-08 19:35
 Hello, I'm trying to multiply a 3x3 with a 3x3 but when i execute it, it just locks up. ```int n = 3; Matrix P = new Matrix(n, n); Matrix Ab_inv = new Matrix(n, n); Matrix Ab_invNew = new Matrix(n, n);   Evaluator Eval_Atinv = (P * Ab_inv ).AssignTo(Ab_invNew);   Eval_Atinv(); // This is where it locks up ``` I was wondering if this was a known bug, and if theres any fixes ? Also is it possible to do something like ```Evaluator Eval_Atinv = (P * Ab_inv).AssignTo(Ab_inv); Eval_Atinv(); ``` Notice its multiplying and setting the same same variable ?
 Re: Matrix Multiplicaion Locking Raul Rangel24-Nov-08 20:13 Raul Rangel 24-Nov-08 20:13
 Performance Dmitri Nesteruk16-Aug-08 10:43 Dmitri Nesteruk 16-Aug-08 10:43
 Sub Matrix strelco25-Sep-07 3:33 strelco 25-Sep-07 3:33
 Generic Efficient Matrix Programming in C# Edberg1-Jun-07 10:02 Edberg 1-Jun-07 10:02
 Re: Generic Efficient Matrix Programming in C# Emanuele Ruffaldi1-Jun-07 10:46 Emanuele Ruffaldi 1-Jun-07 10:46
 Merge image mcavasan10-Mar-07 22:30 mcavasan 10-Mar-07 22:30
 How can I compute matrix determinant? Shokker1220-Nov-06 0:25 Shokker12 20-Nov-06 0:25
 Re: How can I compute matrix determinant? Edberg1-Jun-07 8:52 Edberg 1-Jun-07 8:52
 Hanging code JahBreeze26-Apr-06 12:57 JahBreeze 26-Apr-06 12:57
 Re: Hanging code Raul Rangel24-Nov-08 16:25 Raul Rangel 24-Nov-08 16:25
 Using COMObject from Matlab to C# nelson_yan@bezeqint.net30-Apr-05 23:57 nelson_yan@bezeqint.net 30-Apr-05 23:57
 Solution from matlab's website FinnHawk3-Nov-05 7:23 FinnHawk 3-Nov-05 7:23
 Inverse operator Pål Berg16-Nov-04 1:30 Pål Berg 16-Nov-04 1:30
 Re: Inverse operator Edberg1-Jun-07 8:50 Edberg 1-Jun-07 8:50
 InvalidProgramException Tata Taufik1-Sep-04 16:45 Tata Taufik 1-Sep-04 16:45
 Re: InvalidProgramException Chris Hawkins.12-Dec-04 11:23 Chris Hawkins. 12-Dec-04 11:23
 I can't execute it Darien_Chiba29-Jul-04 17:23 Darien_Chiba 29-Jul-04 17:23
 Very nice, going on Jonathan de Halleux12-Feb-04 23:40 Jonathan de Halleux 12-Feb-04 23:40
 Re: Very nice, going on Sean L. Palmer25-Jan-09 11:56 Sean L. Palmer 25-Jan-09 11:56
 How about performance MATLAB ? khomsak9-Jan-04 18:34 khomsak 9-Jan-04 18:34
 Re: How about performance MATLAB ? Emanuele Ruffaldi10-Jan-04 0:54 Emanuele Ruffaldi 10-Jan-04 0:54
 Error in Assign.cs source file csava5-Jan-04 5:28 csava 5-Jan-04 5:28
 Re: Error in Assign.cs source file meener7771-Jan-06 18:57 meener777 1-Jan-06 18:57
 Re: Error in Assign.cs source file JahBreeze26-Apr-06 13:03 JahBreeze 26-Apr-06 13:03
 How i can import functions of matlab dll in C#? ifafa7-Nov-03 8:35 ifafa 7-Nov-03 8:35
 Re: How i can import functions of matlab dll in C#? Emanuele Ruffaldi7-Nov-03 10:32 Emanuele Ruffaldi 7-Nov-03 10:32
 Re: How i can import functions of matlab dll in C#? Emanuele Ruffaldi23-Nov-03 3:04 Emanuele Ruffaldi 23-Nov-03 3:04
 Solution from matlab's website FinnHawk3-Nov-05 7:22 FinnHawk 3-Nov-05 7:22
 33 Anonymous8-Oct-03 1:53 Anonymous 8-Oct-03 1:53
 very impressed! but missing main project Frank Hileman12-Dec-02 13:09 Frank Hileman 12-Dec-02 13:09
 Re: very impressed! but missing main project pititaly14-Dec-02 5:25 pititaly 14-Dec-02 5:25
 Re: very impressed! but missing main project DaveOH19-Nov-03 18:43 DaveOH 19-Nov-03 18:43
 Stupid question? (Still having math at school...) Lars [Large] Werner12-Dec-02 7:53 Lars [Large] Werner 12-Dec-02 7:53
 Re: Stupid question? (Still having math at school...) pititaly12-Dec-02 8:18 pititaly 12-Dec-02 8:18
 Re: Stupid question? (Still having math at school...) Lars [Large] Werner12-Dec-02 10:40 Lars [Large] Werner 12-Dec-02 10:40
 Re: Stupid question? (Still having math at school...) pititaly12-Dec-02 11:22 pititaly 12-Dec-02 11:22
 Re: Stupid question? (Still having math at school...) Lars [Large] Werner12-Dec-02 23:57 Lars [Large] Werner 12-Dec-02 23:57
 Re: Stupid question? (Still having math at school...) Justin Turney12-Dec-02 15:33 Justin Turney 12-Dec-02 15:33
 Re: Stupid question? (Still having math at school...) Joaquín M López Muñoz12-Dec-02 20:25 Joaquín M López Muñoz 12-Dec-02 20:25
 Re: Stupid question? (Still having math at school...) Daniel Turini12-Dec-02 20:29 Daniel Turini 12-Dec-02 20:29
 Re: Stupid question? (Still having math at school...) Ikke14-Dec-02 2:06 Ikke 14-Dec-02 2:06
 Re: Stupid question? (Still having math at school...) Lars [Large] Werner14-Dec-02 4:27 Lars [Large] Werner 14-Dec-02 4:27
 Last Visit: 31-Dec-99 18:00     Last Update: 31-Aug-15 21:22 Refresh 12 Next »