CodeProject: DotNetMatrix: Simple Matrix Library for .NET. Free source code and programming help

The classes provided with this article give you a basic linear algebra package for .NET. It provides user-level C# classes for constructing and manipulating real, dense matrices. It is meant to provide sufficient functionality for routine problems, packaged in a way that is natural and understandable to non-experts. That said, it is simply a port of a public domain Java matrix library, called JAMA. Check the JAMA web site for more information.

Background

Currently, I work for a small GIS company here in Japan developing GIS components for application developers. Coordinate Transformation and therefore Affine Transformation is a very basic part of our development efforts. Recently, I was assigned a task of designing and implementing a new GIS system for the .NET framework with the ability to easily port to Java and other frameworks. I decided to make maximum use of matrix-based affine transformation, which is also a requirement of the OpenGIS Coordinate Transformation Specifications.

Then, I discovered that the Matrix class provided as part of the GDI+ in the .NET implements the affine transformations in a manner different from standard specifications. In short, while the standard 2D Coordinate System Affine Transformation Matrix is define as

The effect is that most affine transformations with the GDI+ Matrix class will not conform to standard or specifications. For instance, by standard (and mathematically) anti-clockwise (or counter-clockwise) rotations are considered positive but must be negative when using the GDI+ classes.
To solve this problem, I decided to implement a standard affine transformation matrix and found this small Java general matrix library - JAMA. The classes presented with this article are ported from the JAMA with .NET specific improvements like operator overloading etc.

In a different article, I will discuss and present my affine transformation matrix, which should give the same effect (same properties and methods as the GDI+ Matrix class). The rest of the article is extracted from the JAMA documentation, and I will refer to the library as DotNetMatrix (the namespace).

Capabilities

DotNetMatrix is comprised of six C# classes: GeneralMatrix, CholeskyDecomposition, LUDecomposition, QRDecomposition, SingularValueDecomposition and EigenvalueDecomposition.

The GeneralMatrix class provides the fundamental operations of numerical linear algebra. Various constructors create Matrices from two dimensional arrays of double precision floating point numbers. Various gets and sets (properties) provide access to submatrices and matrix elements. The basic arithmetic operations include matrix addition and multiplication, matrix norms and selected element-by-element array operations. A convenient matrix print method is also included.

Five fundamental matrix decompositions, which consist of pairs or triples of matrices, permutation vectors, and the like, produce results in five decomposition classes. These decompositions are accessed by the GeneralMatrix class to compute solutions of simultaneous linear equations, determinants, inverses and other matrix functions. The five decompositions are

The design of DotNetMatrix represents a compromise between the need for pure and elegant object-oriented design and the need to enable high performance implementations.

Using the Code

*Summary of DotNetMatrix Capabilities*
Object Manipulation	constructors set elements get elements copy clone
Elementary Operations	addition subtraction multiplication scalar multiplication element-wise multiplication element-wise division unary minus transpose norm
Decompositions	Cholesky LU QR SVD symmetric eigenvalue nonsymmetric eigenvalue
Equation Solution	nonsingular systems least squares
Derived Quantities	condition number determinant rank inverse pseudoinverse

The following simple example solves a 3x3 linear system Ax=b and computes the norm of the residual.

Points of Interest

There is an empty ISerializable interface implementation in all the classes. I do not know if serialization support is really needed, I could complete the implementations if needed.

License

This is basically a port of public domain codes, so the result is in the public domain.

History

Conclusion

Matrices are valuable mathematical tools which can be used to solve very complex problems, and coordinate transformation is no exception. The classes presented here can help you implement your own affine transformation classes to get a bit far than the .NET framework offers you. Also, since the classes are completely C# unlike the GDI+ classes, which are wrappers for a COM implementation, there is no interop overhead.

In my next article, I will present an affine transformation matrix class based on this library. Suggestions for improvements are welcomed.

You must Sign In to use this message board.


Per page

FirstPrevNext

Delet Rows or Colums ?

muhaa222

9:43 15 Jun '09

Hi I use your Libary. it's very very nice code but i have still a problem i want to delet da row or a col ? i didn#t have found a soulution for this problem cya muhaa
Sign In·View Thread·PermaLink

How to solve a system of 3 linear equations in 3 vars, x,y,z

andeezlemoyo

13:35 2 Dec '08

Hi Paul, first, thank you for sharing your work. Here's a sample linear system I want to solve: 49x - 7y + z = 6 25x - 5y + z = 4 9x - 3y + z = 2 How can I solve this with GeneralMatrix? I would greatly appreciate it if you could give me some sample code. Thanks, Andy
Sign In·View Thread·PermaLink

Re: How to solve a system of 3 linear equations in 3 vars, x,y,z

muhaa222

9:37 15 Jun '09

