public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer
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:
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:
|
Modifier and Type | Field and Description |
---|---|
private double[] |
beta
Coefficients of the Householder transforms vectors.
|
private double |
costRelativeTolerance
Desired relative error in the sum of squares.
|
private double[] |
diagR
Diagonal elements of the R matrix in the Q.R.
|
private 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 double |
orthoTolerance
Desired max cosine on the orthogonality between the function vector
and the columns of the jacobian.
|
private double |
parRelativeTolerance
Desired relative error in the approximate solution parameters.
|
private int[] |
permutation
Columns permutation array.
|
private double |
qrRankingThreshold
Threshold for QR ranking.
|
private int |
rank
Rank of the jacobian matrix.
|
private int |
solvedCols
Number of solved point.
|
checker, cols, cost, DEFAULT_MAX_ITERATIONS, jacobian, objective, point, residuals, residualsWeights, rows, targetValues, wjacobian, wresiduals
Constructor and Description |
---|
LevenbergMarquardtOptimizer()
Build an optimizer for least squares problems.
|
Modifier and Type | Method and Description |
---|---|
private void |
determineLMDirection(double[] qy,
double[] diag,
double[] lmDiag,
double[] work)
Solve a*x = b and d*x = 0 in the least squares sense.
|
private void |
determineLMParameter(double[] qy,
double delta,
double[] diag,
double[] work1,
double[] work2,
double[] work3)
Determine the Levenberg-Marquardt parameter.
|
protected VectorialPointValuePair |
doOptimize()
Perform the bulk of optimization algorithm.
|
private void |
qrDecomposition()
Decompose a matrix A as A.P = Q.R using Householder transforms.
|
private void |
qTy(double[] y)
Compute the product Qt.y for some Q.R.
|
void |
setCostRelativeTolerance(double costRelativeTolerance)
Set the desired relative error in the sum of squares.
|
void |
setInitialStepBoundFactor(double initialStepBoundFactor)
Set the positive input variable used in determining the initial step bound.
|
void |
setOrthoTolerance(double orthoTolerance)
Set the desired max cosine on the orthogonality.
|
void |
setParRelativeTolerance(double parRelativeTolerance)
Set the desired relative error in the approximate solution parameters.
|
void |
setQRRankingThreshold(double threshold)
Set the desired threshold for QR ranking.
|
getChiSquare, getConvergenceChecker, getCovariances, getEvaluations, getIterations, getJacobianEvaluations, getMaxEvaluations, getMaxIterations, getRMS, guessParametersErrors, incrementIterationsCounter, optimize, setConvergenceChecker, setMaxEvaluations, setMaxIterations, updateJacobian, updateResidualsAndCost
private int solvedCols
private double[] diagR
private double[] jacNorm
private double[] beta
private int[] permutation
private int rank
private double lmPar
private double[] lmDir
private double initialStepBoundFactor
private double costRelativeTolerance
private double parRelativeTolerance
private double orthoTolerance
private double qrRankingThreshold
public LevenbergMarquardtOptimizer()
The default values for the algorithm settings are:
vectorial convergence checker
: nullinitial step bound factor
: 100.0maximal iterations
: 1000cost relative tolerance
: 1.0e-10parameters relative tolerance
: 1.0e-10orthogonality tolerance
: 1.0e-10QR ranking threshold
: MathUtils.SAFE_MIN
These default values may be overridden after construction. If the vectorial convergence checker
is set to a non-null value, it
will be used instead of the cost relative tolerance
and parameters relative tolerance
settings.
public void setInitialStepBoundFactor(double initialStepBoundFactor)
initialStepBoundFactor
- initial step bound factorpublic void setCostRelativeTolerance(double costRelativeTolerance)
This setting is used only if the vectorial
convergence checker
is set to null.
costRelativeTolerance
- desired relative error in the sum of squarespublic void setParRelativeTolerance(double parRelativeTolerance)
This setting is used only if the vectorial
convergence checker
is set to null.
parRelativeTolerance
- desired relative error
in the approximate solution parameterspublic void setOrthoTolerance(double orthoTolerance)
This setting is always used, regardless of the vectorial convergence checker
being null or non-null.
orthoTolerance
- desired max cosine on the orthogonality
between the function vector and the columns of the jacobianpublic void setQRRankingThreshold(double threshold)
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
- threshold for QR rankingprotected VectorialPointValuePair doOptimize() throws FunctionEvaluationException, OptimizationException, java.lang.IllegalArgumentException
doOptimize
in class AbstractLeastSquaresOptimizer
FunctionEvaluationException
- if the objective function throws one during
the searchOptimizationException
- if the algorithm failed to convergejava.lang.IllegalArgumentException
- if the start point dimension is wrongprivate void determineLMParameter(double[] qy, double delta, double[] diag, double[] work1, double[] work2, double[] work3)
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:
Luc Maisonobe did the Java translation.
qy
- array containing qTydelta
- upper bound on the euclidean norm of diagR * lmDirdiag
- diagonal matrixwork1
- work arraywork2
- work arraywork3
- work arrayprivate void determineLMDirection(double[] qy, double[] diag, double[] lmDiag, double[] work)
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:
Luc Maisonobe did the Java translation.
qy
- array containing qTydiag
- diagonal matrixlmDiag
- diagonal elements associated with lmDirwork
- work arrayprivate void qrDecomposition() throws OptimizationException
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.uktwe use k non-unit vectors such that:
Hk = I - betakvk.vktwhere 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.
OptimizationException
- if the decomposition cannot be performedprivate void qTy(double[] y)
y
- vector to multiply (will be overwritten with the result)Copyright (c) 2003-2013 Apache Software Foundation