Curve Fitting via SVD

March 7, 2011

I used this method to delete a gradient effect from uneven lighting of some of my book scans.

Here is a simple method for curve fitting various polynomials. This works for non linear equations as long as the actual coeficients are linear.

For example you can NOT fit the following type of function:

Y={a_1cos(a_2x)+a_3x}

Note that the coeficients a_2 is an argument to a nonlinear function. related link

The specific one done here is a bicubic polynomial expressed by the following.

We are going to fit the following curve:

P(x,y) ={a_00+a_01y+a_02y^2+a_03y^3+a_10x+a_11xy+a_12xy^2+a_13xy^3+a_20x^2+a_21x^2y+a_22x^2y^2+a_23x^2y^3+a_30x^3+a_31x^3y+a_32x^3y^2+a_33x^3y^3}

Or more compactly:

P(x,y)=sum{i=0}{3}{sum{j=0}{3}{a_ij}{x^i}{y^j}}

 import cern.colt.matrix.linalg.*;
import cern.colt.matrix.*;
import cern.colt.matrix.impl.*;

import edu.umbc.cs.maple.utils.*;

 public DoubleMatrix2D calcFit(ArrayList<Point> measures)
 {
    int cols = 16; // number of coeficients
    int rows = measures.size();  // size of your dataset

    DenseDoubleMatrix2D M = new DenseDoubleMatrix2D(rows,cols);
    DenseDoubleMatrix2D P = new DenseDoubleMatrix2D(rows,1);

    for (int k=0;k<measures.size();k++)
    {
       Point measure = measures.get(k);  // measures from your dataset
       P.set(k,0,measure.getValue());    // Point contains values P , x and y

       double x = measure.getX();
       double y = measure.getY();

       // populate the matrix with the appropriate polynomial values
       // from the bicubic above
       int col=0;
       for (int i=0;i<4;i++)
         for (int j=0;j<4;j++)
         {
            double cal = pow(x,i)*pow(y,j);
            M.set(k,col++,cal);
         }
    }

    // solve SVD
    SingularValueDecomposition svd = new SingularValueDecomposition(M);
    DoubleMatrix2D U = svd.getU();
    DoubleMatrix2D V = svd.getV();
    DoubleMatrix2D S = svd.getS();

    // S is the diagonal matrix
    // Invert by replacing its diagonal elements by their reciprocals
    for (int i=0;i<16;i++)
    {
       double val = 1.0/S.get(i,i);

       if (S.get(i,i) == 0)  // avoid numeric underflow
         val = 0.0;

       S.set(i,i, val);
    }

    Algebra alg = new Algebra();

    // X pseudoinverse
    DoubleMatrix2D VS = alg.mult(V,S);

    DoubleMatrix2D Minv = alg.mult(VS,alg.transpose(U));

    // matrix has coeficients
    DoubleMatrix2D C = alg.mult(Minv, P);

   // If you would like the closeness of the fit use the following

   DoubleMatrix2D eP = alg.mult(M,C); // estimated values

   ColtUtils utils = new ColtUtils();
   DoubleMatrix2D dP = utils.minus(P,eP); // difference between estimated and observed values
  double sum = utils.dotproduct(utils.getcol(dP,0),utils.getcol(dP,0));

  System.out.println("sum = " + sum);  // how close was our fit?

  return(C);
}

In practice you may need to normalize your measurements if they are large to avoid numeric under and/or over flow.

Tags: ,

This entry was posted on March 7, 2011 at 7:02 pm and is filed under Uncategorized. You can follow any responses to this entry through the RSS 2.0 feed.

Comments are closed.