Add your own alternative version
Stats
393.1K views 11.5K downloads 123 bookmarked
Posted
12 Jan 2004

Comments and Discussions


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:
/// <summary>Return the diagonal matrix of singular values</summary><br />
/// <returns> S<br />
/// </returns><br />
virtual public GeneralMatrix S<br />
{<br />
get<br />
{<br />
GeneralMatrix X = new GeneralMatrix(n, n); // this should be (m, n)<br />
double[][] S = X.Array;<br />
for (int i = 0; i < n; i++) // this loop should go till m not n<br />
{<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.





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






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.

