1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.math.optimization.general; 18 19 import java.util.Arrays; 20 21 import org.apache.commons.math.FunctionEvaluationException; 22 import org.apache.commons.math.exception.util.LocalizedFormats; 23 import org.apache.commons.math.optimization.OptimizationException; 24 import org.apache.commons.math.optimization.VectorialPointValuePair; 25 import org.apache.commons.math.util.FastMath; 26 import org.apache.commons.math.util.MathUtils; 27 28 29 /** 30 * This class solves a least squares problem using the Levenberg-Marquardt algorithm. 31 * 32 * <p>This implementation <em>should</em> work even for over-determined systems 33 * (i.e. systems having more point than equations). Over-determined systems 34 * are solved by ignoring the point which have the smallest impact according 35 * to their jacobian column norm. Only the rank of the matrix and some loop bounds 36 * are changed to implement this.</p> 37 * 38 * <p>The resolution engine is a simple translation of the MINPACK <a 39 * href="http://www.netlib.org/minpack/lmder.f">lmder</a> routine with minor 40 * changes. The changes include the over-determined resolution, the use of 41 * inherited convergence checker and the Q.R. decomposition which has been 42 * rewritten following the algorithm described in the 43 * P. Lascaux and R. Theodor book <i>Analyse numérique matricielle 44 * appliquée à l'art de l'ingénieur</i>, Masson 1986.</p> 45 * <p>The authors of the original fortran version are: 46 * <ul> 47 * <li>Argonne National Laboratory. MINPACK project. March 1980</li> 48 * <li>Burton S. Garbow</li> 49 * <li>Kenneth E. Hillstrom</li> 50 * <li>Jorge J. More</li> 51 * </ul> 52 * The redistribution policy for MINPACK is available <a 53 * href="http://www.netlib.org/minpack/disclaimer">here</a>, for convenience, it 54 * is reproduced below.</p> 55 * 56 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> 57 * <tr><td> 58 * Minpack Copyright Notice (1999) University of Chicago. 59 * All rights reserved 60 * </td></tr> 61 * <tr><td> 62 * Redistribution and use in source and binary forms, with or without 63 * modification, are permitted provided that the following conditions 64 * are met: 65 * <ol> 66 * <li>Redistributions of source code must retain the above copyright 67 * notice, this list of conditions and the following disclaimer.</li> 68 * <li>Redistributions in binary form must reproduce the above 69 * copyright notice, this list of conditions and the following 70 * disclaimer in the documentation and/or other materials provided 71 * with the distribution.</li> 72 * <li>The end-user documentation included with the redistribution, if any, 73 * must include the following acknowledgment: 74 * <code>This product includes software developed by the University of 75 * Chicago, as Operator of Argonne National Laboratory.</code> 76 * Alternately, this acknowledgment may appear in the software itself, 77 * if and wherever such third-party acknowledgments normally appear.</li> 78 * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" 79 * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE 80 * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND 81 * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR 82 * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES 83 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE 84 * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY 85 * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR 86 * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF 87 * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) 88 * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION 89 * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL 90 * BE CORRECTED.</strong></li> 91 * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT 92 * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF 93 * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, 94 * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF 95 * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF 96 * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER 97 * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT 98 * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, 99 * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE 100 * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li> 101 * <ol></td></tr> 102 * </table> 103 * @version $Revision: 1073272 $ $Date: 2011-02-22 10:22:25 +0100 (mar. 22 fvr. 2011) $ 104 * @since 2.0 105 * 106 */ 107 public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer { 108 109 /** Number of solved point. */ 110 private int solvedCols; 111 112 /** Diagonal elements of the R matrix in the Q.R. decomposition. */ 113 private double[] diagR; 114 115 /** Norms of the columns of the jacobian matrix. */ 116 private double[] jacNorm; 117 118 /** Coefficients of the Householder transforms vectors. */ 119 private double[] beta; 120 121 /** Columns permutation array. */ 122 private int[] permutation; 123 124 /** Rank of the jacobian matrix. */ 125 private int rank; 126 127 /** Levenberg-Marquardt parameter. */ 128 private double lmPar; 129 130 /** Parameters evolution direction associated with lmPar. */ 131 private double[] lmDir; 132 133 /** Positive input variable used in determining the initial step bound. */ 134 private double initialStepBoundFactor; 135 136 /** Desired relative error in the sum of squares. */ 137 private double costRelativeTolerance; 138 139 /** Desired relative error in the approximate solution parameters. */ 140 private double parRelativeTolerance; 141 142 /** Desired max cosine on the orthogonality between the function vector 143 * and the columns of the jacobian. */ 144 private double orthoTolerance; 145 146 /** Threshold for QR ranking. */ 147 private double qrRankingThreshold; 148 149 /** 150 * Build an optimizer for least squares problems. 151 * <p>The default values for the algorithm settings are: 152 * <ul> 153 * <li>{@link #setConvergenceChecker(VectorialConvergenceChecker) vectorial convergence checker}: null</li> 154 * <li>{@link #setInitialStepBoundFactor(double) initial step bound factor}: 100.0</li> 155 * <li>{@link #setMaxIterations(int) maximal iterations}: 1000</li> 156 * <li>{@link #setCostRelativeTolerance(double) cost relative tolerance}: 1.0e-10</li> 157 * <li>{@link #setParRelativeTolerance(double) parameters relative tolerance}: 1.0e-10</li> 158 * <li>{@link #setOrthoTolerance(double) orthogonality tolerance}: 1.0e-10</li> 159 * <li>{@link #setQRRankingThreshold(double) QR ranking threshold}: {@link MathUtils#SAFE_MIN}</li> 160 * </ul> 161 * </p> 162 * <p>These default values may be overridden after construction. If the {@link 163 * #setConvergenceChecker vectorial convergence checker} is set to a non-null value, it 164 * will be used instead of the {@link #setCostRelativeTolerance cost relative tolerance} 165 * and {@link #setParRelativeTolerance parameters relative tolerance} settings. 166 */ 167 public LevenbergMarquardtOptimizer() { 168 169 // set up the superclass with a default max cost evaluations setting 170 setMaxIterations(1000); 171 172 // default values for the tuning parameters 173 setConvergenceChecker(null); 174 setInitialStepBoundFactor(100.0); 175 setCostRelativeTolerance(1.0e-10); 176 setParRelativeTolerance(1.0e-10); 177 setOrthoTolerance(1.0e-10); 178 setQRRankingThreshold(MathUtils.SAFE_MIN); 179 180 } 181 182 /** 183 * Set the positive input variable used in determining the initial step bound. 184 * This bound is set to the product of initialStepBoundFactor and the euclidean 185 * norm of diag*x if nonzero, or else to initialStepBoundFactor itself. In most 186 * cases factor should lie in the interval (0.1, 100.0). 100.0 is a generally 187 * recommended value. 188 * 189 * @param initialStepBoundFactor initial step bound factor 190 */ 191 public void setInitialStepBoundFactor(double initialStepBoundFactor) { 192 this.initialStepBoundFactor = initialStepBoundFactor; 193 } 194 195 /** 196 * Set the desired relative error in the sum of squares. 197 * <p>This setting is used only if the {@link #setConvergenceChecker vectorial 198 * convergence checker} is set to null.</p> 199 * @param costRelativeTolerance desired relative error in the sum of squares 200 */ 201 public void setCostRelativeTolerance(double costRelativeTolerance) { 202 this.costRelativeTolerance = costRelativeTolerance; 203 } 204 205 /** 206 * Set the desired relative error in the approximate solution parameters. 207 * <p>This setting is used only if the {@link #setConvergenceChecker vectorial 208 * convergence checker} is set to null.</p> 209 * @param parRelativeTolerance desired relative error 210 * in the approximate solution parameters 211 */ 212 public void setParRelativeTolerance(double parRelativeTolerance) { 213 this.parRelativeTolerance = parRelativeTolerance; 214 } 215 216 /** 217 * Set the desired max cosine on the orthogonality. 218 * <p>This setting is always used, regardless of the {@link #setConvergenceChecker 219 * vectorial convergence checker} being null or non-null.</p> 220 * @param orthoTolerance desired max cosine on the orthogonality 221 * between the function vector and the columns of the jacobian 222 */ 223 public void setOrthoTolerance(double orthoTolerance) { 224 this.orthoTolerance = orthoTolerance; 225 } 226 227 /** 228 * Set the desired threshold for QR ranking. 229 * <p> 230 * If the squared norm of a column vector is smaller or equal to this threshold 231 * during QR decomposition, it is considered to be a zero vector and hence the 232 * rank of the matrix is reduced. 233 * </p> 234 * @param threshold threshold for QR ranking 235 * @since 2.2 236 */ 237 public void setQRRankingThreshold(final double threshold) { 238 this.qrRankingThreshold = threshold; 239 } 240 241 /** {@inheritDoc} */ 242 @Override 243 protected VectorialPointValuePair doOptimize() 244 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException { 245 246 // arrays shared with the other private methods 247 solvedCols = Math.min(rows, cols); 248 diagR = new double[cols]; 249 jacNorm = new double[cols]; 250 beta = new double[cols]; 251 permutation = new int[cols]; 252 lmDir = new double[cols]; 253 254 // local point 255 double delta = 0; 256 double xNorm = 0; 257 double[] diag = new double[cols]; 258 double[] oldX = new double[cols]; 259 double[] oldRes = new double[rows]; 260 double[] oldObj = new double[rows]; 261 double[] qtf = new double[rows]; 262 double[] work1 = new double[cols]; 263 double[] work2 = new double[cols]; 264 double[] work3 = new double[cols]; 265 266 // evaluate the function at the starting point and calculate its norm 267 updateResidualsAndCost(); 268 269 // outer loop 270 lmPar = 0; 271 boolean firstIteration = true; 272 VectorialPointValuePair current = new VectorialPointValuePair(point, objective); 273 while (true) { 274 for (int i=0;i<rows;i++) { 275 qtf[i]=wresiduals[i]; 276 } 277 incrementIterationsCounter(); 278 279 // compute the Q.R. decomposition of the jacobian matrix 280 VectorialPointValuePair previous = current; 281 updateJacobian(); 282 qrDecomposition(); 283 284 // compute Qt.res 285 qTy(qtf); 286 // now we don't need Q anymore, 287 // so let jacobian contain the R matrix with its diagonal elements 288 for (int k = 0; k < solvedCols; ++k) { 289 int pk = permutation[k]; 290 wjacobian[k][pk] = diagR[pk]; 291 } 292 293 if (firstIteration) { 294 295 // scale the point according to the norms of the columns 296 // of the initial jacobian 297 xNorm = 0; 298 for (int k = 0; k < cols; ++k) { 299 double dk = jacNorm[k]; 300 if (dk == 0) { 301 dk = 1.0; 302 } 303 double xk = dk * point[k]; 304 xNorm += xk * xk; 305 diag[k] = dk; 306 } 307 xNorm = FastMath.sqrt(xNorm); 308 309 // initialize the step bound delta 310 delta = (xNorm == 0) ? initialStepBoundFactor : (initialStepBoundFactor * xNorm); 311 312 } 313 314 // check orthogonality between function vector and jacobian columns 315 double maxCosine = 0; 316 if (cost != 0) { 317 for (int j = 0; j < solvedCols; ++j) { 318 int pj = permutation[j]; 319 double s = jacNorm[pj]; 320 if (s != 0) { 321 double sum = 0; 322 for (int i = 0; i <= j; ++i) { 323 sum += wjacobian[i][pj] * qtf[i]; 324 } 325 maxCosine = FastMath.max(maxCosine, FastMath.abs(sum) / (s * cost)); 326 } 327 } 328 } 329 if (maxCosine <= orthoTolerance) { 330 // convergence has been reached 331 updateResidualsAndCost(); 332 current = new VectorialPointValuePair(point, objective); 333 return current; 334 } 335 336 // rescale if necessary 337 for (int j = 0; j < cols; ++j) { 338 diag[j] = FastMath.max(diag[j], jacNorm[j]); 339 } 340 341 // inner loop 342 for (double ratio = 0; ratio < 1.0e-4;) { 343 344 // save the state 345 for (int j = 0; j < solvedCols; ++j) { 346 int pj = permutation[j]; 347 oldX[pj] = point[pj]; 348 } 349 double previousCost = cost; 350 double[] tmpVec = residuals; 351 residuals = oldRes; 352 oldRes = tmpVec; 353 tmpVec = objective; 354 objective = oldObj; 355 oldObj = tmpVec; 356 357 // determine the Levenberg-Marquardt parameter 358 determineLMParameter(qtf, delta, diag, work1, work2, work3); 359 360 // compute the new point and the norm of the evolution direction 361 double lmNorm = 0; 362 for (int j = 0; j < solvedCols; ++j) { 363 int pj = permutation[j]; 364 lmDir[pj] = -lmDir[pj]; 365 point[pj] = oldX[pj] + lmDir[pj]; 366 double s = diag[pj] * lmDir[pj]; 367 lmNorm += s * s; 368 } 369 lmNorm = FastMath.sqrt(lmNorm); 370 // on the first iteration, adjust the initial step bound. 371 if (firstIteration) { 372 delta = FastMath.min(delta, lmNorm); 373 } 374 375 // evaluate the function at x + p and calculate its norm 376 updateResidualsAndCost(); 377 378 // compute the scaled actual reduction 379 double actRed = -1.0; 380 if (0.1 * cost < previousCost) { 381 double r = cost / previousCost; 382 actRed = 1.0 - r * r; 383 } 384 385 // compute the scaled predicted reduction 386 // and the scaled directional derivative 387 for (int j = 0; j < solvedCols; ++j) { 388 int pj = permutation[j]; 389 double dirJ = lmDir[pj]; 390 work1[j] = 0; 391 for (int i = 0; i <= j; ++i) { 392 work1[i] += wjacobian[i][pj] * dirJ; 393 } 394 } 395 double coeff1 = 0; 396 for (int j = 0; j < solvedCols; ++j) { 397 coeff1 += work1[j] * work1[j]; 398 } 399 double pc2 = previousCost * previousCost; 400 coeff1 = coeff1 / pc2; 401 double coeff2 = lmPar * lmNorm * lmNorm / pc2; 402 double preRed = coeff1 + 2 * coeff2; 403 double dirDer = -(coeff1 + coeff2); 404 405 // ratio of the actual to the predicted reduction 406 ratio = (preRed == 0) ? 0 : (actRed / preRed); 407 408 // update the step bound 409 if (ratio <= 0.25) { 410 double tmp = 411 (actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5; 412 if ((0.1 * cost >= previousCost) || (tmp < 0.1)) { 413 tmp = 0.1; 414 } 415 delta = tmp * FastMath.min(delta, 10.0 * lmNorm); 416 lmPar /= tmp; 417 } else if ((lmPar == 0) || (ratio >= 0.75)) { 418 delta = 2 * lmNorm; 419 lmPar *= 0.5; 420 } 421 422 // test for successful iteration. 423 if (ratio >= 1.0e-4) { 424 // successful iteration, update the norm 425 firstIteration = false; 426 xNorm = 0; 427 for (int k = 0; k < cols; ++k) { 428 double xK = diag[k] * point[k]; 429 xNorm += xK * xK; 430 } 431 xNorm = FastMath.sqrt(xNorm); 432 current = new VectorialPointValuePair(point, objective); 433 434 // tests for convergence. 435 if (checker != null) { 436 // we use the vectorial convergence checker 437 if (checker.converged(getIterations(), previous, current)) { 438 return current; 439 } 440 } 441 } else { 442 // failed iteration, reset the previous values 443 cost = previousCost; 444 for (int j = 0; j < solvedCols; ++j) { 445 int pj = permutation[j]; 446 point[pj] = oldX[pj]; 447 } 448 tmpVec = residuals; 449 residuals = oldRes; 450 oldRes = tmpVec; 451 tmpVec = objective; 452 objective = oldObj; 453 oldObj = tmpVec; 454 } 455 if (checker==null) { 456 if (((FastMath.abs(actRed) <= costRelativeTolerance) && 457 (preRed <= costRelativeTolerance) && 458 (ratio <= 2.0)) || 459 (delta <= parRelativeTolerance * xNorm)) { 460 return current; 461 } 462 } 463 // tests for termination and stringent tolerances 464 // (2.2204e-16 is the machine epsilon for IEEE754) 465 if ((FastMath.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) { 466 throw new OptimizationException(LocalizedFormats.TOO_SMALL_COST_RELATIVE_TOLERANCE, 467 costRelativeTolerance); 468 } else if (delta <= 2.2204e-16 * xNorm) { 469 throw new OptimizationException(LocalizedFormats.TOO_SMALL_PARAMETERS_RELATIVE_TOLERANCE, 470 parRelativeTolerance); 471 } else if (maxCosine <= 2.2204e-16) { 472 throw new OptimizationException(LocalizedFormats.TOO_SMALL_ORTHOGONALITY_TOLERANCE, 473 orthoTolerance); 474 } 475 476 } 477 478 } 479 480 } 481 482 /** 483 * Determine the Levenberg-Marquardt parameter. 484 * <p>This implementation is a translation in Java of the MINPACK 485 * <a href="http://www.netlib.org/minpack/lmpar.f">lmpar</a> 486 * routine.</p> 487 * <p>This method sets the lmPar and lmDir attributes.</p> 488 * <p>The authors of the original fortran function are:</p> 489 * <ul> 490 * <li>Argonne National Laboratory. MINPACK project. March 1980</li> 491 * <li>Burton S. Garbow</li> 492 * <li>Kenneth E. Hillstrom</li> 493 * <li>Jorge J. More</li> 494 * </ul> 495 * <p>Luc Maisonobe did the Java translation.</p> 496 * 497 * @param qy array containing qTy 498 * @param delta upper bound on the euclidean norm of diagR * lmDir 499 * @param diag diagonal matrix 500 * @param work1 work array 501 * @param work2 work array 502 * @param work3 work array 503 */ 504 private void determineLMParameter(double[] qy, double delta, double[] diag, 505 double[] work1, double[] work2, double[] work3) { 506 507 // compute and store in x the gauss-newton direction, if the 508 // jacobian is rank-deficient, obtain a least squares solution 509 for (int j = 0; j < rank; ++j) { 510 lmDir[permutation[j]] = qy[j]; 511 } 512 for (int j = rank; j < cols; ++j) { 513 lmDir[permutation[j]] = 0; 514 } 515 for (int k = rank - 1; k >= 0; --k) { 516 int pk = permutation[k]; 517 double ypk = lmDir[pk] / diagR[pk]; 518 for (int i = 0; i < k; ++i) { 519 lmDir[permutation[i]] -= ypk * wjacobian[i][pk]; 520 } 521 lmDir[pk] = ypk; 522 } 523 524 // evaluate the function at the origin, and test 525 // for acceptance of the Gauss-Newton direction 526 double dxNorm = 0; 527 for (int j = 0; j < solvedCols; ++j) { 528 int pj = permutation[j]; 529 double s = diag[pj] * lmDir[pj]; 530 work1[pj] = s; 531 dxNorm += s * s; 532 } 533 dxNorm = FastMath.sqrt(dxNorm); 534 double fp = dxNorm - delta; 535 if (fp <= 0.1 * delta) { 536 lmPar = 0; 537 return; 538 } 539 540 // if the jacobian is not rank deficient, the Newton step provides 541 // a lower bound, parl, for the zero of the function, 542 // otherwise set this bound to zero 543 double sum2; 544 double parl = 0; 545 if (rank == solvedCols) { 546 for (int j = 0; j < solvedCols; ++j) { 547 int pj = permutation[j]; 548 work1[pj] *= diag[pj] / dxNorm; 549 } 550 sum2 = 0; 551 for (int j = 0; j < solvedCols; ++j) { 552 int pj = permutation[j]; 553 double sum = 0; 554 for (int i = 0; i < j; ++i) { 555 sum += wjacobian[i][pj] * work1[permutation[i]]; 556 } 557 double s = (work1[pj] - sum) / diagR[pj]; 558 work1[pj] = s; 559 sum2 += s * s; 560 } 561 parl = fp / (delta * sum2); 562 } 563 564 // calculate an upper bound, paru, for the zero of the function 565 sum2 = 0; 566 for (int j = 0; j < solvedCols; ++j) { 567 int pj = permutation[j]; 568 double sum = 0; 569 for (int i = 0; i <= j; ++i) { 570 sum += wjacobian[i][pj] * qy[i]; 571 } 572 sum /= diag[pj]; 573 sum2 += sum * sum; 574 } 575 double gNorm = FastMath.sqrt(sum2); 576 double paru = gNorm / delta; 577 if (paru == 0) { 578 // 2.2251e-308 is the smallest positive real for IEE754 579 paru = 2.2251e-308 / FastMath.min(delta, 0.1); 580 } 581 582 // if the input par lies outside of the interval (parl,paru), 583 // set par to the closer endpoint 584 lmPar = FastMath.min(paru, FastMath.max(lmPar, parl)); 585 if (lmPar == 0) { 586 lmPar = gNorm / dxNorm; 587 } 588 589 for (int countdown = 10; countdown >= 0; --countdown) { 590 591 // evaluate the function at the current value of lmPar 592 if (lmPar == 0) { 593 lmPar = FastMath.max(2.2251e-308, 0.001 * paru); 594 } 595 double sPar = FastMath.sqrt(lmPar); 596 for (int j = 0; j < solvedCols; ++j) { 597 int pj = permutation[j]; 598 work1[pj] = sPar * diag[pj]; 599 } 600 determineLMDirection(qy, work1, work2, work3); 601 602 dxNorm = 0; 603 for (int j = 0; j < solvedCols; ++j) { 604 int pj = permutation[j]; 605 double s = diag[pj] * lmDir[pj]; 606 work3[pj] = s; 607 dxNorm += s * s; 608 } 609 dxNorm = FastMath.sqrt(dxNorm); 610 double previousFP = fp; 611 fp = dxNorm - delta; 612 613 // if the function is small enough, accept the current value 614 // of lmPar, also test for the exceptional cases where parl is zero 615 if ((FastMath.abs(fp) <= 0.1 * delta) || 616 ((parl == 0) && (fp <= previousFP) && (previousFP < 0))) { 617 return; 618 } 619 620 // compute the Newton correction 621 for (int j = 0; j < solvedCols; ++j) { 622 int pj = permutation[j]; 623 work1[pj] = work3[pj] * diag[pj] / dxNorm; 624 } 625 for (int j = 0; j < solvedCols; ++j) { 626 int pj = permutation[j]; 627 work1[pj] /= work2[j]; 628 double tmp = work1[pj]; 629 for (int i = j + 1; i < solvedCols; ++i) { 630 work1[permutation[i]] -= wjacobian[i][pj] * tmp; 631 } 632 } 633 sum2 = 0; 634 for (int j = 0; j < solvedCols; ++j) { 635 double s = work1[permutation[j]]; 636 sum2 += s * s; 637 } 638 double correction = fp / (delta * sum2); 639 640 // depending on the sign of the function, update parl or paru. 641 if (fp > 0) { 642 parl = FastMath.max(parl, lmPar); 643 } else if (fp < 0) { 644 paru = FastMath.min(paru, lmPar); 645 } 646 647 // compute an improved estimate for lmPar 648 lmPar = FastMath.max(parl, lmPar + correction); 649 650 } 651 } 652 653 /** 654 * Solve a*x = b and d*x = 0 in the least squares sense. 655 * <p>This implementation is a translation in Java of the MINPACK 656 * <a href="http://www.netlib.org/minpack/qrsolv.f">qrsolv</a> 657 * routine.</p> 658 * <p>This method sets the lmDir and lmDiag attributes.</p> 659 * <p>The authors of the original fortran function are:</p> 660 * <ul> 661 * <li>Argonne National Laboratory. MINPACK project. March 1980</li> 662 * <li>Burton S. Garbow</li> 663 * <li>Kenneth E. Hillstrom</li> 664 * <li>Jorge J. More</li> 665 * </ul> 666 * <p>Luc Maisonobe did the Java translation.</p> 667 * 668 * @param qy array containing qTy 669 * @param diag diagonal matrix 670 * @param lmDiag diagonal elements associated with lmDir 671 * @param work work array 672 */ 673 private void determineLMDirection(double[] qy, double[] diag, 674 double[] lmDiag, double[] work) { 675 676 // copy R and Qty to preserve input and initialize s 677 // in particular, save the diagonal elements of R in lmDir 678 for (int j = 0; j < solvedCols; ++j) { 679 int pj = permutation[j]; 680 for (int i = j + 1; i < solvedCols; ++i) { 681 wjacobian[i][pj] = wjacobian[j][permutation[i]]; 682 } 683 lmDir[j] = diagR[pj]; 684 work[j] = qy[j]; 685 } 686 687 // eliminate the diagonal matrix d using a Givens rotation 688 for (int j = 0; j < solvedCols; ++j) { 689 690 // prepare the row of d to be eliminated, locating the 691 // diagonal element using p from the Q.R. factorization 692 int pj = permutation[j]; 693 double dpj = diag[pj]; 694 if (dpj != 0) { 695 Arrays.fill(lmDiag, j + 1, lmDiag.length, 0); 696 } 697 lmDiag[j] = dpj; 698 699 // the transformations to eliminate the row of d 700 // modify only a single element of Qty 701 // beyond the first n, which is initially zero. 702 double qtbpj = 0; 703 for (int k = j; k < solvedCols; ++k) { 704 int pk = permutation[k]; 705 706 // determine a Givens rotation which eliminates the 707 // appropriate element in the current row of d 708 if (lmDiag[k] != 0) { 709 710 final double sin; 711 final double cos; 712 double rkk = wjacobian[k][pk]; 713 if (FastMath.abs(rkk) < FastMath.abs(lmDiag[k])) { 714 final double cotan = rkk / lmDiag[k]; 715 sin = 1.0 / FastMath.sqrt(1.0 + cotan * cotan); 716 cos = sin * cotan; 717 } else { 718 final double tan = lmDiag[k] / rkk; 719 cos = 1.0 / FastMath.sqrt(1.0 + tan * tan); 720 sin = cos * tan; 721 } 722 723 // compute the modified diagonal element of R and 724 // the modified element of (Qty,0) 725 wjacobian[k][pk] = cos * rkk + sin * lmDiag[k]; 726 final double temp = cos * work[k] + sin * qtbpj; 727 qtbpj = -sin * work[k] + cos * qtbpj; 728 work[k] = temp; 729 730 // accumulate the tranformation in the row of s 731 for (int i = k + 1; i < solvedCols; ++i) { 732 double rik = wjacobian[i][pk]; 733 final double temp2 = cos * rik + sin * lmDiag[i]; 734 lmDiag[i] = -sin * rik + cos * lmDiag[i]; 735 wjacobian[i][pk] = temp2; 736 } 737 738 } 739 } 740 741 // store the diagonal element of s and restore 742 // the corresponding diagonal element of R 743 lmDiag[j] = wjacobian[j][permutation[j]]; 744 wjacobian[j][permutation[j]] = lmDir[j]; 745 746 } 747 748 // solve the triangular system for z, if the system is 749 // singular, then obtain a least squares solution 750 int nSing = solvedCols; 751 for (int j = 0; j < solvedCols; ++j) { 752 if ((lmDiag[j] == 0) && (nSing == solvedCols)) { 753 nSing = j; 754 } 755 if (nSing < solvedCols) { 756 work[j] = 0; 757 } 758 } 759 if (nSing > 0) { 760 for (int j = nSing - 1; j >= 0; --j) { 761 int pj = permutation[j]; 762 double sum = 0; 763 for (int i = j + 1; i < nSing; ++i) { 764 sum += wjacobian[i][pj] * work[i]; 765 } 766 work[j] = (work[j] - sum) / lmDiag[j]; 767 } 768 } 769 770 // permute the components of z back to components of lmDir 771 for (int j = 0; j < lmDir.length; ++j) { 772 lmDir[permutation[j]] = work[j]; 773 } 774 775 } 776 777 /** 778 * Decompose a matrix A as A.P = Q.R using Householder transforms. 779 * <p>As suggested in the P. Lascaux and R. Theodor book 780 * <i>Analyse numérique matricielle appliquée à 781 * l'art de l'ingénieur</i> (Masson, 1986), instead of representing 782 * the Householder transforms with u<sub>k</sub> unit vectors such that: 783 * <pre> 784 * H<sub>k</sub> = I - 2u<sub>k</sub>.u<sub>k</sub><sup>t</sup> 785 * </pre> 786 * we use <sub>k</sub> non-unit vectors such that: 787 * <pre> 788 * H<sub>k</sub> = I - beta<sub>k</sub>v<sub>k</sub>.v<sub>k</sub><sup>t</sup> 789 * </pre> 790 * where v<sub>k</sub> = a<sub>k</sub> - alpha<sub>k</sub> e<sub>k</sub>. 791 * The beta<sub>k</sub> coefficients are provided upon exit as recomputing 792 * them from the v<sub>k</sub> vectors would be costly.</p> 793 * <p>This decomposition handles rank deficient cases since the tranformations 794 * are performed in non-increasing columns norms order thanks to columns 795 * pivoting. The diagonal elements of the R matrix are therefore also in 796 * non-increasing absolute values order.</p> 797 * @exception OptimizationException if the decomposition cannot be performed 798 */ 799 private void qrDecomposition() throws OptimizationException { 800 801 // initializations 802 for (int k = 0; k < cols; ++k) { 803 permutation[k] = k; 804 double norm2 = 0; 805 for (int i = 0; i < wjacobian.length; ++i) { 806 double akk = wjacobian[i][k]; 807 norm2 += akk * akk; 808 } 809 jacNorm[k] = FastMath.sqrt(norm2); 810 } 811 812 // transform the matrix column after column 813 for (int k = 0; k < cols; ++k) { 814 815 // select the column with the greatest norm on active components 816 int nextColumn = -1; 817 double ak2 = Double.NEGATIVE_INFINITY; 818 for (int i = k; i < cols; ++i) { 819 double norm2 = 0; 820 for (int j = k; j < wjacobian.length; ++j) { 821 double aki = wjacobian[j][permutation[i]]; 822 norm2 += aki * aki; 823 } 824 if (Double.isInfinite(norm2) || Double.isNaN(norm2)) { 825 throw new OptimizationException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN, 826 rows, cols); 827 } 828 if (norm2 > ak2) { 829 nextColumn = i; 830 ak2 = norm2; 831 } 832 } 833 if (ak2 <= qrRankingThreshold) { 834 rank = k; 835 return; 836 } 837 int pk = permutation[nextColumn]; 838 permutation[nextColumn] = permutation[k]; 839 permutation[k] = pk; 840 841 // choose alpha such that Hk.u = alpha ek 842 double akk = wjacobian[k][pk]; 843 double alpha = (akk > 0) ? -FastMath.sqrt(ak2) : FastMath.sqrt(ak2); 844 double betak = 1.0 / (ak2 - akk * alpha); 845 beta[pk] = betak; 846 847 // transform the current column 848 diagR[pk] = alpha; 849 wjacobian[k][pk] -= alpha; 850 851 // transform the remaining columns 852 for (int dk = cols - 1 - k; dk > 0; --dk) { 853 double gamma = 0; 854 for (int j = k; j < wjacobian.length; ++j) { 855 gamma += wjacobian[j][pk] * wjacobian[j][permutation[k + dk]]; 856 } 857 gamma *= betak; 858 for (int j = k; j < wjacobian.length; ++j) { 859 wjacobian[j][permutation[k + dk]] -= gamma * wjacobian[j][pk]; 860 } 861 } 862 863 } 864 865 rank = solvedCols; 866 867 } 868 869 /** 870 * Compute the product Qt.y for some Q.R. decomposition. 871 * 872 * @param y vector to multiply (will be overwritten with the result) 873 */ 874 private void qTy(double[] y) { 875 for (int k = 0; k < cols; ++k) { 876 int pk = permutation[k]; 877 double gamma = 0; 878 for (int i = k; i < rows; ++i) { 879 gamma += wjacobian[i][pk] * y[i]; 880 } 881 gamma *= beta[pk]; 882 for (int i = k; i < rows; ++i) { 883 y[i] -= gamma * wjacobian[i][pk]; 884 } 885 } 886 } 887 888 } 889
