Apache Commons Math 3.6.1 · org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer

This class solves a least-squares problem using the Levenberg-Marquardt algorithm.
Constraints are not supported: the call to optimize will throw MathUnsupportedOperationException if bounds are passed to it.

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.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
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.
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.
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.
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.

Field Summary

Modifier and Type	Field and Description
private double[]	beta Coefficients of the Householder transforms vectors.
private final double	costRelativeTolerance Desired relative error in the sum of squares.
private double[]	diagR Diagonal elements of the R matrix in the Q.R. decomposition.
private final double	initialStepBoundFactor Positive input variable used in determining the initial step bound.
private double[]	jacNorm Norms of the columns of the jacobian matrix.
private double[]	lmDir Parameters evolution direction associated with lmPar.
private double	lmPar Levenberg-Marquardt parameter.
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 int[]	permutation Columns permutation array.
private final double	qrRankingThreshold Threshold for QR ranking.
private int	rank Rank of the jacobian matrix.
private int	solvedCols Number of solved point.
private static final double	TWO_EPS Twice the "epsilon machine".
private double[][]	weightedJacobian Weighted Jacobian.
private double[]	weightedResidual Weighted residuals.

Constructor Summary

Access	Constructor and Description
public	LevenbergMarquardtOptimizer() Build an optimizer for least squares problems with default values for all the tuning parameters (see the `other contructor`.
public	LevenbergMarquardtOptimizer(ConvergenceChecker<PointVectorValuePair> Convergence checker. checker) Constructor that allows the specification of a custom convergence checker.
public	LevenbergMarquardtOptimizer(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 `initialStepBoundFactor` itself. In most cases factor should lie in the interval `(0.1, 100.0)`. `100` is a generally recommended value. initialStepBoundFactor, ConvergenceChecker<PointVectorValuePair> Convergence checker. checker, double Desired relative error in the sum of squares. costRelativeTolerance, double Desired relative error in the approximate solution parameters. parRelativeTolerance, double Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian. orthoTolerance, 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. threshold) Constructor that allows the specification of a custom convergence checker, in addition to the standard ones.
public	LevenbergMarquardtOptimizer(double Desired relative error in the sum of squares. costRelativeTolerance, double Desired relative error in the approximate solution parameters. parRelativeTolerance, double Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian. orthoTolerance) Build an optimizer for least squares problems with default values for some of the tuning parameters (see the `other contructor`.
public	LevenbergMarquardtOptimizer(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 `initialStepBoundFactor` itself. In most cases factor should lie in the interval `(0.1, 100.0)`. `100` is a generally recommended value. initialStepBoundFactor, double Desired relative error in the sum of squares. costRelativeTolerance, double Desired relative error in the approximate solution parameters. parRelativeTolerance, double Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian. orthoTolerance, 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. threshold) The arguments control the behaviour of the default convergence checking procedure.

Method Summary

Modifier and Type	Method and Description
private void	checkParameters() Hides org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer.checkParameters. Check parameters consistency.
private void	determineLMDirection(double[] array containing qTy qy, double[] diagonal matrix diag, double[] diagonal elements associated with lmDir lmDiag, double[] work array work) Solve ax = b and dx = 0 in the least squares sense.
private void	determineLMParameter(double[] array containing qTy qy, double upper bound on the euclidean norm of diagR * lmDir delta, double[] diagonal matrix diag, double[] work array work1, double[] work array work2, double[] work array work3) Determine the Levenberg-Marquardt parameter.
protected PointVectorValuePair	doOptimize() Implements abstract org.apache.commons.math3.optim.BaseOptimizer.doOptimize. Performs the bulk of the optimization algorithm.
private void	qrDecomposition(RealMatrix Weighted Jacobian matrix at the current point. jacobian) 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) Compute the product Qt.y for some Q.R. decomposition.

Inherited from org.apache.commons.math3.optim.nonlinear.vector.jacobian.AbstractLeastSquaresOptimizer:: computeCost computeCovariances computeResiduals computeSigma computeWeightedJacobian getChiSquare getRMS getWeightSquareRoot optimize parseOptimizationData setCost

Field Detail

beta	back to summary
private double[] beta Coefficients of the Householder transforms vectors.

costRelativeTolerance	back to summary
private final double costRelativeTolerance Desired relative error in the sum of squares.

diagR	back to summary
private double[] diagR Diagonal elements of the R matrix in the Q.R. decomposition.

initialStepBoundFactor	back to summary
private final double initialStepBoundFactor Positive input variable used in determining the initial step bound.

jacNorm	back to summary
private double[] jacNorm Norms of the columns of the jacobian matrix.

lmDir	back to summary
private double[] lmDir Parameters evolution direction associated with lmPar.

lmPar	back to summary
private double lmPar Levenberg-Marquardt parameter.

orthoTolerance	back to summary
private final double orthoTolerance Desired max cosine on the orthogonality between the function vector and the columns of the jacobian.

parRelativeTolerance	back to summary
private final double parRelativeTolerance Desired relative error in the approximate solution parameters.

permutation	back to summary
private int[] permutation Columns permutation array.

qrRankingThreshold	back to summary
private final double qrRankingThreshold Threshold for QR ranking.

rank	back to summary
private int rank Rank of the jacobian matrix.

solvedCols	back to summary
private int solvedCols Number of solved point.

TWO_EPS	back to summary
private static final double TWO_EPS Twice the "epsilon machine".

weightedJacobian	back to summary
private double[][] weightedJacobian Weighted Jacobian.

weightedResidual	back to summary
private double[] weightedResidual Weighted residuals.

Constructor Detail

LevenbergMarquardtOptimizer back to summary

LevenbergMarquardtOptimizer	back to summary
public LevenbergMarquardtOptimizer() Build an optimizer for least squares problems with default values for all the tuning parameters (see the `other contructor`. 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`

public LevenbergMarquardtOptimizer()

Build an optimizer for least squares problems with default values for all the tuning parameters (see the other contructor. 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 back to summary

LevenbergMarquardtOptimizer	back to summary
public LevenbergMarquardtOptimizer(ConvergenceChecker<PointVectorValuePair> checker) Constructor that allows the specification of a custom convergence checker. Note that all the usual convergence checks will be disabled. 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` Parameters checker:ConvergenceChecker<PointVectorValuePair> Convergence checker.

public LevenbergMarquardtOptimizer(ConvergenceChecker<PointVectorValuePair> checker)

Constructor that allows the specification of a custom convergence checker. Note that all the usual convergence checks will be disabled. 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

Parameters

checker:ConvergenceChecker<PointVectorValuePair>: Convergence checker.

LevenbergMarquardtOptimizer back to summary

LevenbergMarquardtOptimizer	back to summary
public LevenbergMarquardtOptimizer(double initialStepBoundFactor, ConvergenceChecker<PointVectorValuePair> checker, double costRelativeTolerance, double parRelativeTolerance, double orthoTolerance, double threshold) Constructor that allows the specification of a custom convergence checker, in addition to the standard ones. Parameters initialStepBoundFactor: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 `initialStepBoundFactor` itself. In most cases factor should lie in the interval `(0.1, 100.0)`. `100` is a generally recommended value. checker:ConvergenceChecker<PointVectorValuePair> Convergence checker. costRelativeTolerance:double Desired relative error in the sum of squares. parRelativeTolerance:double Desired relative error in the approximate solution parameters. orthoTolerance:double Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian. threshold: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.

public LevenbergMarquardtOptimizer(double initialStepBoundFactor, ConvergenceChecker<PointVectorValuePair> checker, double costRelativeTolerance, double parRelativeTolerance, double orthoTolerance, double threshold)

Constructor that allows the specification of a custom convergence checker, in addition to the standard ones.

Parameters

initialStepBoundFactor: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 initialStepBoundFactor itself. In most cases factor should lie in the interval (0.1, 100.0). 100 is a generally recommended value.
checker:ConvergenceChecker<PointVectorValuePair>: Convergence checker.
costRelativeTolerance:double: Desired relative error in the sum of squares.
parRelativeTolerance:double: Desired relative error in the approximate solution parameters.
orthoTolerance:double: Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian.
threshold: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.

LevenbergMarquardtOptimizer back to summary

LevenbergMarquardtOptimizer	back to summary
public LevenbergMarquardtOptimizer(double costRelativeTolerance, double parRelativeTolerance, double orthoTolerance) Build an optimizer for least squares problems with default values for some of the tuning parameters (see the `other contructor`. The default values for the algorithm settings are: Initial step bound factor}: 100 QR ranking threshold}: `Precision#SAFE_MIN` Parameters costRelativeTolerance:double Desired relative error in the sum of squares. parRelativeTolerance:double Desired relative error in the approximate solution parameters. orthoTolerance:double Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian.

public LevenbergMarquardtOptimizer(double costRelativeTolerance, double parRelativeTolerance, double orthoTolerance)

Build an optimizer for least squares problems with default values for some of the tuning parameters (see the other contructor. The default values for the algorithm settings are:

Initial step bound factor}: 100
QR ranking threshold}: Precision#SAFE_MIN

