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