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