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

Comments and Discussions



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.





Here you go, Paul:
QR(A,0)
ans =
5.1405 2.0226
0.1348 1.3817
0.2697 0.0936
0.4045 0.1138
0.5394 0.3193
0.6742 0.5257
Best,
Ralf





I hope you see this is the result the DotNetMatrix gives, it is the "economic size" results as documented by the DotNetMatrix.
You will really not need the full size matrix Q that the Matlab gives in qr(A),
but I will see how to add that support too in the next update, if really required.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Hi Paul,
The alglib.net code worked well enough for me, and so that is what I wound up using. I was not able to use the DotNetMatrix "economic size" results code because I was duplicating another algorithm that uses the full Q and R (whether that's really necessary, I'm not in a position to dig into right now). I certainly do not want to create more work for you. Unless others are asking for the full Q, I would not bother. I would not add this feature in lieu of say, providing some more fully documented examples of the code that is already in there.
I know I was not entirely sure how to properly use the code and what it could do when it came down to the details. I seems to me that any improvements in the documentation for the existing features (e.g., more examples) would be high yield.
Best wishes,
Ralf





Hi. In the CholeskyDecomposition constructor, a flag isspd is set. I have an example where this flag turned out to give me the wrong value.
When looking in the code, I can see a line like this:
isspd = isspd & (A[k][j] == A[j][k]);
since we are dealing with double's, this is a very dangerous way of checking for equality. Dangerous enough, at least, to make the code not work for my example.
I always recommend to use something like this:
isspd = isspd & CloseEnoughToEqual(A[k][j], A[j][k]);
where
CloseEnoughToEqual(a,b)
{
return Math.Abs(ab) < SOMESMALLTHRESHOULD
}
SOMESMALLTHRESHOULD could then be 0.0000001 or less depending on precision. In my case, the values were 0.38490014055502542 versus 0.38490014055503252 making the code believe those two values weren't suppose to be the same
Beside this, I have found the library very useful, and well working.
Claws





You are right, when dealing with double, precision and its effect on robustness is a great factor.
In any case, this library needs a lot of updates and I hope to do that soon.
With love,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Hmm Does this Class handle or is able to handle the PenroseMoore Pseudo Inverse ?
I tried to modify it... however seems there are too many other functions that cannot handle the singular matrixs
also the SVD can it handle singular matrix ?





Hi, how could I calculate the matrix LEFT division
X = A\B
where A is af nonsquare matrix, with your library?
Thanks





In document, it said that GeneralMatrix.Inverse Method will "inverse(A) if A is square, pseudoinverse otherwise."
But I just got an error "Matrix is rank deficient"
The matrix data is
0 0 4 0 4
0 0 0 8 8
0 10 10 10 10
1 1 1 1 1
And if I try to SVD the matrix, the matrix has some problem.
The matrix from GetU() has ColumnDimension = 5, but the array A inside the atrix is 4 x 4.
Is this a bug?





as the author mentioned, the row_dimension must >= column_dimension
in the raw matrix.





public virtual GeneralMatrix Inverse()
{
if (this.m == this.n)
{
return Solve(Identity(m, m));
}
else if (this.m >= this.n)
{
SingularValueDecomposition SVD = this.SVD();
GeneralMatrix U = SVD.GetU();
GeneralMatrix S = SVD.S;
GeneralMatrix V = SVD.GetV();
GeneralMatrix TranS = S.Transpose();
return V*(TranS*S).Inverse()*TranS*U.Transpose();
}
else
{
SingularValueDecomposition SVD = this.Transpose().SVD();
GeneralMatrix U = SVD.GetU();
GeneralMatrix S = SVD.S;
GeneralMatrix V = SVD.GetV();
GeneralMatrix TranS = S.Transpose();
return (V*(TranS*S).Inverse()*TranS*U.Transpose()).Transpose();
}
}





Hi
I use your code and i have following problem.
VB.NET project and dll includes your classes with DotNetMatrix.
When I use GeneralMatrix objects in For ...Next there is OutOfMemory Exception in your source in differnt places at each time : more in constructors and sometimes in other methods. Stack memory and heap are not overflowed.
Exception was raised in line with "new double[][]" code
Can you help with solution problem? Maybe problem isn't in your sources...





Please post a simple code sample to examine.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





I am having the same problem but because I want to compute the SVD of a 500x65536 matrix, which is causing an Out of Memory Exception allocating the U matrix (65536x65536, double precision).
Do you have any smart implementation where the hard disk is used in place of the main memory to compute the SVD of really huge matrices?
Thanks,
Sebastian.





