Implicit Weights

The weights are no longer implicit in LeastSquaresProblem.Evaluation. They are
already included in the computed residuals and Jacobian.

GN and LM multiplied the residuals by the weights immediately, so that was easy
to remove.

Created an AbstractEvaluation class which handles the derived quantitied (cost,
rms, covariance,...) and two implementations. UnweightedEvaluation uses the
straight forward formulas. DenseWeightedEvaluation delegates to an Evaluation
and multiples the residuals and Jacobian by the square root of the weight matrix
before returning them. Allowed me to remove the reference to the full weight
matrix.

Author: wardev

src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractEvaluation.java

@@ -0,0 +1,70 @@
+package org.apache.commons.math3.fitting.leastsquares;
+
+import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation;
+import org.apache.commons.math3.linear.ArrayRealVector;
+import org.apache.commons.math3.linear.DecompositionSolver;
+import org.apache.commons.math3.linear.QRDecomposition;
+import org.apache.commons.math3.linear.RealMatrix;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * An implementation of {@link Evaluation} that is designed for extension. All of the
+ * methods implemented here use the methods that are left unimplemented.
+ * <p/>
+ * TODO cache results?
+ *
+ * @version $Id$
+ */
+abstract class AbstractEvaluation implements Evaluation {
+
+    /** number of observations */
+    private final int observationSize;
+
+    /**
+     * Constructor.
+     *
+     * @param observationSize the number of observation. Needed for {@link
+     *                        #computeRMS()}.
+     */
+    AbstractEvaluation(final int observationSize) {
+        this.observationSize = observationSize;
+    }
+
+    /** {@inheritDoc} */
+    public double[][] computeCovariances(double threshold) {
+        // Set up the Jacobian.
+        final RealMatrix j = this.computeJacobian();
+
+        // Compute transpose(J)J.
+        final RealMatrix jTj = j.transpose().multiply(j);
+
+        // Compute the covariances matrix.
+        final DecompositionSolver solver
+                = new QRDecomposition(jTj, threshold).getSolver();
+        return solver.getInverse().getData();
+    }
+
+    /** {@inheritDoc} */
+    public double[] computeSigma(double covarianceSingularityThreshold) {
+        final double[][] cov = this.computeCovariances(covarianceSingularityThreshold);
+        final int nC = cov.length;
+        final double[] sig = new double[nC];
+        for (int i = 0; i < nC; ++i) {
+            sig[i] = FastMath.sqrt(cov[i][i]);
+        }
+        return sig;
+    }
+
+    /** {@inheritDoc} */
+    public double computeRMS() {
+        final double cost = this.computeCost();
+        return FastMath.sqrt(cost * cost / this.observationSize);
+    }
+
+    /** {@inheritDoc} */
+    public double computeCost() {
+        final ArrayRealVector r = new ArrayRealVector(this.computeResiduals());
+        return FastMath.sqrt(r.dotProduct(r));
+    }
+
+}

src/main/java/org/apache/commons/math3/fitting/leastsquares/DenseWeightedEvaluation.java

@@ -0,0 +1,55 @@
+package org.apache.commons.math3.fitting.leastsquares;
+
+import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation;
+import org.apache.commons.math3.linear.RealMatrix;
+
+/**
+ * Applies a dense weight matrix to an evaluation.
+ *
+ * @version $Id$
+ */
+class DenseWeightedEvaluation extends AbstractEvaluation {
+
+    /** the unweighted evaluation */
+    private final Evaluation unweighted;
+    /** reference to the weight square root matrix */
+    private final RealMatrix weightSqrt;
+
+    /**
+     * Create a weighted evaluation from an unweighted one.
+     *
+     * @param unweighted the evalutation before weights are applied
+     * @param weightSqrt the matrix square root of the weight matrix
+     */
+    DenseWeightedEvaluation(final Evaluation unweighted,
+                            final RealMatrix weightSqrt) {
+        // weight square root is square, nR=nC=number of observations
+        super(weightSqrt.getColumnDimension());
+        this.unweighted = unweighted;
+        this.weightSqrt = weightSqrt;
+    }
+
+    /* apply weights */
+
+    /** {@inheritDoc} */
+    public RealMatrix computeJacobian() {
+        return weightSqrt.multiply(this.unweighted.computeJacobian());
+    }
+
+    /** {@inheritDoc} */
+    public double[] computeResiduals() {
+        return this.weightSqrt.operate(this.unweighted.computeResiduals());
+    }
+
+    /* delegate */
+
+    /** {@inheritDoc} */
+    public double[] getPoint() {
+        return unweighted.getPoint();
+    }
+
+    /** {@inheritDoc} */
+    public double[] computeValue() {
+        return unweighted.computeValue();
+    }
+}

