LevenbergMarquardtEstimator.java

Index Score
org.apache.commons.math.estimation
Commons Math

View: Reasons, Metrics, Source Code

These are the metrics that contribute to the Enerjy Score for this file, ranked by impact. So the metrics listed at the top influence the score to a greater extent that the metrics listed at the bottom.

MetricDescription
LOOPSNumber of loops
EXEC_COMMENTSComments in executable code
LINE_COMMENTNumber of line comments
COMPARISONSNumber of comparison operators
LOGICAL_LINESNumber of statements
SIZESize of the file in bytes
COMMENTSComment lines
DOC_COMMENTNumber of javadoc comment lines
OPERANDSNumber of operands
JAVA0078JAVA0078 Floating point values compared with ==
JAVA0076JAVA0076 Use of magic number
PROGRAM_LENGTHHalstead program length
OPERATORSNumber of operators
LINESNumber of lines in the source file
CYCLOMATICCyclomatic complexity
BLOCKSNumber of blocks
ELOCEffective lines of code
LOCLines of code
JAVA0034JAVA0034 Missing braces in if statement
DECL_COMMENTSComments in declarations
PROGRAM_VOCABHalstead program vocabulary
JAVA0177JAVA0177 Variable declaration missing initializer
UNIQUE_OPERANDSNumber of unique operands
UNIQUE_OPERATORSNumber of unique operators
JAVA0117JAVA0117 Missing javadoc: method 'method'
JAVA0119JAVA0119 Control variable changed within body of for loop
NEST_DEPTHMaximum nesting depth
PROGRAM_VOLUMEHalstead program volume
JAVA0126JAVA0126 Method declares unchecked exception in throws
JAVA0110JAVA0110 Incorrect javadoc: no @return tag
JAVA0123JAVA0123 Use all three components of for loop
JAVA0100JAVA0100 Class contains N non-final fields (maximum: M)
JAVA0108JAVA0108 Incorrect javadoc: no @param tag for 'parameter'
FUNCTIONSNumber of function declarations
JAVA0145JAVA0145 Tab character used in source file
/* * 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.estimation; import java.io.Serializable; import java.util.Arrays; /** * This class solves a least squares problem. * * <p>This implementation <em>should</em> work even for over-determined systems * (i.e. systems having more variables than equations). Over-determined systems * are solved by ignoring the variables 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. This feature has undergone only basic testing * for now and should still be considered experimental.</p> * * <p>The resolution engine is a simple translation of the MINPACK <a * href="http://www.netlib.org/minpack/lmder.f">lmder</a> routine with minor * changes. The changes include the over-determined resolution and the Q.R. * decomposition which has been rewritten following the algorithm described in the * P. Lascaux and R. Theodor book <i>Analyse num&eacute;rique matricielle * appliqu&eacute;e &agrave; l'art de l'ing&eacute;nieur</i>, Masson 1986. The * redistribution policy for MINPACK is available <a * href="http://www.netlib.org/minpack/disclaimer">here</a>, for convenience, it * is reproduced below.</p> * * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> * <tr><td> * Minpack Copyright Notice (1999) University of Chicago. * All rights reserved * </td></tr> * <tr><td> * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * <ol> * <li>Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer.</li> * <li>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.</li> * <li>The end-user documentation included with the redistribution, if any, * must include the following acknowledgment: * <code>This product includes software developed by the University of * Chicago, as Operator of Argonne National Laboratory.</code> * Alternately, this acknowledgment may appear in the software itself, * if and wherever such third-party acknowledgments normally appear.</li> * <li><strong>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.</strong></li> * <li><strong>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.</strong></li> * <ol></td></tr> * </table> * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran) * @author Burton S. Garbow (original fortran) * @author Kenneth E. Hillstrom (original fortran) * @author Jorge J. More (original fortran) * @version $Revision: 640204 $ $Date: 2008-03-23 09:36:03 -0400 (Sun, 23 Mar 2008) $ * @since 1.2 * */ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Serializable { /** * Build an estimator for least squares problems. * <p>The default values for the algorithm settings are: * <ul> * <li>{@link #setInitialStepBoundFactor initial step bound factor}: 100.0</li> * <li>{@link #setMaxCostEval maximal cost evaluations}: 1000</li> * <li>{@link #setCostRelativeTolerance cost relative tolerance}: 1.0e-10</li> * <li>{@link #setParRelativeTolerance parameters relative tolerance}: 1.0e-10</li> * <li>{@link #setOrthoTolerance orthogonality tolerance}: 1.0e-10</li> * </ul> * </p> */ public LevenbergMarquardtEstimator() { // set up the superclass with a default max cost evaluations setting setMaxCostEval(1000); // default values for the tuning parameters setInitialStepBoundFactor(100.0); setCostRelativeTolerance(1.0e-10); setParRelativeTolerance(1.0e-10); setOrthoTolerance(1.0e-10); } /** * 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 * @see #estimate */ public void setInitialStepBoundFactor(double initialStepBoundFactor) { this.initialStepBoundFactor = initialStepBoundFactor; } /** * Set the desired relative error in the sum of squares. * * @param costRelativeTolerance desired relative error in the sum of squares * @see #estimate */ public void setCostRelativeTolerance(double costRelativeTolerance) { this.costRelativeTolerance = costRelativeTolerance; } /** * Set the desired relative error in the approximate solution parameters. * * @param parRelativeTolerance desired relative error * in the approximate solution parameters * @see #estimate */ public void setParRelativeTolerance(double parRelativeTolerance) { this.parRelativeTolerance = parRelativeTolerance; } /** * Set the desired max cosine on the orthogonality. * * @param orthoTolerance desired max cosine on the orthogonality * between the function vector and the columns of the jacobian * @see #estimate */ public void setOrthoTolerance(double orthoTolerance) { this.orthoTolerance = orthoTolerance; } /** * Solve an estimation problem using the Levenberg-Marquardt algorithm. * <p>The algorithm used is a modified Levenberg-Marquardt one, based * on the MINPACK <a href="http://www.netlib.org/minpack/lmder.f">lmder</a> * routine. The algorithm settings must have been set up before this method * is called with the {@link #setInitialStepBoundFactor}, * {@link #setMaxCostEval}, {@link #setCostRelativeTolerance}, * {@link #setParRelativeTolerance} and {@link #setOrthoTolerance} methods. * If these methods have not been called, the default values set up by the * {@link #LevenbergMarquardtEstimator() constructor} will be used.</p> * <p>The authors of the original fortran function are:</p> * <ul> * <li>Argonne National Laboratory. MINPACK project. March 1980</li> * <li>Burton S. Garbow</li> * <li>Kenneth E. Hillstrom</li> * <li>Jorge J. More</li> * </ul> * <p>Luc Maisonobe did the Java translation.</p> * * @param problem estimation problem to solve * @exception EstimationException if convergence cannot be * reached with the specified algorithm settings or if there are more variables * than equations * @see #setInitialStepBoundFactor * @see #setCostRelativeTolerance * @see #setParRelativeTolerance * @see #setOrthoTolerance */ public void estimate(EstimationProblem problem) throws EstimationException { initializeEstimate(problem); // 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 variables double delta = 0, xNorm = 0; double[] diag = new double[cols]; double[] oldX = new double[cols]; double[] oldRes = 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; while (true) { // compute the Q.R. decomposition of the jacobian matrix updateJacobian(); qrDecomposition(); // compute Qt.res qTy(residuals); // now we don't need Q anymore, // so let jacobian contain the R matrix with its diagonal elements for (int k = 0; k < solvedCols; ++k) { int pk = permutation[k]; jacobian[k * cols + pk] = diagR[pk]; } if (firstIteration) { // scale the variables according to the norms of the columns // of the initial jacobian xNorm = 0; for (int k = 0; k < cols; ++k) { double dk = jacNorm[k]; if (dk == 0) { dk = 1.0; } double xk = dk * parameters[k].getEstimate(); xNorm += xk * xk; diag[k] = dk; } xNorm = Math.sqrt(xNorm); // initialize the step bound delta delta = (xNorm == 0) ? initialStepBoundFactor : (initialStepBoundFactor * xNorm); } // check orthogonality between function vector and jacobian columns double maxCosine = 0; if (cost != 0) { for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double s = jacNorm[pj]; if (s != 0) { double sum = 0; for (int i = 0, index = pj; i <= j; ++i, index += cols) { sum += jacobian[index] * residuals[i]; } maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost)); } } } if (maxCosine <= orthoTolerance) { return; } // rescale if necessary for (int j = 0; j < cols; ++j) { diag[j] = Math.max(diag[j], jacNorm[j]); } // inner loop for (double ratio = 0; ratio < 1.0e-4;) { // save the state for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; oldX[pj] = parameters[pj].getEstimate(); } double previousCost = cost; double[] tmpVec = residuals; residuals = oldRes; oldRes = tmpVec; // determine the Levenberg-Marquardt parameter determineLMParameter(oldRes, delta, diag, work1, work2, work3); // compute the new point and the norm of the evolution direction double lmNorm = 0; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; lmDir[pj] = -lmDir[pj]; parameters[pj].setEstimate(oldX[pj] + lmDir[pj]); double s = diag[pj] * lmDir[pj]; lmNorm += s * s; } lmNorm = Math.sqrt(lmNorm); // on the first iteration, adjust the initial step bound. if (firstIteration) { delta = Math.min(delta, lmNorm); } // evaluate the function at x + p and calculate its norm updateResidualsAndCost(); // compute the scaled actual reduction double actRed = -1.0; if (0.1 * cost < previousCost) { double r = cost / previousCost; actRed = 1.0 - r * r; } // compute the scaled predicted reduction // and the scaled directional derivative for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double dirJ = lmDir[pj]; work1[j] = 0; for (int i = 0, index = pj; i <= j; ++i, index += cols) { work1[i] += jacobian[index] * dirJ; } } double coeff1 = 0; for (int j = 0; j < solvedCols; ++j) { coeff1 += work1[j] * work1[j]; } double pc2 = previousCost * previousCost; coeff1 = coeff1 / pc2; double coeff2 = lmPar * lmNorm * lmNorm / pc2; double preRed = coeff1 + 2 * coeff2; double dirDer = -(coeff1 + coeff2); // ratio of the actual to the predicted reduction ratio = (preRed == 0) ? 0 : (actRed / preRed); // update the step bound if (ratio <= 0.25) { double tmp = (actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5; if ((0.1 * cost >= previousCost) || (tmp < 0.1)) { tmp = 0.1; } delta = tmp * Math.min(delta, 10.0 * lmNorm); lmPar /= tmp; } else if ((lmPar == 0) || (ratio >= 0.75)) { delta = 2 * lmNorm; lmPar *= 0.5; } // test for successful iteration. if (ratio >= 1.0e-4) { // successful iteration, update the norm firstIteration = false; xNorm = 0; for (int k = 0; k < cols; ++k) { double xK = diag[k] * parameters[k].getEstimate(); xNorm += xK * xK; } xNorm = Math.sqrt(xNorm); } else { // failed iteration, reset the previous values cost = previousCost; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; parameters[pj].setEstimate(oldX[pj]); } tmpVec = residuals; residuals = oldRes; oldRes = tmpVec; } // tests for convergence. if (((Math.abs(actRed) <= costRelativeTolerance) && (preRed <= costRelativeTolerance) && (ratio <= 2.0)) || (delta <= parRelativeTolerance * xNorm)) { return; } // tests for termination and stringent tolerances // (2.2204e-16 is the machine epsilon for IEEE754) if ((Math.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) { throw new EstimationException("cost relative tolerance is too small ({0})," + " no further reduction in the" + " sum of squares is possible", new Object[] { new Double(costRelativeTolerance) }); } else if (delta <= 2.2204e-16 * xNorm) { throw new EstimationException("parameters relative tolerance is too small" + " ({0}), no further improvement in" + " the approximate solution is possible", new Object[] { new Double(parRelativeTolerance) }); } else if (maxCosine <= 2.2204e-16) { throw new EstimationException("orthogonality tolerance is too small ({0})," + " solution is orthogonal to the jacobian", new Object[] { new Double(orthoTolerance) }); } } } } /** * Determine the Levenberg-Marquardt parameter. * <p>This implementation is a translation in Java of the MINPACK * <a href="http://www.netlib.org/minpack/lmpar.f">lmpar</a> * routine.</p> * <p>This method sets the lmPar and lmDir attributes.</p> * <p>The authors of the original fortran function are:</p> * <ul> * <li>Argonne National Laboratory. MINPACK project. March 1980</li> * <li>Burton S. Garbow</li> * <li>Kenneth E. Hillstrom</li> * <li>Jorge J. More</li> * </ul> * <p>Luc Maisonobe did the Java translation.</p> * * @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, index = pk; i < k; ++i, index += cols) { lmDir[permutation[i]] -= ypk * jacobian[index]; } 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 = Math.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, 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, index = pj; i < j; ++i, index += cols) { sum += jacobian[index] * 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, index = pj; i <= j; ++i, index += cols) { sum += jacobian[index] * qy[i]; } sum /= diag[pj]; sum2 += sum * sum; } double gNorm = Math.sqrt(sum2); double paru = gNorm / delta; if (paru == 0) { // 2.2251e-308 is the smallest positive real for IEE754 paru = 2.2251e-308 / Math.min(delta, 0.1); } // if the input par lies outside of the interval (parl,paru), // set par to the closer endpoint lmPar = Math.min(paru, Math.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 = Math.max(2.2251e-308, 0.001 * paru); } double sPar = Math.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 = Math.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 ((Math.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]] -= jacobian[i * cols + 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 = Math.max(parl, lmPar); } else if (fp < 0) { paru = Math.min(paru, lmPar); } // compute an improved estimate for lmPar lmPar = Math.max(parl, lmPar + correction); } } /** * Solve a*x = b and d*x = 0 in the least squares sense. * <p>This implementation is a translation in Java of the MINPACK * <a href="http://www.netlib.org/minpack/qrsolv.f">qrsolv</a> * routine.</p> * <p>This method sets the lmDir and lmDiag attributes.</p> * <p>The authors of the original fortran function are:</p> * <ul> * <li>Argonne National Laboratory. MINPACK project. March 1980</li> * <li>Burton S. Garbow</li> * <li>Kenneth E. Hillstrom</li> * <li>Jorge J. More</li> * </ul> * <p>Luc Maisonobe did the Java translation.</p> * * @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) { jacobian[i * cols + pj] = jacobian[j * cols + 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) { double sin, cos; double rkk = jacobian[k * cols + pk]; if (Math.abs(rkk) < Math.abs(lmDiag[k])) { double cotan = rkk / lmDiag[k]; sin = 1.0 / Math.sqrt(1.0 + cotan * cotan); cos = sin * cotan; } else { double tan = lmDiag[k] / rkk; cos = 1.0 / Math.sqrt(1.0 + tan * tan); sin = cos * tan; } // compute the modified diagonal element of R and // the modified element of (Qty,0) jacobian[k * cols + pk] = cos * rkk + sin * lmDiag[k]; 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 = jacobian[i * cols + pk]; temp = cos * rik + sin * lmDiag[i]; lmDiag[i] = -sin * rik + cos * lmDiag[i]; jacobian[i * cols + pk] = temp; } } } // store the diagonal element of s and restore // the corresponding diagonal element of R int index = j * cols + permutation[j]; lmDiag[j] = jacobian[index]; jacobian[index] = 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 += jacobian[i * cols + 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. * <p>As suggested in the P. Lascaux and R. Theodor book * <i>Analyse num&eacute;rique matricielle appliqu&eacute;e &agrave; * l'art de l'ing&eacute;nieur</i> (Masson, 1986), instead of representing * the Householder transforms with u<sub>k</sub> unit vectors such that: * <pre> * H<sub>k</sub> = I - 2u<sub>k</sub>.u<sub>k</sub><sup>t</sup> * </pre> * we use <sub>k</sub> non-unit vectors such that: * <pre> * H<sub>k</sub> = I - beta<sub>k</sub>v<sub>k</sub>.v<sub>k</sub><sup>t</sup> * </pre> * where v<sub>k</sub> = a<sub>k</sub> - alpha<sub>k</sub> e<sub>k</sub>. * The beta<sub>k</sub> coefficients are provided upon exit as recomputing * them from the v<sub>k</sub> vectors would be costly.</p> * <p>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.</p> * @exception EstimationException if the decomposition cannot be performed */ private void qrDecomposition() throws EstimationException { // initializations for (int k = 0; k < cols; ++k) { permutation[k] = k; double norm2 = 0; for (int index = k; index < jacobian.length; index += cols) { double akk = jacobian[index]; norm2 += akk * akk; } jacNorm[k] = Math.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; int iDiag = k * cols + permutation[i]; for (int index = iDiag; index < jacobian.length; index += cols) { double aki = jacobian[index]; norm2 += aki * aki; } if (Double.isInfinite(norm2) || Double.isNaN(norm2)) { throw new EstimationException("unable to perform Q.R decomposition on the {0}x{1} jacobian matrix", new Object[] { new Integer(rows), new Integer(cols) }); } if (norm2 > ak2) { nextColumn = i; ak2 = norm2; } } if (ak2 == 0) { rank = k; return; } int pk = permutation[nextColumn]; permutation[nextColumn] = permutation[k]; permutation[k] = pk; // choose alpha such that Hk.u = alpha ek int kDiag = k * cols + pk; double akk = jacobian[kDiag]; double alpha = (akk > 0) ? -Math.sqrt(ak2) : Math.sqrt(ak2); double betak = 1.0 / (ak2 - akk * alpha); beta[pk] = betak; // transform the current column diagR[pk] = alpha; jacobian[kDiag] -= alpha; // transform the remaining columns for (int dk = cols - 1 - k; dk > 0; --dk) { int dkp = permutation[k + dk] - pk; double gamma = 0; for (int index = kDiag; index < jacobian.length; index += cols) { gamma += jacobian[index] * jacobian[index + dkp]; } gamma *= betak; for (int index = kDiag; index < jacobian.length; index += cols) { jacobian[index + dkp] -= gamma * jacobian[index]; } } } 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]; int kDiag = k * cols + pk; double gamma = 0; for (int i = k, index = kDiag; i < rows; ++i, index += cols) { gamma += jacobian[index] * y[i]; } gamma *= beta[pk]; for (int i = k, index = kDiag; i < rows; ++i, index += cols) { y[i] -= gamma * jacobian[index]; } } } /** Number of solved variables. */ private int solvedCols; /** Diagonal elements of the R matrix in the Q.R. decomposition. */ private double[] diagR; /** Norms of the columns of the jacobian matrix. */ private double[] jacNorm; /** Coefficients of the Householder transforms vectors. */ private double[] beta; /** Columns permutation array. */ private int[] permutation; /** Rank of the jacobian matrix. */ private int rank; /** Levenberg-Marquardt parameter. */ private double lmPar; /** Parameters evolution direction associated with lmPar. */ private double[] lmDir; /** Positive input variable used in determining the initial step bound. */ private double initialStepBoundFactor; /** Desired relative error in the sum of squares. */ private double costRelativeTolerance; /** Desired relative error in the approximate solution parameters. */ private double parRelativeTolerance; /** Desired max cosine on the orthogonality between the function vector * and the columns of the jacobian. */ private double orthoTolerance; /** Serializable version identifier */ private static final long serialVersionUID = -5705952631533171019L; }

The table below shows all metrics for LevenbergMarquardtEstimator.java.

MetricValueDescription
BLOCKS95.00Number of blocks
BLOCK_COMMENT16.00Number of block comment lines
COMMENTS318.00Comment lines
COMMENT_DENSITY 0.87Comment density
COMPARISONS104.00Number of comparison operators
CYCLOMATIC106.00Cyclomatic complexity
DECL_COMMENTS25.00Comments in declarations
DOC_COMMENT229.00Number of javadoc comment lines
ELOC367.00Effective lines of code
EXEC_COMMENTS55.00Comments in executable code
EXITS22.00Procedure exits
FUNCTIONS10.00Number of function declarations
HALSTEAD_DIFFICULTY140.39Halstead difficulty
HALSTEAD_EFFORT 0.00Halstead effort
INTERFACE_COMPLEXITY38.00Interface complexity
JAVA0001 0.00JAVA0001 Package name does not contain only lower case letters
JAVA0002 0.00JAVA0002 Package name does not begin with a top level domain name or country code
JAVA0003 0.00JAVA0003 Minimize use of on-demand (.*) imports
JAVA0004 0.00JAVA0004 Unnecessary import from java.lang
JAVA0005 0.00JAVA0005 Imports not in specified order
JAVA0006 0.00JAVA0006 Empty finally block
JAVA0007 0.00JAVA0007 Should not declare public field
JAVA0008 0.00JAVA0008 Empty catch block
JAVA0009 0.00JAVA0009 Protected member in final class
JAVA0010 0.00JAVA0010 Non-instantiable class does not contain a non-private static member
JAVA0011 0.00JAVA0011 Abstract class does not contain an abstract method
JAVA0012 0.00JAVA0012 Non-constructor method with same name as declaring class
JAVA0013 0.00JAVA0013 Non-blank final field is not static
JAVA0014 0.00JAVA0014 Class with only static members has non-private constructor
JAVA0015 0.00JAVA0015 Package class contains public nested type
JAVA0016 0.00JAVA0016 Abstract class contains public constructor
JAVA0017 0.00JAVA0017 Class name does not have required form
JAVA0018 0.00JAVA0018 Method name does not have required form
JAVA0019 0.00JAVA0019 Interface name does not have required form
JAVA0020 0.00JAVA0020 Field name does not have required form
JAVA0021 0.00JAVA0021 Interface method name does not have required form
JAVA0022 0.00JAVA0022 Static final field name does not have required form
JAVA0023 0.00JAVA0023 Empty finalize method
JAVA0024 0.00JAVA0024 Empty class
JAVA0025 0.00JAVA0025 Method override is empty
JAVA0026 0.00JAVA0026 Finalize method with parameters
JAVA0029 0.00JAVA0029 Private method not used
JAVA0030 0.00JAVA0030 Private field not used
JAVA0031 0.00JAVA0031 Case statement not properly closed
JAVA0032 0.00JAVA0032 Switch statement missing default
JAVA0033 0.00JAVA0033 default: not last case in switch statement
JAVA0034 0.00JAVA0034 Missing braces in if statement
JAVA0035 0.00JAVA0035 Missing braces in for statement
JAVA0036 0.00JAVA0036 Missing braces in while statement
JAVA0038 0.00JAVA0038 Non-case label in switch statement
JAVA0039 0.00JAVA0039 Break statement with label
JAVA0040 0.00JAVA0040 Switch statement contains N cases (maximum: M)
JAVA0041 0.00JAVA0041 Nested synchronized block
JAVA0042 0.00JAVA0042 Empty synchronized statement
JAVA0043 0.00JAVA0043 Inner class does not use outer class
JAVA0044 0.00JAVA0044 Serializable class with no instance variables
JAVA0045 0.00JAVA0045 Serializable class with only transient fields
JAVA0046 0.00JAVA0046 Name of class not derived from Exception ends with 'Exception'
JAVA0047 0.00JAVA0047 Serializable class derives from invalid base class
JAVA0048 0.00JAVA0048 Name of class derived from Exception does not end with 'Exception'
JAVA0049 1.00JAVA0049 Nested block at depth N (maximum: M)
JAVA0050 0.00JAVA0050 Class derives from java.lang.Error
JAVA0051 0.00JAVA0051 Class derives from java.lang.RuntimeException
JAVA0052 0.00JAVA0052 Class derives from java.lang.Throwable
JAVA0053 0.00JAVA0053 Unused label
JAVA0054 0.00JAVA0054 Inheritance depth N exceeds maximum M
JAVA0055 0.00JAVA0055 Class should be interface
JAVA0056 0.00JAVA0056 Unnecessary abstract modifier for interface or annotation
JAVA0057 0.00JAVA0057 Unnecessary default constructor
JAVA0058 0.00JAVA0058 Constructor calls super()
JAVA0059 0.00JAVA0059 Method override only calls super()
JAVA0061 0.00JAVA0061 Inaccessible member in anonymous class
JAVA0062 0.00JAVA0062 Public class missing public member or protected constructor
JAVA0063 0.00JAVA0063 Identifier name should not contain '$'
JAVA0064 1.00JAVA0064 N variations of identifier name (maximum: M)
JAVA0065 0.00JAVA0065 Unnecessary final modifier for method in final class
JAVA0066 0.00JAVA0066 Unnecessary modifier for interface nested type
JAVA0067 0.00JAVA0067 Array descriptor on identifier name
JAVA0068 0.00JAVA0068 Modifiers not declared in recommended order
JAVA0071 0.00JAVA0071 Strings compared with ==
JAVA0073 0.00JAVA0073 Integer division in floating-point context
JAVA0074 0.00JAVA0074 Use of Object.notify()
JAVA0075 0.00JAVA0075 Method parameter hides field
JAVA007622.00JAVA0076 Use of magic number
JAVA0077 0.00JAVA0077 Private field not used in declaring class
JAVA007814.00JAVA0078 Floating point values compared with ==
JAVA0079 0.00JAVA0079 Use of instance to reference static member
JAVA0080 0.00JAVA0080 Import declaration not used
JAVA0081 0.00JAVA0081 Boolean literal in comparison
JAVA0082 0.00JAVA0082 Unnecessary widening cast
JAVA0083 0.00JAVA0083 Unnecessary instanceof test
JAVA0084 1.00JAVA0084 Should use compound assignment operator
JAVA0085 0.00JAVA0085 Use of sun.* class
JAVA0087 0.00JAVA0087 Use of Thread.sleep()
JAVA0089 0.00JAVA0089 Use of restricted package
JAVA0092 0.00JAVA0092 Use of restricted type
JAVA0093 0.00JAVA0093 Redundant assignment
JAVA0094 0.00JAVA0094 Field hides a superclass field
JAVA0095 0.00JAVA0095 Uninitialized private field
JAVA0096 0.00JAVA0096 Field in nested class hides outer field
JAVA0098 0.00JAVA0098 Minimize use of implicit field initializers
JAVA0100 1.00JAVA0100 Class contains N non-final fields (maximum: M)
JAVA0101 0.00JAVA0101 Unnecessary modifier for field in interface
JAVA0102 0.00JAVA0102 Last statement in finalize() not super.finalize()
JAVA0103 0.00JAVA0103 Explicit call to finalize()
JAVA0104 0.00JAVA0104 finalize() only calls super.finalize()
JAVA0105 0.00JAVA0105 Duplicate import declaration
JAVA0106 0.00JAVA0106 Unnecessary import from current package
JAVA0108 0.00JAVA0108 Incorrect javadoc: no @param tag for 'parameter'
JAVA0109 0.00JAVA0109 Incorrect javadoc: no parameter 'parameter'
JAVA0110 0.00JAVA0110 Incorrect javadoc: no @return tag
JAVA0111 0.00JAVA0111 Incorrect javadoc: @return tag for void method
JAVA0112 0.00JAVA0112 Incorrect javadoc: no exception 'exception' in throws
JAVA0113 0.00JAVA0113 Incorrect javadoc: no @author tag
JAVA0114 0.00JAVA0114 Incorrect javadoc: no @version tag
JAVA0115 0.00JAVA0115 Incorrect javadoc: no @throws or @exception tag for 'exception'
JAVA0116 0.00JAVA0116 Missing javadoc: field 'field'
JAVA0117 0.00JAVA0117 Missing javadoc: method 'method'
JAVA0118 0.00JAVA0118 Missing javadoc: type 'type'
JAVA0119 1.00JAVA0119 Control variable changed within body of for loop
JAVA0123 1.00JAVA0123 Use all three components of for loop
JAVA0125 0.00JAVA0125 Continue statement with label
JAVA0126 0.00JAVA0126 Method declares unchecked exception in throws
JAVA0128 0.00JAVA0128 Public constructor in non-public class
JAVA0130 0.00JAVA0130 Non-static method does not use instance fields
JAVA0131 0.00JAVA0131 Compatible method does not override base
JAVA0132 0.00JAVA0132 Method overload with compatible signature
JAVA0133 0.00JAVA0133 Non-synchronized method overrides synchronized method
JAVA0135 0.00JAVA0135 Only one of Object.equals and Object.hashCode defined: missing 'method'
JAVA0136 0.00JAVA0136 N methods defined in class (maximum: M)
JAVA0137 0.00JAVA0137 Non-abstract class missing constructor
JAVA0138 1.00JAVA0138 N parameters defined for method (maximum: M)
JAVA0139 0.00JAVA0139 Definition of main other than public static void main(java.lang.String[])
JAVA0141 0.00JAVA0141 Unnecessary modifier for method in interface
JAVA0143 0.00JAVA0143 Synchronized method
JAVA0144 0.00JAVA0144 Line exceeds maximum M characters
JAVA0145 0.00JAVA0145 Tab character used in source file
JAVA0150 0.00JAVA0150 java.lang.Error (or subclass) thrown
JAVA0153 0.00JAVA0153 Inefficient conversion of integer to string
JAVA0159 0.00JAVA0159 Inefficient conversion of string to integer
JAVA0160 0.00JAVA0160 Method does not throw specified exception
JAVA0161 0.00JAVA0161 Conditional wait() not in loop
JAVA0163 0.00JAVA0163 Empty statement
JAVA0165 0.00JAVA0165 Conflicting return statement in finally block
JAVA0166 0.00JAVA0166 Generic exception caught
JAVA0167 0.00JAVA0167 ThreadDeath not rethrown
JAVA0169 0.00JAVA0169 Unnecessary catch block: exception 'exception'
JAVA0170 0.00JAVA0170 Caught exception not derived from java.lang.Exception
JAVA0171 0.00JAVA0171 Unused local variable
JAVA0173 0.00JAVA0173 Unused method parameter
JAVA0174 0.00JAVA0174 Assigned local variable never used
JAVA0175 0.00JAVA0175 Successive assignment to variable
JAVA0176 0.00JAVA0176 Local variable name does not have required form
JAVA0177 3.00JAVA0177 Variable declaration missing initializer
JAVA0179 0.00JAVA0179 Local variable hides visible field
JAVA0233 0.00JAVA0233 Definition of serialVersionUID other than 'private static final long serialVersionUID'
JAVA0234 0.00JAVA0234 Class is Serializable but does not define serialVersionUID
JAVA0235 0.00JAVA0235 Class defines serialVersionUID but does not implement Serializable
JAVA0236 0.00JAVA0236 Attempt to clone an object which does not implement Cloneable
JAVA0237 0.00JAVA0237 Class implements Cloneable but does not have public clone method
JAVA0238 0.00JAVA0238 Clone method does not call super.clone()
JAVA0239 0.00JAVA0239 Class declares 'readObject' or 'writeObject' but does not implement Serializable
JAVA0240 0.00JAVA0240 Serializable class which declares readObject or writeObject but not both
JAVA0241 0.00JAVA0241 'readObject' or 'writeObject' should be declared private in Serializable class
JAVA0242 0.00JAVA0242 Transient field in non-Serializable class
JAVA0243 0.00JAVA0243 'readResolve' or 'writeReplace' should be declared private or protected
JAVA0244 0.00JAVA0244 Field or method name in subclass differs only by case from inherited field or method
JAVA0245 0.00JAVA0245 JUnit TestCase with non-trivial constructor
JAVA0246 0.00JAVA0246 JUnit assertXXX statement missing message parameter
JAVA0247 0.00JAVA0247 JUnit 'setUp()' and 'tearDown()' should call super method
JAVA0248 0.00JAVA0248 JUnit method 'setUp' or 'tearDown' with incorrect signature
JAVA0249 0.00JAVA0249 JUnit TestCase 'suite()' should be declared static
JAVA0250 0.00JAVA0250 JUnit TestCase declares testXXX method with incorrect signature
JAVA0251 0.00JAVA0251 Use '%n' for line breaks in printf/format for platform independence
JAVA0252 0.00JAVA0252 'enum' is a Java 1.5 reserved word
JAVA0253 0.00JAVA0253 Not all enum constants consumed in switch statement
JAVA0254 0.00JAVA0254 Use enhanced for loop construct instead of Iterator
JAVA0255 0.00JAVA0255 Result of method invocation not used
JAVA0256 0.00JAVA0256 Assignment of external collection/array to field
JAVA0257 0.00JAVA0257 Use of 'Constant Interface' anti-pattern
JAVA0258 0.00JAVA0258 Implement Iterable for foreach compatibility
JAVA0259 0.00JAVA0259 Return of collection/array field
JAVA0260 0.00JAVA0260 Use 'enum' instead of Enumerated Type pattern
JAVA0261 0.00JAVA0261 Use specialized Enum collection types
JAVA0262 0.00JAVA0262 Use of char in integer context
JAVA0263 0.00JAVA0263 Long literal ends with 'l' instead of 'L'
JAVA0264 0.00JAVA0264 Integer math in long context - check for overflow
JAVA0265 0.00JAVA0265 Use of Throwable.printStackTrace()
JAVA0266 0.00JAVA0266 Use of System.out
JAVA0267 0.00JAVA0267 Use of System.err
JAVA0269 0.00JAVA0269 Contents of StringBuffer never used
JAVA0270 0.00JAVA0270 Use Java 5.0 enhanced for loop construct to iterate over all elements in an array
JAVA0271 0.00JAVA0271 Minimize use of on-demand (.*) static imports
JAVA0272 0.00JAVA0272 Thread.run() called
JAVA0273 0.00JAVA0273 Non-final derivative of Thread calls start() in constructor
JAVA0274 0.00JAVA0274 Serializable class has a synchronized readObject()
JAVA0275 0.00JAVA0275 Serializable class has a synchronized writeObject() and no other synchronized methods
JAVA0276 0.00JAVA0276 Unnecessary use of String constructor
JAVA0277 0.00JAVA0277 Iterator.next() implementation does not throw NoSuchElementException
JAVA0278 0.00JAVA0278 Unnecessary use of Boolean constructor
JAVA0279 0.00JAVA0279 Serialization method readObject or readObjectNoData calls an overridable method
JAVA0280 0.00JAVA0280 IllegalMonitorStateException caught
JAVA0281 0.00JAVA0281 Iterator.next() not called in loop
JAVA0282 0.00JAVA0282 Call to Iterator.next() in loop which does not test Iterator.hasNext()
JAVA0283 0.00JAVA0283 Control variable not updated in loop body
JAVA0284 0.00JAVA0284 Explicit garbage collection
JAVA0285 0.00JAVA0285 Dereference of potentially null variable
JAVA0286 0.00JAVA0286 Dereference of null variable
JAVA0287 0.00JAVA0287 Unnecessary null check
JAVA0288 0.00JAVA0288 Inconsistent null check
LINES871.00Number of lines in the source file
LINE_COMMENT73.00Number of line comments
LOC457.00Lines of code
LOGICAL_LINES354.00Number of statements
LOOPS51.00Number of loops
NEST_DEPTH 6.00Maximum nesting depth
OPERANDS1312.00Number of operands
OPERATORS2366.00Number of operators
PARAMS16.00Number of formal parameter declarations
PROGRAM_LENGTH3678.00Halstead program length
PROGRAM_VOCAB312.00Halstead program vocabulary
PROGRAM_VOLUME 0.00Halstead program volume
RETURNS22.00Number of return points from functions
SIZE31312.00Size of the file in bytes
UNIQUE_OPERANDS257.00Number of unique operands
UNIQUE_OPERATORS55.00Number of unique operators
WHITESPACE96.00Number of whitespace lines