Add your own alternative version
Stats
337.8K views 10.8K downloads 120 bookmarked
Posted
12 Jan 2004

Comments and Discussions



The author has described this library as a port from the JAMA library. That library is stated by NIST to be preliminary and is written in Java. After testing in the latest Java NetBeans IDE I have discovered that, unlike C#, Java does not recognize the double[,] (matrix) syntax. Consequently, the NIST/MatLab authors were stuck with the 'jagged array' construct, double[][]. The explanation below as to 'why double[][]' makes no sense to me at all, nor do I think the issue is any way significantly related to GC or performance. It's simply that the author of this current C# library would have to have rewritten the whole port using double[,], a daunting chore to say the least.
Here are a couple of relevant code snippets that might help. The first is my own, the second from
Converting jagged array to 2D array C#  Stack Overflow[^]
<pre>
double[][] D3jj = { new double[] { 1.0, 4.0, 7.0, 10.0, 12.0, 14.0, 8.2 }, new double[] { 2.0, 5.0, 8.0, 11.0, 14.0 }, new double[] { 3.0, 6.0, 9.0, 12.0 } };
PrintJaggedArray(D3jj);
/* Output:
1 4 7 10 12 14 8.2
2 5 8 11 14
3 6 9 12
*/
static void PrintJaggedArray(double[][] array)
{
int nrows = array.GetLength(0);
for (int i = 0; i < nrows; i++)
{
for (int j = 0; j < array[i].GetLength(0); j++) // array[i].GetLength(0) gives i th row length
{
Console.Write(array[i][j] + " ");
} Console.WriteLine(" ");
} Console.WriteLine(" ");
}
static T[,] To2D<t>(T[][] source)
{
try
{
int FirstDim = source.Length;
int SecondDim = source.GroupBy(row => row.Length).Single().Key; // throws InvalidOperationException if source is not rectangular
var result = new T[FirstDim, SecondDim];
for (int i = 0; i < FirstDim; ++i)
for (int j = 0; j < SecondDim; ++j)
result[i, j] = source[i][j];
return result;
}
catch (InvalidOperationException)
{
throw new InvalidOperationException("The given jagged array is not rectangular.");
}
}
Now this is not an insignificant decision on the part of NIST to go with Java jagged array syntax. Had they selected C# in the first place, the implementation of complex matrices becomes almost trivial.
Using System.Numerics;
Complex[,] Cx = new Complex[3, 3];
Cx[0, 0] = new Complex(1.0, 1.0); Cx[0, 1] = new Complex(2.0, 1.0); Cx[0, 2] = new Complex(1.0, 1.0);
Cx[1, 0] = new Complex(2.0, 2.0); Cx[1, 1] = new Complex(1.0, 2.0); Cx[1, 2] = new Complex(2.0, 3.0);
Cx[2, 0] = new Complex(3.0, 1.0); Cx[2, 1] = new Complex(4.0, 1.0); Cx[2, 2] = new Complex(3.0, 1.0);
PrintComplexMatrix(Cx);
static void PrintComplexMatrix(Complex[,] Cx)
{
for (int i = 0; i < Cx.GetLength(0); i++)
{
for (int j = 0; j < Cx.GetLength(1); j++)
{
Console.Write(Cx[i, j] + " ");
} Console.WriteLine(" ");
} Console.WriteLine(" ");
}
modified 15Aug16 18:25pm.





For example, I have an array A declared as:
double[,] A = new double[6,3]
How do I create a GeneralMatrix from A?





Hi All, I am newbie here ,i really beginner in C# . So my problem is how can I use EigenvalueDecomposition class for finding the eigenvalue of my matrix I hope u can also provide the code example Thanx u
Peace!!





hOw aBout sUpporting cOmplex nUmbers iN tHe nEw vErsion? ?





I've been looking for something like this for a while now, and this seems to fit my needs wonderfully.





Any one can Help How to use DotNetMatrix to find LU Decomposition. And to return L, U,P . Please Help..