Hi I think the solve is GeneralMatrix A like 49 -7 1 25 -5 1 9 -3 1 and GeneralMatrix B like 6 4 2 GeneralMatrix N=A.Transpose()A; GeneralMatrix n=A.Transpose()B; GeneralMatrix x=N.Inverse()*n x.GetElement(0,0) is your x x.GetElement(1,0) is your y x.GetElement(2,0) is your z i use this for solving overdetermined equation but is still work for your problem cya muhaa
Sign In·View Thread·PermaLink

QR Decomposition - Possible Bug

LachlanGro

0:40 10 Nov '08

Hi,

I have a 3x3 matrix, V:
V =
100 10 1
625 25 1
900 30 1

In MatLab I do the QR and get:
>> [Q,R] = qr(V,0)
Q =
-0.0909 0.8789 0.4683
-0.5680 0.3405 -0.7493
-0.8180 -0.3341 0.4683
R =
1.0e+003 *
-1.1003 -0.0396 -0.0015
0 0.0073 0.0009
0 0 0.0002

However when using the DotNetMatrix I get:
Matrix Q =
-0.0908856214131766 0.878886544558686 -0.468292905790847
-0.568035133832354 0.340488733969466 0.749268649265355
-0.817970592718589 -0.334104570207539 -0.468292905790847

While the values are good, the right hand column is all a negative of the MatLab results.

Any suggestions?

Regards
Lachlan.

Re: QR Decomposition - Possible Bug

Paul Selormey

16:13 10 Nov '08

LachlanGro wrote:
While the values are good, the right hand column is all a negative of the MatLab results.

True, I have verified it. I have also used the example matrix in the wikipedia http://en.wikipedia.org/wiki/QR_decomposition[^] and the signs do not correspond.

LachlanGro wrote:
Any suggestions?

I am currently looking into it, I do not have any immediate idea of the source of the problem.

Best regards,
Paul.

Jesus Christ is LOVE! Please tell somebody.

Re: QR Decomposition - Possible Bug

LachlanGro

16:48 10 Nov '08

Hi Paul, I posted a similar question on the Jama discussion list and apparently there is no bug, its just that the solution is "not unique". While I believe this is the case, I could not find an indication of what would be unique. Either way I pass exactly the same data into other libraries (such as BlueBit, Matlab) and the results are fine. Regards, Lachlan.
Sign In·View Thread·PermaLink	5.00/5

Re: QR Decomposition - Possible Bug

Paul Selormey

18:08 10 Nov '08

LachlanGro wrote:
I posted a similar question on the Jama discussion list and apparently there is no bug, its just that the solution is "not unique".

Thanks for the information.

LachlanGro wrote:
While I believe this is the case, I could not find an indication of what would be unique.

Computations based on the results are valid, also the Q'Q = I, so that may be the case.

LachlanGro wrote:
Either way I pass exactly the same data into other libraries (such as BlueBit, Matlab) and the results are fine.

Currently, I do not have any math package so I was only looking for books and online examples to compare. I will look into the algorithm used to see if we could get the same results as Matlab to prevent any doubt.

Best regards,
Paul.

Jesus Christ is LOVE! Please tell somebody.

5.00/5

Bad Matrix class

Kuryn

18:35 6 Aug '08

Why make a class virtual? That causes a performance issue. You should never make a class virtual.
Sign In·View Thread·PermaLink

Re: Bad Matrix class

Paul Selormey

19:03 6 Aug '08

Kuryn wrote: You should never make a class virtual. Hmmm, I do not know about this new rule in town. However, you have the source codes, so please do whatever you like with it - it is a public domain license. Best regards, Paul. Jesus Christ is LOVE! Please tell somebody.
Sign In·View Thread·PermaLink	5.00/5

Re: Bad Matrix class

Kuryn

6:24 9 Aug '08

Well you should make a class virtual when needed. But not a mathematical class, because nothing changes different. Like in a Vector, a dotProduct is a dotProduct.
Sign In·View Thread·PermaLink	2.00/5

Pseudoinverse

FueledByIgnorance

22:32 7 Jul '08

Thanks for this library, it is very useful!

I have a question: I arrive at some point with a Matrix that is singular, and I want to perform an Inverse on it.

GeneralMatrix matrixInvertQ = matrixMultiplyN.Inverse();
...

public virtual Matrix Inverse()
{
     return Solve(Identity(m, m));
}

The Inverse method calls Solve, which instantiates a LUDecomposition object and calls its Solve method. There it tests for singularity, and I get the "Matrix is singular" exception.

I have very limited knowledge in mathematics, so bear with me. At the point where the matrix is considered singular, shouldn't a pseudoinverse by performed?
In the comments for the Inverse method, it says "inverse(A) if A is square, pseudoinverse otherwise".
Is there a part missing from implementation?

Thanks!

5.00/5