Parameters

costRelativeTolerance:double: Desired relative error in the sum of squares.
parRelativeTolerance:double: Desired relative error in the approximate solution parameters.
orthoTolerance:double: Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian.

LevenbergMarquardtOptimizer back to summary

LevenbergMarquardtOptimizer	back to summary
public LevenbergMarquardtOptimizer(double initialStepBoundFactor, double costRelativeTolerance, double parRelativeTolerance, double orthoTolerance, double threshold) The arguments control the behaviour of the default convergence checking procedure. Additional criteria can defined through the setting of a `ConvergenceChecker`. Parameters initialStepBoundFactor: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 `initialStepBoundFactor` itself. In most cases factor should lie in the interval `(0.1, 100.0)`. `100` is a generally recommended value. costRelativeTolerance:double Desired relative error in the sum of squares. parRelativeTolerance:double Desired relative error in the approximate solution parameters. orthoTolerance:double Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian. threshold: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.

public LevenbergMarquardtOptimizer(double initialStepBoundFactor, double costRelativeTolerance, double parRelativeTolerance, double orthoTolerance, double threshold)

The arguments control the behaviour of the default convergence checking procedure. Additional criteria can defined through the setting of a ConvergenceChecker.

Parameters

initialStepBoundFactor: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 initialStepBoundFactor itself. In most cases factor should lie in the interval (0.1, 100.0). 100 is a generally recommended value.
costRelativeTolerance:double: Desired relative error in the sum of squares.
parRelativeTolerance:double: Desired relative error in the approximate solution parameters.
orthoTolerance:double: Desired max cosine on the orthogonality between the function vector and the columns of the Jacobian.
threshold: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.

