Example Polynomial Fitting - Revision history

Peter at 15:30, 7 July 2021

2021-07-07T15:30:10Z

Peter at 15:29, 7 July 2021

2021-07-07T15:29:51Z

Peter at 19:49, 18 May 2017

2017-05-18T19:49:52Z

Peter at 00:59, 10 August 2015

2015-08-10T00:59:49Z

Peter at 02:42, 1 April 2015

2015-04-01T02:42:29Z

Peter: Created page with "In this example it is shown how EJML can be used to fit a polynomial of arbitrary degree to a set of data. The key concepts shown here are; 1) how to create a linear using Li..."

2015-03-22T05:34:19Z

Created page with "In this example it is shown how EJML can be used to fit a polynomial of arbitrary degree to a set of data. The key concepts shown here are; 1) how to create a linear using Li..."

New page

In this example it is shown how EJML can be used to fit a polynomial of arbitrary degree to a set of data. The key concepts shown here are; 1) how to create a linear using LinearSolverFactory, 2) use an adjustable linear solver, 3) and effective matrix reshaping. This is all done using the procedural interface.

First a best fit polynomial is fit to a set of data and then a outliers are removed from the observation set and the coefficients recomputed. Outliers are removed efficiently using an adjustable solver that does not resolve the whole system again.

Example on GitHub:
* [https://github.com/lessthanoptimal/ejml/blob/master/examples/src/org/ejml/example/PolynomialFit.java PolynomialFit]

= PolynomialFit Example Code =

<syntaxhighlight lang="java">
/**
* <p>
* This example demonstrates how a polynomial can be fit to a set of data. This is done by
* using a least squares solver that is adjustable. By using an adjustable solver elements
* can be inexpensively removed and the coefficients recomputed. This is much less expensive
* than resolving the whole system from scratch.
* </p>
* <p>
* The following is demonstrated:<br>
* <ol>
* <li>Creating a solver using LinearSolverFactory</li>
* <li>Using an adjustable solver</li>
* <li>reshaping</li>
* </ol>
* @author Peter Abeles
*/
public class PolynomialFit {

// Vandermonde matrix
DenseMatrix64F A;
// matrix containing computed polynomial coefficients
DenseMatrix64F coef;
// observation matrix
DenseMatrix64F y;

// solver used to compute
AdjustableLinearSolver solver;

/**
* Constructor.
*
* @param degree The polynomial's degree which is to be fit to the observations.
*/
public PolynomialFit( int degree ) {
coef = new DenseMatrix64F(degree+1,1);
A = new DenseMatrix64F(1,degree+1);
y = new DenseMatrix64F(1,1);

// create a solver that allows elements to be added or removed efficiently
solver = LinearSolverFactory.adjustable();
}

/**
* Returns the computed coefficients
*
* @return polynomial coefficients that best fit the data.
*/
public double[] getCoef() {
return coef.data;
}

/**
* Computes the best fit set of polynomial coefficients to the provided observations.
*
* @param samplePoints where the observations were sampled.
* @param observations A set of observations.
*/
public void fit( double samplePoints[] , double[] observations ) {
// Create a copy of the observations and put it into a matrix
y.reshape(observations.length,1,false);
System.arraycopy(observations,0, y.data,0,observations.length);

// reshape the matrix to avoid unnecessarily declaring new memory
// save values is set to false since its old values don't matter
A.reshape(y.numRows, coef.numRows,false);

// set up the A matrix
for( int i = 0; i < observations.length; i++ ) {

double obs = 1;

for( int j = 0; j < coef.numRows; j++ ) {
A.set(i,j,obs);
obs *= samplePoints[i];
}
}

// process the A matrix and see if it failed
if( !solver.setA(A) )
throw new RuntimeException("Solver failed");

// solver the the coefficients
solver.solve(y,coef);
}

/**
* Removes the observation that fits the model the worst and recomputes the coefficients.
* This is done efficiently by using an adjustable solver. Often times the elements with
* the largest errors are outliers and not part of the system being modeled. By removing them
* a more accurate set of coefficients can be computed.
*/
public void removeWorstFit() {
// find the observation with the most error
int worstIndex=-1;
double worstError = -1;

for( int i = 0; i < y.numRows; i++ ) {
double predictedObs = 0;

for( int j = 0; j < coef.numRows; j++ ) {
predictedObs += A.get(i,j)*coef.get(j,0);
}

double error = Math.abs(predictedObs- y.get(i,0));

if( error > worstError ) {
worstError = error;
worstIndex = i;
}
}

// nothing left to remove, so just return
if( worstIndex == -1 )
return;

// remove that observation
removeObservation(worstIndex);

// update A
solver.removeRowFromA(worstIndex);

// solve for the parameters again
solver.solve(y,coef);
}

/**
* Removes an element from the observation matrix.
*
* @param index which element is to be removed
*/
private void removeObservation( int index ) {
final int N = y.numRows-1;
final double d[] = y.data;

// shift
for( int i = index; i < N; i++ ) {
d[i] = d[i+1];
}
y.numRows--;
}
}
</syntaxhighlight>

@@ Line 5: / Line 5: @@
 External Resources:
 * [https://github.com/lessthanoptimal/ejml/blob/v0.41/examples/src/org/ejml/example/PolynomialFit.java PolynomialFit.java source code]
 = PolynomialFit Example Code =

@@ Line 4: / Line 4: @@
 External Resources:
-* [https://github.com/lessthanoptimal/ejml/blob/v0.27/examples/src/org/ejml/example/PolynomialFit.java PolynomialFit.java source code]
+* [https://github.com/lessthanoptimal/ejml/blob/v0.31/examples/src/org/ejml/example/PolynomialFit.java PolynomialFit.java source code]
 * <disqus>Discuss this example</disqus>
@@ Line 29: / Line 29: @@
      // Vandermonde matrix
-     DenseMatrix64F A;
+     DMatrixRMaj A;
      // matrix containing computed polynomial coefficients
-     DenseMatrix64F coef;
+     DMatrixRMaj coef;
      // observation matrix
-     DenseMatrix64F y;
+     DMatrixRMaj y;
      // solver used to compute
-     AdjustableLinearSolver solver;
+     AdjustableLinearSolver_DDRM solver;
      /**
@@ Line 44: / Line 44: @@
       */
      public PolynomialFit( int degree ) {
-         coef = new DenseMatrix64F(degree+1,1);
+         coef = new DMatrixRMaj(degree+1,1);
-         A = new DenseMatrix64F(1,degree+1);
+         A = new DMatrixRMaj(1,degree+1);
-         y = new DenseMatrix64F(1,1);
+         y = new DMatrixRMaj(1,1);
          // create a solver that allows elements to be added or removed efficiently
-         solver = LinearSolverFactory.adjustable();
+         solver = LinearSolverFactory_DDRM.adjustable();
+     }

@@ Line 3: / Line 3: @@
 First a best fit polynomial is fit to a set of data and then a outliers are removed from the observation set and the coefficients recomputed.  Outliers are removed efficiently using an adjustable solver that does not resolve the whole system again.
-Example on GitHub:
+External Resources:
-* [https://github.com/lessthanoptimal/ejml/blob/v0.27/examples/src/org/ejml/example/PolynomialFit.java PolynomialFit]
+* [https://github.com/lessthanoptimal/ejml/blob/v0.27/examples/src/org/ejml/example/PolynomialFit.java PolynomialFit.java source code]
 = PolynomialFit Example Code =

@@ Line 4: / Line 4: @@
 External Resources:
-* [https://github.com/lessthanoptimal/ejml/blob/v0.31/examples/src/org/ejml/example/PolynomialFit.java PolynomialFit.java source code]
+* [https://github.com/lessthanoptimal/ejml/blob/v0.41/examples/src/org/ejml/example/PolynomialFit.java PolynomialFit.java source code]
 * <disqus>Discuss this example</disqus>
@@ Line 24: / Line 24: @@
   *  <li>reshaping</li>
   * </ol>
   * @author Peter Abeles
   */
 public class PolynomialFit {
      // Vandermonde matrix
      DMatrixRMaj A;
@@ Line 44: / Line 44: @@
       */
      public PolynomialFit( int degree ) {
-         coef = new DMatrixRMaj(degree+1,1);
+         coef = new DMatrixRMaj(degree + 1, 1);
-         A = new DMatrixRMaj(1,degree+1);
+         A = new DMatrixRMaj(1, degree + 1);
-         y = new DMatrixRMaj(1,1);
+         y = new DMatrixRMaj(1, 1);
          // create a solver that allows elements to be added or removed efficiently
@@ Line 67: / Line 67: @@
       * @param observations A set of observations.
       */
-     public void fit( double samplePoints[] , double[] observations ) {
+     public void fit( double[] samplePoints, double[] observations ) {
          // Create a copy of the observations and put it into a matrix
-         y.reshape(observations.length,1,false);
+         y.reshape(observations.length, 1, false);
-         System.arraycopy(observations,0, y.data,0,observations.length);
+         System.arraycopy(observations, 0, y.data, 0, observations.length);
          // reshape the matrix to avoid unnecessarily declaring new memory
          // save values is set to false since its old values don't matter
-         A.reshape(y.numRows, coef.numRows,false);
+         A.reshape(y.numRows, coef.numRows, false);
          // set up the A matrix
-         for( int i = 0; i < observations.length; i++ ) {
+         for (int i = 0; i < observations.length; i++) {
              double obs = 1;
-             for( int j = 0; j < coef.numRows; j++ ) {
+             for (int j = 0; j < coef.numRows; j++) {
-                 A.set(i,j,obs);
+                 A.set(i, j, obs);
                  obs *= samplePoints[i];
+             }
@@ Line 88: / Line 88: @@
          // process the A matrix and see if it failed
-         if( !solver.setA(A) )
+         if (!solver.setA(A))
              throw new RuntimeException("Solver failed");
          // solver the the coefficients
-         solver.solve(y,coef);
+         solver.solve(y, coef);
+     }
@@ Line 103: / Line 103: @@
      public void removeWorstFit() {
          // find the observation with the most error
-         int worstIndex=-1;
+         int worstIndex = -1;
          double worstError = -1;
-         for( int i = 0; i < y.numRows; i++ ) {
+         for (int i = 0; i < y.numRows; i++) {
              double predictedObs = 0;
-             for( int j = 0; j < coef.numRows; j++ ) {
+             for (int j = 0; j < coef.numRows; j++) {
-                 predictedObs += A.get(i,j)*coef.get(j,0);
+                 predictedObs += A.get(i, j)*coef.get(j, 0);
+             }
-             double error = Math.abs(predictedObs- y.get(i,0));
+             double error = Math.abs(predictedObs - y.get(i, 0));
-             if( error > worstError ) {
+             if (error > worstError) {
                  worstError = error;
                  worstIndex = i;
@@ Line 122: / Line 122: @@
          // nothing left to remove, so just return
-         if( worstIndex == -1 )
+         if (worstIndex == -1)
              return;
@@ Line 132: / Line 132: @@
          // solve for the parameters again
-         solver.solve(y,coef);
+         solver.solve(y, coef);
+     }
@@ Line 141: / Line 141: @@
       */
      private void removeObservation( int index ) {
-         final int N = y.numRows-1;
+         final int N = y.numRows - 1;
-         final double d[] = y.data;
+         final double[] d = y.data;
          // shift
-         for( int i = index; i < N; i++ ) {
+         for (int i = index; i < N; i++) {
-             d[i] = d[i+1];
+             d[i] = d[i + 1];
+         }
          y.numRows--;