11,925,292 members (59,846 online)
Add your own
alternative version

187K views
14.2K downloads
66 bookmarked
Posted

# C# Matrix Library

, 28 Nov 2007 CPOL
 Rate this:
Please Sign up or sign in to vote.
A C# library for basic numerical linear algebra.

## Introduction

CSML - C# Matrix Library - is a compact and lightweight package for numerical linear algebra. Many matrix operations known from Matlab, Scilab and Co. are implemented.

### Version Info

Version 0.91 - Last update: November 29, 2007.

### Remark

Make sure to return to this article once in a while for updates. A project of this size a is big thing for one man to handle. Scilab (free Matlab clone), for instance, has been created by an academic consortium, and Matlab's creator, Mathworks, is a well-fed enterprise.

Bug fixes will always come as updates, and if you have a look at the history below, you will see there have been quit a number of, ehm, improvements in a short period of time.

## What it can do

The core of the library, the `Matrix` class, includes over 90 methods for matrix operations such as multiplication, summation, exponentiation, and solving linear equation systems; matrix manipulations such as concatenation, insertion, transpose, inverse, flipping, symmetrizing, insertion, and extraction; for matrix computations such as determinant, trace, permanent, norm (Frobenius, Euclidian, maximum norm, taxi norm, p-norm), condition number; for matrix decompositions such as LU factorization, Gram-Schmidtian orthogonalization, and Cholesky factorization.

Now the entire library has been updated to work with complex arithmetic. A real matrix M is to be considered a special complex matrix, where M.Re() == M, e.i., the imaginary part of each entry equals 0.

The project is nearly entirely PDF-documented; most methods are also illustrated with examples. Difficult algorithms like the computation of the inverse are explained on a mathematical level as well. If vitally necessary, complexity classes of certain algorithms are noted.

## What it cannot do

At this time, there is no implementation of

• polynomials, interpolation, integration (*);
• eigenvectors (*);

(*) - Working on it. Any contribution is welcome. A complex numbers library has been issued by me here on CodeProject, look for CompLib. A library for polynomials has been released as well, PolyLib.

## How to use CSML

Two general ways for using the code:

1. Include CSML.dll as a reference to your project and type "`using CSML;`" on top of your source file
2. Include classes Matrix.cs and Complex.cs to your project and adjust namespaces

Let us see an example. Say, we want to compute the determinant and the inverse of the 2 by 2 matrix [1, 1; 1, 2]:

```Matrix M = new Matrix("1,1;1,2"); // init
Complex det = M.Determinant(); // det = 1
Matrix Minv = M.Inverse(); // Minv = [2, -1; -1, 1]
string buf = M.ToString(); // for outputting M in a multiline textbox```

For details of implementation and usage, refer to the documentation.

## Points of interest

This project is issued without license and warranty, it should be considered a gift to the developer's community in general and to this page in particular - most of my programming knowledge is based on free code, on examples and tutorials submitted without the greed for money. I am now in the happy position not having to turn anything into bucks: this is the result.

## A word about speed issues

None of the algorithms presented here is optimized for speed. The idea behind this project is to demonstrate many common and well-known algorithms for matrices in a straightforward and understandable way.

Possible optimization ideas are:

• port the methods you need to C++
• many small methods should be declared inline
• use double arrays instead of a matrix data type; C++ allows for dimension extensions during runtime via pointer arithmetic

## A word about Eigenvalues

The `Eigenvalues()` method in its current state uses basic QR iteration based upon Gram-Schmidtian orthogonalization. This implies two things:

1. for now, computation of Eigenvalues is possible only for real matrices
2. if a matrix has multiple Eigenvalues or complex Eigenvalues, the method works, but junk may be returned

In fact, I had thoroughly satisfying results only for triangular matrices and symmetric positive definite matrices.

These problems mirror the difficulties buried under the Eigenvalue problem. Since the Eigenvalues of a matrix A are defined as the roots of the characteristic polynomial:

`p(L) = det(A-L*id)`

