You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2013/07/30 17:04:23 UTC
svn commit: r1508481 [1/5] - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math3/fitting/leastsquares/
main/java/org/apache/commons/math3/optim/
test/java/org/apache/commons/math3/fitting/leastsquares/
test/resources/org/apache/common...
Author: erans
Date: Tue Jul 30 15:04:22 2013
New Revision: 1508481
URL: http://svn.apache.org/r1508481
Log:
MATH-1008
Created package "o.a.c.m.fitting.leastsquares" to contain a modified version
of the contents of (to-be-deprecated) "o.a.c.m.optim.nonlinear.vector", the
main purpose being to provide a new "fluent" API (cf. "withXxx" methods).
Along the way, class "LevenbergMarquardtOptimizer" has been further cleaned
up (e.g. removing protected fields and deprecated methods and using local
variables instead of instance fields).
An additional constructor was necessary in "BaseOptimizer".
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithConvergenceChecker.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxEvaluations.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxIterations.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithModelAndJacobian.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithStartPoint.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithTarget.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithWeight.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerAbstractTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTestValidation.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleProblem.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleVectorial.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/MinpackTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomCirclePointGenerator.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomStraightLinePointGenerator.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDataset.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDatasetFactory.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/StraightLineProblem.java (with props)
commons/proper/math/trunk/src/test/resources/org/apache/commons/math3/fitting/
commons/proper/math/trunk/src/test/resources/org/apache/commons/math3/fitting/leastsquares/
commons/proper/math/trunk/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Hahn1.dat
commons/proper/math/trunk/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Kirby2.dat
commons/proper/math/trunk/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Lanczos1.dat
commons/proper/math/trunk/src/test/resources/org/apache/commons/math3/fitting/leastsquares/MGH17.dat
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optim/BaseOptimizer.java
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java?rev=1508481&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java Tue Jul 30 15:04:22 2013
@@ -0,0 +1,362 @@
+/*
+ * 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.DimensionMismatchException;
+import org.apache.commons.math3.exception.MathUnsupportedOperationException;
+import org.apache.commons.math3.analysis.MultivariateVectorFunction;
+import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
+import org.apache.commons.math3.linear.ArrayRealVector;
+import org.apache.commons.math3.linear.RealMatrix;
+import org.apache.commons.math3.linear.DiagonalMatrix;
+import org.apache.commons.math3.linear.DecompositionSolver;
+import org.apache.commons.math3.linear.MatrixUtils;
+import org.apache.commons.math3.linear.QRDecomposition;
+import org.apache.commons.math3.linear.EigenDecomposition;
+import org.apache.commons.math3.optim.ConvergenceChecker;
+import org.apache.commons.math3.optim.BaseOptimizer;
+import org.apache.commons.math3.optim.PointVectorValuePair;
+import org.apache.commons.math3.optim.OptimizationData;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * Base class for implementing least-squares optimizers.
+ * It provides methods for error estimation.
+ *
+ * @version $Id$
+ * @since 3.3
+ */
+public abstract class AbstractLeastSquaresOptimizer
+ extends BaseOptimizer<PointVectorValuePair> {
+ /** Target values for the model function at optimum. */
+ private final double[] target;
+ /** Weight matrix. */
+ private final RealMatrix weight;
+ /** Model function. */
+ private final MultivariateVectorFunction model;
+ /** Jacobian of the model function. */
+ private final MultivariateMatrixFunction jacobian;
+ /** Square-root of the weight matrix. */
+ private final RealMatrix weightSqrt;
+ /** Initial guess. */
+ private final double[] start;
+
+ /**
+ * @param target Observations.
+ * @param weight Weight of the observations.
+ * For performance, no defensive copy is performed.
+ * @param weightSqrt Square-root of the {@code weight} matrix.
+ * If {@code null}, it will be computed; otherwise it is the caller's
+ * responsibility that {@code weight} and {@code weightSqrt} are
+ * consistent.
+ * No defensive copy is performed.
+ * @param model ModelFunction.
+ * @param jacobian Jacobian of the model function.
+ * @param checker Convergence checker.
+ * @param start Initial guess.
+ * @param maxEval Maximum number of evaluations of the model
+ * function.
+ * @param maxIter Maximum number of iterations.
+ */
+ protected AbstractLeastSquaresOptimizer(double[] target,
+ RealMatrix weight,
+ RealMatrix weightSqrt,
+ MultivariateVectorFunction model,
+ MultivariateMatrixFunction jacobian,
+ ConvergenceChecker<PointVectorValuePair> checker,
+ double[] start,
+ int maxEval,
+ int maxIter) {
+ super(checker, maxEval, maxIter);
+
+ this.target = target;
+ this.weight = weight;
+ this.model = model;
+ this.jacobian = jacobian;
+ this.start = start;
+
+ this.weightSqrt = weightSqrt == null ?
+ (weight == null ?
+ null : squareRoot(weight)) : weightSqrt;
+ }
+
+ /**
+ * Gets the target values.
+ *
+ * @return the target values.
+ */
+ public double[] getTarget() {
+ return target == null ? null : target.clone();
+ }
+
+ /**
+ * Gets the initial guess.
+ *
+ * @return the initial guess values.
+ */
+ public double[] getStart() {
+ return start == null ? null : start.clone();
+ }
+
+ /**
+ * Gets the square-root of the weight matrix.
+ *
+ * @return the square-root of the weight matrix.
+ */
+ public RealMatrix getWeightSquareRoot() {
+ return weightSqrt == null ? null : weightSqrt.copy();
+ }
+
+ /**
+ * Gets the model function.
+ *
+ * @return the model function.
+ */
+ public MultivariateVectorFunction getModel() {
+ return model;
+ }
+
+ /**
+ * Gets the model function's Jacobian.
+ *
+ * @return the Jacobian.
+ */
+ public MultivariateMatrixFunction getJacobian() {
+ return jacobian;
+ }
+
+ /**
+ * 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 params Model parameters.
+ * @param threshold Singularity threshold.
+ * @return the covariance matrix.
+ * @throws org.apache.commons.math3.linear.SingularMatrixException
+ * if the covariance matrix cannot be computed (singular problem).
+ */
+ public double[][] computeCovariances(double[] params,
+ double threshold) {
+ // Set up the Jacobian.
+ final RealMatrix j = computeWeightedJacobian(params);
+
+ // Compute transpose(J)J.
+ final RealMatrix jTj = j.transpose().multiply(j);
+
+ // Compute the covariances matrix.
+ final DecompositionSolver solver
+ = new QRDecomposition(jTj, threshold).getSolver();
+ return solver.getInverse().getData();
+ }
+
+ /**
+ * Computes 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 params Model parameters.
+ * @param covarianceSingularityThreshold Singularity threshold (see
+ * {@link #computeCovariances(double[],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.
+ */
+ public double[] computeSigma(double[] params,
+ double covarianceSingularityThreshold) {
+ final int nC = params.length;
+ final double[] sig = new double[nC];
+ final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
+ for (int i = 0; i < nC; ++i) {
+ sig[i] = FastMath.sqrt(cov[i][i]);
+ }
+ return sig;
+ }
+
+ /**
+ * Gets the weight matrix of the observations.
+ *
+ * @return the weight matrix.
+ */
+ public RealMatrix getWeight() {
+ return weight.copy();
+ }
+
+ /**
+ * Computes the normalized cost.
+ * It is the square-root of the sum of squared of the residuals, divided
+ * by the number of measurements.
+ *
+ * @param params Model function parameters.
+ * @return the cost.
+ */
+ public double computeRMS(double[] params) {
+ final double cost = computeCost(computeResiduals(getModel().value(params)));
+ return FastMath.sqrt(cost * cost / target.length);
+ }
+
+ /**
+ * Calling this method will raise an exception.
+ *
+ * @param optData Obsolete.
+ * @return nothing.
+ * @throws MathUnsupportedOperationException if called.
+ * @deprecated Do not use this method.
+ */
+ @Deprecated
+ @Override
+ public PointVectorValuePair optimize(OptimizationData... optData)
+ throws MathUnsupportedOperationException {
+ throw new MathUnsupportedOperationException();
+ }
+
+ /**
+ * Gets a reference to the corresponding field.
+ * Altering it could jeopardize the consistency of this class.
+ *
+ * @return the reference.
+ */
+ protected double[] getTargetInternal() {
+ return target;
+ }
+
+ /**
+ * Gets a reference to the corresponding field.
+ * Altering it could jeopardize the consistency of this class.
+ *
+ * @return the reference.
+ */
+ protected RealMatrix getWeightInternal() {
+ return weight;
+ }
+
+ /**
+ * Gets a reference to the corresponding field.
+ * Altering it could jeopardize the consistency of this class.
+ *
+ * @return the reference.
+ */
+ protected RealMatrix getWeightSquareRootInternal() {
+ return weightSqrt;
+ }
+
+ /**
+ * Computes the objective function value.
+ * This method <em>must</em> be called by subclasses to enforce the
+ * evaluation counter limit.
+ *
+ * @param params Point at which the objective function must be evaluated.
+ * @return the objective function value at the specified point.
+ * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
+ * if the maximal number of evaluations (of the model vector function) is
+ * exceeded.
+ */
+ protected double[] computeObjectiveValue(double[] params) {
+ super.incrementEvaluationCount();
+ return model.value(params);
+ }
+
+ /**
+ * Computes the weighted Jacobian matrix.
+ *
+ * @param params Model parameters at which to compute the Jacobian.
+ * @return the weighted Jacobian: W<sup>1/2</sup> J.
+ * @throws DimensionMismatchException if the Jacobian dimension does not
+ * match problem dimension.
+ */
+ protected RealMatrix computeWeightedJacobian(double[] params) {
+ return weightSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params)));
+ }
+
+ /**
+ * Computes the Jacobian matrix.
+ *
+ * @param params Point at which the Jacobian must be evaluated.
+ * @return the Jacobian at the specified point.
+ */
+ protected double[][] computeJacobian(final double[] params) {
+ return jacobian.value(params);
+ }
+
+ /**
+ * Computes the cost.
+ *
+ * @param residuals Residuals.
+ * @return the cost.
+ * @see #computeResiduals(double[])
+ */
+ protected double computeCost(double[] residuals) {
+ final ArrayRealVector r = new ArrayRealVector(residuals);
+ return FastMath.sqrt(r.dotProduct(weight.operate(r)));
+ }
+
+ /**
+ * Computes the 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.
+ *
+ * @param objectiveValue Value of the the objective function. This is
+ * the value returned from a call to
+ * {@link #computeObjectiveValue(double[]) computeObjectiveValue}
+ * (whose array argument contains the model parameters).
+ * @return the residuals.
+ * @throws DimensionMismatchException if {@code params} has a wrong
+ * length.
+ */
+ protected double[] computeResiduals(double[] objectiveValue) {
+ if (objectiveValue.length != target.length) {
+ throw new DimensionMismatchException(target.length,
+ objectiveValue.length);
+ }
+
+ final double[] residuals = new double[target.length];
+ for (int i = 0; i < target.length; i++) {
+ residuals[i] = target[i] - objectiveValue[i];
+ }
+
+ return residuals;
+ }
+
+ /**
+ * Computes the square-root of the weight matrix.
+ *
+ * @param m Symmetric, positive-definite (weight) matrix.
+ * @return the square-root of the weight matrix.
+ */
+ private RealMatrix squareRoot(RealMatrix m) {
+ if (m instanceof DiagonalMatrix) {
+ final int dim = m.getRowDimension();
+ final RealMatrix sqrtM = new DiagonalMatrix(dim);
+ for (int i = 0; i < dim; i++) {
+ sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
+ }
+ return sqrtM;
+ } else {
+ final EigenDecomposition dec = new EigenDecomposition(m);
+ return dec.getSquareRoot();
+ }
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java
------------------------------------------------------------------------------
svn:keywords = Id Revision
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java?rev=1508481&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java Tue Jul 30 15:04:22 2013
@@ -0,0 +1,325 @@
+/*
+ * 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.MultivariateVectorFunction;
+import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.ConvergenceException;
+import org.apache.commons.math3.exception.NullArgumentException;
+import org.apache.commons.math3.exception.MathInternalError;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.linear.ArrayRealVector;
+import org.apache.commons.math3.linear.BlockRealMatrix;
+import org.apache.commons.math3.linear.DecompositionSolver;
+import org.apache.commons.math3.linear.LUDecomposition;
+import org.apache.commons.math3.linear.QRDecomposition;
+import org.apache.commons.math3.linear.RealMatrix;
+import org.apache.commons.math3.linear.SingularMatrixException;
+import org.apache.commons.math3.optim.ConvergenceChecker;
+import org.apache.commons.math3.optim.PointVectorValuePair;
+
+/**
+ * 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
+ * QR decomposition can be used to solve the normal equations. LU decomposition
+ * is faster but QR decomposition is more robust for difficult problems.
+ * </p>
+ *
+ * @version $Id$
+ * @since 3.3
+ *
+ */
+public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer
+ implements WithTarget<GaussNewtonOptimizer>,
+ WithWeight<GaussNewtonOptimizer>,
+ WithModelAndJacobian<GaussNewtonOptimizer>,
+ WithConvergenceChecker<GaussNewtonOptimizer>,
+ WithStartPoint<GaussNewtonOptimizer>,
+ WithMaxIterations<GaussNewtonOptimizer>,
+ WithMaxEvaluations<GaussNewtonOptimizer> {
+ /** Indicator for using LU decomposition. */
+ private final boolean useLU;
+
+ /**
+ * Constructor called by the various {@code withXxx} methods.
+ *
+ * @param target Observations.
+ * @param weight Weight of the observations.
+ * For performance, no defensive copy is performed.
+ * @param weightSqrt Square-root of the {@code weight} matrix.
+ * If {@code null}, it will be computed; otherwise it is the caller's
+ * responsibility that {@code weight} and {@code weightSqrt} are
+ * consistent.
+ * No defensive copy is performed.
+ * @param model ModelFunction.
+ * @param jacobian Jacobian of the model function.
+ * @param checker Convergence checker.
+ * @param start Initial guess.
+ * @param maxEval Maximum number of evaluations of the model
+ * function.
+ * @param maxIter Maximum number of iterations.
+ * @param useLU Whether to use LU decomposition.
+ */
+ private GaussNewtonOptimizer(double[] target,
+ RealMatrix weight,
+ RealMatrix weightSqrt,
+ MultivariateVectorFunction model,
+ MultivariateMatrixFunction jacobian,
+ ConvergenceChecker<PointVectorValuePair> checker,
+ double[] start,
+ int maxEval,
+ int maxIter,
+ boolean useLU) {
+ super(target, weight, weightSqrt, model, jacobian, checker, start, maxEval, maxIter);
+
+ this.useLU = useLU;
+ }
+
+ /**
+ * Creates a bare-bones instance.
+ * Several calls to {@code withXxx} methods are necessary to obtain
+ * an object with all necessary fields set to sensible values.
+ * <br/>
+ * The default for the algorithm is to solve the normal equations
+ * using LU decomposition.
+ *
+ * @return an instance of this class.
+ */
+ public static GaussNewtonOptimizer create() {
+ return new GaussNewtonOptimizer(null, null, null, null, null, null, null,
+ 0, 0, true);
+ }
+
+ /** {@inheritDoc} */
+ public GaussNewtonOptimizer withTarget(double[] target) {
+ return new GaussNewtonOptimizer(target,
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ useLU);
+ }
+
+ /** {@inheritDoc} */
+ public GaussNewtonOptimizer withWeight(RealMatrix weight) {
+ return new GaussNewtonOptimizer(getTargetInternal(),
+ weight,
+ null,
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ useLU);
+ }
+
+ /** {@inheritDoc} */
+ public GaussNewtonOptimizer withModelAndJacobian(MultivariateVectorFunction model,
+ MultivariateMatrixFunction jacobian) {
+ return new GaussNewtonOptimizer(getTargetInternal(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ model,
+ jacobian,
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ useLU);
+ }
+
+ /** {@inheritDoc} */
+ public GaussNewtonOptimizer withConvergenceChecker(ConvergenceChecker<PointVectorValuePair> checker) {
+ return new GaussNewtonOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ checker,
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ useLU);
+ }
+
+ /** {@inheritDoc} */
+ public GaussNewtonOptimizer withStartPoint(double[] start) {
+ return new GaussNewtonOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ start,
+ getMaxEvaluations(),
+ getMaxIterations(),
+ useLU);
+ }
+
+ /** {@inheritDoc} */
+ public GaussNewtonOptimizer withMaxIterations(int maxIter) {
+ return new GaussNewtonOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ maxIter,
+ useLU);
+ }
+
+ /** {@inheritDoc} */
+ public GaussNewtonOptimizer withMaxEvaluations(int maxEval) {
+ return new GaussNewtonOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ maxEval,
+ getMaxIterations(),
+ useLU);
+ }
+
+ /**
+ * Creates a new instance.
+ *
+ * @param withLU Whether to use LU decomposition.
+ * @return a new instance with all fields identical to this instance except
+ * for the givens arguments.
+ */
+ public GaussNewtonOptimizer withLU(boolean withLU) {
+ return new GaussNewtonOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ withLU);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public PointVectorValuePair doOptimize() {
+ final ConvergenceChecker<PointVectorValuePair> checker
+ = getConvergenceChecker();
+
+ // Computation will be useless without a checker (see "for-loop").
+ if (checker == null) {
+ throw new NullArgumentException();
+ }
+
+ final double[] targetValues = getTarget();
+ final int nR = targetValues.length; // Number of observed data.
+
+ final RealMatrix weightMatrix = getWeight();
+ if (weightMatrix.getRowDimension() != nR) {
+ throw new DimensionMismatchException(weightMatrix.getRowDimension(), nR);
+ }
+ if (weightMatrix.getColumnDimension() != nR) {
+ throw new DimensionMismatchException(weightMatrix.getColumnDimension(), nR);
+ }
+
+ // Diagonal of the weight matrix.
+ final double[] residualsWeights = new double[nR];
+ for (int i = 0; i < nR; i++) {
+ residualsWeights[i] = weightMatrix.getEntry(i, i);
+ }
+
+ final double[] currentPoint = getStart();
+ final int nC = currentPoint.length;
+
+ // iterate until convergence is reached
+ PointVectorValuePair current = null;
+ for (boolean converged = false; !converged;) {
+ incrementIterationCount();
+
+ // evaluate the objective function and its jacobian
+ PointVectorValuePair previous = current;
+ // Value of the objective function at "currentPoint".
+ final double[] currentObjective = computeObjectiveValue(currentPoint);
+ final double[] currentResiduals = computeResiduals(currentObjective);
+ final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
+ current = new PointVectorValuePair(currentPoint, currentObjective);
+
+ // build the linear problem
+ final double[] b = new double[nC];
+ final double[][] a = new double[nC][nC];
+ for (int i = 0; i < nR; ++i) {
+
+ final double[] grad = weightedJacobian.getRow(i);
+ final double weight = residualsWeights[i];
+ final double residual = currentResiduals[i];
+
+ // compute the normal equation
+ final double wr = weight * residual;
+ for (int j = 0; j < nC; ++j) {
+ b[j] += wr * grad[j];
+ }
+
+ // build the contribution matrix for measurement i
+ for (int k = 0; k < nC; ++k) {
+ double[] ak = a[k];
+ double wgk = weight * grad[k];
+ for (int l = 0; l < nC; ++l) {
+ ak[l] += wgk * grad[l];
+ }
+ }
+ }
+
+ // Check convergence.
+ if (previous != null) {
+ converged = checker.converged(getIterations(), previous, current);
+ if (converged) {
+ return current;
+ }
+ }
+
+ try {
+ // solve the linearized least squares problem
+ RealMatrix mA = new BlockRealMatrix(a);
+ DecompositionSolver solver = useLU ?
+ new LUDecomposition(mA).getSolver() :
+ new QRDecomposition(mA).getSolver();
+ final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
+ // update the estimated parameters
+ for (int i = 0; i < nC; ++i) {
+ currentPoint[i] += dX[i];
+ }
+ } catch (SingularMatrixException e) {
+ throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
+ }
+ }
+ // Must never happen.
+ throw new MathInternalError();
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java
------------------------------------------------------------------------------
svn:keywords = Id Revision
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java?rev=1508481&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java Tue Jul 30 15:04:22 2013
@@ -0,0 +1,1172 @@
+/*
+ * 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.analysis.MultivariateVectorFunction;
+import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
+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.optim.PointVectorValuePair;
+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>
+ *
+ * @version $Id$
+ * @since 2.0
+ */
+public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer
+ implements WithTarget<LevenbergMarquardtOptimizer>,
+ WithWeight<LevenbergMarquardtOptimizer>,
+ WithModelAndJacobian<LevenbergMarquardtOptimizer>,
+ WithConvergenceChecker<LevenbergMarquardtOptimizer>,
+ WithStartPoint<LevenbergMarquardtOptimizer>,
+ WithMaxIterations<LevenbergMarquardtOptimizer>,
+ WithMaxEvaluations<LevenbergMarquardtOptimizer> {
+ /** Twice the "epsilon machine". */
+ private static final double TWO_EPS = 2 * Precision.EPSILON;
+ /** 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;
+ /** Levenberg-Marquardt parameter. */
+ private double lmPar;
+ /** Parameters evolution direction associated with lmPar. */
+ private double[] lmDir;
+
+ /**
+ * Constructor called by the various {@code withXxx} methods.
+ *
+ * @param target Observations.
+ * @param weight Weight of the observations.
+ * For performance, no defensive copy is performed.
+ * @param weightSqrt Square-root of the {@code weight} matrix.
+ * If {@code null}, it will be computed; otherwise it is the caller's
+ * responsibility that {@code weight} and {@code weightSqrt} are
+ * consistent.
+ * No defensive copy is performed.
+ * @param model ModelFunction.
+ * @param jacobian Jacobian of the model function.
+ * @param checker Convergence checker.
+ * @param start Initial guess.
+ * @param maxEval Maximum number of evaluations of the model
+ * function.
+ * @param maxIter Maximum number of iterations.
+ * @param initialStepBoundFactor 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 initialStepBoundFactor}
+ * itself. In most cases factor should lie in the interval
+ * {@code (0.1, 100.0)}. {@code 100} is a generally recommended value.
+ * @param costRelativeTolerance Desired relative error in the sum of
+ * squares.
+ * @param parRelativeTolerance Desired relative error in the approximate
+ * solution parameters.
+ * @param orthoTolerance Desired max cosine on the orthogonality between
+ * the function vector and the columns of the Jacobian.
+ * @param threshold 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.
+ */
+ private LevenbergMarquardtOptimizer(double[] target,
+ RealMatrix weight,
+ RealMatrix weightSqrt,
+ MultivariateVectorFunction model,
+ MultivariateMatrixFunction jacobian,
+ ConvergenceChecker<PointVectorValuePair> checker,
+ double[] start,
+ int maxEval,
+ int maxIter,
+ double initialStepBoundFactor,
+ double costRelativeTolerance,
+ double parRelativeTolerance,
+ double orthoTolerance,
+ double threshold) {
+ super(target, weight, weightSqrt, model, jacobian, checker, start, maxEval, maxIter);
+
+ this.initialStepBoundFactor = initialStepBoundFactor;
+ this.costRelativeTolerance = costRelativeTolerance;
+ this.parRelativeTolerance = parRelativeTolerance;
+ this.orthoTolerance = orthoTolerance;
+ this.qrRankingThreshold = threshold;
+ }
+
+ /**
+ * Creates a bare-bones instance.
+ * Several calls to {@code withXxx} methods are necessary to obtain
+ * an object with all necessary fields set to sensible values.
+ * <br/>
+ * 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>
+ *
+ * @return an instance of this class.
+ */
+ public static LevenbergMarquardtOptimizer create() {
+ return new LevenbergMarquardtOptimizer(null, null, null, null, null, null, null,
+ 0, 0, 100, 1e-10, 1e-10, 1e-10,
+ Precision.SAFE_MIN);
+ }
+
+ /** {@inheritDoc} */
+ public LevenbergMarquardtOptimizer withTarget(double[] target) {
+ return new LevenbergMarquardtOptimizer(target,
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ initialStepBoundFactor,
+ costRelativeTolerance,
+ parRelativeTolerance,
+ orthoTolerance,
+ qrRankingThreshold);
+ }
+
+ /** {@inheritDoc} */
+ public LevenbergMarquardtOptimizer withWeight(RealMatrix weight) {
+ return new LevenbergMarquardtOptimizer(getTargetInternal(),
+ weight,
+ null,
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ initialStepBoundFactor,
+ costRelativeTolerance,
+ parRelativeTolerance,
+ orthoTolerance,
+ qrRankingThreshold);
+ }
+
+ /** {@inheritDoc} */
+ public LevenbergMarquardtOptimizer withModelAndJacobian(MultivariateVectorFunction model,
+ MultivariateMatrixFunction jacobian) {
+ return new LevenbergMarquardtOptimizer(getTargetInternal(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ model,
+ jacobian,
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ initialStepBoundFactor,
+ costRelativeTolerance,
+ parRelativeTolerance,
+ orthoTolerance,
+ qrRankingThreshold);
+ }
+
+ /** {@inheritDoc} */
+ public LevenbergMarquardtOptimizer withConvergenceChecker(ConvergenceChecker<PointVectorValuePair> checker) {
+ return new LevenbergMarquardtOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ checker,
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ initialStepBoundFactor,
+ costRelativeTolerance,
+ parRelativeTolerance,
+ orthoTolerance,
+ qrRankingThreshold);
+ }
+
+ /** {@inheritDoc} */
+ public LevenbergMarquardtOptimizer withStartPoint(double[] start) {
+ return new LevenbergMarquardtOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ start,
+ getMaxEvaluations(),
+ getMaxIterations(),
+ initialStepBoundFactor,
+ costRelativeTolerance,
+ parRelativeTolerance,
+ orthoTolerance,
+ qrRankingThreshold);
+ }
+
+ /** {@inheritDoc} */
+ public LevenbergMarquardtOptimizer withMaxIterations(int maxIter) {
+ return new LevenbergMarquardtOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ maxIter,
+ initialStepBoundFactor,
+ costRelativeTolerance,
+ parRelativeTolerance,
+ orthoTolerance,
+ qrRankingThreshold);
+ }
+
+ /** {@inheritDoc} */
+ public LevenbergMarquardtOptimizer withMaxEvaluations(int maxEval) {
+ return new LevenbergMarquardtOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ maxEval,
+ getMaxIterations(),
+ initialStepBoundFactor,
+ costRelativeTolerance,
+ parRelativeTolerance,
+ orthoTolerance,
+ qrRankingThreshold);
+ }
+
+ /**
+ * Creates a new instance.
+ *
+ * @param initStepBoundFactor 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 initialStepBoundFactor}
+ * itself. In most cases factor should lie in the interval
+ * {@code (0.1, 100.0)}. {@code 100} is a generally recommended value.
+ * @param costRelTol Desired relative error in the sum of squares.
+ * @param parRelTol Desired relative error in the approximate solution
+ * parameters.
+ * @param orthoTol Desired max cosine on the orthogonality between
+ * the function vector and the columns of the Jacobian.
+ * @param threshold 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 with all fields identical to this instance except
+ * for the givens arguments.
+ */
+ public LevenbergMarquardtOptimizer withTuningParameters(double initStepBoundFactor,
+ double costRelTol,
+ double parRelTol,
+ double orthoTol,
+ double threshold) {
+ return new LevenbergMarquardtOptimizer(getTarget(),
+ getWeightInternal(),
+ getWeightSquareRootInternal(),
+ getModel(),
+ getJacobian(),
+ getConvergenceChecker(),
+ getStart(),
+ getMaxEvaluations(),
+ getMaxIterations(),
+ initStepBoundFactor,
+ costRelTol,
+ parRelTol,
+ orthoTol,
+ threshold);
+ }
+
+ /**
+ * Gets the value of a tuning parameter.
+ * @see #withTuningParameters(double,double,double,double,double)
+ *
+ * @return the parameter's value.
+ */
+ public double getInitialStepBoundFactor() {
+ return initialStepBoundFactor;
+ }
+
+ /**
+ * Gets the value of a tuning parameter.
+ * @see #withTuningParameters(double,double,double,double,double)
+ *
+ * @return the parameter's value.
+ */
+ public double getCostRelativeTolerance() {
+ return costRelativeTolerance;
+ }
+
+ /**
+ * Gets the value of a tuning parameter.
+ * @see #withTuningParameters(double,double,double,double,double)
+ *
+ * @return the parameter's value.
+ */
+ public double getParRelativeTolerance() {
+ return parRelativeTolerance;
+ }
+
+ /**
+ * Gets the value of a tuning parameter.
+ * @see #withTuningParameters(double,double,double,double,double)
+ *
+ * @return the parameter's value.
+ */
+ public double getOrthoTolerance() {
+ return orthoTolerance;
+ }
+
+ /**
+ * Gets the value of a tuning parameter.
+ * @see #withTuningParameters(double,double,double,double,double)
+ *
+ * @return the parameter's value.
+ */
+ public double getRankingThreshold() {
+ return qrRankingThreshold;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected PointVectorValuePair doOptimize() {
+ final int nR = getTarget().length; // Number of observed data.
+ final double[] currentPoint = getStart();
+ final int nC = currentPoint.length; // Number of parameters.
+
+ // arrays shared with the other private methods
+ final int solvedCols = FastMath.min(nR, nC);
+ lmDir = new double[nC];
+ 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[] oldObj = new double[nR];
+ double[] qtf = new double[nR];
+ double[] work1 = new double[nC];
+ double[] work2 = new double[nC];
+ double[] work3 = new double[nC];
+
+ final RealMatrix weightMatrixSqrt = getWeightSquareRoot();
+
+ // Evaluate the function at the starting point and calculate its norm.
+ double[] currentObjective = computeObjectiveValue(currentPoint);
+ double[] currentResiduals = computeResiduals(currentObjective);
+ PointVectorValuePair current = new PointVectorValuePair(currentPoint, currentObjective);
+ double currentCost = computeCost(currentResiduals);
+
+ // Outer loop.
+ boolean firstIteration = true;
+ final ConvergenceChecker<PointVectorValuePair> checker = getConvergenceChecker();
+ while (true) {
+ incrementIterationCount();
+
+ final PointVectorValuePair previous = current;
+
+ // QR decomposition of the jacobian matrix
+ final InternalData internalData = qrDecomposition(computeWeightedJacobian(currentPoint),
+ solvedCols);
+ final double[][] weightedJacobian = internalData.weightedJacobian;
+ final int[] permutation = internalData.permutation;
+ final double[] diagR = internalData.diagR;
+ final double[] jacNorm = internalData.jacNorm;
+
+ double[] weightedResidual = weightMatrixSqrt.operate(currentResiduals);
+ 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 current;
+ }
+
+ // 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;
+ tmpVec = currentObjective;
+ currentObjective = oldObj;
+ oldObj = tmpVec;
+
+ // determine the Levenberg-Marquardt parameter
+ determineLMParameter(qtf, delta, diag,
+ internalData, solvedCols,
+ work1, work2, work3);
+
+ // 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.
+ currentObjective = computeObjectiveValue(currentPoint);
+ currentResiduals = computeResiduals(currentObjective);
+ current = new PointVectorValuePair(currentPoint, currentObjective);
+ currentCost = computeCost(currentResiduals);
+
+ // 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 = 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(getIterations(), previous, current)) {
+ return current;
+ }
+ } 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;
+ tmpVec = currentObjective;
+ currentObjective = oldObj;
+ oldObj = tmpVec;
+ // Reset "current" to previous values.
+ current = new PointVectorValuePair(currentPoint, currentObjective);
+ }
+
+ // Default convergence criteria.
+ if ((FastMath.abs(actRed) <= costRelativeTolerance &&
+ preRed <= costRelativeTolerance &&
+ ratio <= 2.0) ||
+ delta <= parRelativeTolerance * xNorm) {
+ return current;
+ }
+
+ // 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. */
+ final double[][] weightedJacobian;
+ /** Columns permutation array. */
+ final int[] permutation;
+ /** Rank of the Jacobian matrix. */
+ final int rank;
+ /** Diagonal elements of the R matrix in the QR decomposition. */
+ final double[] diagR;
+ /** Norms of the columns of the jacobian matrix. */
+ final double[] jacNorm;
+ /** Coefficients of the Householder transforms vectors. */
+ 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
+ */
+ private void determineLMParameter(double[] qy, double delta, double[] diag,
+ InternalData internalData, int solvedCols,
+ double[] work1, double[] work2, double[] work3) {
+ 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;
+ }
+
+ // 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);
+
+ 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;
+ }
+
+ // 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;
+ }
+
+ /**
+ * 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
+ */
+ private void determineLMDirection(double[] qy, double[] diag,
+ double[] lmDiag,
+ InternalData internalData,
+ int solvedCols,
+ double[] work) {
+ 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];
+ }
+ }
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java
------------------------------------------------------------------------------
svn:keywords = Id Revision