Method Detail

checkParameters back to summary

checkParameters	back to summary
private void checkParameters() Hides org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer.checkParameters. Doc from org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer.checkParameters. Check parameters consistency. Exceptions MathUnsupportedOperationException: if bounds were passed to the `optimize` method.

private void checkParameters()

Hides org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer.checkParameters.

Doc from org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer.checkParameters.

Check parameters consistency.

Exceptions

MathUnsupportedOperationException:: if bounds were passed to the optimize method.

determineLMDirection	back to summary
private void determineLMDirection(double[] qy, double[] diag, double[] lmDiag, double[] work) Solve ax = b and dx = 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 work:double[] work array

determineLMDirection

back to summary

private void determineLMDirection(double[] qy, double[] diag, double[] lmDiag, double[] work)

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
work:double[]: work array

determineLMParameter	back to summary
private void determineLMParameter(double[] qy, double delta, double[] diag, double[] work1, double[] work2, double[] work3) Determine 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 work1:double[] work array work2:double[] work array work3:double[] work array

determineLMParameter

back to summary

private void determineLMParameter(double[] qy, double delta, double[] diag, double[] work1, double[] work2, double[] work3)

Determine 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
work1:double[]: work array
work2:double[]: work array
work3:double[]: work array

doOptimize	back to summary
protected PointVectorValuePair doOptimize() Implements abstract org.apache.commons.math3.optim.BaseOptimizer.doOptimize. Doc from org.apache.commons.math3.optim.BaseOptimizer.doOptimize. Performs the bulk of the optimization algorithm. Returns:PointVectorValuePair the point/value pair giving the optimal value of the objective function. Annotations @Override

doOptimize

back to summary

protected PointVectorValuePair doOptimize()

Implements abstract org.apache.commons.math3.optim.BaseOptimizer.doOptimize.

Doc from org.apache.commons.math3.optim.BaseOptimizer.doOptimize.

Performs the bulk of the optimization algorithm.

Returns:PointVectorValuePair

the point/value pair giving the optimal value of the objective function.

Annotations

@Override

qrDecomposition	back to summary
private void qrDecomposition(RealMatrix jacobian) 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 u_k unit vectors such that: H_k = I - 2u_k.u_k^t we use _k non-unit vectors such that: H_k = I - beta_kv_k.v_k^t where v_k = a_k - alpha_k e_k. The beta_k coefficients are provided upon exit as recomputing them from the v_k 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. Exceptions ConvergenceException: if the decomposition cannot be performed

qrDecomposition

back to summary

private void qrDecomposition(RealMatrix jacobian) 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 u_k unit vectors such that:

H_k = I - 2u_k.u_k^t

we use _k non-unit vectors such that:

H_k = I - beta_kv_k.v_k^t

where v_k = a_k - alpha_k e_k. The beta_k coefficients are provided upon exit as recomputing them from the v_k 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.

Exceptions

ConvergenceException:: if the decomposition cannot be performed

qTy	back to summary
private void qTy(double[] y) Compute the product Qt.y for some Q.R. decomposition. Parameters y:double[] vector to multiply (will be overwritten with the result)

qTy

back to summary

private void qTy(double[] y)

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

Parameters

y:double[]: vector to multiply (will be overwritten with the result)

public Class LevenbergMarquardtOptimizer

Field Summary

Constructor Summary

Method Summary

Field Detail

Constructor Detail

Method Detail