computation is mathematically equivalent to the computation of a determinant and the n roots of p. (Well, this would work for all matrices with any distribution of Eigenvalues, but it is numerically the worst idea, since Weierstrass iteration (compare the `Roots()` method in PolyLib) is badly conditioned for polynomials with roots not being well-separated.)

Therefore, I am working on a canonical double-shift QR iteration based upon Givens rotations. That is the way Matlab's function `eig` works.

That I am forced to talk at length about Eigenvalues, although there are so many other difficult computations implemented, reveals to me that this problem is one of the deepest and most bothersome in basic numerical linear algebra.

## History

Conjugate gradient method (`SolveCG()` is implemented but buggy); I'm going to issue an example project showing the usage and basic functionality of CSML one of these days.

### Update November 29, 2007

• Bug in matrix addition/subtraction (found by Raven123) and bug in horizontal/vertical concatenation fixed (found by alexei_s).

### Update July 4, 2007

• Added Hessenberg-Householder reduction, Householder reflection, and Givens rotation; `HessenbergHouseholder()` works nicely; `QRGivens()` and `QRIterationHessenberg()` are still buggy.
• Alightly modified constructors for vector initialization.
• Bug in `InverseLeverrier()` fixed; this black box method should now be standard for matrix inversion (`Inverse()` is much slower for general matrices, but fast for special matrices being orthogonal, unitary, or diagonal).
• `ToString(string format)` method in Complex.cs beautified; multiplication of double and complex values slightly changed (no bugs there, but inconsistencies).

### Update July 3, 2007 #2

• Major bug in `Arg()` fixed (thanks Petr Stanislav!); this affected `Log()`, `Pow()`, and `Sqrt()`.

### Update July 3, 2007

1. Minor bug within the matrix access routine fixed, making two `try-catch` blocks obsolete and increasing speed enormously.
2. Defined matrix exponentiation with negative exponent k as exponentiation of the inverse with (-k); (3) major bug within the `Insert()` method fixed.
3. Old matrix class (without complex number support) deleted.

### Update July 2, 2007

• A rainy Monday in Bavaria... Instead of learning for my exams, I updated and finished the documentation and added the hyperbolic functions `Sinh`, `Cosh`, `Tanh`, `Coth`, `Sech`, `Csch` to the `Complex` class.
• `SolveSafe` fixed, renamed to `Solve` and made static. This method is now nearly equivalent to the Matlab backslash operator. Fixation of `SolveSafe` (using LU decomposition with column pivoting) suddenly made the old `Solve` method superfluous, which has been removed.
• `Size` method removed and replaced with read-only properties `RowCount` and `ColumnCount`.
• `IsVector` method renamed with `VectorLength` (tests if a matrix is a vector, e.i. n by 1 or 1 by n).

### Update June 29, 2007

• Finally complex numbers and complex matrices; any matrix consists of complex entries now.

New in this version:

• Major bug in `IsPermutation()` removed (now: permutation matrices are precisely the involuntary 0-1 matrices; refer to doc).
• Tests for complex matrices implemented: `IsUnitary()`, `IsReal()`, `IsImaginary()`, `IsHermitian()`.
• Test added: `IsZeroOneMatrix()`.
• `KroneckerDelta()` implemented.
• Gram-Schmidtian orthogonalization.
• QR iteration based on Gram-Schmidt.
• Function for Eigenvalues (slow!).
• General method `Insert()` to insert a submatrix at a certain position.
• General method `Extract()` to copy a submatrix from a given matrix.
• `Conjugate()`, `ConjugateTranpose()`.
• `BlockMatrix()` to build a matrix from four given matrices.
• `ChessboardMatrix()` to construct matrices with interchanging 0-1 entries.
• `RandomZeroOne()`
• `Random()` extended to generate not just double, but also integer matrices.
• `Re()`, `Im()` to extract the real/imaginary parts of a matrix.
• `SolveSafe()` marked as buggy.
• `Eigenvector()`, `SolveCG()`, `QRHessenbergHouseholder()`, `HouseholderVector()` implemented but yet (correctly) marked as buggy.
• Documentation updated (but not finished).
• Faster method for matrix inversion: `InverseLeverrier()`, using the Leverrier method (RTFM).

