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