Hello, I have tried to use your sample code, but VS 2010 gives me this error:
"Array initializers can only be used in a variable or field initializer. Try using a new expression instead."
VS 2010 simply does not like the way the array is initialized:
double[][] array = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
I have tried to change it to:
double[,] array = {{1,2,3},{4,5,6},{7,8,10}};
which is correct way for VS 2010, but then GeneralMatrix(array) does not like this form.
How shall one initialize an array so that it will work with this probably awesome tool?
Many thanks for any hints.
hugo





I also encounter the very same problem with VS2003.





double[][] avals = {new double[]{1.0, 4.0, 7.0, 10.0}, new double[]{2.0, 5.0, 8.0, 11.0}, new double[]{3.0, 6.0, 9.0, 12.0}};
look at the examples in the sources and here http://msdn.microsoft.com/enus/library/sf3awfty(v=vs.71).aspx
Please correct this in the example on this page
modified 6Nov12 20:34pm.






I created two functions to get the Max and Min Values in a given matrix. If this is something you want to add to the project let me know and I will provide the code. Its just a couple of loops and a compare statement... Not sure how to contribute on Code Project, yet...





Hi Paul,
Great job. I just want to know if you have some updated version of theses codes?
Is there a C++ version available?
Thanks a lot.





Hi,
I used this dll really effectively in a research project. I like to make my work accessible to other researchers free. But with the dependency of Dll, I dont know whether I can do it or not?
Is there an copy right issue if I distribute your DLL with my dll?
Thank you
Tharindu





Hello Tharindu,
As you can see in the License section, this is a public domain code,
so you can anything you like with it.
You do not need to have it as a DLL, you can directly use the source codes in your projects.
The source code is very old now, you can clean it and keep anything changes you make, it is
a public domain stuff, just maintain any license notifications in the source code itself.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Thank you for the clarification.





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





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





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





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.





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.





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.





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.





The third column of Q is eigenvector, a unit direction vector, negative value is OK, just denotes the opposite direction.





Why make a class virtual? That causes a performance issue. You should never make a class virtual.





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.





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.





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!





Hello,
Singularity and pseudoinverse are different.
To understand the pseudoinverse, see the following
1. In PlanetMath[^]
2. Example[^]
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





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:
<code>
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;
</code>
The bug is reported here:
http://cio.nist.gov/esd/emaildir/lists/jama/msg00430.html





Thank you, I will take a look.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





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:
<br />
m = Arg.RowDimension;<br />
n = Arg.ColumnDimension;<br />
int nu = System.Math.Min(m, n);<br />
U = new double[m][];<br />
for (int i = 0; i < m; i++)<br />
{<br />
U[i] = new double[<big>nu</big>];<br />
}<br />
Shouldn't it be like this:
m = Arg.RowDimension;<br />
n = Arg.ColumnDimension;<br />
int nu = System.Math.Min(m, n);<br />
s = new double[System.Math.Min(m + 1, n)];<br />
U = new double[m][];<br />
for (int i = 0; i < m; i++)<br />
{<br />
U[i] = new double[<big>m</big>];<br />
}<br />





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





Hi there,
I encountered two problems when solving underdetermined equation systems:
The first problem is an ArrayIndexOutOfBoundsException when trying to get the SMatrix. I think the problem is the wrong size of S, because of A[mxn] = U[mxm]*S[mxn]*V[nxn]' the dimensions should be mbyn but in the code they are nbyn:
virtual public GeneralMatrix S<br />
{<br />
get<br />
{<br />
GeneralMatrix X = new GeneralMatrix(n, n);
double[][] S = X.Array;<br />
for (int i = 0; i < n; i++)
{<br />
for (int j = 0; j < n; j++)<br />
{<br />
S[i][j] = 0.0;<br />
}<br />
S[i][i] = this.s[i];<br />
}<br />
return X;<br />
}<br />
}
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)];<br />
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 svdfunction of Matlab. The singular values and the Umatrix are exactly the same, but in the Vmatrix the last nm 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.
modified on Friday, August 29, 2008 10:16 AM





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 4by5 matrix), the S matrix from the SVD object is a 5by5 matrix, whereas the correct result should be a 4by5 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 nonsquare matrices.
Thanks!
Martin