### Update June 10, 2007

• Major multiplication bug removed; in "`Matrix operator *(Matrix A, double x)`".
• Test for definiteness added (works for symmetric matrices only); needs beta-testing!
• Test for symmetric positive definiteness (SPD) added and used within `Cholesky()`; Cholesky decomposition is possible only for SPD matrices.
• Undocumented basic graph algorithms: BFS (broad-first search), DFS (depth-first search), and Floyd (shortest paths between all vertices of a graph).

## License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

## About the Author

 Germany
No Biography provided

## Comments and Discussions

 First PrevNext
 How to Compile CSML with Another Project Member 17199711-Sep-15 6:45 Member 1719971 1-Sep-15 6:45
 operator [-] Benn Eifert3-Jul-15 12:49 Benn Eifert 3-Jul-15 12:49
 Matrix left division Marek Bar12-Jul-14 23:23 Marek Bar 12-Jul-14 23:23
 Great library, but code duplication. KommuSoft19-Jun-14 4:46 KommuSoft 19-Jun-14 4:46
 The cofactor function of a matrix Babagana Zannah7-May-14 3:53 Babagana Zannah 7-May-14 3:53
 Possible Bug in HorizontalConcat Bedir Yilmaz27-Apr-14 10:47 Bedir Yilmaz 27-Apr-14 10:47
 Constructor (String to Complex) wieschoo5-Apr-13 9:34 wieschoo 5-Apr-13 9:34
 Big bug found for matrix operator Jaxam8-Mar-13 8:11 Jaxam 8-Mar-13 8:11
 Re: Big bug found for matrix operator Jaxam8-Mar-13 9:53 Jaxam 8-Mar-13 9:53
 My vote of 5 eshark1119-Oct-11 14:06 eshark11 19-Oct-11 14:06
 My vote of 1 tonatoz17-Mar-11 23:46 tonatoz 17-Mar-11 23:46
 Re: My vote of 1 Member 778524129-Mar-11 4:02 Member 7785241 29-Mar-11 4:02
 Matrix Multiplication nikolay-202-Mar-11 16:22 nikolay-20 2-Mar-11 16:22
 Re: Matrix Multiplication Member 778524129-Mar-11 4:03 Member 7785241 29-Mar-11 4:03
 Very Useful jorge.aqp.2323-Oct-10 9:22 jorge.aqp.23 23-Oct-10 9:22
 хуйня gizmonder14-Sep-10 14:24 gizmonder 14-Sep-10 14:24
 My vote of 1 gizmonder14-Sep-10 14:24 gizmonder 14-Sep-10 14:24
 Solution to Bug affecting Inverse/Solve/Determinant- easy code fix eshark119-Aug-10 19:16 eshark11 9-Aug-10 19:16
 Great library, Hanzzoid. Hope you return to this project some day. In the meantime there is a bug which causes a "Warning: Matrix close to singular." exception, which we found during matrix.Inverse() and would also affect Solve and Determinant. The solution is easy, and requires moving two lines down a little. All users of this library should do this. 1. Open the Matrix.cs file 2. Find this exception at line 1935 in LUSafe() ```if (this[j, j] == 0) throw new DivideByZeroException("Warning: Matrix close to singular."); ``` 3. Move those two lines down by 14 lines so they are inside the loop at 1949 ```for (int i = j + 1; i <= n; i++) { this[i, j] = this[i, j] / this[j, j]; } ``` 4. The result should look like this ```for (int i = j + 1; i <= n; i++) { if (this[j, j] == 0) throw new DivideByZeroException("Warning: Matrix close to singular."); this[i, j] = this[i, j] / this[j, j]; } ``` This is for the code version 0.9 (Last Update: July 4, 2007) Hope this helps someone, Regards, David Pierson
 Compatibility Raymas9-Aug-10 18:39 Raymas 9-Aug-10 18:39
 Re: Compatibility eshark119-Aug-10 19:20 eshark11 9-Aug-10 19:20
 Re: Compatibility Raymas13-Aug-10 2:39 Raymas 13-Aug-10 2:39
 My vote of 5 InvaderXive30-Jun-10 14:48 InvaderXive 30-Jun-10 14:48
 Retrieving double element from Matrix [modified] holyfetzer23-May-10 1:32 holyfetzer 23-May-10 1:32
 Re: Retrieving double element from Matrix InvaderXive30-Jun-10 14:50 InvaderXive 30-Jun-10 14:50
 accessing matrix values! mingas28-Feb-10 9:50 mingas 28-Feb-10 9:50
 Re: accessing matrix values! Member 778524128-Mar-11 5:34 Member 7785241 28-Mar-11 5:34
 retrieval of a specific value from the matrix Innayat Ullah8-Feb-10 3:43 Innayat Ullah 8-Feb-10 3:43
 Re: retrieval of a specific value from the matrix yopieadrianto26-Feb-10 19:12 yopieadrianto 26-Feb-10 19:12
 Latest developments cornelius172929-Oct-09 8:04 cornelius1729 29-Oct-09 8:04
 Bug in Matrix Inverse Rajnikant Sharma BYU17-May-09 6:35 Rajnikant Sharma BYU 17-May-09 6:35
 Re: Bug in Matrix Inverse safee ullah2-Sep-09 18:32 safee ullah 2-Sep-09 18:32
 Re: Bug in Matrix Inverse emre unsal17-Nov-09 4:11 emre unsal 17-Nov-09 4:11
 Re: Bug in Matrix Inverse eshark119-Aug-10 19:22 eshark11 9-Aug-10 19:22
 3 warrnings detected chen dingjun17-Dec-08 13:15 chen dingjun 17-Dec-08 13:15
 HorizontalConcat seems to have a bug Matthias Zeintlinger3-Dec-08 6:50 Matthias Zeintlinger 3-Dec-08 6:50
 Re: HorizontalConcat seems to have a bug Håkan MacLean5-Apr-09 6:56 Håkan MacLean 5-Apr-09 6:56
 bug in Inverse()-Method arthur_dent428-Apr-08 7:35 arthur_dent42 8-Apr-08 7:35
 Inverse Matrix of Array get13-Dec-07 2:45 get 13-Dec-07 2:45
 strange operator consequences? Raven12327-Nov-07 6:00 Raven123 27-Nov-07 6:00
 Re: strange operator consequences? hanzzoid29-Nov-07 0:05 hanzzoid 29-Nov-07 0:05
 Re: strange operator consequences? yopieadrianto26-Feb-10 19:16 yopieadrianto 26-Feb-10 19:16
 Great zhengdong jin23-Nov-07 5:18 zhengdong jin 23-Nov-07 5:18
 Re: Great Rajnikant Sharma BYU28-Dec-08 15:56 Rajnikant Sharma BYU 28-Dec-08 15:56
 HorizontalConcat() and VerticalConcat() do not work alexei_c21-Oct-07 14:51 alexei_c 21-Oct-07 14:51
 Re: HorizontalConcat() and VerticalConcat() do not work hanzzoid29-Nov-07 0:24 hanzzoid 29-Nov-07 0:24
 Inverse of m x n matrix? joselvra29-Sep-07 13:56 joselvra 29-Sep-07 13:56
 Re: Inverse of m x n matrix? mvasiliev18-Oct-07 22:42 mvasiliev 18-Oct-07 22:42
 Re: Inverse of m x n matrix? hanzzoid29-Nov-07 0:09 hanzzoid 29-Nov-07 0:09
 exception thrown slady.net28-Aug-07 18:05 slady.net 28-Aug-07 18:05
 operations matiasjaure28-Aug-07 4:55 matiasjaure 28-Aug-07 4:55
 Last Visit: 31-Dec-99 19:00     Last Update: 26-Nov-15 16:20 Refresh 12 Next »

General    News    Suggestion    Question    Bug    Answer    Joke    Praise    Rant    Admin

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Terms of Use | Mobile
Web03 | 2.8.151126.1 | Last Updated 29 Nov 2007
Article Copyright 2007 by hanzzoid
Everything else Copyright © CodeProject, 1999-2015
Layout: fixed | fluid