diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/fitting/leastsquares')
15 files changed, 2877 insertions, 0 deletions
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 0000000..b164380 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractEvaluation.java @@ -0,0 +1,87 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +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.linear.RealVector; +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? + * + * @since 3.3 + */ +public abstract class AbstractEvaluation implements Evaluation { + + /** number of observations */ + private final int observationSize; + + /** + * Constructor. + * + * @param observationSize the number of observation. Needed for {@link + * #getRMS()}. + */ + AbstractEvaluation(final int observationSize) { + this.observationSize = observationSize; + } + + /** {@inheritDoc} */ + public RealMatrix getCovariances(double threshold) { + // Set up the Jacobian. + final RealMatrix j = this.getJacobian(); + + // 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(); + } + + /** {@inheritDoc} */ + public RealVector getSigma(double covarianceSingularityThreshold) { + final RealMatrix cov = this.getCovariances(covarianceSingularityThreshold); + final int nC = cov.getColumnDimension(); + final RealVector sig = new ArrayRealVector(nC); + for (int i = 0; i < nC; ++i) { + sig.setEntry(i, FastMath.sqrt(cov.getEntry(i,i))); + } + return sig; + } + + /** {@inheritDoc} */ + public double getRMS() { + final double cost = this.getCost(); + return FastMath.sqrt(cost * cost / this.observationSize); + } + + /** {@inheritDoc} */ + public double getCost() { + final ArrayRealVector r = new ArrayRealVector(this.getResiduals()); + 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 0000000..89f5f1f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/DenseWeightedEvaluation.java @@ -0,0 +1,68 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; + +/** + * Applies a dense weight matrix to an evaluation. + * + * @since 3.3 + */ +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 getJacobian() { + return weightSqrt.multiply(this.unweighted.getJacobian()); + } + + /** {@inheritDoc} */ + public RealVector getResiduals() { + return this.weightSqrt.operate(this.unweighted.getResiduals()); + } + + /* delegate */ + + /** {@inheritDoc} */ + public RealVector getPoint() { + return unweighted.getPoint(); + } + +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/EvaluationRmsChecker.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/EvaluationRmsChecker.java new file mode 100644 index 0000000..ceb5988 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/EvaluationRmsChecker.java @@ -0,0 +1,75 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.util.Precision; + +/** + * Check if an optimization has converged based on the change in computed RMS. + * + * @since 3.4 + */ +public class EvaluationRmsChecker implements ConvergenceChecker<Evaluation> { + + /** relative tolerance for comparisons. */ + private final double relTol; + /** absolute tolerance for comparisons. */ + private final double absTol; + + /** + * Create a convergence checker for the RMS with the same relative and absolute + * tolerance. + * + * <p>Convenience constructor for when the relative and absolute tolerances are the + * same. Same as {@code new EvaluationRmsChecker(tol, tol)}. + * + * @param tol the relative and absolute tolerance. + * @see #EvaluationRmsChecker(double, double) + */ + public EvaluationRmsChecker(final double tol) { + this(tol, tol); + } + + /** + * Create a convergence checker for the RMS with a relative and absolute tolerance. + * + * <p>The optimization has converged when the RMS of consecutive evaluations are equal + * to within the given relative tolerance or absolute tolerance. + * + * @param relTol the relative tolerance. + * @param absTol the absolute tolerance. + * @see Precision#equals(double, double, double) + * @see Precision#equalsWithRelativeTolerance(double, double, double) + */ + public EvaluationRmsChecker(final double relTol, final double absTol) { + this.relTol = relTol; + this.absTol = absTol; + } + + /** {@inheritDoc} */ + public boolean converged(final int iteration, + final Evaluation previous, + final Evaluation current) { + final double prevRms = previous.getRMS(); + final double currRms = current.getRMS(); + return Precision.equals(prevRms, currRms, this.absTol) || + Precision.equalsWithRelativeTolerance(prevRms, currRms, this.relTol); + } + +} 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 new file mode 100644 index 0000000..8157706 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java @@ -0,0 +1,299 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.CholeskyDecomposition; +import org.apache.commons.math3.linear.LUDecomposition; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException; +import org.apache.commons.math3.linear.QRDecomposition; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.apache.commons.math3.linear.SingularMatrixException; +import org.apache.commons.math3.linear.SingularValueDecomposition; +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.util.Incrementor; +import org.apache.commons.math3.util.Pair; + +/** + * Gauss-Newton least-squares solver. + * <p> This class solve a least-square problem by + * solving the normal equations of the linearized problem at each iteration. Either LU + * decomposition or Cholesky decomposition can be used to solve the normal equations, + * or QR decomposition or SVD decomposition can be used to solve the linear system. LU + * decomposition is faster but QR decomposition is more robust for difficult problems, + * and SVD can compute a solution for rank-deficient problems. + * </p> + * + * @since 3.3 + */ +public class GaussNewtonOptimizer implements LeastSquaresOptimizer { + + /** The decomposition algorithm to use to solve the normal equations. */ + //TODO move to linear package and expand options? + public enum Decomposition { + /** + * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and + * using the {@link LUDecomposition}. + * + * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the + * normal matrix and n<sup>3</sup>/3 operations (m > n) to solve the system using + * the LU decomposition. </p> + */ + LU { + @Override + protected RealVector solve(final RealMatrix jacobian, + final RealVector residuals) { + try { + final Pair<RealMatrix, RealVector> normalEquation = + computeNormalMatrix(jacobian, residuals); + final RealMatrix normal = normalEquation.getFirst(); + final RealVector jTr = normalEquation.getSecond(); + return new LUDecomposition(normal, SINGULARITY_THRESHOLD) + .getSolver() + .solve(jTr); + } catch (SingularMatrixException e) { + throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); + } + } + }, + /** + * Solve the linear least squares problem (Jx=r) using the {@link + * QRDecomposition}. + * + * <p> Theoretically this method takes mn<sup>2</sup> - n<sup>3</sup>/3 operations + * (m > n) and has better numerical accuracy than any method that forms the normal + * equations. </p> + */ + QR { + @Override + protected RealVector solve(final RealMatrix jacobian, + final RealVector residuals) { + try { + return new QRDecomposition(jacobian, SINGULARITY_THRESHOLD) + .getSolver() + .solve(residuals); + } catch (SingularMatrixException e) { + throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); + } + } + }, + /** + * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and + * using the {@link CholeskyDecomposition}. + * + * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the + * normal matrix and n<sup>3</sup>/6 operations (m > n) to solve the system using + * the Cholesky decomposition. </p> + */ + CHOLESKY { + @Override + protected RealVector solve(final RealMatrix jacobian, + final RealVector residuals) { + try { + final Pair<RealMatrix, RealVector> normalEquation = + computeNormalMatrix(jacobian, residuals); + final RealMatrix normal = normalEquation.getFirst(); + final RealVector jTr = normalEquation.getSecond(); + return new CholeskyDecomposition( + normal, SINGULARITY_THRESHOLD, SINGULARITY_THRESHOLD) + .getSolver() + .solve(jTr); + } catch (NonPositiveDefiniteMatrixException e) { + throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); + } + } + }, + /** + * Solve the linear least squares problem using the {@link + * SingularValueDecomposition}. + * + * <p> This method is slower, but can provide a solution for rank deficient and + * nearly singular systems. + */ + SVD { + @Override + protected RealVector solve(final RealMatrix jacobian, + final RealVector residuals) { + return new SingularValueDecomposition(jacobian) + .getSolver() + .solve(residuals); + } + }; + + /** + * Solve the linear least squares problem Jx=r. + * + * @param jacobian the Jacobian matrix, J. the number of rows >= the number or + * columns. + * @param residuals the computed residuals, r. + * @return the solution x, to the linear least squares problem Jx=r. + * @throws ConvergenceException if the matrix properties (e.g. singular) do not + * permit a solution. + */ + protected abstract RealVector solve(RealMatrix jacobian, + RealVector residuals); + } + + /** + * The singularity threshold for matrix decompositions. Determines when a {@link + * ConvergenceException} is thrown. The current value was the default value for {@link + * LUDecomposition}. + */ + private static final double SINGULARITY_THRESHOLD = 1e-11; + + /** Indicator for using LU decomposition. */ + private final Decomposition decomposition; + + /** + * Creates a Gauss Newton optimizer. + * <p/> + * The default for the algorithm is to solve the normal equations using QR + * decomposition. + */ + public GaussNewtonOptimizer() { + this(Decomposition.QR); + } + + /** + * Create a Gauss Newton optimizer that uses the given decomposition algorithm to + * solve the normal equations. + * + * @param decomposition the {@link Decomposition} algorithm. + */ + public GaussNewtonOptimizer(final Decomposition decomposition) { + this.decomposition = decomposition; + } + + /** + * Get the matrix decomposition algorithm used to solve the normal equations. + * + * @return the matrix {@link Decomposition} algoritm. + */ + public Decomposition getDecomposition() { + return this.decomposition; + } + + /** + * Configure the decomposition algorithm. + * + * @param newDecomposition the {@link Decomposition} algorithm to use. + * @return a new instance. + */ + public GaussNewtonOptimizer withDecomposition(final Decomposition newDecomposition) { + return new GaussNewtonOptimizer(newDecomposition); + } + + /** {@inheritDoc} */ + public Optimum optimize(final LeastSquaresProblem lsp) { + //create local evaluation and iteration counts + final Incrementor evaluationCounter = lsp.getEvaluationCounter(); + final Incrementor iterationCounter = lsp.getIterationCounter(); + final ConvergenceChecker<Evaluation> checker + = lsp.getConvergenceChecker(); + + // Computation will be useless without a checker (see "for-loop"). + if (checker == null) { + throw new NullArgumentException(); + } + + RealVector currentPoint = lsp.getStart(); + + // iterate until convergence is reached + Evaluation current = null; + while (true) { + iterationCounter.incrementCount(); + + // evaluate the objective function and its jacobian + Evaluation previous = current; + // Value of the objective function at "currentPoint". + evaluationCounter.incrementCount(); + current = lsp.evaluate(currentPoint); + final RealVector currentResiduals = current.getResiduals(); + final RealMatrix weightedJacobian = current.getJacobian(); + currentPoint = current.getPoint(); + + // Check convergence. + if (previous != null && + checker.converged(iterationCounter.getCount(), previous, current)) { + return new OptimumImpl(current, + evaluationCounter.getCount(), + iterationCounter.getCount()); + } + + // solve the linearized least squares problem + final RealVector dX = this.decomposition.solve(weightedJacobian, currentResiduals); + // update the estimated parameters + currentPoint = currentPoint.add(dX); + } + } + + /** {@inheritDoc} */ + @Override + public String toString() { + return "GaussNewtonOptimizer{" + + "decomposition=" + decomposition + + '}'; + } + + /** + * Compute the normal matrix, J<sup>T</sup>J. + * + * @param jacobian the m by n jacobian matrix, J. Input. + * @param residuals the m by 1 residual vector, r. Input. + * @return the n by n normal matrix and the n by 1 J<sup>Tr vector. + */ + private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian, + final RealVector residuals) { + //since the normal matrix is symmetric, we only need to compute half of it. + final int nR = jacobian.getRowDimension(); + final int nC = jacobian.getColumnDimension(); + //allocate space for return values + final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC); + final RealVector jTr = new ArrayRealVector(nC); + //for each measurement + for (int i = 0; i < nR; ++i) { + //compute JTr for measurement i + for (int j = 0; j < nC; j++) { + jTr.setEntry(j, jTr.getEntry(j) + + residuals.getEntry(i) * jacobian.getEntry(i, j)); + } + + // add the the contribution to the normal matrix for measurement i + for (int k = 0; k < nC; ++k) { + //only compute the upper triangular part + for (int l = k; l < nC; ++l) { + normal.setEntry(k, l, normal.getEntry(k, l) + + jacobian.getEntry(i, k) * jacobian.getEntry(i, l)); + } + } + } + //copy the upper triangular part to the lower triangular part. + for (int i = 0; i < nC; i++) { + for (int j = 0; j < i; j++) { + normal.setEntry(i, j, normal.getEntry(j, i)); + } + } + return new Pair<RealMatrix, RealVector>(normal, jTr); + } + +} 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 0000000..1c09874 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresAdapter.java @@ -0,0 +1,77 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.linear.RealVector; +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.util.Incrementor; + +/** + * An adapter that delegates to another implementation of {@link LeastSquaresProblem}. + * + * @since 3.3 + */ +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 RealVector getStart() { + return problem.getStart(); + } + + /** {@inheritDoc} */ + public int getObservationSize() { + return problem.getObservationSize(); + } + + /** {@inheritDoc} */ + public int getParameterSize() { + return problem.getParameterSize(); + } + + /** {@inheritDoc} + * @param point*/ + public Evaluation evaluate(final RealVector point) { + return problem.evaluate(point); + } + + /** {@inheritDoc} */ + public Incrementor getEvaluationCounter() { + return problem.getEvaluationCounter(); + } + + /** {@inheritDoc} */ + public Incrementor getIterationCounter() { + return problem.getIterationCounter(); + } + + /** {@inheritDoc} */ + public ConvergenceChecker<Evaluation> getConvergenceChecker() { + return problem.getConvergenceChecker(); + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresBuilder.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresBuilder.java new file mode 100644 index 0000000..7b14b37 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresBuilder.java @@ -0,0 +1,226 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +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.fitting.leastsquares.LeastSquaresProblem.Evaluation; +import org.apache.commons.math3.linear.ArrayRealVector; +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; + +/** + * A mutable builder for {@link LeastSquaresProblem}s. + * + * @see LeastSquaresFactory + * @since 3.3 + */ +public class LeastSquaresBuilder { + + /** max evaluations */ + private int maxEvaluations; + /** max iterations */ + private int maxIterations; + /** convergence checker */ + private ConvergenceChecker<Evaluation> checker; + /** model function */ + private MultivariateJacobianFunction model; + /** observed values */ + private RealVector target; + /** initial guess */ + private RealVector start; + /** weight matrix */ + private RealMatrix weight; + /** + * Lazy evaluation. + * + * @since 3.4 + */ + private boolean lazyEvaluation; + /** Validator. + * + * @since 3.4 + */ + private ParameterValidator paramValidator; + + + /** + * Construct a {@link LeastSquaresProblem} from the data in this builder. + * + * @return a new {@link LeastSquaresProblem}. + */ + public LeastSquaresProblem build() { + return LeastSquaresFactory.create(model, + target, + start, + weight, + checker, + maxEvaluations, + maxIterations, + lazyEvaluation, + paramValidator); + } + + /** + * Configure the max evaluations. + * + * @param newMaxEvaluations the maximum number of evaluations permitted. + * @return this + */ + public LeastSquaresBuilder maxEvaluations(final int newMaxEvaluations) { + this.maxEvaluations = newMaxEvaluations; + return this; + } + + /** + * Configure the max iterations. + * + * @param newMaxIterations the maximum number of iterations permitted. + * @return this + */ + public LeastSquaresBuilder maxIterations(final int newMaxIterations) { + this.maxIterations = newMaxIterations; + return this; + } + + /** + * Configure the convergence checker. + * + * @param newChecker the convergence checker. + * @return this + */ + public LeastSquaresBuilder checker(final ConvergenceChecker<Evaluation> newChecker) { + this.checker = newChecker; + return this; + } + + /** + * Configure the convergence checker. + * <p/> + * This function is an overloaded version of {@link #checker(ConvergenceChecker)}. + * + * @param newChecker the convergence checker. + * @return this + */ + public LeastSquaresBuilder checkerPair(final ConvergenceChecker<PointVectorValuePair> newChecker) { + return this.checker(LeastSquaresFactory.evaluationChecker(newChecker)); + } + + /** + * Configure the model function. + * + * @param value the model function value + * @param jacobian the Jacobian of {@code value} + * @return this + */ + public LeastSquaresBuilder model(final MultivariateVectorFunction value, + final MultivariateMatrixFunction jacobian) { + return model(LeastSquaresFactory.model(value, jacobian)); + } + + /** + * Configure the model function. + * + * @param newModel the model function value and Jacobian + * @return this + */ + public LeastSquaresBuilder model(final MultivariateJacobianFunction newModel) { + this.model = newModel; + return this; + } + + /** + * Configure the observed data. + * + * @param newTarget the observed data. + * @return this + */ + public LeastSquaresBuilder target(final RealVector newTarget) { + this.target = newTarget; + return this; + } + + /** + * Configure the observed data. + * + * @param newTarget the observed data. + * @return this + */ + public LeastSquaresBuilder target(final double[] newTarget) { + return target(new ArrayRealVector(newTarget, false)); + } + + /** + * Configure the initial guess. + * + * @param newStart the initial guess. + * @return this + */ + public LeastSquaresBuilder start(final RealVector newStart) { + this.start = newStart; + return this; + } + + /** + * Configure the initial guess. + * + * @param newStart the initial guess. + * @return this + */ + public LeastSquaresBuilder start(final double[] newStart) { + return start(new ArrayRealVector(newStart, false)); + } + + /** + * Configure the weight matrix. + * + * @param newWeight the weight matrix + * @return this + */ + public LeastSquaresBuilder weight(final RealMatrix newWeight) { + this.weight = newWeight; + return this; + } + + /** + * Configure whether evaluation will be lazy or not. + * + * @param newValue Whether to perform lazy evaluation. + * @return this object. + * + * @since 3.4 + */ + public LeastSquaresBuilder lazyEvaluation(final boolean newValue) { + lazyEvaluation = newValue; + return this; + } + + /** + * Configure the validator of the model parameters. + * + * @param newValidator Parameter validator. + * @return this object. + * + * @since 3.4 + */ + public LeastSquaresBuilder parameterValidator(final ParameterValidator newValidator) { + paramValidator = newValidator; + return this; + } +} 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 new file mode 100644 index 0000000..42cdf89 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresFactory.java @@ -0,0 +1,532 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.ArrayRealVector; +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.AbstractOptimizationProblem; +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; +import org.apache.commons.math3.util.Pair; + +/** + * A Factory for creating {@link LeastSquaresProblem}s. + * + * @since 3.3 + */ +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 (unit weights). + * + * @param model the model function. Produces the computed values. + * @param observed the observed (target) values + * @param start the initial guess. + * @param weight the weight matrix + * @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 + * @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)} + * will defer the evaluation until access to the value is requested. + * @param paramValidator Model parameters validator. + * @return the specified General Least Squares problem. + * + * @since 3.4 + */ + public static LeastSquaresProblem create(final MultivariateJacobianFunction model, + final RealVector observed, + final RealVector start, + final RealMatrix weight, + final ConvergenceChecker<Evaluation> checker, + final int maxEvaluations, + final int maxIterations, + final boolean lazyEvaluation, + final ParameterValidator paramValidator) { + final LeastSquaresProblem p = new LocalLeastSquaresProblem(model, + observed, + start, + checker, + maxEvaluations, + maxIterations, + lazyEvaluation, + paramValidator); + if (weight != null) { + return weightMatrix(p, weight); + } else { + return p; + } + } + + /** + * Create a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem} + * from the given elements. There will be no weights applied (unit weights). + * + * @param model the model function. Produces the computed values. + * @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 MultivariateJacobianFunction model, + final RealVector observed, + final RealVector start, + final ConvergenceChecker<Evaluation> checker, + final int maxEvaluations, + final int maxIterations) { + return create(model, + observed, + start, + null, + checker, + maxEvaluations, + maxIterations, + false, + null); + } + + /** + * Create a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem} + * from the given elements. + * + * @param model the model function. Produces the computed values. + * @param observed the observed (target) values + * @param start the initial guess. + * @param weight the weight matrix + * @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 MultivariateJacobianFunction model, + final RealVector observed, + final RealVector start, + final RealMatrix weight, + final ConvergenceChecker<Evaluation> checker, + final int maxEvaluations, + final int maxIterations) { + return weightMatrix(create(model, + observed, + start, + checker, + maxEvaluations, + maxIterations), + weight); + } + + /** + * Create a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem} + * from the given elements. + * <p> + * This factory method is provided for continuity with previous interfaces. Newer + * applications should use {@link #create(MultivariateJacobianFunction, RealVector, + * RealVector, ConvergenceChecker, int, int)}, or {@link #create(MultivariateJacobianFunction, + * RealVector, RealVector, RealMatrix, ConvergenceChecker, int, int)}. + * + * @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 weight the weight matrix + * @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 RealMatrix weight, + final ConvergenceChecker<Evaluation> checker, + final int maxEvaluations, + final int maxIterations) { + return create(model(model, jacobian), + new ArrayRealVector(observed, false), + new ArrayRealVector(start, false), + weight, + checker, + maxEvaluations, + maxIterations); + } + + /** + * 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) { + /** {@inheritDoc} */ + @Override + public Evaluation evaluate(final RealVector point) { + return new DenseWeightedEvaluation(super.evaluate(point), weightSquareRoot); + } + }; + } + + /** + * Apply a diagonal 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())); + } + + /** + * Count the evaluations of a particular problem. The {@code counter} will be + * incremented every time {@link LeastSquaresProblem#evaluate(RealVector)} is called on + * the <em>returned</em> problem. + * + * @param problem the problem to track. + * @param counter the counter to increment. + * @return a least squares problem that tracks evaluations + */ + public static LeastSquaresProblem countEvaluations(final LeastSquaresProblem problem, + final Incrementor counter) { + return new LeastSquaresAdapter(problem) { + + /** {@inheritDoc} */ + @Override + public Evaluation evaluate(final RealVector point) { + counter.incrementCount(); + return super.evaluate(point); + } + + // Delegate the rest. + }; + } + + /** + * View a convergence checker specified for a {@link PointVectorValuePair} as one + * specified for an {@link Evaluation}. + * + * @param checker the convergence checker to adapt. + * @return a convergence checker that delegates to {@code checker}. + */ + public static ConvergenceChecker<Evaluation> evaluationChecker(final ConvergenceChecker<PointVectorValuePair> checker) { + return new ConvergenceChecker<Evaluation>() { + /** {@inheritDoc} */ + public boolean converged(final int iteration, + final Evaluation previous, + final Evaluation current) { + return checker.converged( + iteration, + new PointVectorValuePair( + previous.getPoint().toArray(), + previous.getResiduals().toArray(), + false), + new PointVectorValuePair( + current.getPoint().toArray(), + current.getResiduals().toArray(), + false) + ); + } + }; + } + + /** + * 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(); + } + } + + /** + * Combine a {@link MultivariateVectorFunction} with a {@link + * MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}. + * + * @param value the vector value function + * @param jacobian the Jacobian function + * @return a function that computes both at the same time + */ + public static MultivariateJacobianFunction model(final MultivariateVectorFunction value, + final MultivariateMatrixFunction jacobian) { + return new LocalValueAndJacobianFunction(value, jacobian); + } + + /** + * Combine a {@link MultivariateVectorFunction} with a {@link + * MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}. + * + * @param value the vector value function + * @param jacobian the Jacobian function + * @return a function that computes both at the same time + */ + private static class LocalValueAndJacobianFunction + implements ValueAndJacobianFunction { + /** Model. */ + private final MultivariateVectorFunction value; + /** Model's Jacobian. */ + private final MultivariateMatrixFunction jacobian; + + /** + * @param value Model function. + * @param jacobian Model's Jacobian function. + */ + LocalValueAndJacobianFunction(final MultivariateVectorFunction value, + final MultivariateMatrixFunction jacobian) { + this.value = value; + this.jacobian = jacobian; + } + + /** {@inheritDoc} */ + public Pair<RealVector, RealMatrix> value(final RealVector point) { + //TODO get array from RealVector without copying? + final double[] p = point.toArray(); + + // Evaluate. + return new Pair<RealVector, RealMatrix>(computeValue(p), + computeJacobian(p)); + } + + /** {@inheritDoc} */ + public RealVector computeValue(final double[] params) { + return new ArrayRealVector(value.value(params), false); + } + + /** {@inheritDoc} */ + public RealMatrix computeJacobian(final double[] params) { + return new Array2DRowRealMatrix(jacobian.value(params), false); + } + } + + + /** + * A private, "field" immutable (not "real" immutable) implementation of {@link + * LeastSquaresProblem}. + * @since 3.3 + */ + private static class LocalLeastSquaresProblem + extends AbstractOptimizationProblem<Evaluation> + implements LeastSquaresProblem { + + /** Target values for the model function at optimum. */ + private final RealVector target; + /** Model function. */ + private final MultivariateJacobianFunction model; + /** Initial guess. */ + private final RealVector start; + /** Whether to use lazy evaluation. */ + private final boolean lazyEvaluation; + /** Model parameters validator. */ + private final ParameterValidator paramValidator; + + /** + * Create a {@link LeastSquaresProblem} from the given data. + * + * @param model the model function + * @param target the observed data + * @param start the initial guess + * @param checker the convergence checker + * @param maxEvaluations the allowed evaluations + * @param maxIterations the allowed iterations + * @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)} + * will defer the evaluation until access to the value is requested. + * @param paramValidator Model parameters validator. + */ + LocalLeastSquaresProblem(final MultivariateJacobianFunction model, + final RealVector target, + final RealVector start, + final ConvergenceChecker<Evaluation> checker, + final int maxEvaluations, + final int maxIterations, + final boolean lazyEvaluation, + final ParameterValidator paramValidator) { + super(maxEvaluations, maxIterations, checker); + this.target = target; + this.model = model; + this.start = start; + this.lazyEvaluation = lazyEvaluation; + this.paramValidator = paramValidator; + + if (lazyEvaluation && + !(model instanceof ValueAndJacobianFunction)) { + // Lazy evaluation requires that value and Jacobian + // can be computed separately. + throw new MathIllegalStateException(LocalizedFormats.INVALID_IMPLEMENTATION, + model.getClass().getName()); + } + } + + /** {@inheritDoc} */ + public int getObservationSize() { + return target.getDimension(); + } + + /** {@inheritDoc} */ + public int getParameterSize() { + return start.getDimension(); + } + + /** {@inheritDoc} */ + public RealVector getStart() { + return start == null ? null : start.copy(); + } + + /** {@inheritDoc} */ + public Evaluation evaluate(final RealVector point) { + // Copy so optimizer can change point without changing our instance. + final RealVector p = paramValidator == null ? + point.copy() : + paramValidator.validate(point.copy()); + + if (lazyEvaluation) { + return new LazyUnweightedEvaluation((ValueAndJacobianFunction) model, + target, + p); + } else { + // Evaluate value and jacobian in one function call. + final Pair<RealVector, RealMatrix> value = model.value(p); + return new UnweightedEvaluation(value.getFirst(), + value.getSecond(), + target, + p); + } + } + + /** + * Container with the model evaluation at a particular point. + */ + private static class UnweightedEvaluation extends AbstractEvaluation { + /** Point of evaluation. */ + private final RealVector point; + /** Derivative at point. */ + private final RealMatrix jacobian; + /** Computed residuals. */ + private final RealVector residuals; + + /** + * Create an {@link Evaluation} with no weights. + * + * @param values the computed function values + * @param jacobian the computed function Jacobian + * @param target the observed values + * @param point the abscissa + */ + private UnweightedEvaluation(final RealVector values, + final RealMatrix jacobian, + final RealVector target, + final RealVector point) { + super(target.getDimension()); + this.jacobian = jacobian; + this.point = point; + this.residuals = target.subtract(values); + } + + /** {@inheritDoc} */ + public RealMatrix getJacobian() { + return jacobian; + } + + /** {@inheritDoc} */ + public RealVector getPoint() { + return point; + } + + /** {@inheritDoc} */ + public RealVector getResiduals() { + return residuals; + } + } + + /** + * Container with the model <em>lazy</em> evaluation at a particular point. + */ + private static class LazyUnweightedEvaluation extends AbstractEvaluation { + /** Point of evaluation. */ + private final RealVector point; + /** Model and Jacobian functions. */ + private final ValueAndJacobianFunction model; + /** Target values for the model function at optimum. */ + private final RealVector target; + + /** + * Create an {@link Evaluation} with no weights. + * + * @param model the model function + * @param target the observed values + * @param point the abscissa + */ + private LazyUnweightedEvaluation(final ValueAndJacobianFunction model, + final RealVector target, + final RealVector point) { + super(target.getDimension()); + // Safe to cast as long as we control usage of this class. + this.model = model; + this.point = point; + this.target = target; + } + + /** {@inheritDoc} */ + public RealMatrix getJacobian() { + return model.computeJacobian(point.toArray()); + } + + /** {@inheritDoc} */ + public RealVector getPoint() { + return point; + } + + /** {@inheritDoc} */ + public RealVector getResiduals() { + return target.subtract(model.computeValue(point.toArray())); + } + } + } +} + diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresOptimizer.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresOptimizer.java new file mode 100644 index 0000000..50d5b8a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresOptimizer.java @@ -0,0 +1,62 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +/** + * An algorithm that can be applied to a non-linear least squares problem. + * + * @since 3.3 + */ +public interface LeastSquaresOptimizer { + + /** + * Solve the non-linear least squares problem. + * + * + * @param leastSquaresProblem the problem definition, including model function and + * convergence criteria. + * @return The optimum. + */ + Optimum optimize(LeastSquaresProblem leastSquaresProblem); + + /** + * The optimum found by the optimizer. This object contains the point, its value, and + * some metadata. + */ + //TODO Solution? + interface Optimum extends LeastSquaresProblem.Evaluation { + + /** + * Get the number of times the model was evaluated in order to produce this + * optimum. + * + * @return the number of model (objective) function evaluations + */ + int getEvaluations(); + + /** + * Get the number of times the algorithm iterated in order to produce this + * optimum. In general least squares it is common to have one {@link + * #getEvaluations() evaluation} per iterations. + * + * @return the number of iterations + */ + int getIterations(); + + } + +} 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 new file mode 100644 index 0000000..097ff81 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblem.java @@ -0,0 +1,156 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.apache.commons.math3.optim.OptimizationProblem; + +/** + * The data necessary to define a non-linear least squares problem. + * <p> + * Includes the observed values, computed model function, and + * convergence/divergence criteria. Weights are implicit in {@link + * Evaluation#getResiduals()} and {@link Evaluation#getJacobian()}. + * </p> + * <p> + * Instances are typically either created progressively using a {@link + * LeastSquaresBuilder builder} or created at once using a {@link LeastSquaresFactory + * factory}. + * </p> + * @see LeastSquaresBuilder + * @see LeastSquaresFactory + * @see LeastSquaresAdapter + * + * @since 3.3 + */ +public interface LeastSquaresProblem extends OptimizationProblem<LeastSquaresProblem.Evaluation> { + + /** + * Gets the initial guess. + * + * @return the initial guess values. + */ + RealVector getStart(); + + /** + * Get the number of observations (rows in the Jacobian) in this problem. + * + * @return the number of scalar observations + */ + int getObservationSize(); + + /** + * Get the number of parameters (columns in the Jacobian) in this problem. + * + * @return the number of scalar parameters + */ + int getParameterSize(); + + /** + * Evaluate the model at the specified point. + * + * + * @param point the parameter values. + * @return the model's value and derivative at the given point. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximal number of evaluations (of the model vector function) is + * exceeded. + */ + Evaluation evaluate(RealVector point); + + /** + * An evaluation of a {@link LeastSquaresProblem} at a particular point. This class + * also computes several quantities derived from the value and its Jacobian. + */ + public interface Evaluation { + + /** + * Get the covariance matrix of the optimized parameters. <br/> Note that this + * operation involves the inversion of the <code>J<sup>T</sup>J</code> matrix, + * where {@code J} is the Jacobian matrix. The {@code threshold} parameter is a + * way for the caller to specify that the result of this computation should be + * considered meaningless, and thus trigger an exception. + * + * + * @param threshold Singularity threshold. + * @return the covariance matrix. + * @throws org.apache.commons.math3.linear.SingularMatrixException + * if the covariance matrix cannot be computed (singular problem). + */ + RealMatrix getCovariances(double threshold); + + /** + * Get an estimate of the standard deviation of the parameters. The returned + * values are the square root of the diagonal coefficients of the covariance + * matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]} is the optimized + * value of the {@code i}-th parameter, and {@code C} is the covariance matrix. + * + * + * @param covarianceSingularityThreshold Singularity threshold (see {@link + * #getCovariances(double) computeCovariances}). + * @return an estimate of the standard deviation of the optimized parameters + * @throws org.apache.commons.math3.linear.SingularMatrixException + * if the covariance matrix cannot be computed. + */ + RealVector getSigma(double covarianceSingularityThreshold); + + /** + * Get the normalized cost. It is the square-root of the sum of squared of + * the residuals, divided by the number of measurements. + * + * @return the cost. + */ + double getRMS(); + + /** + * Get the weighted Jacobian matrix. + * + * @return the weighted Jacobian: W<sup>1/2</sup> J. + * @throws org.apache.commons.math3.exception.DimensionMismatchException + * if the Jacobian dimension does not match problem dimension. + */ + RealMatrix getJacobian(); + + /** + * Get the cost. + * + * @return the cost. + * @see #getResiduals() + */ + double getCost(); + + /** + * Get the weighted residuals. The residual is the difference between the + * observed (target) values and the model (objective function) value. There is one + * residual for each element of the vector-valued function. The raw residuals are + * then multiplied by the square root of the weight matrix. + * + * @return the weighted residuals: W<sup>1/2</sup> K. + * @throws org.apache.commons.math3.exception.DimensionMismatchException + * if the residuals have the wrong length. + */ + RealVector getResiduals(); + + /** + * Get the abscissa (independent variables) of this evaluation. + * + * @return the point provided to {@link #evaluate(RealVector)}. + */ + RealVector getPoint(); + } +} 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 new file mode 100644 index 0000000..358d240 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java @@ -0,0 +1,1042 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import java.util.Arrays; + +import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.util.Incrementor; +import org.apache.commons.math3.util.Precision; +import org.apache.commons.math3.util.FastMath; + + +/** + * This class solves a least-squares problem using the Levenberg-Marquardt + * algorithm. + * + * <p>This implementation <em>should</em> work even for over-determined systems + * (i.e. systems having more point than equations). Over-determined systems + * are solved by ignoring the point which have the smallest impact according + * to their jacobian column norm. Only the rank of the matrix and some loop bounds + * are changed to implement this.</p> + * + * <p>The resolution engine is a simple translation of the MINPACK <a + * href="http://www.netlib.org/minpack/lmder.f">lmder</a> routine with minor + * changes. The changes include the over-determined resolution, the use of + * inherited convergence checker and the Q.R. decomposition which has been + * rewritten following the algorithm described in the + * P. Lascaux and R. Theodor book <i>Analyse numérique matricielle + * appliquée à l'art de l'ingénieur</i>, Masson 1986.</p> + * <p>The authors of the original fortran version are: + * <ul> + * <li>Argonne National Laboratory. MINPACK project. March 1980</li> + * <li>Burton S. Garbow</li> + * <li>Kenneth E. Hillstrom</li> + * <li>Jorge J. More</li> + * </ul> + * The redistribution policy for MINPACK is available <a + * href="http://www.netlib.org/minpack/disclaimer">here</a>, for convenience, it + * is reproduced below.</p> + * + * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> + * <tr><td> + * Minpack Copyright Notice (1999) University of Chicago. + * All rights reserved + * </td></tr> + * <tr><td> + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * <ol> + * <li>Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer.</li> + * <li>Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following + * disclaimer in the documentation and/or other materials provided + * with the distribution.</li> + * <li>The end-user documentation included with the redistribution, if any, + * must include the following acknowledgment: + * <code>This product includes software developed by the University of + * Chicago, as Operator of Argonne National Laboratory.</code> + * Alternately, this acknowledgment may appear in the software itself, + * if and wherever such third-party acknowledgments normally appear.</li> + * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" + * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE + * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND + * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE + * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY + * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR + * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF + * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) + * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION + * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL + * BE CORRECTED.</strong></li> + * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT + * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF + * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, + * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF + * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF + * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER + * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT + * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, + * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE + * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li> + * <ol></td></tr> + * </table> + * + * @since 3.3 + */ +public class LevenbergMarquardtOptimizer implements LeastSquaresOptimizer { + + /** Twice the "epsilon machine". */ + private static final double TWO_EPS = 2 * Precision.EPSILON; + + /* configuration parameters */ + /** Positive input variable used in determining the initial step bound. */ + private final double initialStepBoundFactor; + /** Desired relative error in the sum of squares. */ + private final double costRelativeTolerance; + /** Desired relative error in the approximate solution parameters. */ + private final double parRelativeTolerance; + /** Desired max cosine on the orthogonality between the function vector + * and the columns of the jacobian. */ + private final double orthoTolerance; + /** Threshold for QR ranking. */ + private final double qrRankingThreshold; + + /** Default constructor. + * <p> + * The default values for the algorithm settings are: + * <ul> + * <li>Initial step bound factor: 100</li> + * <li>Cost relative tolerance: 1e-10</li> + * <li>Parameters relative tolerance: 1e-10</li> + * <li>Orthogonality tolerance: 1e-10</li> + * <li>QR ranking threshold: {@link Precision#SAFE_MIN}</li> + * </ul> + **/ + public LevenbergMarquardtOptimizer() { + this(100, 1e-10, 1e-10, 1e-10, Precision.SAFE_MIN); + } + + /** + * Construct an instance with all parameters specified. + * + * @param initialStepBoundFactor initial step bound factor + * @param costRelativeTolerance cost relative tolerance + * @param parRelativeTolerance parameters relative tolerance + * @param orthoTolerance orthogonality tolerance + * @param qrRankingThreshold threshold in the QR decomposition. Columns with a 2 + * norm less than this threshold are considered to be + * all 0s. + */ + public LevenbergMarquardtOptimizer( + final double initialStepBoundFactor, + final double costRelativeTolerance, + final double parRelativeTolerance, + final double orthoTolerance, + final double qrRankingThreshold) { + this.initialStepBoundFactor = initialStepBoundFactor; + this.costRelativeTolerance = costRelativeTolerance; + this.parRelativeTolerance = parRelativeTolerance; + this.orthoTolerance = orthoTolerance; + this.qrRankingThreshold = qrRankingThreshold; + } + + /** + * @param newInitialStepBoundFactor Positive input variable used in + * determining the initial step bound. This bound is set to the + * product of initialStepBoundFactor and the euclidean norm of + * {@code diag * x} if non-zero, or else to {@code newInitialStepBoundFactor} + * itself. In most cases factor should lie in the interval + * {@code (0.1, 100.0)}. {@code 100} is a generally recommended value. + * of the matrix is reduced. + * @return a new instance. + */ + public LevenbergMarquardtOptimizer withInitialStepBoundFactor(double newInitialStepBoundFactor) { + return new LevenbergMarquardtOptimizer( + newInitialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** + * @param newCostRelativeTolerance Desired relative error in the sum of squares. + * @return a new instance. + */ + public LevenbergMarquardtOptimizer withCostRelativeTolerance(double newCostRelativeTolerance) { + return new LevenbergMarquardtOptimizer( + initialStepBoundFactor, + newCostRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** + * @param newParRelativeTolerance Desired relative error in the approximate solution + * parameters. + * @return a new instance. + */ + public LevenbergMarquardtOptimizer withParameterRelativeTolerance(double newParRelativeTolerance) { + return new LevenbergMarquardtOptimizer( + initialStepBoundFactor, + costRelativeTolerance, + newParRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** + * Modifies the given parameter. + * + * @param newOrthoTolerance Desired max cosine on the orthogonality between + * the function vector and the columns of the Jacobian. + * @return a new instance. + */ + public LevenbergMarquardtOptimizer withOrthoTolerance(double newOrthoTolerance) { + return new LevenbergMarquardtOptimizer( + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + newOrthoTolerance, + qrRankingThreshold); + } + + /** + * @param newQRRankingThreshold Desired threshold for QR ranking. + * If the squared norm of a column vector is smaller or equal to this + * threshold during QR decomposition, it is considered to be a zero vector + * and hence the rank of the matrix is reduced. + * @return a new instance. + */ + public LevenbergMarquardtOptimizer withRankingThreshold(double newQRRankingThreshold) { + return new LevenbergMarquardtOptimizer( + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + newQRRankingThreshold); + } + + /** + * Gets the value of a tuning parameter. + * @see #withInitialStepBoundFactor(double) + * + * @return the parameter's value. + */ + public double getInitialStepBoundFactor() { + return initialStepBoundFactor; + } + + /** + * Gets the value of a tuning parameter. + * @see #withCostRelativeTolerance(double) + * + * @return the parameter's value. + */ + public double getCostRelativeTolerance() { + return costRelativeTolerance; + } + + /** + * Gets the value of a tuning parameter. + * @see #withParameterRelativeTolerance(double) + * + * @return the parameter's value. + */ + public double getParameterRelativeTolerance() { + return parRelativeTolerance; + } + + /** + * Gets the value of a tuning parameter. + * @see #withOrthoTolerance(double) + * + * @return the parameter's value. + */ + public double getOrthoTolerance() { + return orthoTolerance; + } + + /** + * Gets the value of a tuning parameter. + * @see #withRankingThreshold(double) + * + * @return the parameter's value. + */ + public double getRankingThreshold() { + return qrRankingThreshold; + } + + /** {@inheritDoc} */ + public Optimum optimize(final LeastSquaresProblem problem) { + // Pull in relevant data from the problem as locals. + final int nR = problem.getObservationSize(); // Number of observed data. + final int nC = problem.getParameterSize(); // Number of parameters. + // Counters. + final Incrementor iterationCounter = problem.getIterationCounter(); + final Incrementor evaluationCounter = problem.getEvaluationCounter(); + // Convergence criterion. + final ConvergenceChecker<Evaluation> checker = problem.getConvergenceChecker(); + + // arrays shared with the other private methods + final int solvedCols = FastMath.min(nR, nC); + /* Parameters evolution direction associated with lmPar. */ + double[] lmDir = new double[nC]; + /* Levenberg-Marquardt parameter. */ + double lmPar = 0; + + // local point + double delta = 0; + double xNorm = 0; + double[] diag = new double[nC]; + double[] oldX = new double[nC]; + double[] oldRes = new double[nR]; + double[] qtf = new double[nR]; + double[] work1 = new double[nC]; + double[] work2 = new double[nC]; + double[] work3 = new double[nC]; + + + // Evaluate the function at the starting point and calculate its norm. + evaluationCounter.incrementCount(); + //value will be reassigned in the loop + Evaluation current = problem.evaluate(problem.getStart()); + double[] currentResiduals = current.getResiduals().toArray(); + double currentCost = current.getCost(); + double[] currentPoint = current.getPoint().toArray(); + + // Outer loop. + boolean firstIteration = true; + while (true) { + iterationCounter.incrementCount(); + + final Evaluation previous = current; + + // QR decomposition of the jacobian matrix + final InternalData internalData + = qrDecomposition(current.getJacobian(), solvedCols); + final double[][] weightedJacobian = internalData.weightedJacobian; + final int[] permutation = internalData.permutation; + final double[] diagR = internalData.diagR; + final double[] jacNorm = internalData.jacNorm; + + //residuals already have weights applied + double[] weightedResidual = currentResiduals; + for (int i = 0; i < nR; i++) { + qtf[i] = weightedResidual[i]; + } + + // compute Qt.res + qTy(qtf, internalData); + + // now we don't need Q anymore, + // so let jacobian contain the R matrix with its diagonal elements + for (int k = 0; k < solvedCols; ++k) { + int pk = permutation[k]; + weightedJacobian[k][pk] = diagR[pk]; + } + + if (firstIteration) { + // scale the point according to the norms of the columns + // of the initial jacobian + xNorm = 0; + for (int k = 0; k < nC; ++k) { + double dk = jacNorm[k]; + if (dk == 0) { + dk = 1.0; + } + double xk = dk * currentPoint[k]; + xNorm += xk * xk; + diag[k] = dk; + } + xNorm = FastMath.sqrt(xNorm); + + // initialize the step bound delta + delta = (xNorm == 0) ? initialStepBoundFactor : (initialStepBoundFactor * xNorm); + } + + // check orthogonality between function vector and jacobian columns + double maxCosine = 0; + if (currentCost != 0) { + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double s = jacNorm[pj]; + if (s != 0) { + double sum = 0; + for (int i = 0; i <= j; ++i) { + sum += weightedJacobian[i][pj] * qtf[i]; + } + maxCosine = FastMath.max(maxCosine, FastMath.abs(sum) / (s * currentCost)); + } + } + } + if (maxCosine <= orthoTolerance) { + // Convergence has been reached. + return new OptimumImpl( + current, + evaluationCounter.getCount(), + iterationCounter.getCount()); + } + + // rescale if necessary + for (int j = 0; j < nC; ++j) { + diag[j] = FastMath.max(diag[j], jacNorm[j]); + } + + // Inner loop. + for (double ratio = 0; ratio < 1.0e-4;) { + + // save the state + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + oldX[pj] = currentPoint[pj]; + } + final double previousCost = currentCost; + double[] tmpVec = weightedResidual; + weightedResidual = oldRes; + oldRes = tmpVec; + + // determine the Levenberg-Marquardt parameter + lmPar = determineLMParameter(qtf, delta, diag, + internalData, solvedCols, + work1, work2, work3, lmDir, lmPar); + + // compute the new point and the norm of the evolution direction + double lmNorm = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + lmDir[pj] = -lmDir[pj]; + currentPoint[pj] = oldX[pj] + lmDir[pj]; + double s = diag[pj] * lmDir[pj]; + lmNorm += s * s; + } + lmNorm = FastMath.sqrt(lmNorm); + // on the first iteration, adjust the initial step bound. + if (firstIteration) { + delta = FastMath.min(delta, lmNorm); + } + + // Evaluate the function at x + p and calculate its norm. + evaluationCounter.incrementCount(); + current = problem.evaluate(new ArrayRealVector(currentPoint)); + currentResiduals = current.getResiduals().toArray(); + currentCost = current.getCost(); + currentPoint = current.getPoint().toArray(); + + // compute the scaled actual reduction + double actRed = -1.0; + if (0.1 * currentCost < previousCost) { + double r = currentCost / previousCost; + actRed = 1.0 - r * r; + } + + // compute the scaled predicted reduction + // and the scaled directional derivative + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double dirJ = lmDir[pj]; + work1[j] = 0; + for (int i = 0; i <= j; ++i) { + work1[i] += weightedJacobian[i][pj] * dirJ; + } + } + double coeff1 = 0; + for (int j = 0; j < solvedCols; ++j) { + coeff1 += work1[j] * work1[j]; + } + double pc2 = previousCost * previousCost; + coeff1 /= pc2; + double coeff2 = lmPar * lmNorm * lmNorm / pc2; + double preRed = coeff1 + 2 * coeff2; + double dirDer = -(coeff1 + coeff2); + + // ratio of the actual to the predicted reduction + ratio = (preRed == 0) ? 0 : (actRed / preRed); + + // update the step bound + if (ratio <= 0.25) { + double tmp = + (actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5; + if ((0.1 * currentCost >= previousCost) || (tmp < 0.1)) { + tmp = 0.1; + } + delta = tmp * FastMath.min(delta, 10.0 * lmNorm); + lmPar /= tmp; + } else if ((lmPar == 0) || (ratio >= 0.75)) { + delta = 2 * lmNorm; + lmPar *= 0.5; + } + + // test for successful iteration. + if (ratio >= 1.0e-4) { + // successful iteration, update the norm + firstIteration = false; + xNorm = 0; + for (int k = 0; k < nC; ++k) { + double xK = diag[k] * currentPoint[k]; + xNorm += xK * xK; + } + xNorm = FastMath.sqrt(xNorm); + + // tests for convergence. + if (checker != null && checker.converged(iterationCounter.getCount(), previous, current)) { + return new OptimumImpl(current, evaluationCounter.getCount(), iterationCounter.getCount()); + } + } else { + // failed iteration, reset the previous values + currentCost = previousCost; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + currentPoint[pj] = oldX[pj]; + } + tmpVec = weightedResidual; + weightedResidual = oldRes; + oldRes = tmpVec; + // Reset "current" to previous values. + current = previous; + } + + // Default convergence criteria. + if ((FastMath.abs(actRed) <= costRelativeTolerance && + preRed <= costRelativeTolerance && + ratio <= 2.0) || + delta <= parRelativeTolerance * xNorm) { + return new OptimumImpl(current, evaluationCounter.getCount(), iterationCounter.getCount()); + } + + // tests for termination and stringent tolerances + if (FastMath.abs(actRed) <= TWO_EPS && + preRed <= TWO_EPS && + ratio <= 2.0) { + throw new ConvergenceException(LocalizedFormats.TOO_SMALL_COST_RELATIVE_TOLERANCE, + costRelativeTolerance); + } else if (delta <= TWO_EPS * xNorm) { + throw new ConvergenceException(LocalizedFormats.TOO_SMALL_PARAMETERS_RELATIVE_TOLERANCE, + parRelativeTolerance); + } else if (maxCosine <= TWO_EPS) { + throw new ConvergenceException(LocalizedFormats.TOO_SMALL_ORTHOGONALITY_TOLERANCE, + orthoTolerance); + } + } + } + } + + /** + * Holds internal data. + * This structure was created so that all optimizer fields can be "final". + * Code should be further refactored in order to not pass around arguments + * that will modified in-place (cf. "work" arrays). + */ + private static class InternalData { + /** Weighted Jacobian. */ + private final double[][] weightedJacobian; + /** Columns permutation array. */ + private final int[] permutation; + /** Rank of the Jacobian matrix. */ + private final int rank; + /** Diagonal elements of the R matrix in the QR decomposition. */ + private final double[] diagR; + /** Norms of the columns of the jacobian matrix. */ + private final double[] jacNorm; + /** Coefficients of the Householder transforms vectors. */ + private final double[] beta; + + /** + * @param weightedJacobian Weighted Jacobian. + * @param permutation Columns permutation array. + * @param rank Rank of the Jacobian matrix. + * @param diagR Diagonal elements of the R matrix in the QR decomposition. + * @param jacNorm Norms of the columns of the jacobian matrix. + * @param beta Coefficients of the Householder transforms vectors. + */ + InternalData(double[][] weightedJacobian, + int[] permutation, + int rank, + double[] diagR, + double[] jacNorm, + double[] beta) { + this.weightedJacobian = weightedJacobian; + this.permutation = permutation; + this.rank = rank; + this.diagR = diagR; + this.jacNorm = jacNorm; + this.beta = beta; + } + } + + /** + * Determines the Levenberg-Marquardt parameter. + * + * <p>This implementation is a translation in Java of the MINPACK + * <a href="http://www.netlib.org/minpack/lmpar.f">lmpar</a> + * routine.</p> + * <p>This method sets the lmPar and lmDir attributes.</p> + * <p>The authors of the original fortran function are:</p> + * <ul> + * <li>Argonne National Laboratory. MINPACK project. March 1980</li> + * <li>Burton S. Garbow</li> + * <li>Kenneth E. Hillstrom</li> + * <li>Jorge J. More</li> + * </ul> + * <p>Luc Maisonobe did the Java translation.</p> + * + * @param qy Array containing qTy. + * @param delta Upper bound on the euclidean norm of diagR * lmDir. + * @param diag Diagonal matrix. + * @param internalData Data (modified in-place in this method). + * @param solvedCols Number of solved point. + * @param work1 work array + * @param work2 work array + * @param work3 work array + * @param lmDir the "returned" LM direction will be stored in this array. + * @param lmPar the value of the LM parameter from the previous iteration. + * @return the new LM parameter + */ + private double determineLMParameter(double[] qy, double delta, double[] diag, + InternalData internalData, int solvedCols, + double[] work1, double[] work2, double[] work3, + double[] lmDir, double lmPar) { + final double[][] weightedJacobian = internalData.weightedJacobian; + final int[] permutation = internalData.permutation; + final int rank = internalData.rank; + final double[] diagR = internalData.diagR; + + final int nC = weightedJacobian[0].length; + + // compute and store in x the gauss-newton direction, if the + // jacobian is rank-deficient, obtain a least squares solution + for (int j = 0; j < rank; ++j) { + lmDir[permutation[j]] = qy[j]; + } + for (int j = rank; j < nC; ++j) { + lmDir[permutation[j]] = 0; + } + for (int k = rank - 1; k >= 0; --k) { + int pk = permutation[k]; + double ypk = lmDir[pk] / diagR[pk]; + for (int i = 0; i < k; ++i) { + lmDir[permutation[i]] -= ypk * weightedJacobian[i][pk]; + } + lmDir[pk] = ypk; + } + + // evaluate the function at the origin, and test + // for acceptance of the Gauss-Newton direction + double dxNorm = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double s = diag[pj] * lmDir[pj]; + work1[pj] = s; + dxNorm += s * s; + } + dxNorm = FastMath.sqrt(dxNorm); + double fp = dxNorm - delta; + if (fp <= 0.1 * delta) { + lmPar = 0; + return lmPar; + } + + // if the jacobian is not rank deficient, the Newton step provides + // a lower bound, parl, for the zero of the function, + // otherwise set this bound to zero + double sum2; + double parl = 0; + if (rank == solvedCols) { + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] *= diag[pj] / dxNorm; + } + sum2 = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double sum = 0; + for (int i = 0; i < j; ++i) { + sum += weightedJacobian[i][pj] * work1[permutation[i]]; + } + double s = (work1[pj] - sum) / diagR[pj]; + work1[pj] = s; + sum2 += s * s; + } + parl = fp / (delta * sum2); + } + + // calculate an upper bound, paru, for the zero of the function + sum2 = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double sum = 0; + for (int i = 0; i <= j; ++i) { + sum += weightedJacobian[i][pj] * qy[i]; + } + sum /= diag[pj]; + sum2 += sum * sum; + } + double gNorm = FastMath.sqrt(sum2); + double paru = gNorm / delta; + if (paru == 0) { + paru = Precision.SAFE_MIN / FastMath.min(delta, 0.1); + } + + // if the input par lies outside of the interval (parl,paru), + // set par to the closer endpoint + lmPar = FastMath.min(paru, FastMath.max(lmPar, parl)); + if (lmPar == 0) { + lmPar = gNorm / dxNorm; + } + + for (int countdown = 10; countdown >= 0; --countdown) { + + // evaluate the function at the current value of lmPar + if (lmPar == 0) { + lmPar = FastMath.max(Precision.SAFE_MIN, 0.001 * paru); + } + double sPar = FastMath.sqrt(lmPar); + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] = sPar * diag[pj]; + } + determineLMDirection(qy, work1, work2, internalData, solvedCols, work3, lmDir); + + dxNorm = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double s = diag[pj] * lmDir[pj]; + work3[pj] = s; + dxNorm += s * s; + } + dxNorm = FastMath.sqrt(dxNorm); + double previousFP = fp; + fp = dxNorm - delta; + + // if the function is small enough, accept the current value + // of lmPar, also test for the exceptional cases where parl is zero + if (FastMath.abs(fp) <= 0.1 * delta || + (parl == 0 && + fp <= previousFP && + previousFP < 0)) { + return lmPar; + } + + // compute the Newton correction + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] = work3[pj] * diag[pj] / dxNorm; + } + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] /= work2[j]; + double tmp = work1[pj]; + for (int i = j + 1; i < solvedCols; ++i) { + work1[permutation[i]] -= weightedJacobian[i][pj] * tmp; + } + } + sum2 = 0; + for (int j = 0; j < solvedCols; ++j) { + double s = work1[permutation[j]]; + sum2 += s * s; + } + double correction = fp / (delta * sum2); + + // depending on the sign of the function, update parl or paru. + if (fp > 0) { + parl = FastMath.max(parl, lmPar); + } else if (fp < 0) { + paru = FastMath.min(paru, lmPar); + } + + // compute an improved estimate for lmPar + lmPar = FastMath.max(parl, lmPar + correction); + } + + return lmPar; + } + + /** + * Solve a*x = b and d*x = 0 in the least squares sense. + * <p>This implementation is a translation in Java of the MINPACK + * <a href="http://www.netlib.org/minpack/qrsolv.f">qrsolv</a> + * routine.</p> + * <p>This method sets the lmDir and lmDiag attributes.</p> + * <p>The authors of the original fortran function are:</p> + * <ul> + * <li>Argonne National Laboratory. MINPACK project. March 1980</li> + * <li>Burton S. Garbow</li> + * <li>Kenneth E. Hillstrom</li> + * <li>Jorge J. More</li> + * </ul> + * <p>Luc Maisonobe did the Java translation.</p> + * + * @param qy array containing qTy + * @param diag diagonal matrix + * @param lmDiag diagonal elements associated with lmDir + * @param internalData Data (modified in-place in this method). + * @param solvedCols Number of sloved point. + * @param work work array + * @param lmDir the "returned" LM direction is stored in this array + */ + private void determineLMDirection(double[] qy, double[] diag, + double[] lmDiag, + InternalData internalData, + int solvedCols, + double[] work, + double[] lmDir) { + final int[] permutation = internalData.permutation; + final double[][] weightedJacobian = internalData.weightedJacobian; + final double[] diagR = internalData.diagR; + + // copy R and Qty to preserve input and initialize s + // in particular, save the diagonal elements of R in lmDir + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + for (int i = j + 1; i < solvedCols; ++i) { + weightedJacobian[i][pj] = weightedJacobian[j][permutation[i]]; + } + lmDir[j] = diagR[pj]; + work[j] = qy[j]; + } + + // eliminate the diagonal matrix d using a Givens rotation + for (int j = 0; j < solvedCols; ++j) { + + // prepare the row of d to be eliminated, locating the + // diagonal element using p from the Q.R. factorization + int pj = permutation[j]; + double dpj = diag[pj]; + if (dpj != 0) { + Arrays.fill(lmDiag, j + 1, lmDiag.length, 0); + } + lmDiag[j] = dpj; + + // the transformations to eliminate the row of d + // modify only a single element of Qty + // beyond the first n, which is initially zero. + double qtbpj = 0; + for (int k = j; k < solvedCols; ++k) { + int pk = permutation[k]; + + // determine a Givens rotation which eliminates the + // appropriate element in the current row of d + if (lmDiag[k] != 0) { + + final double sin; + final double cos; + double rkk = weightedJacobian[k][pk]; + if (FastMath.abs(rkk) < FastMath.abs(lmDiag[k])) { + final double cotan = rkk / lmDiag[k]; + sin = 1.0 / FastMath.sqrt(1.0 + cotan * cotan); + cos = sin * cotan; + } else { + final double tan = lmDiag[k] / rkk; + cos = 1.0 / FastMath.sqrt(1.0 + tan * tan); + sin = cos * tan; + } + + // compute the modified diagonal element of R and + // the modified element of (Qty,0) + weightedJacobian[k][pk] = cos * rkk + sin * lmDiag[k]; + final double temp = cos * work[k] + sin * qtbpj; + qtbpj = -sin * work[k] + cos * qtbpj; + work[k] = temp; + + // accumulate the tranformation in the row of s + for (int i = k + 1; i < solvedCols; ++i) { + double rik = weightedJacobian[i][pk]; + final double temp2 = cos * rik + sin * lmDiag[i]; + lmDiag[i] = -sin * rik + cos * lmDiag[i]; + weightedJacobian[i][pk] = temp2; + } + } + } + + // store the diagonal element of s and restore + // the corresponding diagonal element of R + lmDiag[j] = weightedJacobian[j][permutation[j]]; + weightedJacobian[j][permutation[j]] = lmDir[j]; + } + + // solve the triangular system for z, if the system is + // singular, then obtain a least squares solution + int nSing = solvedCols; + for (int j = 0; j < solvedCols; ++j) { + if ((lmDiag[j] == 0) && (nSing == solvedCols)) { + nSing = j; + } + if (nSing < solvedCols) { + work[j] = 0; + } + } + if (nSing > 0) { + for (int j = nSing - 1; j >= 0; --j) { + int pj = permutation[j]; + double sum = 0; + for (int i = j + 1; i < nSing; ++i) { + sum += weightedJacobian[i][pj] * work[i]; + } + work[j] = (work[j] - sum) / lmDiag[j]; + } + } + + // permute the components of z back to components of lmDir + for (int j = 0; j < lmDir.length; ++j) { + lmDir[permutation[j]] = work[j]; + } + } + + /** + * Decompose a matrix A as A.P = Q.R using Householder transforms. + * <p>As suggested in the P. Lascaux and R. Theodor book + * <i>Analyse numérique matricielle appliquée à + * l'art de l'ingénieur</i> (Masson, 1986), instead of representing + * the Householder transforms with u<sub>k</sub> unit vectors such that: + * <pre> + * H<sub>k</sub> = I - 2u<sub>k</sub>.u<sub>k</sub><sup>t</sup> + * </pre> + * we use <sub>k</sub> non-unit vectors such that: + * <pre> + * H<sub>k</sub> = I - beta<sub>k</sub>v<sub>k</sub>.v<sub>k</sub><sup>t</sup> + * </pre> + * where v<sub>k</sub> = a<sub>k</sub> - alpha<sub>k</sub> e<sub>k</sub>. + * The beta<sub>k</sub> coefficients are provided upon exit as recomputing + * them from the v<sub>k</sub> vectors would be costly.</p> + * <p>This decomposition handles rank deficient cases since the tranformations + * are performed in non-increasing columns norms order thanks to columns + * pivoting. The diagonal elements of the R matrix are therefore also in + * non-increasing absolute values order.</p> + * + * @param jacobian Weighted Jacobian matrix at the current point. + * @param solvedCols Number of solved point. + * @return data used in other methods of this class. + * @throws ConvergenceException if the decomposition cannot be performed. + */ + private InternalData qrDecomposition(RealMatrix jacobian, + int solvedCols) throws ConvergenceException { + // Code in this class assumes that the weighted Jacobian is -(W^(1/2) J), + // hence the multiplication by -1. + final double[][] weightedJacobian = jacobian.scalarMultiply(-1).getData(); + + final int nR = weightedJacobian.length; + final int nC = weightedJacobian[0].length; + + final int[] permutation = new int[nC]; + final double[] diagR = new double[nC]; + final double[] jacNorm = new double[nC]; + final double[] beta = new double[nC]; + + // initializations + for (int k = 0; k < nC; ++k) { + permutation[k] = k; + double norm2 = 0; + for (int i = 0; i < nR; ++i) { + double akk = weightedJacobian[i][k]; + norm2 += akk * akk; + } + jacNorm[k] = FastMath.sqrt(norm2); + } + + // transform the matrix column after column + for (int k = 0; k < nC; ++k) { + + // select the column with the greatest norm on active components + int nextColumn = -1; + double ak2 = Double.NEGATIVE_INFINITY; + for (int i = k; i < nC; ++i) { + double norm2 = 0; + for (int j = k; j < nR; ++j) { + double aki = weightedJacobian[j][permutation[i]]; + norm2 += aki * aki; + } + if (Double.isInfinite(norm2) || Double.isNaN(norm2)) { + throw new ConvergenceException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN, + nR, nC); + } + if (norm2 > ak2) { + nextColumn = i; + ak2 = norm2; + } + } + if (ak2 <= qrRankingThreshold) { + return new InternalData(weightedJacobian, permutation, k, diagR, jacNorm, beta); + } + int pk = permutation[nextColumn]; + permutation[nextColumn] = permutation[k]; + permutation[k] = pk; + + // choose alpha such that Hk.u = alpha ek + double akk = weightedJacobian[k][pk]; + double alpha = (akk > 0) ? -FastMath.sqrt(ak2) : FastMath.sqrt(ak2); + double betak = 1.0 / (ak2 - akk * alpha); + beta[pk] = betak; + + // transform the current column + diagR[pk] = alpha; + weightedJacobian[k][pk] -= alpha; + + // transform the remaining columns + for (int dk = nC - 1 - k; dk > 0; --dk) { + double gamma = 0; + for (int j = k; j < nR; ++j) { + gamma += weightedJacobian[j][pk] * weightedJacobian[j][permutation[k + dk]]; + } + gamma *= betak; + for (int j = k; j < nR; ++j) { + weightedJacobian[j][permutation[k + dk]] -= gamma * weightedJacobian[j][pk]; + } + } + } + + return new InternalData(weightedJacobian, permutation, solvedCols, diagR, jacNorm, beta); + } + + /** + * Compute the product Qt.y for some Q.R. decomposition. + * + * @param y vector to multiply (will be overwritten with the result) + * @param internalData Data. + */ + private void qTy(double[] y, + InternalData internalData) { + final double[][] weightedJacobian = internalData.weightedJacobian; + final int[] permutation = internalData.permutation; + final double[] beta = internalData.beta; + + final int nR = weightedJacobian.length; + final int nC = weightedJacobian[0].length; + + for (int k = 0; k < nC; ++k) { + int pk = permutation[k]; + double gamma = 0; + for (int i = k; i < nR; ++i) { + gamma += weightedJacobian[i][pk] * y[i]; + } + gamma *= beta[pk]; + for (int i = k; i < nR; ++i) { + y[i] -= gamma * weightedJacobian[i][pk]; + } + } + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/MultivariateJacobianFunction.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/MultivariateJacobianFunction.java new file mode 100644 index 0000000..e673855 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/MultivariateJacobianFunction.java @@ -0,0 +1,39 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.apache.commons.math3.util.Pair; + +/** + * A interface for functions that compute a vector of values and can compute their + * derivatives (Jacobian). + * + * @since 3.3 + */ +public interface MultivariateJacobianFunction { + + /** + * Compute the function value and its Jacobian. + * + * @param point the abscissae + * @return the values and their Jacobian of this vector valued function. + */ + Pair<RealVector, RealMatrix> value(RealVector point); + +} 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 new file mode 100644 index 0000000..698f86c --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/OptimumImpl.java @@ -0,0 +1,97 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum; +import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; + +/** + * A pedantic implementation of {@link Optimum}. + * + * @since 3.3 + */ +class OptimumImpl implements Optimum { + + /** abscissa and ordinate */ + private final Evaluation value; + /** number of evaluations to compute this optimum */ + private final int evaluations; + /** number of iterations to compute this optimum */ + private final int iterations; + + /** + * Construct an optimum from an evaluation and the values of the counters. + * + * @param value the function value + * @param evaluations number of times the function was evaluated + * @param iterations number of iterations of the algorithm + */ + OptimumImpl(final Evaluation value, final int evaluations, final int iterations) { + this.value = value; + this.evaluations = evaluations; + this.iterations = iterations; + } + + /* auto-generated implementations */ + + /** {@inheritDoc} */ + public int getEvaluations() { + return evaluations; + } + + /** {@inheritDoc} */ + public int getIterations() { + return iterations; + } + + /** {@inheritDoc} */ + public RealMatrix getCovariances(double threshold) { + return value.getCovariances(threshold); + } + + /** {@inheritDoc} */ + public RealVector getSigma(double covarianceSingularityThreshold) { + return value.getSigma(covarianceSingularityThreshold); + } + + /** {@inheritDoc} */ + public double getRMS() { + return value.getRMS(); + } + + /** {@inheritDoc} */ + public RealMatrix getJacobian() { + return value.getJacobian(); + } + + /** {@inheritDoc} */ + public double getCost() { + return value.getCost(); + } + + /** {@inheritDoc} */ + public RealVector getResiduals() { + return value.getResiduals(); + } + + /** {@inheritDoc} */ + public RealVector getPoint() { + return value.getPoint(); + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/ParameterValidator.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/ParameterValidator.java new file mode 100644 index 0000000..d5b8529 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/ParameterValidator.java @@ -0,0 +1,34 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.linear.RealVector; + +/** + * Interface for validating a set of model parameters. + * + * @since 3.4 + */ +public interface ParameterValidator { + /** + * Validates the set of parameters. + * + * @param params Input parameters. + * @return the validated values. + */ + RealVector validate(RealVector params); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/ValueAndJacobianFunction.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/ValueAndJacobianFunction.java new file mode 100644 index 0000000..180e328 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/ValueAndJacobianFunction.java @@ -0,0 +1,44 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting.leastsquares; + +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; + +/** + * A interface for functions that compute a vector of values and can compute their + * derivatives (Jacobian). + * + * @since 3.4 + */ +public interface ValueAndJacobianFunction extends MultivariateJacobianFunction { + /** + * Compute the value. + * + * @param params Point. + * @return the value at the given point. + */ + RealVector computeValue(final double[] params); + + /** + * Compute the Jacobian. + * + * @param params Point. + * @return the Jacobian at the given point. + */ + RealMatrix computeJacobian(final double[] params); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java new file mode 100644 index 0000000..98623b5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java @@ -0,0 +1,39 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +/** + * This package provides algorithms that minimize the residuals + * between observations and model values. + * The {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer + * least-squares optimizers} minimize the distance (called + * <em>cost</em> or <em>χ<sup>2</sup></em>) between model and + * observations. + * + * <br/> + * Algorithms in this category need access to a <em>problem</em> + * (represented by a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem + * LeastSquaresProblem}). + * Such a model predicts a set of values which the algorithm tries to match + * with a set of given set of observed values. + * <br/> + * The problem can be created progressively using a {@link + * org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder builder} or it can + * be created at once using a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory + * factory}. + * @since 3.3 + */ +package org.apache.commons.math3.fitting.leastsquares; |