/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math.optimization.general; import java.util.Arrays; import org.apache.commons.math.FunctionEvaluationException; import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.optimization.OptimizationException; import org.apache.commons.math.optimization.VectorialPointValuePair; import org.apache.commons.math.util.FastMath; import org.apache.commons.math.util.MathUtils; /** * 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: *
* 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:
*
|
The default values for the algorithm settings are: *
These default values may be overridden after construction. If the {@link * #setConvergenceChecker vectorial convergence checker} is set to a non-null value, it * will be used instead of the {@link #setCostRelativeTolerance cost relative tolerance} * and {@link #setParRelativeTolerance parameters relative tolerance} settings. */ public LevenbergMarquardtOptimizer() { // set up the superclass with a default max cost evaluations setting setMaxIterations(1000); // default values for the tuning parameters setConvergenceChecker(null); setInitialStepBoundFactor(100.0); setCostRelativeTolerance(1.0e-10); setParRelativeTolerance(1.0e-10); setOrthoTolerance(1.0e-10); setQRRankingThreshold(MathUtils.SAFE_MIN); } /** * 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. * * @param initialStepBoundFactor initial step bound factor */ public void setInitialStepBoundFactor(double initialStepBoundFactor) { this.initialStepBoundFactor = initialStepBoundFactor; } /** * Set the desired relative error in the sum of squares. *
This setting is used only if the {@link #setConvergenceChecker vectorial * convergence checker} is set to null.
* @param costRelativeTolerance desired relative error in the sum of squares */ public void setCostRelativeTolerance(double costRelativeTolerance) { this.costRelativeTolerance = costRelativeTolerance; } /** * Set the desired relative error in the approximate solution parameters. *This setting is used only if the {@link #setConvergenceChecker vectorial * convergence checker} is set to null.
* @param parRelativeTolerance desired relative error * in the approximate solution parameters */ public void setParRelativeTolerance(double parRelativeTolerance) { this.parRelativeTolerance = parRelativeTolerance; } /** * Set the desired max cosine on the orthogonality. *This setting is always used, regardless of the {@link #setConvergenceChecker * vectorial convergence checker} being null or non-null.
* @param orthoTolerance desired max cosine on the orthogonality * between the function vector and the columns of the jacobian */ public void setOrthoTolerance(double orthoTolerance) { this.orthoTolerance = orthoTolerance; } /** * 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. *
* @param threshold threshold for QR ranking * @since 2.2 */ public void setQRRankingThreshold(final double threshold) { this.qrRankingThreshold = threshold; } /** {@inheritDoc} */ @Override protected VectorialPointValuePair doOptimize() throws FunctionEvaluationException, OptimizationException, IllegalArgumentException { // arrays shared with the other private methods solvedCols = Math.min(rows, cols); diagR = new double[cols]; jacNorm = new double[cols]; beta = new double[cols]; permutation = new int[cols]; lmDir = new double[cols]; // local point double delta = 0; double xNorm = 0; double[] diag = new double[cols]; double[] oldX = new double[cols]; double[] oldRes = new double[rows]; double[] oldObj = new double[rows]; double[] qtf = new double[rows]; double[] work1 = new double[cols]; double[] work2 = new double[cols]; double[] work3 = new double[cols]; // evaluate the function at the starting point and calculate its norm updateResidualsAndCost(); // outer loop lmPar = 0; boolean firstIteration = true; VectorialPointValuePair current = new VectorialPointValuePair(point, objective); while (true) { for (int i=0;iThis 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.
* * @param qy array containing qTy * @param delta upper bound on the euclidean norm of diagR * lmDir * @param diag diagonal matrix * @param work1 work array * @param work2 work array * @param work3 work array */ private void determineLMParameter(double[] qy, double delta, double[] diag, double[] work1, double[] work2, double[] work3) { // compute and store in x the gauss-newton direction, if the // jacobian is rank-deficient, obtain a least squares solution for (int j = 0; j < rank; ++j) { lmDir[permutation[j]] = qy[j]; } for (int j = rank; j < cols; ++j) { lmDir[permutation[j]] = 0; } for (int k = rank - 1; k >= 0; --k) { int pk = permutation[k]; double ypk = lmDir[pk] / diagR[pk]; for (int i = 0; i < k; ++i) { lmDir[permutation[i]] -= ypk * wjacobian[i][pk]; } lmDir[pk] = ypk; } // evaluate the function at the origin, and test // for acceptance of the Gauss-Newton direction double dxNorm = 0; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double s = diag[pj] * lmDir[pj]; work1[pj] = s; dxNorm += s * s; } dxNorm = FastMath.sqrt(dxNorm); double fp = dxNorm - delta; if (fp <= 0.1 * delta) { lmPar = 0; return; } // if the jacobian is not rank deficient, the Newton step provides // a lower bound, parl, for the zero of the function, // otherwise set this bound to zero double sum2; double parl = 0; if (rank == solvedCols) { for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; work1[pj] *= diag[pj] / dxNorm; } sum2 = 0; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double sum = 0; for (int i = 0; i < j; ++i) { sum += wjacobian[i][pj] * work1[permutation[i]]; } double s = (work1[pj] - sum) / diagR[pj]; work1[pj] = s; sum2 += s * s; } parl = fp / (delta * sum2); } // calculate an upper bound, paru, for the zero of the function sum2 = 0; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double sum = 0; for (int i = 0; i <= j; ++i) { sum += wjacobian[i][pj] * qy[i]; } sum /= diag[pj]; sum2 += sum * sum; } double gNorm = FastMath.sqrt(sum2); double paru = gNorm / delta; if (paru == 0) { // 2.2251e-308 is the smallest positive real for IEE754 paru = 2.2251e-308 / FastMath.min(delta, 0.1); } // if the input par lies outside of the interval (parl,paru), // set par to the closer endpoint lmPar = FastMath.min(paru, FastMath.max(lmPar, parl)); if (lmPar == 0) { lmPar = gNorm / dxNorm; } for (int countdown = 10; countdown >= 0; --countdown) { // evaluate the function at the current value of lmPar if (lmPar == 0) { lmPar = FastMath.max(2.2251e-308, 0.001 * paru); } double sPar = FastMath.sqrt(lmPar); for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; work1[pj] = sPar * diag[pj]; } determineLMDirection(qy, work1, work2, work3); dxNorm = 0; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double s = diag[pj] * lmDir[pj]; work3[pj] = s; dxNorm += s * s; } dxNorm = FastMath.sqrt(dxNorm); double previousFP = fp; fp = dxNorm - delta; // if the function is small enough, accept the current value // of lmPar, also test for the exceptional cases where parl is zero if ((FastMath.abs(fp) <= 0.1 * delta) || ((parl == 0) && (fp <= previousFP) && (previousFP < 0))) { return; } // compute the Newton correction for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; work1[pj] = work3[pj] * diag[pj] / dxNorm; } for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; work1[pj] /= work2[j]; double tmp = work1[pj]; for (int i = j + 1; i < solvedCols; ++i) { work1[permutation[i]] -= wjacobian[i][pj] * tmp; } } sum2 = 0; for (int j = 0; j < solvedCols; ++j) { double s = work1[permutation[j]]; sum2 += s * s; } double correction = fp / (delta * sum2); // depending on the sign of the function, update parl or paru. if (fp > 0) { parl = FastMath.max(parl, lmPar); } else if (fp < 0) { paru = FastMath.min(paru, lmPar); } // compute an improved estimate for lmPar lmPar = FastMath.max(parl, lmPar + correction); } } /** * 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:
*Luc Maisonobe did the Java translation.
* * @param qy array containing qTy * @param diag diagonal matrix * @param lmDiag diagonal elements associated with lmDir * @param work work array */ private void determineLMDirection(double[] qy, double[] diag, double[] lmDiag, double[] work) { // copy R and Qty to preserve input and initialize s // in particular, save the diagonal elements of R in lmDir for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; for (int i = j + 1; i < solvedCols; ++i) { wjacobian[i][pj] = wjacobian[j][permutation[i]]; } lmDir[j] = diagR[pj]; work[j] = qy[j]; } // eliminate the diagonal matrix d using a Givens rotation for (int j = 0; j < solvedCols; ++j) { // prepare the row of d to be eliminated, locating the // diagonal element using p from the Q.R. factorization int pj = permutation[j]; double dpj = diag[pj]; if (dpj != 0) { Arrays.fill(lmDiag, j + 1, lmDiag.length, 0); } lmDiag[j] = dpj; // the transformations to eliminate the row of d // modify only a single element of Qty // beyond the first n, which is initially zero. double qtbpj = 0; for (int k = j; k < solvedCols; ++k) { int pk = permutation[k]; // determine a Givens rotation which eliminates the // appropriate element in the current row of d if (lmDiag[k] != 0) { final double sin; final double cos; double rkk = wjacobian[k][pk]; if (FastMath.abs(rkk) < FastMath.abs(lmDiag[k])) { final double cotan = rkk / lmDiag[k]; sin = 1.0 / FastMath.sqrt(1.0 + cotan * cotan); cos = sin * cotan; } else { final double tan = lmDiag[k] / rkk; cos = 1.0 / FastMath.sqrt(1.0 + tan * tan); sin = cos * tan; } // compute the modified diagonal element of R and // the modified element of (Qty,0) wjacobian[k][pk] = cos * rkk + sin * lmDiag[k]; final double temp = cos * work[k] + sin * qtbpj; qtbpj = -sin * work[k] + cos * qtbpj; work[k] = temp; // accumulate the tranformation in the row of s for (int i = k + 1; i < solvedCols; ++i) { double rik = wjacobian[i][pk]; final double temp2 = cos * rik + sin * lmDiag[i]; lmDiag[i] = -sin * rik + cos * lmDiag[i]; wjacobian[i][pk] = temp2; } } } // store the diagonal element of s and restore // the corresponding diagonal element of R lmDiag[j] = wjacobian[j][permutation[j]]; wjacobian[j][permutation[j]] = lmDir[j]; } // solve the triangular system for z, if the system is // singular, then obtain a least squares solution int nSing = solvedCols; for (int j = 0; j < solvedCols; ++j) { if ((lmDiag[j] == 0) && (nSing == solvedCols)) { nSing = j; } if (nSing < solvedCols) { work[j] = 0; } } if (nSing > 0) { for (int j = nSing - 1; j >= 0; --j) { int pj = permutation[j]; double sum = 0; for (int i = j + 1; i < nSing; ++i) { sum += wjacobian[i][pj] * work[i]; } work[j] = (work[j] - sum) / lmDiag[j]; } } // permute the components of z back to components of lmDir for (int j = 0; j < lmDir.length; ++j) { lmDir[permutation[j]] = work[j]; } } /** * 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.
* @exception OptimizationException if the decomposition cannot be performed */ private void qrDecomposition() throws OptimizationException { // initializations for (int k = 0; k < cols; ++k) { permutation[k] = k; double norm2 = 0; for (int i = 0; i < wjacobian.length; ++i) { double akk = wjacobian[i][k]; norm2 += akk * akk; } jacNorm[k] = FastMath.sqrt(norm2); } // transform the matrix column after column for (int k = 0; k < cols; ++k) { // select the column with the greatest norm on active components int nextColumn = -1; double ak2 = Double.NEGATIVE_INFINITY; for (int i = k; i < cols; ++i) { double norm2 = 0; for (int j = k; j < wjacobian.length; ++j) { double aki = wjacobian[j][permutation[i]]; norm2 += aki * aki; } if (Double.isInfinite(norm2) || Double.isNaN(norm2)) { throw new OptimizationException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN, rows, cols); } if (norm2 > ak2) { nextColumn = i; ak2 = norm2; } } if (ak2 <= qrRankingThreshold) { rank = k; return; } int pk = permutation[nextColumn]; permutation[nextColumn] = permutation[k]; permutation[k] = pk; // choose alpha such that Hk.u = alpha ek double akk = wjacobian[k][pk]; double alpha = (akk > 0) ? -FastMath.sqrt(ak2) : FastMath.sqrt(ak2); double betak = 1.0 / (ak2 - akk * alpha); beta[pk] = betak; // transform the current column diagR[pk] = alpha; wjacobian[k][pk] -= alpha; // transform the remaining columns for (int dk = cols - 1 - k; dk > 0; --dk) { double gamma = 0; for (int j = k; j < wjacobian.length; ++j) { gamma += wjacobian[j][pk] * wjacobian[j][permutation[k + dk]]; } gamma *= betak; for (int j = k; j < wjacobian.length; ++j) { wjacobian[j][permutation[k + dk]] -= gamma * wjacobian[j][pk]; } } } rank = solvedCols; } /** * Compute the product Qt.y for some Q.R. decomposition. * * @param y vector to multiply (will be overwritten with the result) */ private void qTy(double[] y) { for (int k = 0; k < cols; ++k) { int pk = permutation[k]; double gamma = 0; for (int i = k; i < rows; ++i) { gamma += wjacobian[i][pk] * y[i]; } gamma *= beta[pk]; for (int i = k; i < rows; ++i) { y[i] -= gamma * wjacobian[i][pk]; } } } }