http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresBuilder.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 7b14b37..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresBuilder.java +++ /dev/null @@ -1,226 +0,0 @@ -/* - * 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; - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresFactory.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index d483c26..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresFactory.java +++ /dev/null @@ -1,528 +0,0 @@ -/* - * 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) { - @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) { - - 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>() { - 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())); - } - } - } -} - http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresOptimizer.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 50d5b8a..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresOptimizer.java +++ /dev/null @@ -1,62 +0,0 @@ -/* - * 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(); - - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblem.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 097ff81..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresProblem.java +++ /dev/null @@ -1,156 +0,0 @@ -/* - * 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(); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 358d240..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java +++ /dev/null @@ -1,1042 +0,0 @@ -/* - * 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]; - } - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/MultivariateJacobianFunction.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index e673855..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/MultivariateJacobianFunction.java +++ /dev/null @@ -1,39 +0,0 @@ -/* - * 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); - -}
