Java tutorial
/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math3.fitting.leastsquares; import org.apache.commons.math3.exception.ConvergenceException; import org.apache.commons.math3.exception.NullArgumentException; import org.apache.commons.math3.exception.util.LocalizedFormats; import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.CholeskyDecomposition; import org.apache.commons.math3.linear.LUDecomposition; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException; import org.apache.commons.math3.linear.QRDecomposition; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; import org.apache.commons.math3.linear.SingularMatrixException; import org.apache.commons.math3.linear.SingularValueDecomposition; import org.apache.commons.math3.optim.ConvergenceChecker; import org.apache.commons.math3.util.Incrementor; import org.apache.commons.math3.util.Pair; /** * Gauss-Newton least-squares solver. * <p> This class solve a least-square problem by * solving the normal equations of the linearized problem at each iteration. Either LU * decomposition or Cholesky decomposition can be used to solve the normal equations, * or QR decomposition or SVD decomposition can be used to solve the linear system. LU * decomposition is faster but QR decomposition is more robust for difficult problems, * and SVD can compute a solution for rank-deficient problems. * </p> * * @since 3.3 */ public class GaussNewtonOptimizer implements LeastSquaresOptimizer { /** The decomposition algorithm to use to solve the normal equations. */ //TODO move to linear package and expand options? public static enum Decomposition { /** * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and * using the {@link LUDecomposition}. * * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the * normal matrix and n<sup>3</sup>/3 operations (m > n) to solve the system using * the LU decomposition. </p> */ LU { @Override protected RealVector solve(final RealMatrix jacobian, final RealVector residuals) { try { final Pair<RealMatrix, RealVector> normalEquation = computeNormalMatrix(jacobian, residuals); final RealMatrix normal = normalEquation.getFirst(); final RealVector jTr = normalEquation.getSecond(); return new LUDecomposition(normal, SINGULARITY_THRESHOLD).getSolver().solve(jTr); } catch (SingularMatrixException e) { throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); } } }, /** * Solve the linear least squares problem (Jx=r) using the {@link * QRDecomposition}. * * <p> Theoretically this method takes mn<sup>2</sup> - n<sup>3</sup>/3 operations * (m > n) and has better numerical accuracy than any method that forms the normal * equations. </p> */ QR { @Override protected RealVector solve(final RealMatrix jacobian, final RealVector residuals) { try { return new QRDecomposition(jacobian, SINGULARITY_THRESHOLD).getSolver().solve(residuals); } catch (SingularMatrixException e) { throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); } } }, /** * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and * using the {@link CholeskyDecomposition}. * * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the * normal matrix and n<sup>3</sup>/6 operations (m > n) to solve the system using * the Cholesky decomposition. </p> */ CHOLESKY { @Override protected RealVector solve(final RealMatrix jacobian, final RealVector residuals) { try { final Pair<RealMatrix, RealVector> normalEquation = computeNormalMatrix(jacobian, residuals); final RealMatrix normal = normalEquation.getFirst(); final RealVector jTr = normalEquation.getSecond(); return new CholeskyDecomposition(normal, SINGULARITY_THRESHOLD, SINGULARITY_THRESHOLD) .getSolver().solve(jTr); } catch (NonPositiveDefiniteMatrixException e) { throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); } } }, /** * Solve the linear least squares problem using the {@link * SingularValueDecomposition}. * * <p> This method is slower, but can provide a solution for rank deficient and * nearly singular systems. */ SVD { @Override protected RealVector solve(final RealMatrix jacobian, final RealVector residuals) { return new SingularValueDecomposition(jacobian).getSolver().solve(residuals); } }; /** * Solve the linear least squares problem Jx=r. * * @param jacobian the Jacobian matrix, J. the number of rows >= the number or * columns. * @param residuals the computed residuals, r. * @return the solution x, to the linear least squares problem Jx=r. * @throws ConvergenceException if the matrix properties (e.g. singular) do not * permit a solution. */ protected abstract RealVector solve(RealMatrix jacobian, RealVector residuals); } /** * The singularity threshold for matrix decompositions. Determines when a {@link * ConvergenceException} is thrown. The current value was the default value for {@link * LUDecomposition}. */ private static final double SINGULARITY_THRESHOLD = 1e-11; /** Indicator for using LU decomposition. */ private final Decomposition decomposition; /** * Creates a Gauss Newton optimizer. * <p/> * The default for the algorithm is to solve the normal equations using QR * decomposition. */ public GaussNewtonOptimizer() { this(Decomposition.QR); } /** * Create a Gauss Newton optimizer that uses the given decomposition algorithm to * solve the normal equations. * * @param decomposition the {@link Decomposition} algorithm. */ public GaussNewtonOptimizer(final Decomposition decomposition) { this.decomposition = decomposition; } /** * Get the matrix decomposition algorithm used to solve the normal equations. * * @return the matrix {@link Decomposition} algoritm. */ public Decomposition getDecomposition() { return this.decomposition; } /** * Configure the decomposition algorithm. * * @param newDecomposition the {@link Decomposition} algorithm to use. * @return a new instance. */ public GaussNewtonOptimizer withDecomposition(final Decomposition newDecomposition) { return new GaussNewtonOptimizer(newDecomposition); } /** {@inheritDoc} */ public Optimum optimize(final LeastSquaresProblem lsp) { //create local evaluation and iteration counts final Incrementor evaluationCounter = lsp.getEvaluationCounter(); final Incrementor iterationCounter = lsp.getIterationCounter(); final ConvergenceChecker<Evaluation> checker = lsp.getConvergenceChecker(); // Computation will be useless without a checker (see "for-loop"). if (checker == null) { throw new NullArgumentException(); } RealVector currentPoint = lsp.getStart(); // iterate until convergence is reached Evaluation current = null; while (true) { iterationCounter.incrementCount(); // evaluate the objective function and its jacobian Evaluation previous = current; // Value of the objective function at "currentPoint". evaluationCounter.incrementCount(); current = lsp.evaluate(currentPoint); final RealVector currentResiduals = current.getResiduals(); final RealMatrix weightedJacobian = current.getJacobian(); currentPoint = current.getPoint(); // Check convergence. if (previous != null) { if (checker.converged(iterationCounter.getCount(), previous, current)) { return new OptimumImpl(current, evaluationCounter.getCount(), iterationCounter.getCount()); } } // solve the linearized least squares problem final RealVector dX = this.decomposition.solve(weightedJacobian, currentResiduals); // update the estimated parameters currentPoint = currentPoint.add(dX); } } @Override public String toString() { return "GaussNewtonOptimizer{" + "decomposition=" + decomposition + '}'; } /** * Compute the normal matrix, J<sup>T</sup>J. * * @param jacobian the m by n jacobian matrix, J. Input. * @param residuals the m by 1 residual vector, r. Input. * @return the n by n normal matrix and the n by 1 J<sup>Tr vector. */ private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian, final RealVector residuals) { //since the normal matrix is symmetric, we only need to compute half of it. final int nR = jacobian.getRowDimension(); final int nC = jacobian.getColumnDimension(); //allocate space for return values final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC); final RealVector jTr = new ArrayRealVector(nC); //for each measurement for (int i = 0; i < nR; ++i) { //compute JTr for measurement i for (int j = 0; j < nC; j++) { jTr.setEntry(j, jTr.getEntry(j) + residuals.getEntry(i) * jacobian.getEntry(i, j)); } // add the the contribution to the normal matrix for measurement i for (int k = 0; k < nC; ++k) { //only compute the upper triangular part for (int l = k; l < nC; ++l) { normal.setEntry(k, l, normal.getEntry(k, l) + jacobian.getEntry(i, k) * jacobian.getEntry(i, l)); } } } //copy the upper triangular part to the lower triangular part. for (int i = 0; i < nC; i++) { for (int j = 0; j < i; j++) { normal.setEntry(i, j, normal.getEntry(j, i)); } } return new Pair<RealMatrix, RealVector>(normal, jTr); } }