Polynomial Least Squares Question

Oct 11, 2012 at 4:03 PM
Edited Oct 11, 2012 at 4: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 9:31 AM
Edited Oct 12, 2012 at 9: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 4:56 PM
Edited Oct 12, 2012 at 4: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 Gauss-Jordan 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/c2-1.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 7:38 AM
Edited Oct 13, 2012 at 7: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/INT-APP/SURF-APP-global.html