Top Description Inners Fields Constructors Methods
org.apache.commons.math3.fitting.leastsquares

public Class LevenbergMarquardtOptimizer

extends Object
implements LeastSquaresOptimizer
Class Inheritance
All Implemented Interfaces
org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer
Imports
java.util.Arrays, org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation, org.apache.commons.math3.linear.ArrayRealVector, .RealMatrix, org.apache.commons.math3.exception.ConvergenceException, org.apache.commons.math3.exception.util.LocalizedFormats, org.apache.commons.math3.optim.ConvergenceChecker, org.apache.commons.math3.util.Incrementor, .Precision, .FastMath
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:

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

    Nested and Inner Type Summary

    Modifier and TypeClass and Description
    private static class
    LevenbergMarquardtOptimizer.InternalData

    Holds internal data.

    Field Summary

    Modifier and TypeField and Description
    private final double
    costRelativeTolerance

    Desired relative error in the sum of squares.

    private final double
    initialStepBoundFactor

    Positive input variable used in determining the initial step bound.

    private final double
    orthoTolerance

    Desired max cosine on the orthogonality between the function vector and the columns of the jacobian.

    private final double
    parRelativeTolerance

    Desired relative error in the approximate solution parameters.

    private final double
    qrRankingThreshold

    Threshold for QR ranking.

    private static final double
    TWO_EPS

    Twice the "epsilon machine".

    Constructor Summary

    AccessConstructor and Description
    public
    LevenbergMarquardtOptimizer()

    Default constructor.

    public
    LevenbergMarquardtOptimizer(final double
    initial step bound factor
    initialStepBoundFactor    ,     final double
    cost relative tolerance
    costRelativeTolerance    ,     final double
    parameters relative tolerance
    parRelativeTolerance    ,     final double
    orthogonality tolerance
    orthoTolerance    ,     final double
    threshold in the QR decomposition. Columns with a 2 norm less than this threshold are considered to be all 0s.
    qrRankingThreshold    )

    Construct an instance with all parameters specified.

    Method Summary

    Modifier and TypeMethod and Description
    private void
    determineLMDirection(double[]
    array containing qTy
    qy    ,     double[]
    diagonal matrix
    diag    ,     double[]
    diagonal elements associated with lmDir
    lmDiag    ,     LevenbergMarquardtOptimizer.InternalData
    Data (modified in-place in this method).
    internalData    ,     int
    Number of sloved point.
    solvedCols    ,     double[]
    work array
    work    ,     double[]
    the "returned" LM direction is stored in this array
    lmDir    )

    Solve a*x = b and d*x = 0 in the least squares sense.

    private double

    Returns:

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

    Determines the Levenberg-Marquardt parameter.

    public double

    Returns:

    the parameter's value.
    getCostRelativeTolerance    ()

    Gets the value of a tuning parameter.

    public double

    Returns:

    the parameter's value.
    getInitialStepBoundFactor    ()

    Gets the value of a tuning parameter.

    public double

    Returns:

    the parameter's value.
    getOrthoTolerance    ()

    Gets the value of a tuning parameter.

    public double

    Returns:

    the parameter's value.
    getParameterRelativeTolerance    ()

    Gets the value of a tuning parameter.

    public double

    Returns:

    the parameter's value.
    getRankingThreshold    ()

    Gets the value of a tuning parameter.

    public LeastSquaresOptimizer.Optimum
    optimize(final LeastSquaresProblem
    the problem definition, including model function and convergence criteria.
    problem    )

    Implements org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.optimize.

    Solve the non-linear least squares problem.

    private LevenbergMarquardtOptimizer.InternalData

    Returns:

    data used in other methods of this class.
    qrDecomposition    (RealMatrix
    Weighted Jacobian matrix at the current point.
    jacobian    ,     int
    Number of solved point.
    solvedCols    )

    Decompose a matrix A as A.P = Q.R using Householder transforms.

    private void
    qTy(double[]
    vector to multiply (will be overwritten with the result)
    y    ,     LevenbergMarquardtOptimizer.InternalData
    Data.
    internalData    )

    Compute the product Qt.y for some Q.R. decomposition.

    public LevenbergMarquardtOptimizer

    Returns:

    a new instance.
    withCostRelativeTolerance    (double
    Desired relative error in the sum of squares.
    newCostRelativeTolerance    )

    public LevenbergMarquardtOptimizer

    Returns:

    a new instance.
    withInitialStepBoundFactor    (double
    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.
    newInitialStepBoundFactor    )

    public LevenbergMarquardtOptimizer

    Returns:

    a new instance.
    withOrthoTolerance    (double
    Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian.
    newOrthoTolerance    )

    Modifies the given parameter.

    public LevenbergMarquardtOptimizer

    Returns:

    a new instance.
    withParameterRelativeTolerance    (double
    Desired relative error in the approximate solution parameters.
    newParRelativeTolerance    )

    public LevenbergMarquardtOptimizer

    Returns:

    a new instance.
    withRankingThreshold    (double
    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.
    newQRRankingThreshold    )

    Inherited from java.lang.Object:
    cloneequalsfinalizegetClasshashCodenotifynotifyAlltoStringwaitwaitwait

    Field Detail

    costRelativeToleranceback to summary
    private final double costRelativeTolerance

    Desired relative error in the sum of squares.

    initialStepBoundFactorback to summary
    private final double initialStepBoundFactor

    Positive input variable used in determining the initial step bound.

    orthoToleranceback to summary
    private final double orthoTolerance

    Desired max cosine on the orthogonality between the function vector and the columns of the jacobian.

    parRelativeToleranceback to summary
    private final double parRelativeTolerance

    Desired relative error in the approximate solution parameters.

    qrRankingThresholdback to summary
    private final double qrRankingThreshold

    Threshold for QR ranking.

    TWO_EPSback to summary
    private static final double TWO_EPS

    Twice the "epsilon machine".

    Constructor Detail

    LevenbergMarquardtOptimizerback to summary
    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
    LevenbergMarquardtOptimizerback to summary
    public LevenbergMarquardtOptimizer(final double initialStepBoundFactor, final double costRelativeTolerance, final double parRelativeTolerance, final double orthoTolerance, final double qrRankingThreshold)

    Construct an instance with all parameters specified.

    Parameters
    initialStepBoundFactor:double

    initial step bound factor

    costRelativeTolerance:double

    cost relative tolerance

    parRelativeTolerance:double

    parameters relative tolerance

    orthoTolerance:double

    orthogonality tolerance

    qrRankingThreshold:double

    threshold in the QR decomposition. Columns with a 2 norm less than this threshold are considered to be all 0s.

    Method Detail

    determineLMDirectionback to summary
    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:double[]

    array containing qTy

    diag:double[]

    diagonal matrix

    lmDiag:double[]

    diagonal elements associated with lmDir

    internalData:LevenbergMarquardtOptimizer.InternalData

    Data (modified in-place in this method).

    solvedCols:int

    Number of sloved point.

    work:double[]

    work array

    lmDir:double[]

    the "returned" LM direction is stored in this array

    determineLMParameterback to summary
    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:double[]

    Array containing qTy.

    delta:double

    Upper bound on the euclidean norm of diagR * lmDir.

    diag:double[]

    Diagonal matrix.

    internalData:LevenbergMarquardtOptimizer.InternalData

    Data (modified in-place in this method).

    solvedCols:int

    Number of solved point.

    work1:double[]

    work array

    work2:double[]

    work array

    work3:double[]

    work array

    lmDir:double[]

    the "returned" LM direction will be stored in this array.

    lmPar:double

    the value of the LM parameter from the previous iteration.

    Returns:double

    the new LM parameter

    getCostRelativeToleranceback to summary
    public double getCostRelativeTolerance()

    Gets the value of a tuning parameter.

    Returns:double

    the parameter's value.

    See Also
    withCostRelativeTolerance(double)
    getInitialStepBoundFactorback to summary
    public double getInitialStepBoundFactor()

    Gets the value of a tuning parameter.

    Returns:double

    the parameter's value.

    See Also
    withInitialStepBoundFactor(double)
    getOrthoToleranceback to summary
    public double getOrthoTolerance()

    Gets the value of a tuning parameter.

    Returns:double

    the parameter's value.

    See Also
    withOrthoTolerance(double)
    getParameterRelativeToleranceback to summary
    public double getParameterRelativeTolerance()

    Gets the value of a tuning parameter.

    Returns:double

    the parameter's value.

    See Also
    withParameterRelativeTolerance(double)
    getRankingThresholdback to summary
    public double getRankingThreshold()

    Gets the value of a tuning parameter.

    Returns:double

    the parameter's value.

    See Also
    withRankingThreshold(double)
    optimizeback to summary
    public LeastSquaresOptimizer.Optimum optimize(final LeastSquaresProblem problem)

    Implements org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.optimize.

    Doc from org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.optimize.

    Solve the non-linear least squares problem.

    Parameters
    problem:LeastSquaresProblem

    the problem definition, including model function and convergence criteria.

    Returns:LeastSquaresOptimizer.Optimum

    The optimum.

    qrDecompositionback to summary
    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:RealMatrix

    Weighted Jacobian matrix at the current point.

    solvedCols:int

    Number of solved point.

    Returns:LevenbergMarquardtOptimizer.InternalData

    data used in other methods of this class.

    Exceptions
    ConvergenceException:
    if the decomposition cannot be performed.
    qTyback to summary
    private void qTy(double[] y, LevenbergMarquardtOptimizer.InternalData internalData)

    Compute the product Qt.y for some Q.R. decomposition.

    Parameters
    y:double[]

    vector to multiply (will be overwritten with the result)

    internalData:LevenbergMarquardtOptimizer.InternalData

    Data.

    withCostRelativeToleranceback to summary
    public LevenbergMarquardtOptimizer withCostRelativeTolerance(double newCostRelativeTolerance)
    Parameters
    newCostRelativeTolerance:double

    Desired relative error in the sum of squares.

    Returns:LevenbergMarquardtOptimizer

    a new instance.

    withInitialStepBoundFactorback to summary
    public LevenbergMarquardtOptimizer withInitialStepBoundFactor(double newInitialStepBoundFactor)
    Parameters
    newInitialStepBoundFactor:double

    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:LevenbergMarquardtOptimizer

    a new instance.

    withOrthoToleranceback to summary
    public LevenbergMarquardtOptimizer withOrthoTolerance(double newOrthoTolerance)

    Modifies the given parameter.

    Parameters
    newOrthoTolerance:double

    Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian.

    Returns:LevenbergMarquardtOptimizer

    a new instance.

    withParameterRelativeToleranceback to summary
    public LevenbergMarquardtOptimizer withParameterRelativeTolerance(double newParRelativeTolerance)
    Parameters
    newParRelativeTolerance:double

    Desired relative error in the approximate solution parameters.

    Returns:LevenbergMarquardtOptimizer

    a new instance.

    withRankingThresholdback to summary
    public LevenbergMarquardtOptimizer withRankingThreshold(double newQRRankingThreshold)
    Parameters
    newQRRankingThreshold:double

    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:LevenbergMarquardtOptimizer

    a new instance.

    org.apache.commons.math3.fitting.leastsquares back to summary

    private Class LevenbergMarquardtOptimizer.InternalData

    extends Object
    Class Inheritance
    • java.lang.Object
    • org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer.InternalData
    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).

    Field Summary

    Modifier and TypeField and Description
    private final double[]
    beta

    Coefficients of the Householder transforms vectors.

    private final double[]
    diagR

    Diagonal elements of the R matrix in the QR decomposition.

    private final double[]
    jacNorm

    Norms of the columns of the jacobian matrix.

    private final int[]
    permutation

    Columns permutation array.

    private final int
    rank

    Rank of the Jacobian matrix.

    private final double[][]
    weightedJacobian

    Weighted Jacobian.

    Constructor Summary

    AccessConstructor and Description
    pack-priv
    InternalData(double[][]
    Weighted Jacobian.
    weightedJacobian    ,     int[]
    Columns permutation array.
    permutation    ,     int
    Rank of the Jacobian matrix.
    rank    ,     double[]
    Diagonal elements of the R matrix in the QR decomposition.
    diagR    ,     double[]
    Norms of the columns of the jacobian matrix.
    jacNorm    ,     double[]
    Coefficients of the Householder transforms vectors.
    beta    )

    Method Summary

    Inherited from java.lang.Object:
    cloneequalsfinalizegetClasshashCodenotifynotifyAlltoStringwaitwaitwait

    Field Detail

    betaback to summary
    private final double[] beta

    Coefficients of the Householder transforms vectors.

    diagRback to summary
    private final double[] diagR

    Diagonal elements of the R matrix in the QR decomposition.

    jacNormback to summary
    private final double[] jacNorm

    Norms of the columns of the jacobian matrix.

    permutationback to summary
    private final int[] permutation

    Columns permutation array.

    rankback to summary
    private final int rank

    Rank of the Jacobian matrix.

    weightedJacobianback to summary
    private final double[][] weightedJacobian

    Weighted Jacobian.

    Constructor Detail

    InternalDataback to summary
    pack-priv InternalData(double[][] weightedJacobian, int[] permutation, int rank, double[] diagR, double[] jacNorm, double[] beta)
    Parameters
    weightedJacobian:double[][]

    Weighted Jacobian.

    permutation:int[]

    Columns permutation array.

    rank:int

    Rank of the Jacobian matrix.

    diagR:double[]

    Diagonal elements of the R matrix in the QR decomposition.

    jacNorm:double[]

    Norms of the columns of the jacobian matrix.

    beta:double[]

    Coefficients of the Householder transforms vectors.