senrique wrote:
(65536x65536, double precision).
Wow! this is a simple library for very simple matrix operations.
senrique wrote:
Do you have any smart implementation where the hard disk is used in place of the main memory to compute the SVD of really huge matrices?
Currently no, in fact I never thought many would have interest in this matrix stuff. To handle huge matrix will involve defining complex data structures, we have to take on this when the generics support in .NET 2.0 is released.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Great Matrix code. Very useful.
Just one thing: I can't seem to find out how to print the matrix's, even though its said "Methods for reading and printing matrices are also included." Maybe I'm being stupid or maybe these methods haven't been writen.
Anyway. I added my own, very simple ToString() method. I'll post it here as it might save someone a bit of time. Just paste it into the GeneralMatrix.cs file and compile.
/// <summary>
/// Converts the array into a string.
/// "First row \n Second row \n Last row \n"
///
/// <returns>String of Array Elements
public override string ToString()
{
string elementsString = "";
for(int j=0; j < m; j++)
{
for(int i=0; i < n; i++)
elementsString += A[j][i].ToString() + " ";
elementsString += "\n";
}
return elementsString;
}
Mark Sugrue
Machine Vision Research Group
Royal Holloway, Uni of London
m.sugrue@rhul.ac.uk





seabhcan wrote:
t one thing:...
Just ignore it, it is a mistake in the documentations  sorry.
seabhcan wrote:
Anyway. I added my own, very simple ToString() method.
Use StringBuilder, it is more efficient for such work.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





thanks for your help, it is save me a lot of time...





Hi,
Is this library can resolve a least square non linear problem ??
Which algorithm is implemented for the least square (LevenbergMarquardt, Gauss, ...) ???
In advance, thanks.
See you,
Philippe





Sorry for the delay in responding.
The simple answer is no. This is a simple matrix library for simple matrix operations. Even though you can extend it for more complex algebra, it is not designed for such works.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Sorry, I deleted your last mail you directly sent to me. But, I have read it and answer is:
If I have such library, I will not let you ask
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Great job! I just have a look at the code. I was wondering why the matrices are represented by a double[][] array instead of a double[ , ] array. Since, we will be dealing only with rectangular structures, I don't really see the advantage of using double[][]. Additionnaly, the use of double[][] implies one additionnal indirection to access to any element. That's not going to make a big change, but since the overall time of any operation on the matrix is going to be almost proportionnal to the element accessing time, it could still make a difference.
Also, there are details bothering me. Why implementing the IDisposable and destructor interfaces for GeneralMatrix ? As far, I have seen, it's not required. Remember that implementing IDisposable is really going to slow down the garbadge collection process (idem for destructor).
Idem for ISerializable, the attribute [Serializable] is itself sufficient. So why implementing the ISerializable interface ?





Hello,
Thanks for the information.
Joannes Vermorel wrote:
I was wondering why the matrices are represented by a double[][] array instead of a double[ , ] array.
The DotNetMatrix is basically a port of the JAMA which also uses the double[][]. I will look into this and see how best to correct it  thanks.
Joannes Vermorel wrote:
Remember that implementing IDisposable is really going to slow down the garbadge collection process (idem for destructor).
I have not seen that information anywhere, any link? Every .NET class has a destructor implemented as Finalize, the Dipose and its GC.SuppressFinalize(true) reduces the overhead of the GC calling the Finalizer. This is the picture I have so far from the docs.
Joannes Vermorel wrote:
Idem for ISerializable, the attribute [Serializable] is itself sufficient. So why implementing the ISerializable interface
I have stated this in the article and as you can see the ISerializable implementation is empty. I do not know if we could better control the serialization thru the ISerialization than the [Serializable], still looking it the issue.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Paul Selormey wrote:
I have not seen that information anywhere, any link? Every .NET class has a destructor implemented as Finalize, the Dipose and its GC.SuppressFinalize(true) reduces the overhead of the GC calling the Finalizer. This is the picture I have so far from the docs.
Well, basically, IDisposable means for the GC, "caution, this object is special, do not treat this object like the others", thus implying overhead (sorry I do not have any reference in mind, but try Google). With .Net, unless you are currently handling physical ressources, basically you never need IDisposable nor the destructor. The GC will take care of destroying the objects in an efficient fashion.
Paul Selormey wrote:
I have stated this in the article and as you can see the ISerializable implementation is empty. I do not know if we could better control the serialization thru the ISerialization than the [Serializable], still looking it the issue.
The question is "Why do you want to control the serialization?". Your objects are very simple. Just remove the ISerializable interfaces and leave the [Serializable] attribute. The BinarySerializer will do everything for you. With .Net, ISerializable is basically needed only for versionning issue.