I do confirm you results my friend the sign of resulted U and V are messed up





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






sergiomad wrote: Is available a method for getting the KarhunenLoève vectors?
No, KarhunenLoè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.





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





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.





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.





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/leastsquares.c[^] (some adaptation required since this is C code, but it should be easy to adapt).
All the best,
Ralf





lukner wrote: (it didn't do the QR decomposition as I expected ... I never got the expected Q matrix for some reason).
I am currently trying to update this matrix, please can you give me some examples of the problems you have with the Qmatrix? (may be a matrix and expected Q).
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





(reply is in the Dynamic size thread below by accident ...





This library is great! But it seems that it's lack the capability of handling dynamic size matrix. Sometime we won't know what the matrix size would be as it would grow/shrink with time, so it will be better if the library can support dynamic size.





Currently you can shink by using the GetMatrix method, which allows you to retrieve a submatrix.
I will add a general method, say Resize(int newRows, int newColumns), to support your request in the next version.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Hi Paul,
Wiki has an example you can use:
http://en.wikipedia.org/wiki/QR_decomposition
I used the following:
A =
0 1.0000
0.6931 1.0000
1.3863 1.0000
2.0794 1.0000
2.7726 1.0000
3.4657 1.0000
Correct answer from MATLAB:
>> [Q,R] = qr(A)
Q =
0 0.7237 0.3248 0.3381 0.3513 0.3646
0.1348 0.5264 0.1776 0.1156 0.4088 0.7020
0.2697 0.3290 0.8999 0.0761 0.0520 0.0280
0.4045 0.1316 0.1163 0.8451 0.1936 0.2322
0.5394 0.0658 0.1325 0.2338 0.6649 0.4365
0.6742 0.2632 0.1487 0.3127 0.4767 0.3593
R =
5.1405 2.0226
0 1.3817
0 0
0 0
0 0
0 0
DotNetMatrix doesn't work correctly because it only gives the first two columns of Q.
I did this to use the DotNetMatrix:
// Create a jagged "array" based on A suitable for Method 1
// Create a b vector for use in the QR decomposition of A
double[][] array = new double[A.GetLength(0)  1][];
double[][] bMeth1 = new double[A.GetLength(0)  1][];
(assign array & bMeth1)
// Create a GeneralMatrix AMeth1
GeneralMatrix AMeth1 = new GeneralMatrix(array);
//// Perform the QR decomposition
QRDecomposition QRofA = new QRDecomposition(AMeth1);
GeneralMatrix QRsoln = QRofA.Solve(new GeneralMatrix(bMeth1));
GeneralMatrix QMeth1 = QRofA.Q;
GeneralMatrix RMeth1 = QRofA.R;
Best regards,
Ralf





Thanks for the information. I used these to test the classes
for the next update.
BTW, you have posted to a wrong thread.
Also, this is the easy way to do the QRD
class Program
{
static void Main(string[] args)
{
double[] values = { 0, 0.6931, 1.3863, 2.0794, 2.7726, 3.4657,
1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000 };
GeneralMatrix matrix = new GeneralMatrix(values, 6);
Print("The Matrix A", matrix);
QRDecomposition qr = matrix.QRD();
Print("The Matrix Q", qr.Q);
Print("The Matrix R", qr.R);
}
static void Print(string caption, GeneralMatrix matrix)
{
Console.WriteLine(caption);
int m = matrix.RowDimension;
int n = matrix.ColumnDimension;
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
double dE = matrix.GetElement(i, j);
if (dE >= 0)
{
Console.Write(" {0:f4} ", dE);
}
else
{
Console.Write("{0:f4} ", dE);
}
}
Console.WriteLine();
}
Console.WriteLine();
}
}
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Oops, sorry Paul about posting in the wrong place; I was in a bit of a hurry this morning. Anyway, you got the information. Glad to see that you are generously providing this.
Sincerely,
Ralf





Please try "QR(A,0)" in matlab and see the results.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.







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.

