diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractEvaluation.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractEvaluation.java new file mode 100644 index 000000000..d0d820c62 --- /dev/null +++ b/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. + *

+ * 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)); + } + +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/DenseWeightedEvaluation.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/DenseWeightedEvaluation.java new file mode 100644 index 000000000..ca6a4adeb --- /dev/null +++ b/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(); + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java index 84e8e538d..dddaed529 100644 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java @@ -104,17 +104,10 @@ public class GaussNewtonOptimizer implements LeastSquaresOptimizer { 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 class GaussNewtonOptimizer implements LeastSquaresOptimizer { 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 class GaussNewtonOptimizer implements LeastSquaresOptimizer { 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]; } } } diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresAdapter.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresAdapter.java new file mode 100644 index 000000000..3e88d035c --- /dev/null +++ b/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 getConvergenceChecker() { + return problem.getConvergenceChecker(); + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresFactory.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresFactory.java index 7872788fb..67bf660a4 100644 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresFactory.java +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresFactory.java @@ -2,14 +2,57 @@ package org.apache.commons.math3.fitting.leastsquares; import org.apache.commons.math3.analysis.MultivariateMatrixFunction; import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.apache.commons.math3.linear.EigenDecomposition; import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; import org.apache.commons.math3.optim.ConvergenceChecker; import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.util.FastMath; import org.apache.commons.math3.util.Incrementor; -/** @author Evan Ward */ +/** + * A Factory for creating {@link LeastSquaresProblem}s. + * + * @version $Id$ + */ public class LeastSquaresFactory { + /** Prevent instantiation. */ + private LeastSquaresFactory() { + } + + /** + * Create a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem} + * from the given elements. There will be no weights applied (Identity weights). + * + * @param model the model function. Produces the computed values. + * @param jacobian the jacobian of the model with respect to the parameters + * @param observed the observed (target) values + * @param start the initial guess. + * @param checker convergence checker + * @param maxEvaluations the maximum number of times to evaluate the model + * @param maxIterations the maximum number to times to iterate in the algorithm + * @return the specified General Least Squares problem. + */ + public static LeastSquaresProblem create(final MultivariateVectorFunction model, + final MultivariateMatrixFunction jacobian, + final double[] observed, + final double[] start, + final ConvergenceChecker checker, + final int maxEvaluations, + final int maxIterations) { + return new LeastSquaresProblemImpl( + maxEvaluations, + maxIterations, + checker, + observed, + model, + jacobian, + start + ); + } + /** * Create a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem} * from the given elements. @@ -32,16 +75,50 @@ public class LeastSquaresFactory { final ConvergenceChecker checker, final int maxEvaluations, final int maxIterations) { - return new LeastSquaresProblemImpl( - maxEvaluations, - maxIterations, - checker, - observed, - weight, - model, - jacobian, - start - ); + return weightMatrix( + create( + model, + jacobian, + observed, + start, + checker, + maxEvaluations, + maxIterations + ), + weight); + } + + /** + * Apply a dense weight matrix to the {@link LeastSquaresProblem}. + * + * @param problem the unweighted problem + * @param weights the matrix of weights + * @return a new {@link LeastSquaresProblem} with the weights applied. The original + * {@code problem} is not modified. + */ + public static LeastSquaresProblem weightMatrix(final LeastSquaresProblem problem, + final RealMatrix weights) { + final RealMatrix weightSquareRoot = squareRoot(weights); + return new LeastSquaresAdapter(problem) { + @Override + public Evaluation evaluate(final double[] point) { + return new DenseWeightedEvaluation(super.evaluate(point), weightSquareRoot); + } + }; + } + + /** + * Apply a diagon weight matrix to the {@link LeastSquaresProblem}. + * + * @param problem the unweighted problem + * @param weights the diagonal of the weight matrix + * @return a new {@link LeastSquaresProblem} with the weights applied. The original + * {@code problem} is not modified. + */ + public static LeastSquaresProblem weightDiagonal(final LeastSquaresProblem problem, + final RealVector weights) { + //TODO more efficient implementation + return weightMatrix(problem, new DiagonalMatrix(weights.toArray())); } /** @@ -55,48 +132,36 @@ public class LeastSquaresFactory { */ public static LeastSquaresProblem countEvaluations(final LeastSquaresProblem problem, final Incrementor counter) { - //TODO adapter? - return new LeastSquaresProblem() { + return new LeastSquaresAdapter(problem) { - public Evaluation evaluate(double[] point) { + public Evaluation evaluate(final double[] point) { counter.incrementCount(); - return problem.evaluate(point); + return super.evaluate(point); } /* delegate the rest */ - public double[] getStart() { - return problem.getStart(); - } - - public int getObservationSize() { - return problem.getObservationSize(); - } - - public int getParameterSize() { - return problem.getParameterSize(); - } - - public RealMatrix getWeight() { - return problem.getWeight(); - } - - public RealMatrix getWeightSquareRoot() { - return problem.getWeightSquareRoot(); - } - - public Incrementor getEvaluationCounter() { - return problem.getEvaluationCounter(); - } - - public Incrementor getIterationCounter() { - return problem.getIterationCounter(); - } - - public ConvergenceChecker getConvergenceChecker() { - return problem.getConvergenceChecker(); - } }; } + /** + * Computes the square-root of the weight matrix. + * + * @param m Symmetric, positive-definite (weight) matrix. + * @return the square-root of the weight matrix. + */ + private static RealMatrix squareRoot(final RealMatrix m) { + if (m instanceof DiagonalMatrix) { + final int dim = m.getRowDimension(); + final RealMatrix sqrtM = new DiagonalMatrix(dim); + for (int i = 0; i < dim; i++) { + sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i))); + } + return sqrtM; + } else { + final EigenDecomposition dec = new EigenDecomposition(m); + return dec.getSquareRoot(); + } + } } + diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblem.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblem.java index 08b5ac7c5..a8d8a497c 100644 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblem.java +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblem.java @@ -39,22 +39,6 @@ public interface LeastSquaresProblem extends OptimizationProblem - * TODO Is it possible to leave out this method and have the weights implicit in the - * {@link Evaluation}? - * - * @return the weight matrix. - */ - RealMatrix getWeight(); - - /**Get the square root of the weight matrix. - * TODO delete this method - * @return the square root of the weight matrix - */ - RealMatrix getWeightSquareRoot(); - public interface Evaluation { /** @@ -111,14 +95,7 @@ public interface LeastSquaresProblem extends OptimizationProblem1/2 K. * @throws DimensionMismatchException if {@code params} has a wrong length. */ double[] computeResiduals(); - /** * Get the abscissa (independent variables) of this evaluation. * diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblemImpl.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblemImpl.java index 5e11fab88..3cb49d39d 100644 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblemImpl.java +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblemImpl.java @@ -19,17 +19,11 @@ package org.apache.commons.math3.fitting.leastsquares; import org.apache.commons.math3.analysis.MultivariateMatrixFunction; import org.apache.commons.math3.analysis.MultivariateVectorFunction; import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.DecompositionSolver; -import org.apache.commons.math3.linear.DiagonalMatrix; -import org.apache.commons.math3.linear.EigenDecomposition; import org.apache.commons.math3.linear.MatrixUtils; -import org.apache.commons.math3.linear.QRDecomposition; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.optim.AbstractOptimizationProblem; import org.apache.commons.math3.optim.ConvergenceChecker; import org.apache.commons.math3.optim.PointVectorValuePair; -import org.apache.commons.math3.util.FastMath; /** * A private, "field" immutable (not "real" immutable) implementation of {@link @@ -44,14 +38,10 @@ class LeastSquaresProblemImpl /** Target values for the model function at optimum. */ private double[] target; - /** Weight matrix. */ - private RealMatrix weight; /** Model function. */ private MultivariateVectorFunction model; /** Jacobian of the model function. */ private MultivariateMatrixFunction jacobian; - /** Square-root of the weight matrix. */ - private RealMatrix weightSqrt; /** Initial guess. */ private double[] start; @@ -59,16 +49,13 @@ class LeastSquaresProblemImpl final int maxIterations, final ConvergenceChecker checker, final double[] target, - final RealMatrix weight, final MultivariateVectorFunction model, final MultivariateMatrixFunction jacobian, final double[] start) { super(maxEvaluations, maxIterations, checker); this.target = target; - this.weight = weight; this.model = model; this.jacobian = jacobian; - this.weightSqrt = squareRoot(weight); this.start = start; } @@ -93,44 +80,11 @@ class LeastSquaresProblemImpl return start == null ? null : start.clone(); } - /** - * Gets the square-root of the weight matrix. - * - * @return the square-root of the weight matrix. - */ - public RealMatrix getWeightSquareRoot() { - return weightSqrt == null ? null : weightSqrt.copy(); - } - - /** - * Gets the model function. - * - * @return the model function. - */ - public MultivariateVectorFunction getModel() { - return model; - } - - /** - * Gets the model function's Jacobian. - * - * @return the Jacobian. - */ - public MultivariateMatrixFunction getJacobian() { - return jacobian; - } - - public RealMatrix getWeight() { - return weight.copy(); - } - public Evaluation evaluate(final double[] point) { //TODO evaluate value and jacobian in one function call - return new EvaluationImpl( + return new UnweightedEvaluation( this.model.value(point), this.jacobian.value(point), - this.weight, - this.weightSqrt, this.target, point); } @@ -140,7 +94,7 @@ class LeastSquaresProblemImpl *

* TODO revisit lazy evaluation */ - private static class EvaluationImpl implements Evaluation { + private static class UnweightedEvaluation extends AbstractEvaluation { /** the point of evaluation */ private final double[] point; @@ -148,70 +102,32 @@ class LeastSquaresProblemImpl private final double[] values; /** deriviative at point */ private final double[][] jacobian; - /* references to data defined by the least squares problem. Not modified. - * Could be a reference to the problem. - */ - private final RealMatrix weight; - private final RealMatrix weightSqrt; + /** reference to the observed values */ private final double[] target; - private EvaluationImpl(final double[] values, - final double[][] jacobian, - final RealMatrix weight, - final RealMatrix weightSqrt, - final double[] target, - final double[] point) { + private UnweightedEvaluation(final double[] values, + final double[][] jacobian, + final double[] target, + final double[] point) { + super(target.length); this.values = values; this.jacobian = jacobian; - this.weight = weight; - this.weightSqrt = weightSqrt; this.target = target; this.point = point; } - public double[][] computeCovariances(double threshold) { - // Set up the Jacobian. - final RealMatrix j = computeWeightedJacobian(); - - // 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(); - } - - public double[] computeSigma(double covarianceSingularityThreshold) { - final double[][] cov = 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; - } - - public double computeRMS() { - final double cost = computeCost(); - return FastMath.sqrt(cost * cost / target.length); - } public double[] computeValue() { return this.values; } - public RealMatrix computeWeightedJacobian() { - return weightSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian())); + public RealMatrix computeJacobian() { + return MatrixUtils.createRealMatrix(this.jacobian); } - public double[][] computeJacobian() { - return this.jacobian; - } - public double computeCost() { - final ArrayRealVector r = new ArrayRealVector(computeResiduals()); - return FastMath.sqrt(r.dotProduct(weight.operate(r))); + public double[] getPoint() { + return this.point; } public double[] computeResiduals() { @@ -229,29 +145,6 @@ class LeastSquaresProblemImpl return residuals; } - public double[] getPoint() { - //TODO copy? - return this.point; - } } - /** - * Computes the square-root of the weight matrix. - * - * @param m Symmetric, positive-definite (weight) matrix. - * @return the square-root of the weight matrix. - */ - private RealMatrix squareRoot(RealMatrix m) { - if (m instanceof DiagonalMatrix) { - final int dim = m.getRowDimension(); - final RealMatrix sqrtM = new DiagonalMatrix(dim); - for (int i = 0; i < dim; i++) { - sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i))); - } - return sqrtM; - } else { - final EigenDecomposition dec = new EigenDecomposition(m); - return dec.getSquareRoot(); - } - } } diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java index 59855cb6b..13b7831c9 100644 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java @@ -311,7 +311,6 @@ public class LevenbergMarquardtOptimizer implements LeastSquaresOptimizer { //convergence criterion final ConvergenceChecker checker = problem.getConvergenceChecker(); - final RealMatrix weightMatrixSqrt = problem.getWeightSquareRoot(); // arrays shared with the other private methods final int solvedCols = FastMath.min(nR, nC); @@ -350,13 +349,14 @@ public class LevenbergMarquardtOptimizer implements LeastSquaresOptimizer { // QR decomposition of the jacobian matrix final InternalData internalData - = qrDecomposition(value.computeWeightedJacobian(), solvedCols); + = qrDecomposition(value.computeJacobian(), solvedCols); final double[][] weightedJacobian = internalData.weightedJacobian; final int[] permutation = internalData.permutation; final double[] diagR = internalData.diagR; final double[] jacNorm = internalData.jacNorm; - double[] weightedResidual = weightMatrixSqrt.operate(currentResiduals); + //residuals already have weights applied + double[] weightedResidual = currentResiduals; for (int i = 0; i < nR; i++) { qtf[i] = weightedResidual[i]; } diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/OptimumImpl.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/OptimumImpl.java index b0e14013d..589ed3b1d 100644 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/OptimumImpl.java +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/OptimumImpl.java @@ -64,12 +64,7 @@ class OptimumImpl implements Optimum { } /** {@inheritDoc} */ - public RealMatrix computeWeightedJacobian() { - return value.computeWeightedJacobian(); - } - - /** {@inheritDoc} */ - public double[][] computeJacobian() { + public RealMatrix computeJacobian() { return value.computeJacobian(); }