This project has moved and is read-only. For the latest updates, please go here.

Polynomial Least Squares Question

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

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?


~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();


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

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"

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.





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