GaussNewtonOptimizer.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.legacy.fitting.leastsquares;

  19. import org.apache.commons.math4.legacy.exception.ConvergenceException;
  20. import org.apache.commons.math4.legacy.exception.NullArgumentException;
  21. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  22. import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem.Evaluation;
  23. import org.apache.commons.math4.legacy.linear.ArrayRealVector;
  24. import org.apache.commons.math4.legacy.linear.CholeskyDecomposition;
  25. import org.apache.commons.math4.legacy.linear.LUDecomposition;
  26. import org.apache.commons.math4.legacy.linear.MatrixUtils;
  27. import org.apache.commons.math4.legacy.linear.NonPositiveDefiniteMatrixException;
  28. import org.apache.commons.math4.legacy.linear.QRDecomposition;
  29. import org.apache.commons.math4.legacy.linear.RealMatrix;
  30. import org.apache.commons.math4.legacy.linear.RealVector;
  31. import org.apache.commons.math4.legacy.linear.SingularMatrixException;
  32. import org.apache.commons.math4.legacy.linear.SingularValueDecomposition;
  33. import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
  34. import org.apache.commons.math4.legacy.core.IntegerSequence;
  35. import org.apache.commons.math4.legacy.core.Pair;

  37. /**
  38.  * Gauss-Newton least-squares solver.
  39.  * <p> This class solve a least-square problem by
  40.  * solving the normal equations of the linearized problem at each iteration. Either LU
  41.  * decomposition or Cholesky decomposition can be used to solve the normal equations,
  42.  * or QR decomposition or SVD decomposition can be used to solve the linear system. LU
  43.  * decomposition is faster but QR decomposition is more robust for difficult problems,
  44.  * and SVD can compute a solution for rank-deficient problems.
  45.  * </p>
  46.  *
  47.  * @since 3.3
  48.  */
  49. public class GaussNewtonOptimizer implements LeastSquaresOptimizer {

  51.     /** The decomposition algorithm to use to solve the normal equations. */
  52.     //TODO move to linear package and expand options?
  53.     public enum Decomposition {
  54.         /**
  55.          * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and
  56.          * using the {@link LUDecomposition}.
  57.          *
  58.          * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the
  59.          * normal matrix and n<sup>3</sup>/3 operations (m &gt; n) to solve the system using
  60.          * the LU decomposition. </p>
  61.          */
  62.         LU {
  63.             @Override
  64.             protected RealVector solve(final RealMatrix jacobian,
  65.                                        final RealVector residuals) {
  66.                 try {
  67.                     final Pair<RealMatrix, RealVector> normalEquation =
  68.                             computeNormalMatrix(jacobian, residuals);
  69.                     final RealMatrix normal = normalEquation.getFirst();
  70.                     final RealVector jTr = normalEquation.getSecond();
  71.                     return new LUDecomposition(normal, SINGULARITY_THRESHOLD)
  72.                             .getSolver()
  73.                             .solve(jTr);
  74.                 } catch (SingularMatrixException e) {
  75.                     throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
  76.                 }
  77.             }
  78.         },
  79.         /**
  80.          * Solve the linear least squares problem (Jx=r) using the {@link
  81.          * QRDecomposition}.
  82.          *
  83.          * <p> Theoretically this method takes mn<sup>2</sup> - n<sup>3</sup>/3 operations
  84.          * (m &gt; n) and has better numerical accuracy than any method that forms the normal
  85.          * equations. </p>
  86.          */
  87.         QR {
  88.             @Override
  89.             protected RealVector solve(final RealMatrix jacobian,
  90.                                        final RealVector residuals) {
  91.                 try {
  92.                     return new QRDecomposition(jacobian, SINGULARITY_THRESHOLD)
  93.                             .getSolver()
  94.                             .solve(residuals);
  95.                 } catch (SingularMatrixException e) {
  96.                     throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
  97.                 }
  98.             }
  99.         },
  100.         /**
  101.          * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and
  102.          * using the {@link CholeskyDecomposition}.
  103.          *
  104.          * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the
  105.          * normal matrix and n<sup>3</sup>/6 operations (m &gt; n) to solve the system using
  106.          * the Cholesky decomposition. </p>
  107.          */
  108.         CHOLESKY {
  109.             @Override
  110.             protected RealVector solve(final RealMatrix jacobian,
  111.                                        final RealVector residuals) {
  112.                 try {
  113.                     final Pair<RealMatrix, RealVector> normalEquation =
  114.                             computeNormalMatrix(jacobian, residuals);
  115.                     final RealMatrix normal = normalEquation.getFirst();
  116.                     final RealVector jTr = normalEquation.getSecond();
  117.                     return new CholeskyDecomposition(
  118.                             normal, SINGULARITY_THRESHOLD, SINGULARITY_THRESHOLD)
  119.                             .getSolver()
  120.                             .solve(jTr);
  121.                 } catch (NonPositiveDefiniteMatrixException e) {
  122.                     throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
  123.                 }
  124.             }
  125.         },
  126.         /**
  127.          * Solve the linear least squares problem using the {@link
  128.          * SingularValueDecomposition}.
  129.          *
  130.          * <p> This method is slower, but can provide a solution for rank deficient and
  131.          * nearly singular systems.
  132.          */
  133.         SVD {
  134.             @Override
  135.             protected RealVector solve(final RealMatrix jacobian,
  136.                                        final RealVector residuals) {
  137.                 return new SingularValueDecomposition(jacobian)
  138.                         .getSolver()
  139.                         .solve(residuals);
  140.             }
  141.         };

  143.         /**
  144.          * Solve the linear least squares problem Jx=r.
  145.          *
  146.          * @param jacobian  the Jacobian matrix, J. the number of rows &gt;= the number or
  147.          *                  columns.
  148.          * @param residuals the computed residuals, r.
  149.          * @return the solution x, to the linear least squares problem Jx=r.
  150.          * @throws ConvergenceException if the matrix properties (e.g. singular) do not
  151.          *                              permit a solution.
  152.          */
  153.         protected abstract RealVector solve(RealMatrix jacobian,
  154.                                             RealVector residuals);
  155.     }

  157.     /**
  158.      * The singularity threshold for matrix decompositions. Determines when a {@link
  159.      * ConvergenceException} is thrown. The current value was the default value for {@link
  160.      * LUDecomposition}.
  161.      */
  162.     private static final double SINGULARITY_THRESHOLD = 1e-11;

  164.     /** Indicator for using LU decomposition. */
  165.     private final Decomposition decomposition;

  167.     /**
  168.      * Creates a Gauss Newton optimizer.
  169.      * <p>
  170.      * The default for the algorithm is to solve the normal equations using QR
  171.      * decomposition.
  172.      */
  173.     public GaussNewtonOptimizer() {
  174.         this(Decomposition.QR);
  175.     }

  177.     /**
  178.      * Create a Gauss Newton optimizer that uses the given decomposition algorithm to
  179.      * solve the normal equations.
  180.      *
  181.      * @param decomposition the {@link Decomposition} algorithm.
  182.      */
  183.     public GaussNewtonOptimizer(final Decomposition decomposition) {
  184.         this.decomposition = decomposition;
  185.     }

  187.     /**
  188.      * Get the matrix decomposition algorithm used to solve the normal equations.
  189.      *
  190.      * @return the matrix {@link Decomposition} algoritm.
  191.      */
  192.     public Decomposition getDecomposition() {
  193.         return this.decomposition;
  194.     }

  196.     /**
  197.      * Configure the decomposition algorithm.
  198.      *
  199.      * @param newDecomposition the {@link Decomposition} algorithm to use.
  200.      * @return a new instance.
  201.      */
  202.     public GaussNewtonOptimizer withDecomposition(final Decomposition newDecomposition) {
  203.         return new GaussNewtonOptimizer(newDecomposition);
  204.     }

  206.     /** {@inheritDoc} */
  207.     @Override
  208.     public Optimum optimize(final LeastSquaresProblem lsp) {
  209.         //create local evaluation and iteration counts
  210.         final IntegerSequence.Incrementor evaluationCounter = lsp.getEvaluationCounter();
  211.         final IntegerSequence.Incrementor iterationCounter = lsp.getIterationCounter();
  212.         final ConvergenceChecker<Evaluation> checker
  213.                 = lsp.getConvergenceChecker();

  215.         // Computation will be useless without a checker (see "for-loop").
  216.         if (checker == null) {
  217.             throw new NullArgumentException();
  218.         }

  220.         RealVector currentPoint = lsp.getStart();

  222.         // iterate until convergence is reached
  223.         Evaluation current = null;
  224.         while (true) {
  225.             iterationCounter.increment();

  227.             // evaluate the objective function and its jacobian
  228.             Evaluation previous = current;
  229.             // Value of the objective function at "currentPoint".
  230.             evaluationCounter.increment();
  231.             current = lsp.evaluate(currentPoint);
  232.             final RealVector currentResiduals = current.getResiduals();
  233.             final RealMatrix weightedJacobian = current.getJacobian();
  234.             currentPoint = current.getPoint();

  236.             // Check convergence.
  237.             if (previous != null &&
  238.                 checker.converged(iterationCounter.getCount(), previous, current)) {
  239.                 return new OptimumImpl(current,
  240.                                        evaluationCounter.getCount(),
  241.                                        iterationCounter.getCount());
  242.             }

  244.             // solve the linearized least squares problem
  245.             final RealVector dX = this.decomposition.solve(weightedJacobian, currentResiduals);
  246.             // update the estimated parameters
  247.             currentPoint = currentPoint.add(dX);
  248.         }
  249.     }

  251.     /** {@inheritDoc} */
  252.     @Override
  253.     public String toString() {
  254.         return "GaussNewtonOptimizer{" +
  255.                 "decomposition=" + decomposition +
  256.                 '}';
  257.     }

  259.     /**
  260.      * Compute the normal matrix, J<sup>T</sup>J.
  261.      *
  262.      * @param jacobian  the m by n jacobian matrix, J. Input.
  263.      * @param residuals the m by 1 residual vector, r. Input.
  264.      * @return  the n by n normal matrix and  the n by 1 J<sup>Tr</sup> vector.
  265.      */
  266.     private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian,
  267.                                                                     final RealVector residuals) {
  268.         //since the normal matrix is symmetric, we only need to compute half of it.
  269.         final int nR = jacobian.getRowDimension();
  270.         final int nC = jacobian.getColumnDimension();
  271.         //allocate space for return values
  272.         final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC);
  273.         final RealVector jTr = new ArrayRealVector(nC);
  274.         //for each measurement
  275.         for (int i = 0; i < nR; ++i) {
  276.             //compute JTr for measurement i
  277.             for (int j = 0; j < nC; j++) {
  278.                 jTr.setEntry(j, jTr.getEntry(j) +
  279.                         residuals.getEntry(i) * jacobian.getEntry(i, j));
  280.             }

  282.             // add the contribution to the normal matrix for measurement i
  283.             for (int k = 0; k < nC; ++k) {
  284.                 //only compute the upper triangular part
  285.                 for (int l = k; l < nC; ++l) {
  286.                     normal.setEntry(k, l, normal.getEntry(k, l) +
  287.                             jacobian.getEntry(i, k) * jacobian.getEntry(i, l));
  288.                 }
  289.             }
  290.         }
  291.         //copy the upper triangular part to the lower triangular part.
  292.         for (int i = 0; i < nC; i++) {
  293.             for (int j = 0; j < i; j++) {
  294.                 normal.setEntry(i, j, normal.getEntry(j, i));
  295.             }
  296.         }
  297.         return new Pair<>(normal, jTr);
  298.     }
  299. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT