Class LevenbergMarquardtOptimizer

java.lang.Object
org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer
All Implemented Interfaces:
LeastSquaresOptimizer

public class LevenbergMarquardtOptimizer extends Object implements LeastSquaresOptimizer
This class solves a least-squares problem using the Levenberg-Marquardt algorithm.

This implementation should 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.

The resolution engine is a simple translation of the MINPACK lmder 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 Analyse numérique matricielle appliquée à l'art de l'ingénieur, Masson 1986.

The authors of the original fortran version are:

  • Argonne National Laboratory. MINPACK project. March 1980
  • Burton S. Garbow
  • Kenneth E. Hillstrom
  • Jorge J. More
The redistribution policy for MINPACK is available here, for convenience, it is reproduced below.

Minpack Copyright Notice (1999) University of Chicago. All rights reserved
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
  2. 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.
  3. The end-user documentation included with the redistribution, if any, must include the following acknowledgment: This product includes software developed by the University of Chicago, as Operator of Argonne National Laboratory. Alternately, this acknowledgment may appear in the software itself, if and wherever such third-party acknowledgments normally appear.
  4. 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.
  5. 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.
    Since:
    3.3
    • Field Details

      • TWO_EPS

        private static final double TWO_EPS
        Twice the "epsilon machine".
      • initialStepBoundFactor

        private final double initialStepBoundFactor
        Positive input variable used in determining the initial step bound.
      • costRelativeTolerance

        private final double costRelativeTolerance
        Desired relative error in the sum of squares.
      • parRelativeTolerance

        private final double parRelativeTolerance
        Desired relative error in the approximate solution parameters.
      • orthoTolerance

        private final double orthoTolerance
        Desired max cosine on the orthogonality between the function vector and the columns of the jacobian.
      • qrRankingThreshold

        private final double qrRankingThreshold
        Threshold for QR ranking.
    • Constructor Details

      • LevenbergMarquardtOptimizer

        public LevenbergMarquardtOptimizer()
        Default constructor.

        The default values for the algorithm settings are:

        • Initial step bound factor: 100
        • Cost relative tolerance: 1e-10
        • Parameters relative tolerance: 1e-10
        • Orthogonality tolerance: 1e-10
        • QR ranking threshold: Precision.SAFE_MIN
      • LevenbergMarquardtOptimizer

        public LevenbergMarquardtOptimizer(double initialStepBoundFactor, double costRelativeTolerance, double parRelativeTolerance, double orthoTolerance, double qrRankingThreshold)
        Construct an instance with all parameters specified.
        Parameters:
        initialStepBoundFactor - initial step bound factor
        costRelativeTolerance - cost relative tolerance
        parRelativeTolerance - parameters relative tolerance
        orthoTolerance - orthogonality tolerance
        qrRankingThreshold - threshold in the QR decomposition. Columns with a 2 norm less than this threshold are considered to be all 0s.
    • Method Details

      • withInitialStepBoundFactor

        public LevenbergMarquardtOptimizer withInitialStepBoundFactor(double newInitialStepBoundFactor)
        Parameters:
        newInitialStepBoundFactor - Positive input variable used in determining the initial step bound. This bound is set to the product of initialStepBoundFactor and the euclidean norm of diag * x if non-zero, or else to newInitialStepBoundFactor itself. In most cases factor should lie in the interval (0.1, 100.0). 100 is a generally recommended value. of the matrix is reduced.
        Returns:
        a new instance.
      • withCostRelativeTolerance

        public LevenbergMarquardtOptimizer withCostRelativeTolerance(double newCostRelativeTolerance)
        Parameters:
        newCostRelativeTolerance - Desired relative error in the sum of squares.
        Returns:
        a new instance.
      • withParameterRelativeTolerance

        public LevenbergMarquardtOptimizer withParameterRelativeTolerance(double newParRelativeTolerance)
        Parameters:
        newParRelativeTolerance - Desired relative error in the approximate solution parameters.
        Returns:
        a new instance.
      • withOrthoTolerance

        public LevenbergMarquardtOptimizer withOrthoTolerance(double newOrthoTolerance)
        Modifies the given parameter.
        Parameters:
        newOrthoTolerance - Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian.
        Returns:
        a new instance.
      • withRankingThreshold

        public LevenbergMarquardtOptimizer withRankingThreshold(double newQRRankingThreshold)
        Parameters:
        newQRRankingThreshold - Desired threshold for QR ranking. If the squared norm of a column vector is smaller or equal to this threshold during QR decomposition, it is considered to be a zero vector and hence the rank of the matrix is reduced.
        Returns:
        a new instance.
      • getInitialStepBoundFactor

        public double getInitialStepBoundFactor()
        Gets the value of a tuning parameter.
        Returns:
        the parameter's value.
        See Also:
      • getCostRelativeTolerance

        public double getCostRelativeTolerance()
        Gets the value of a tuning parameter.
        Returns:
        the parameter's value.
        See Also:
      • getParameterRelativeTolerance

        public double getParameterRelativeTolerance()
        Gets the value of a tuning parameter.
        Returns:
        the parameter's value.
        See Also:
      • getOrthoTolerance

        public double getOrthoTolerance()
        Gets the value of a tuning parameter.
        Returns:
        the parameter's value.
        See Also:
      • getRankingThreshold

        public double getRankingThreshold()
        Gets the value of a tuning parameter.
        Returns:
        the parameter's value.
        See Also:
      • optimize

        Solve the non-linear least squares problem.
        Specified by:
        optimize in interface LeastSquaresOptimizer
        Parameters:
        problem - the problem definition, including model function and convergence criteria.
        Returns:
        The optimum.
      • determineLMParameter

        private double determineLMParameter(double[] qy, double delta, double[] diag, LevenbergMarquardtOptimizer.InternalData internalData, int solvedCols, double[] work1, double[] work2, double[] work3, double[] lmDir, double lmPar)
        Determines the Levenberg-Marquardt parameter.

        This implementation is a translation in Java of the MINPACK lmpar routine.

        This method sets the lmPar and lmDir attributes.

        The authors of the original fortran function are:

        • Argonne National Laboratory. MINPACK project. March 1980
        • Burton S. Garbow
        • Kenneth E. Hillstrom
        • Jorge J. More

        Luc Maisonobe did the Java translation.

        Parameters:
        qy - Array containing qTy.
        delta - Upper bound on the euclidean norm of diagR * lmDir.
        diag - Diagonal matrix.
        internalData - Data (modified in-place in this method).
        solvedCols - Number of solved point.
        work1 - work array
        work2 - work array
        work3 - work array
        lmDir - the "returned" LM direction will be stored in this array.
        lmPar - the value of the LM parameter from the previous iteration.
        Returns:
        the new LM parameter
      • determineLMDirection

        private void determineLMDirection(double[] qy, double[] diag, double[] lmDiag, LevenbergMarquardtOptimizer.InternalData internalData, int solvedCols, double[] work, double[] lmDir)
        Solve a*x = b and d*x = 0 in the least squares sense.

        This implementation is a translation in Java of the MINPACK qrsolv routine.

        This method sets the lmDir and lmDiag attributes.

        The authors of the original fortran function are:

        • Argonne National Laboratory. MINPACK project. March 1980
        • Burton S. Garbow
        • Kenneth E. Hillstrom
        • Jorge J. More

        Luc Maisonobe did the Java translation.

        Parameters:
        qy - array containing qTy
        diag - diagonal matrix
        lmDiag - diagonal elements associated with lmDir
        internalData - Data (modified in-place in this method).
        solvedCols - Number of sloved point.
        work - work array
        lmDir - the "returned" LM direction is stored in this array
      • qrDecomposition

        private LevenbergMarquardtOptimizer.InternalData qrDecomposition(RealMatrix jacobian, int solvedCols) throws ConvergenceException
        Decompose a matrix A as A.P = Q.R using Householder transforms.

        As suggested in the P. Lascaux and R. Theodor book Analyse numérique matricielle appliquée à l'art de l'ingénieur (Masson, 1986), instead of representing the Householder transforms with uk unit vectors such that:

         Hk = I - 2uk.ukt
         
        we use k non-unit vectors such that:
         Hk = I - betakvk.vkt
         
        where vk = ak - alphak ek. The betak coefficients are provided upon exit as recomputing them from the vk vectors would be costly.

        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.

        Parameters:
        jacobian - Weighted Jacobian matrix at the current point.
        solvedCols - Number of solved point.
        Returns:
        data used in other methods of this class.
        Throws:
        ConvergenceException - if the decomposition cannot be performed.
      • qTy

        private void qTy(double[] y, LevenbergMarquardtOptimizer.InternalData internalData)
        Compute the product Qt.y for some Q.R. decomposition.
        Parameters:
        y - vector to multiply (will be overwritten with the result)
        internalData - Data.