src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java

@@ -104,17 +104,10 @@ public Optimum optimize(final LeastSquaresProblem lsp) {
             throw new NullArgumentException();
         }
 
-        final RealMatrix weightMatrix = lsp.getWeight();
-        final int nR = weightMatrix.getRowDimension(); // Number of observed data.
-
-        // Diagonal of the weight matrix.
-        final double[] residualsWeights = new double[nR];
-        for (int i = 0; i < nR; i++) {
-            residualsWeights[i] = weightMatrix.getEntry(i, i);
-        }
+        final int nR = lsp.getObservationSize(); // Number of observed data.
+        final int nC = lsp.getParameterSize();
 
         final double[] currentPoint = lsp.getStart();
-        final int nC = currentPoint.length;
 
         // iterate until convergence is reached
         PointVectorValuePair current = null;
@@ -128,7 +121,7 @@ public Optimum optimize(final LeastSquaresProblem lsp) {
             final Evaluation value = lsp.evaluate(currentPoint);
             final double[] currentObjective = value.computeValue();
             final double[] currentResiduals = value.computeResiduals();
-            final RealMatrix weightedJacobian = value.computeWeightedJacobian();
+            final RealMatrix weightedJacobian = value.computeJacobian();
             current = new PointVectorValuePair(currentPoint, currentObjective);
 
             // build the linear problem
@@ -137,21 +130,20 @@ public Optimum optimize(final LeastSquaresProblem lsp) {
             for (int i = 0; i < nR; ++i) {
 
                 final double[] grad = weightedJacobian.getRow(i);
-                final double weight = residualsWeights[i];
                 final double residual = currentResiduals[i];
 
                 // compute the normal equation
-                final double wr = weight * residual;
+                //residual is already weighted
                 for (int j = 0; j < nC; ++j) {
-                    b[j] += wr * grad[j];
+                    b[j] += residual * grad[j];
                 }
 
                 // build the contribution matrix for measurement i
                 for (int k = 0; k < nC; ++k) {
                     double[] ak = a[k];
-                    double wgk = weight * grad[k];
+                    //Jacobian/gradient is already weighted
                     for (int l = 0; l < nC; ++l) {
-                        ak[l] += wgk * grad[l];
+                        ak[l] += grad[k] * grad[l];
                     }
                 }
             }

src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresAdapter.java

@@ -0,0 +1,60 @@
+package org.apache.commons.math3.fitting.leastsquares;
+
+import org.apache.commons.math3.optim.ConvergenceChecker;
+import org.apache.commons.math3.optim.PointVectorValuePair;
+import org.apache.commons.math3.util.Incrementor;
+
+/**
+ * An adapter that delegates to another implementation of {@link LeastSquaresProblem}.
+ *
+ * @version $Id$
+ */
+public class LeastSquaresAdapter implements LeastSquaresProblem {
+
+    /** the delegate problem */
+    private final LeastSquaresProblem problem;
+
+    /**
+     * Delegate the {@link LeastSquaresProblem} interface to the given implementation.
+     *
+     * @param problem the delegate
+     */
+    public LeastSquaresAdapter(final LeastSquaresProblem problem) {
+        this.problem = problem;
+    }
+
+    /** {@inheritDoc} */
+    public double[] getStart() {
+        return problem.getStart();
+    }
+
+    /** {@inheritDoc} */
+    public int getObservationSize() {
+        return problem.getObservationSize();
+    }
+
+    /** {@inheritDoc} */
+    public int getParameterSize() {
+        return problem.getParameterSize();
+    }
+
+    /** {@inheritDoc} */
+    public Evaluation evaluate(final double[] point) {
+        return problem.evaluate(point);
+    }
+
+    /** {@inheritDoc} */
+    public Incrementor getEvaluationCounter() {
+        return problem.getEvaluationCounter();
+    }
+
+    /** {@inheritDoc} */
+    public Incrementor getIterationCounter() {
+        return problem.getIterationCounter();
+    }
+
+    /** {@inheritDoc} */
+    public ConvergenceChecker<PointVectorValuePair> getConvergenceChecker() {
+        return problem.getConvergenceChecker();
+    }
+}