Re: Pseudoinverse

Paul Selormey

7:34 10 Jul '08

Hello, Singularity and pseudo-inverse are different. To understand the pseudo-inverse, see the following 1. In PlanetMath[^] 2. Example[^] Best regards, Paul. Jesus Christ is LOVE! Please tell somebody.
Sign In·View Thread·PermaLink

Bug in Singular Value Decomposition

ajwo

9:48 19 Nov '07

Hi:

There is a bug in the svd decomposition, inherited from Jama. The problem is with the V matrix, when n>m. The correction is:


if (wantv)
{
    for (int k = n - 1; k >= 0; k--)
    {
        if ((k < nrt) & (e[k] != 0.0))
	{
	    //for (int j = k + 1; j < nu; j++) error
            for (int j = k + 1; j < n; j++) //correction
            {
	       double t = 0;

The bug is reported here:
http://cio.nist.gov/esd/emaildir/lists/jama/msg00430.html

Re: Bug in Singular Value Decomposition

Paul Selormey

11:02 19 Nov '07

Thank you, I will take a look. Best regards, Paul. Jesus Christ is LOVE! Please tell somebody.
Sign In·View Thread·PermaLink

Re: Bug in Singular Value Decomposition

iotha

21:42 3 Dec '07

Thank you!

I'm no math guy but I have a question upon another thing (may be a "bug"?):

The SVD is A[mxn] = U[mxm]*S[mxn]*V[nxn]^T - so U should be a square matrix. But it is defined here as:


			m = Arg.RowDimension;
			n = Arg.ColumnDimension;
			int nu = System.Math.Min(m, n);
			U = new double[m][];
			for (int i = 0; i < m; i++)
			{
				U[i] = new double[nu];
			}

Shouldn't it be like this:

			m = Arg.RowDimension;
			n = Arg.ColumnDimension;
			int nu = System.Math.Min(m, n);
			s = new double[System.Math.Min(m + 1, n)];
			U = new double[m][];
			for (int i = 0; i < m; i++)
			{
				U[i] = new double[m];
			}

Re: Bug in Singular Value Decomposition

celisdelorena

5:42 27 Jun '08

Hello everybody, I am new using matrix in C# and I would like to thank and congratulate to the author of General Matrix. This is a very usefull and easy class to work with it. Currently I am working with the SVD descomposition to obtain the PseudoInverse... But it doesn't work. The code blocks during the SVD calculation. Does anybody have the same problem? Thanks in advance
Sign In·View Thread·PermaLink	2.67/5

Re: Bug in Singular Value Decomposition when solving under-determined equation systems [modified]

raz0rblade

4:52 29 Aug '08

Hi there,
I encountered two problems when solving under-determined equation systems:
The first problem is an ArrayIndexOutOfBoundsException when trying to get the S-Matrix. I think the problem is the wrong size of S, because of A[mxn] = U[mxm]*S[mxn]*V[nxn]' the dimensions should be m-by-n but in the code they are n-by-n:

/// Return the diagonal matrix of singular values

///      S
/// 
virtual public GeneralMatrix S
{
	get
	{
		GeneralMatrix X = new GeneralMatrix(n, n); // this should be (m, n)
		double[][] S = X.Array;
		for (int i = 0; i < n; i++) // this loop should go till m not n
		{
			for (int j = 0; j < n; j++)
			{
				S[i][j] = 0.0;
			}
			S[i][i] = this.s[i];
		}
		return X;
	}
}

Then there is a problem with the dimension of the singular values vector in the constructor of SingularValueDecomposition:

s = new double[System.Math.Min(m + 1, n)];

For example, my equation system has six equations and eight variables, so m=6 and n=8. Then s will be the size of seven and there will be a singular value that doesn't really exists.

In addition to this, I get different results when using the svd-function of Matlab. The singular values and the U-matrix are exactly the same, but in the V-matrix the last n-m coloumns have different values, the others have the same like in Matlab. But I didn't find a suitable place in the code for that problem.

Please correct me, if I understood something mathematically wrong. Smile

modified on Friday, August 29, 2008 10:16 AM

Re: Bug in Singular Value Decomposition

mbreidt

5:16 9 Apr '09

Not sure this discussion board is still active, but I recently discovered this library and found it extremely useful, so first of all thanks very much for making it available!

Unfortunately, I have also discovered problems with the SVD computation, and wrong dimensions of the resulting matrices.

If you do a test with the example provided on Wikipedia[^], this library produces wrong results as far as I can tell (FWIW, Matlab returns the same values given on Wikipedia).

After running the example (an SVD on a 4-by-5 matrix), the S matrix from the SVD object is a 5-by-5 matrix, whereas the correct result should be a 4-by-5 matrix. Also, the entries of the U and V matrix seem to be in the wrong order and/or have the wrong sign.

Also, the singular values vector returned by the SVD object has five elements in this example, but, according to Wikipedia

"An m × n matrix M has at least one and at most p = min(m,n) distinct singular values."

so there should be a maximum of four singular values, not five.

If anyone is still following this board, I'd very much like to hear whether anyone can confirm my findings and/or has a bug fix!
I'd really like to use this library for my work, but only if works with non-square matrices.

Thanks!
Martin

Problem with eigenvalues

Cukaricn

7:07 30 Oct '07

I have problem with relation between eigenvalues and eigenvectors. My eigenvectors do not correspond with their eigenvalues. Instead they are mixed (for exaple eigenvalue a1 has vector v3, instead of v1). All eigenvalues i have are negative so can that be problem? Thanks Cukaric
Sign In·View Thread·PermaLink

Karhunen-Loève

sergiomad

15:32 25 Jun '07

Is available a method for getting the Karhunen-Loève vectors? I'm trying to translate this code in Java http://astro.u-strasbg.fr/~fmurtagh/mda-sw/java/PCAcorr.java[^], that uses JAMA, but I can't find methods like "times()" Thank you
Sign In·View Thread·PermaLink

Re: Karhunen-Loève

Paul Selormey

18:37 25 Jun '07

sergiomad wrote:
Is available a method for getting the Karhunen-Loève vectors?

No, Karhunen-Loève is not directly related to matrix.

sergiomad wrote:
but I can't find methods like "times()"

This port uses the .NET naming, it should be "Multiply(...)"

Best regards,
Paul.

Jesus Christ is LOVE! Please tell somebody.

Is it a bug in the constructor of GeneralMatrix?

danielmark

11:09 23 May '07

public GeneralMatrix(double[] vals, int m)
		{
			this.m = m;
			n = (m != 0?vals.Length / m:0);
			if (m * n != vals.Length)
			{
				throw new System.ArgumentException("Array length must be a multiple of m.");
			}
			A = new double[m][];
			for (int i = 0; i < m; i++)
			{
				A[i] = new double[n];
			}
			for (int i = 0; i < m; i++)
			{
				for (int j = 0; j < n; j++)
				{
					A[i][j] = vals[i + j * m];
				}
			}
		}

I think the statement

A[i][j] = vals[i + j * m];

should be revised as

A[i][j] = vals[j + i * m];

Re: Is it a bug in the constructor of GeneralMatrix?

Paul Selormey

19:16 25 Jun '07

Sorry for the late response.
Well, not a bug. If you see the documentation of that constructor, vals is said to be packed by columns.

I think you needed the "packed by row" array filling, I kept the original implementations. I will update this library soon and give the user the options
to choose between the two with an extra parameter to avoid confusion.

Best regards,
Paul.

-- modified at 11:32 Tuesday 17th July, 2007

Jesus Christ is LOVE! Please tell somebody.

Any sample code for a newbie

tgpritchard

15:00 28 Apr '07

Paul,
I am new to DotNet programming and I'm working on a project to determine whether a series of values over time are trending upward, down or steady. I'm working with performance data from the Windows OS.
I believe the simplest approach will be to run a least squares analysis on the data (40 points, each from a successive day). I would like to use your DotNetMatrix library but am a little out of my league. Do you have any script samples/instructions for implementing the library and especially for a least squares analyis.

Thanks for your help,
Tom

P.S. Love your message postscript, my brother.

1.00/5

Re: Any sample code for a newbie

lukner

13:46 20 Jun '07

Tom,

You have many options for least squares, although I'm unsure how to use this particular library for the job (it didn't do the QR decomposition as I expected ... I never got the expected Q matrix for some reason).

Anyway, here is a least squares solution using QR decomposition (sometimes called QR Factorization). Feel free to look up least squares on Wiki as well. Also, I give links to a much easier way of doing least squares at the end of this message.

http://www.alkires.com/teaching/ee103/Rec8_LLSAndQRFactorization.htm[^]

You can get QR decomposition and linear equation solvers which actually do work (they seem to be based on LAPACK) from here:

http://www.alglib.net/matrixops/general/qr.php routine, which seems to have come from LAPACK CGEQR

and

http://www.alglib.net/linequations/general/lu.php

Both of the links above have code samples and test code that seems to work.

**** However, you may want something much simpler, such as ... ****

http://www.codeguru.com/forum/showpost.php?p=1139483&postcount=2[^]

or

http://www.pdc.kth.se/training/Tutor/MPI/Pt2pt.Lab/least-squares.c[^] (some adaptation required since this is C code, but it should be easy to adapt).

All the best,
--Ralf

Last Visit: 19:07 21 Nov '09 Last Update: 19:07 21 Nov '09

12 3 4 Next »

Introduction

Background

Capabilities

Using the Code

Points of Interest

License

History

Conclusion