Oct 11, 2012 at 3:03 PM
Edited Oct 11, 2012 at 3:13 PM

Hello,
I am trying to setup a least squares filter, but whenever I invert the desnse matrix, it becomes full of NANs. I slapped a debugger on Mathdotnet and found that
line1113ish in Algorithms\LinearAlgebra\ManagedLinearAlgebraProvider.Double.cs
b[k + (j * order)] /= a[korder];
has the a[korder] set to zero and is causing these NANs. I believe "a" is getting set during LU factorization inside the DenseLU class. Please let me know what I am doing wrong? Is this just a case of hitting the limits of double during factorization?
~Thanks,
~Matthew Johnson
//THis code will repro the issue I am seeing
//using polynomial least squares Xa = A
//Populate X
// X = 1 x0 y0 x0^2 x0y0 y0^2 x0^3 x0^2y0 x0y0^2 y0^3
// 1 x1 y1 x1^2 x1y1 y1^2 x1^3 x1^2y1 x1y1^2 y1^3
// : : : : : : : : : :
// row n
int index = 0;
for (int n = 0; n < 35; ++n, ++index)
{
X[n, 0] = 1; // 1
X[n, 1] = index; //x0
X[n, 2] = index; //y0
X[n, 3] = (double)Math.Pow(index, 2); //x0^2
X[n, 4] = (double)Math.Pow(index, 2); //x0y0
X[n, 5] = (double)Math.Pow(index, 2); //y0^2
X[n, 6] = (double)Math.Pow(index, 3); //x0^3
X[n, 7] = (double)Math.Pow(index, 2) * index; //x0^2y0
X[n, 8] = (double)Math.Pow(index, 2) * index; //x0y0^2
X[n, 9] = (double)Math.Pow(index, 3); //y0^3
}
// solving for a
//a = (XTX)^(1)*XTA (XT is the transpose of the polynomial coeff).
var C = X.Transpose();
C = C.Multiply(X);
C = C.Inverse();


Coordinator
Oct 12, 2012 at 8:31 AM
Edited Oct 12, 2012 at 8:48 AM

Without actually looking at the numbers: The matrix inverse is numerically a very unreliable way to solve least squares problems. Have you tried using the QR decomposition instead?
b = inv(XTX)*XTy
Hence: XTXb = XTy b = (XTX).QR().Solve(XTy)
or (if I'm not mistaken) simply: b = X.QR().Solve(y)
Let me know whether this still yields NaN's


Oct 12, 2012 at 3:56 PM
Edited Oct 12, 2012 at 3:58 PM

Thanks for the direction:
The QR decomposition also gave me NANs. However, I went through the whole list of functions in ExtentionMethods.cs and found the only one that did not give me NANs was:
C = C.Evd().Solve(X.Transpose());
Except for some Linear Algebra taken over 10 years ago in school, I am not really qualified to determine what issues I may run into using the EVD decomp process. Since I am more interested in accuracy over speed (and I have been trying to do some research
on the subject). This has led me to a few more questions:
Am I correct in assuming the GaussJordan elimination method for getting the inverse is more "accurate", or should I not be concerned with EVD?
Here is the research that was the basis of the above statement:
"The method’s principal strength is that it is as stable as any other direct method, perhaps even a bit more stable when
full pivoting is used"
http://courses.cms.caltech.edu/cs171/c21.pdf
Is there a better (in terms of accuracy) method for interpolating a 2D surface (aka Bspline ect) using mathdotnet (or is this something needed in mathdotnet)? All of the examples show 1D fitting, 2D did not seem to be an option.
~Thanks,
~Matt


Oct 13, 2012 at 6:38 AM
Edited Oct 13, 2012 at 6:38 AM

I think this may be the answer I was looking for (avoiding matrix math altogether). Instead of using a Bspline surface I think any interpolation method would work using this links column first then column interpolated point to row interpolation method
http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/INTAPP/SURFAPPglobal.html

