org.apache.commons.math.optimization.general

## Class LevenbergMarquardtOptimizer

- java.lang.Object
- - org.apache.commons.math.optimization.general.AbstractLeastSquaresOptimizer
  - - org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer

All Implemented Interfaces:
  DifferentiableMultivariateVectorialOptimizer

    public class LevenbergMarquardtOptimizer
    extends AbstractLeastSquaresOptimizer

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:
  2.0
Version:
  $Revision: 1073272 $ $Date: 2011-02-22 10:22:25 +0100 (mar. 22 févr. 2011) $

### Field Summary

Fields
`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.

### Fields inherited from class org.apache.commons.math.optimization.general.AbstractLeastSquaresOptimizer

`checker, cols, cost, DEFAULT_MAX_ITERATIONS, jacobian, objective, point, residuals, residualsWeights, rows, targetValues, wjacobian, wresiduals`

### Constructor Summary

Constructors
Constructor and Description
`LevenbergMarquardtOptimizer()` Build an optimizer for least squares problems.

### Method Summary

Methods
`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.

### Methods inherited from class org.apache.commons.math.optimization.general.AbstractLeastSquaresOptimizer

`getChiSquare, getConvergenceChecker, getCovariances, getEvaluations, getIterations, getJacobianEvaluations, getMaxEvaluations, getMaxIterations, getRMS, guessParametersErrors, incrementIterationsCounter, optimize, setConvergenceChecker, setMaxEvaluations, setMaxIterations, updateJacobian, updateResidualsAndCost`

### Methods inherited from class java.lang.Object

`clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait`

### Field Detail

#### solvedCols

    private int solvedCols

Number of solved point.

#### diagR

    private double[] diagR

Diagonal elements of the R matrix in the Q.R. decomposition.

#### jacNorm

    private double[] jacNorm

Norms of the columns of the jacobian matrix.

#### beta

    private double[] beta

Coefficients of the Householder transforms vectors.

#### permutation

    private int[] permutation

Columns permutation array.

#### rank

    private int rank

Rank of the jacobian matrix.

#### lmPar

    private double lmPar

Levenberg-Marquardt parameter.

#### lmDir

    private double[] lmDir

Parameters evolution direction associated with lmPar.

#### initialStepBoundFactor

    private double initialStepBoundFactor

Positive input variable used in determining the initial step bound.

#### costRelativeTolerance

    private double costRelativeTolerance

Desired relative error in the sum of squares.

#### parRelativeTolerance

    private double parRelativeTolerance

Desired relative error in the approximate solution parameters.

#### orthoTolerance

    private double orthoTolerance

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

#### qrRankingThreshold

    private double qrRankingThreshold

Threshold for QR ranking.

### Constructor Detail

#### LevenbergMarquardtOptimizer

    public LevenbergMarquardtOptimizer()

Build an optimizer for least squares problems.

The default values for the algorithm settings are:

- `vectorial convergence checker`: null
- `initial step bound factor`: 100.0
- `maximal iterations`: 1000
- `cost relative tolerance`: 1.0e-10
- `parameters relative tolerance`: 1.0e-10
- `orthogonality tolerance`: 1.0e-10
- `QR 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.

### Method Detail

#### setInitialStepBoundFactor

    public void setInitialStepBoundFactor(double initialStepBoundFactor)

Set the 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 nonzero, or else to initialStepBoundFactor itself. In most cases factor should lie in the interval (0.1, 100.0). 100.0 is a generally recommended value.

Parameters:
  `initialStepBoundFactor` - initial step bound factor

#### setCostRelativeTolerance

    public void setCostRelativeTolerance(double costRelativeTolerance)

Set the desired relative error in the sum of squares.

This setting is used only if the `vectorial convergence checker` is set to null.

Parameters:
  `costRelativeTolerance` - desired relative error in the sum of squares

#### setParRelativeTolerance

    public void setParRelativeTolerance(double parRelativeTolerance)

Set the desired relative error in the approximate solution parameters.

This setting is used only if the `vectorial convergence checker` is set to null.

Parameters:
  `parRelativeTolerance` - desired relative error in the approximate solution parameters

#### setOrthoTolerance

    public void setOrthoTolerance(double orthoTolerance)

Set the desired max cosine on the orthogonality.

This setting is always used, regardless of the `vectorial convergence checker` being null or non-null.

Parameters:
  `orthoTolerance` - desired max cosine on the orthogonality between the function vector and the columns of the jacobian

#### setQRRankingThreshold

    public void setQRRankingThreshold(double threshold)

Set the 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.

Parameters:
  `threshold` - threshold for QR ranking
Since:
  2.2

#### doOptimize

    protected VectorialPointValuePair doOptimize()
                                          throws FunctionEvaluationException,
                                                 OptimizationException,
                                                 IllegalArgumentException

Perform the bulk of optimization algorithm.

Specified by:
  `doOptimize` in class `AbstractLeastSquaresOptimizer`
Returns:
  the point/value pair giving the optimal value for objective function
Throws:
  `FunctionEvaluationException` - if the objective function throws one during the search
  `OptimizationException` - if the algorithm failed to converge
  `IllegalArgumentException` - if the start point dimension is wrong

#### determineLMParameter

    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` - array containing qTy
  `delta` - upper bound on the euclidean norm of diagR * lmDir
  `diag` - diagonal matrix
  `work1` - work array
  `work2` - work array
  `work3` - work array

#### determineLMDirection

    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` - array containing qTy
  `diag` - diagonal matrix
  `lmDiag` - diagonal elements associated with lmDir
  `work` - work array

#### qrDecomposition

    private void qrDecomposition()
                          throws OptimizationException

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.

Throws:
  `OptimizationException` - if the decomposition cannot be performed

#### qTy

    private void qTy(double[] y)

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

Parameters:
  `y` - vector to multiply (will be overwritten with the result)