Joannes Vermorel wrote:
Well, basically, IDisposable means for the GC
Well, untill I point to an authoritative source, I do not think it is an information to rely on. For all I know, GC does not know IDisposable. Also, since the actual purpose of here is to implement Affine Transformation, I will be dealing with the system matrix, which implements IDisposable, somehow.
Joannes Vermorel wrote:
Why do you want to control the serialization?".
Efficiency, if there is. As I have said, it is an issue I am still looking at.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





IDisposable is typically implemented when you use unmanaged resources. See Cleaning Up Unmanaged Resources[^] in the .NET Framework SDK (how's that for authoritative?) for a discussion of what it is, when to use it, and how to tie together your Dispose and Finalize (destructor in C# and MC++) implementations.
Microsoft MVP, Visual C#
My Articles





Here is explanation of use double[][]…
I understand performance implications of double[][], but KarhunenLoeve transform often requires sorting of Eigen vectors in order of modulus of Eigen values, so vectors with biggest eigen values are on the top of the matrix (this frequently used in data dimensiality reduction). It will be tough to sort double[,], since it will require copying of row elements over to a new index, where in double[][] t is just switching the pointer to array of double.
Sincerely,
Alex Cherkasov
http://www.xpidea.com





Hi Paul,
very interesting tutorial !
Any target date planned for the affine transformation ?
Jan.





Thanks.
Never knew others are really interested. I will post the codes by next week Feb. 6.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Hi Paul,
ehm,...curious as when the aditional code will be posted.
Assuming you ment Feb. 6 2004 instead of 2005
Jan.





Hello Jan,
I understand and please accept my apology. I am working on related stuff and will post all as soon as possible. The Affine Transformation is just a little part of the package.
At least it will not be 2005.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





No problem





Hi Paul,
been a while since i read the Message Board, but i see i am just in time,...
"At least it will not be 2005."
Ehm,... hope you have fast fingers
Later,
Jan.





I believe that the GDI+ Matrix class provides several methods for Affine transforms.





Hi StenSat,
yes correct.
I am looking for a method which calculates the parameters between sourcepoints and destinationpoints.





We needed a replacement for the GDI+ Matrix class in our graphics system, VG.net. We found that for a small matrix, 3x3 affine, which is internally 2x3 floats for vector graphics, a struct provides much better performance than a class, even if you pass the matrices around as method arguments often. Furthermore, you get value semantics for your operator overloading without having to hit the heap or bother the garbage collector. Our idea was validated recently when we noticed Avalon does the same thing. Not quite as fast as C++ but close!
If you are interested, get a VG.net beta (see sig), and look at the FMatrix class. We use it internally instead of System.Drawing.Matrix. On rendering, we must use the System.Drawing class, so we do, but we avoid it as much as possible, as it uses interop and is very slow.
Regards,
Frank Hileman
check out VG.net: www.vgdotnet.com
An animated vector graphics system integrated in VS.net





Thanks for the information. I will look into it.
Frank Hileman wrote:
If you are interested, get a VG.net beta (see sig), and look at the FMatrix class.
Just submitted the beta form. Nice gradient stuff you have in there, I will love your article on the use of GDI+ gradient classes
Frank Hileman wrote:
We use it internally instead of System.Drawing.Matrix. On rendering, we must use the System.Drawing class, so we do, but we avoid it as much as possible, as it uses interop and is very slow.
Good to know that I am not alone on this. I hope to post an article on my affine transformation matrix soon. Hope you will take a look.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Hello,
great work. I have the same type of library in VB on Code project, but not nearly as extended. Also it doesn't define a matrix type, but is just a sealed class that provides shared methods for working with 1 and 2 D arrays. I do have a multidimensional nonlinear solver though, based on the NelderMead simplex algorithm as found in 'Numerical Recipes in C'. It is literally translated from C, and is unreadable, but if you just transform the in and output data formats it should work .... This great for data and model fitting.
Cheers





Thanks, just found your article on the algebra. Sorry, I did not look into the VB.NET section before posting this work.
I will read your article more closely to see how to use some of the ideas and if possible codes there.
Thanks for the contribution and the suggestions.
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.






Hmmmm, never knew TeXnicCenter fans are here too
Best regards,
Paul.
Jesus Christ is LOVE! Please tell somebody.





Hi paul, it is nice to see someone taking care of the dirty work (implementing a matrix library).
Jonathan de Halleux.
www.dotnetwiki.org







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.

