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.optimization.OptimizationException;
23 import org.apache.commons.math.optimization.VectorialPointValuePair;
24
25
26 /**
27 * This class solves a least squares problem using the Levenberg-Marquardt algorithm.
28 *
29 * <p>This implementation <em>should</em> work even for over-determined systems
30 * (i.e. systems having more point than equations). Over-determined systems
31 * are solved by ignoring the point 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 * @version $Revision: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 2009) $
100 * @since 2.0
101 *
102 */
103 public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
104
105 /** Number of solved point. */
106 private int solvedCols;
107
108 /** Diagonal elements of the R matrix in the Q.R. decomposition. */
109 private double[] diagR;
110
111 /** Norms of the columns of the jacobian matrix. */
112 private double[] jacNorm;
113
114 /** Coefficients of the Householder transforms vectors. */
115 private double[] beta;
116
117 /** Columns permutation array. */
118 private int[] permutation;
119
120 /** Rank of the jacobian matrix. */
121 private int rank;
122
123 /** Levenberg-Marquardt parameter. */
124 private double lmPar;
125
126 /** Parameters evolution direction associated with lmPar. */
127 private double[] lmDir;
128
129 /** Positive input variable used in determining the initial step bound. */
130 private double initialStepBoundFactor;
131
132 /** Desired relative error in the sum of squares. */
133 private double costRelativeTolerance;
134
135 /** Desired relative error in the approximate solution parameters. */
136 private double parRelativeTolerance;
137
138 /** Desired max cosine on the orthogonality between the function vector
139 * and the columns of the jacobian. */
140 private double orthoTolerance;
141
142 /**
143 * Build an optimizer for least squares problems.
144 * <p>The default values for the algorithm settings are:
145 * <ul>
146 * <li>{@link #setInitialStepBoundFactor initial step bound factor}: 100.0</li>
147 * <li>{@link #setMaxIterations maximal iterations}: 1000</li>
148 * <li>{@link #setCostRelativeTolerance cost relative tolerance}: 1.0e-10</li>
149 * <li>{@link #setParRelativeTolerance parameters relative tolerance}: 1.0e-10</li>
150 * <li>{@link #setOrthoTolerance orthogonality tolerance}: 1.0e-10</li>
151 * </ul>
152 * </p>
153 */
154 public LevenbergMarquardtOptimizer() {
155
156 // set up the superclass with a default max cost evaluations setting
157 setMaxIterations(1000);
158
159 // default values for the tuning parameters
160 setInitialStepBoundFactor(100.0);
161 setCostRelativeTolerance(1.0e-10);
162 setParRelativeTolerance(1.0e-10);
163 setOrthoTolerance(1.0e-10);
164
165 }
166
167 /**
168 * Set the positive input variable used in determining the initial step bound.
169 * This bound is set to the product of initialStepBoundFactor and the euclidean
170 * norm of diag*x if nonzero, or else to initialStepBoundFactor itself. In most
171 * cases factor should lie in the interval (0.1, 100.0). 100.0 is a generally
172 * recommended value.
173 *
174 * @param initialStepBoundFactor initial step bound factor
175 */
176 public void setInitialStepBoundFactor(double initialStepBoundFactor) {
177 this.initialStepBoundFactor = initialStepBoundFactor;
178 }
179
180 /**
181 * Set the desired relative error in the sum of squares.
182 *
183 * @param costRelativeTolerance desired relative error in the sum of squares
184 */
185 public void setCostRelativeTolerance(double costRelativeTolerance) {
186 this.costRelativeTolerance = costRelativeTolerance;
187 }
188
189 /**
190 * Set the desired relative error in the approximate solution parameters.
191 *
192 * @param parRelativeTolerance desired relative error
193 * in the approximate solution parameters
194 */
195 public void setParRelativeTolerance(double parRelativeTolerance) {
196 this.parRelativeTolerance = parRelativeTolerance;
197 }
198
199 /**
200 * Set the desired max cosine on the orthogonality.
201 *
202 * @param orthoTolerance desired max cosine on the orthogonality
203 * between the function vector and the columns of the jacobian
204 */
205 public void setOrthoTolerance(double orthoTolerance) {
206 this.orthoTolerance = orthoTolerance;
207 }
208
209 /** {@inheritDoc} */
210 @Override
211 protected VectorialPointValuePair doOptimize()
212 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
213
214 // arrays shared with the other private methods
215 solvedCols = Math.min(rows, cols);
216 diagR = new double[cols];
217 jacNorm = new double[cols];
218 beta = new double[cols];
219 permutation = new int[cols];
220 lmDir = new double[cols];
221
222 // local point
223 double delta = 0;
224 double xNorm = 0;
225 double[] diag = new double[cols];
226 double[] oldX = new double[cols];
227 double[] oldRes = new double[rows];
228 double[] work1 = new double[cols];
229 double[] work2 = new double[cols];
230 double[] work3 = new double[cols];
231
232 // evaluate the function at the starting point and calculate its norm
233 updateResidualsAndCost();
234
235 // outer loop
236 lmPar = 0;
237 boolean firstIteration = true;
238 while (true) {
239
240 incrementIterationsCounter();
241
242 // compute the Q.R. decomposition of the jacobian matrix
243 updateJacobian();
244 qrDecomposition();
245
246 // compute Qt.res
247 qTy(residuals);
248
249 // now we don't need Q anymore,
250 // so let jacobian contain the R matrix with its diagonal elements
251 for (int k = 0; k < solvedCols; ++k) {
252 int pk = permutation[k];
253 jacobian[k][pk] = diagR[pk];
254 }
255
256 if (firstIteration) {
257
258 // scale the point according to the norms of the columns
259 // of the initial jacobian
260 xNorm = 0;
261 for (int k = 0; k < cols; ++k) {
262 double dk = jacNorm[k];
263 if (dk == 0) {
264 dk = 1.0;
265 }
266 double xk = dk * point[k];
267 xNorm += xk * xk;
268 diag[k] = dk;
269 }
270 xNorm = Math.sqrt(xNorm);
271
272 // initialize the step bound delta
273 delta = (xNorm == 0) ? initialStepBoundFactor : (initialStepBoundFactor * xNorm);
274
275 }
276
277 // check orthogonality between function vector and jacobian columns
278 double maxCosine = 0;
279 if (cost != 0) {
280 for (int j = 0; j < solvedCols; ++j) {
281 int pj = permutation[j];
282 double s = jacNorm[pj];
283 if (s != 0) {
284 double sum = 0;
285 for (int i = 0; i <= j; ++i) {
286 sum += jacobian[i][pj] * residuals[i];
287 }
288 maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost));
289 }
290 }
291 }
292 if (maxCosine <= orthoTolerance) {
293 // convergence has been reached
294 return new VectorialPointValuePair(point, objective);
295 }
296
297 // rescale if necessary
298 for (int j = 0; j < cols; ++j) {
299 diag[j] = Math.max(diag[j], jacNorm[j]);
300 }
301
302 // inner loop
303 for (double ratio = 0; ratio < 1.0e-4;) {
304
305 // save the state
306 for (int j = 0; j < solvedCols; ++j) {
307 int pj = permutation[j];
308 oldX[pj] = point[pj];
309 }
310 double previousCost = cost;
311 double[] tmpVec = residuals;
312 residuals = oldRes;
313 oldRes = tmpVec;
314
315 // determine the Levenberg-Marquardt parameter
316 determineLMParameter(oldRes, delta, diag, work1, work2, work3);
317
318 // compute the new point and the norm of the evolution direction
319 double lmNorm = 0;
320 for (int j = 0; j < solvedCols; ++j) {
321 int pj = permutation[j];
322 lmDir[pj] = -lmDir[pj];
323 point[pj] = oldX[pj] + lmDir[pj];
324 double s = diag[pj] * lmDir[pj];
325 lmNorm += s * s;
326 }
327 lmNorm = Math.sqrt(lmNorm);
328
329 // on the first iteration, adjust the initial step bound.
330 if (firstIteration) {
331 delta = Math.min(delta, lmNorm);
332 }
333
334 // evaluate the function at x + p and calculate its norm
335 updateResidualsAndCost();
336
337 // compute the scaled actual reduction
338 double actRed = -1.0;
339 if (0.1 * cost < previousCost) {
340 double r = cost / previousCost;
341 actRed = 1.0 - r * r;
342 }
343
344 // compute the scaled predicted reduction
345 // and the scaled directional derivative
346 for (int j = 0; j < solvedCols; ++j) {
347 int pj = permutation[j];
348 double dirJ = lmDir[pj];
349 work1[j] = 0;
350 for (int i = 0; i <= j; ++i) {
351 work1[i] += jacobian[i][pj] * dirJ;
352 }
353 }
354 double coeff1 = 0;
355 for (int j = 0; j < solvedCols; ++j) {
356 coeff1 += work1[j] * work1[j];
357 }
358 double pc2 = previousCost * previousCost;
359 coeff1 = coeff1 / pc2;
360 double coeff2 = lmPar * lmNorm * lmNorm / pc2;
361 double preRed = coeff1 + 2 * coeff2;
362 double dirDer = -(coeff1 + coeff2);
363
364 // ratio of the actual to the predicted reduction
365 ratio = (preRed == 0) ? 0 : (actRed / preRed);
366
367 // update the step bound
368 if (ratio <= 0.25) {
369 double tmp =
370 (actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5;
371 if ((0.1 * cost >= previousCost) || (tmp < 0.1)) {
372 tmp = 0.1;
373 }
374 delta = tmp * Math.min(delta, 10.0 * lmNorm);
375 lmPar /= tmp;
376 } else if ((lmPar == 0) || (ratio >= 0.75)) {
377 delta = 2 * lmNorm;
378 lmPar *= 0.5;
379 }
380
381 // test for successful iteration.
382 if (ratio >= 1.0e-4) {
383 // successful iteration, update the norm
384 firstIteration = false;
385 xNorm = 0;
386 for (int k = 0; k < cols; ++k) {
387 double xK = diag[k] * point[k];
388 xNorm += xK * xK;
389 }
390 xNorm = Math.sqrt(xNorm);
391 } else {
392 // failed iteration, reset the previous values
393 cost = previousCost;
394 for (int j = 0; j < solvedCols; ++j) {
395 int pj = permutation[j];
396 point[pj] = oldX[pj];
397 }
398 tmpVec = residuals;
399 residuals = oldRes;
400 oldRes = tmpVec;
401 }
402
403 // tests for convergence.
404 if (((Math.abs(actRed) <= costRelativeTolerance) &&
405 (preRed <= costRelativeTolerance) &&
406 (ratio <= 2.0)) ||
407 (delta <= parRelativeTolerance * xNorm)) {
408 return new VectorialPointValuePair(point, objective);
409 }
410
411 // tests for termination and stringent tolerances
412 // (2.2204e-16 is the machine epsilon for IEEE754)
413 if ((Math.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) {
414 throw new OptimizationException("cost relative tolerance is too small ({0})," +
415 " no further reduction in the" +
416 " sum of squares is possible",
417 costRelativeTolerance);
418 } else if (delta <= 2.2204e-16 * xNorm) {
419 throw new OptimizationException("parameters relative tolerance is too small" +
420 " ({0}), no further improvement in" +
421 " the approximate solution is possible",
422 parRelativeTolerance);
423 } else if (maxCosine <= 2.2204e-16) {
424 throw new OptimizationException("orthogonality tolerance is too small ({0})," +
425 " solution is orthogonal to the jacobian",
426 orthoTolerance);
427 }
428
429 }
430
431 }
432
433 }
434
435 /**
436 * Determine the Levenberg-Marquardt parameter.
437 * <p>This implementation is a translation in Java of the MINPACK
438 * <a href="http://www.netlib.org/minpack/lmpar.f">lmpar</a>
439 * routine.</p>
440 * <p>This method sets the lmPar and lmDir attributes.</p>
441 * <p>The authors of the original fortran function are:</p>
442 * <ul>
443 * <li>Argonne National Laboratory. MINPACK project. March 1980</li>
444 * <li>Burton S. Garbow</li>
445 * <li>Kenneth E. Hillstrom</li>
446 * <li>Jorge J. More</li>
447 * </ul>
448 * <p>Luc Maisonobe did the Java translation.</p>
449 *
450 * @param qy array containing qTy
451 * @param delta upper bound on the euclidean norm of diagR * lmDir
452 * @param diag diagonal matrix
453 * @param work1 work array
454 * @param work2 work array
455 * @param work3 work array
456 */
457 private void determineLMParameter(double[] qy, double delta, double[] diag,
458 double[] work1, double[] work2, double[] work3) {
459
460 // compute and store in x the gauss-newton direction, if the
461 // jacobian is rank-deficient, obtain a least squares solution
462 for (int j = 0; j < rank; ++j) {
463 lmDir[permutation[j]] = qy[j];
464 }
465 for (int j = rank; j < cols; ++j) {
466 lmDir[permutation[j]] = 0;
467 }
468 for (int k = rank - 1; k >= 0; --k) {
469 int pk = permutation[k];
470 double ypk = lmDir[pk] / diagR[pk];
471 for (int i = 0; i < k; ++i) {
472 lmDir[permutation[i]] -= ypk * jacobian[i][pk];
473 }
474 lmDir[pk] = ypk;
475 }
476
477 // evaluate the function at the origin, and test
478 // for acceptance of the Gauss-Newton direction
479 double dxNorm = 0;
480 for (int j = 0; j < solvedCols; ++j) {
481 int pj = permutation[j];
482 double s = diag[pj] * lmDir[pj];
483 work1[pj] = s;
484 dxNorm += s * s;
485 }
486 dxNorm = Math.sqrt(dxNorm);
487 double fp = dxNorm - delta;
488 if (fp <= 0.1 * delta) {
489 lmPar = 0;
490 return;
491 }
492
493 // if the jacobian is not rank deficient, the Newton step provides
494 // a lower bound, parl, for the zero of the function,
495 // otherwise set this bound to zero
496 double sum2;
497 double parl = 0;
498 if (rank == solvedCols) {
499 for (int j = 0; j < solvedCols; ++j) {
500 int pj = permutation[j];
501 work1[pj] *= diag[pj] / dxNorm;
502 }
503 sum2 = 0;
504 for (int j = 0; j < solvedCols; ++j) {
505 int pj = permutation[j];
506 double sum = 0;
507 for (int i = 0; i < j; ++i) {
508 sum += jacobian[i][pj] * work1[permutation[i]];
509 }
510 double s = (work1[pj] - sum) / diagR[pj];
511 work1[pj] = s;
512 sum2 += s * s;
513 }
514 parl = fp / (delta * sum2);
515 }
516
517 // calculate an upper bound, paru, for the zero of the function
518 sum2 = 0;
519 for (int j = 0; j < solvedCols; ++j) {
520 int pj = permutation[j];
521 double sum = 0;
522 for (int i = 0; i <= j; ++i) {
523 sum += jacobian[i][pj] * qy[i];
524 }
525 sum /= diag[pj];
526 sum2 += sum * sum;
527 }
528 double gNorm = Math.sqrt(sum2);
529 double paru = gNorm / delta;
530 if (paru == 0) {
531 // 2.2251e-308 is the smallest positive real for IEE754
532 paru = 2.2251e-308 / Math.min(delta, 0.1);
533 }
534
535 // if the input par lies outside of the interval (parl,paru),
536 // set par to the closer endpoint
537 lmPar = Math.min(paru, Math.max(lmPar, parl));
538 if (lmPar == 0) {
539 lmPar = gNorm / dxNorm;
540 }
541
542 for (int countdown = 10; countdown >= 0; --countdown) {
543
544 // evaluate the function at the current value of lmPar
545 if (lmPar == 0) {
546 lmPar = Math.max(2.2251e-308, 0.001 * paru);
547 }
548 double sPar = Math.sqrt(lmPar);
549 for (int j = 0; j < solvedCols; ++j) {
550 int pj = permutation[j];
551 work1[pj] = sPar * diag[pj];
552 }
553 determineLMDirection(qy, work1, work2, work3);
554
555 dxNorm = 0;
556 for (int j = 0; j < solvedCols; ++j) {
557 int pj = permutation[j];
558 double s = diag[pj] * lmDir[pj];
559 work3[pj] = s;
560 dxNorm += s * s;
561 }
562 dxNorm = Math.sqrt(dxNorm);
563 double previousFP = fp;
564 fp = dxNorm - delta;
565
566 // if the function is small enough, accept the current value
567 // of lmPar, also test for the exceptional cases where parl is zero
568 if ((Math.abs(fp) <= 0.1 * delta) ||
569 ((parl == 0) && (fp <= previousFP) && (previousFP < 0))) {
570 return;
571 }
572
573 // compute the Newton correction
574 for (int j = 0; j < solvedCols; ++j) {
575 int pj = permutation[j];
576 work1[pj] = work3[pj] * diag[pj] / dxNorm;
577 }
578 for (int j = 0; j < solvedCols; ++j) {
579 int pj = permutation[j];
580 work1[pj] /= work2[j];
581 double tmp = work1[pj];
582 for (int i = j + 1; i < solvedCols; ++i) {
583 work1[permutation[i]] -= jacobian[i][pj] * tmp;
584 }
585 }
586 sum2 = 0;
587 for (int j = 0; j < solvedCols; ++j) {
588 double s = work1[permutation[j]];
589 sum2 += s * s;
590 }
591 double correction = fp / (delta * sum2);
592
593 // depending on the sign of the function, update parl or paru.
594 if (fp > 0) {
595 parl = Math.max(parl, lmPar);
596 } else if (fp < 0) {
597 paru = Math.min(paru, lmPar);
598 }
599
600 // compute an improved estimate for lmPar
601 lmPar = Math.max(parl, lmPar + correction);
602
603 }
604 }
605
606 /**
607 * Solve a*x = b and d*x = 0 in the least squares sense.
608 * <p>This implementation is a translation in Java of the MINPACK
609 * <a href="http://www.netlib.org/minpack/qrsolv.f">qrsolv</a>
610 * routine.</p>
611 * <p>This method sets the lmDir and lmDiag attributes.</p>
612 * <p>The authors of the original fortran function are:</p>
613 * <ul>
614 * <li>Argonne National Laboratory. MINPACK project. March 1980</li>
615 * <li>Burton S. Garbow</li>
616 * <li>Kenneth E. Hillstrom</li>
617 * <li>Jorge J. More</li>
618 * </ul>
619 * <p>Luc Maisonobe did the Java translation.</p>
620 *
621 * @param qy array containing qTy
622 * @param diag diagonal matrix
623 * @param lmDiag diagonal elements associated with lmDir
624 * @param work work array
625 */
626 private void determineLMDirection(double[] qy, double[] diag,
627 double[] lmDiag, double[] work) {
628
629 // copy R and Qty to preserve input and initialize s
630 // in particular, save the diagonal elements of R in lmDir
631 for (int j = 0; j < solvedCols; ++j) {
632 int pj = permutation[j];
633 for (int i = j + 1; i < solvedCols; ++i) {
634 jacobian[i][pj] = jacobian[j][permutation[i]];
635 }
636 lmDir[j] = diagR[pj];
637 work[j] = qy[j];
638 }
639
640 // eliminate the diagonal matrix d using a Givens rotation
641 for (int j = 0; j < solvedCols; ++j) {
642
643 // prepare the row of d to be eliminated, locating the
644 // diagonal element using p from the Q.R. factorization
645 int pj = permutation[j];
646 double dpj = diag[pj];
647 if (dpj != 0) {
648 Arrays.fill(lmDiag, j + 1, lmDiag.length, 0);
649 }
650 lmDiag[j] = dpj;
651
652 // the transformations to eliminate the row of d
653 // modify only a single element of Qty
654 // beyond the first n, which is initially zero.
655 double qtbpj = 0;
656 for (int k = j; k < solvedCols; ++k) {
657 int pk = permutation[k];
658
659 // determine a Givens rotation which eliminates the
660 // appropriate element in the current row of d
661 if (lmDiag[k] != 0) {
662
663 final double sin;
664 final double cos;
665 double rkk = jacobian[k][pk];
666 if (Math.abs(rkk) < Math.abs(lmDiag[k])) {
667 final double cotan = rkk / lmDiag[k];
668 sin = 1.0 / Math.sqrt(1.0 + cotan * cotan);
669 cos = sin * cotan;
670 } else {
671 final double tan = lmDiag[k] / rkk;
672 cos = 1.0 / Math.sqrt(1.0 + tan * tan);
673 sin = cos * tan;
674 }
675
676 // compute the modified diagonal element of R and
677 // the modified element of (Qty,0)
678 jacobian[k][pk] = cos * rkk + sin * lmDiag[k];
679 final double temp = cos * work[k] + sin * qtbpj;
680 qtbpj = -sin * work[k] + cos * qtbpj;
681 work[k] = temp;
682
683 // accumulate the tranformation in the row of s
684 for (int i = k + 1; i < solvedCols; ++i) {
685 double rik = jacobian[i][pk];
686 final double temp2 = cos * rik + sin * lmDiag[i];
687 lmDiag[i] = -sin * rik + cos * lmDiag[i];
688 jacobian[i][pk] = temp2;
689 }
690
691 }
692 }
693
694 // store the diagonal element of s and restore
695 // the corresponding diagonal element of R
696 lmDiag[j] = jacobian[j][permutation[j]];
697 jacobian[j][permutation[j]] = lmDir[j];
698
699 }
700
701 // solve the triangular system for z, if the system is
702 // singular, then obtain a least squares solution
703 int nSing = solvedCols;
704 for (int j = 0; j < solvedCols; ++j) {
705 if ((lmDiag[j] == 0) && (nSing == solvedCols)) {
706 nSing = j;
707 }
708 if (nSing < solvedCols) {
709 work[j] = 0;
710 }
711 }
712 if (nSing > 0) {
713 for (int j = nSing - 1; j >= 0; --j) {
714 int pj = permutation[j];
715 double sum = 0;
716 for (int i = j + 1; i < nSing; ++i) {
717 sum += jacobian[i][pj] * work[i];
718 }
719 work[j] = (work[j] - sum) / lmDiag[j];
720 }
721 }
722
723 // permute the components of z back to components of lmDir
724 for (int j = 0; j < lmDir.length; ++j) {
725 lmDir[permutation[j]] = work[j];
726 }
727
728 }
729
730 /**
731 * Decompose a matrix A as A.P = Q.R using Householder transforms.
732 * <p>As suggested in the P. Lascaux and R. Theodor book
733 * <i>Analyse numérique matricielle appliquée à
734 * l'art de l'ingénieur</i> (Masson, 1986), instead of representing
735 * the Householder transforms with u<sub>k</sub> unit vectors such that:
736 * <pre>
737 * H<sub>k</sub> = I - 2u<sub>k</sub>.u<sub>k</sub><sup>t</sup>
738 * </pre>
739 * we use <sub>k</sub> non-unit vectors such that:
740 * <pre>
741 * H<sub>k</sub> = I - beta<sub>k</sub>v<sub>k</sub>.v<sub>k</sub><sup>t</sup>
742 * </pre>
743 * where v<sub>k</sub> = a<sub>k</sub> - alpha<sub>k</sub> e<sub>k</sub>.
744 * The beta<sub>k</sub> coefficients are provided upon exit as recomputing
745 * them from the v<sub>k</sub> vectors would be costly.</p>
746 * <p>This decomposition handles rank deficient cases since the tranformations
747 * are performed in non-increasing columns norms order thanks to columns
748 * pivoting. The diagonal elements of the R matrix are therefore also in
749 * non-increasing absolute values order.</p>
750 * @exception OptimizationException if the decomposition cannot be performed
751 */
752 private void qrDecomposition() throws OptimizationException {
753
754 // initializations
755 for (int k = 0; k < cols; ++k) {
756 permutation[k] = k;
757 double norm2 = 0;
758 for (int i = 0; i < jacobian.length; ++i) {
759 double akk = jacobian[i][k];
760 norm2 += akk * akk;
761 }
762 jacNorm[k] = Math.sqrt(norm2);
763 }
764
765 // transform the matrix column after column
766 for (int k = 0; k < cols; ++k) {
767
768 // select the column with the greatest norm on active components
769 int nextColumn = -1;
770 double ak2 = Double.NEGATIVE_INFINITY;
771 for (int i = k; i < cols; ++i) {
772 double norm2 = 0;
773 for (int j = k; j < jacobian.length; ++j) {
774 double aki = jacobian[j][permutation[i]];
775 norm2 += aki * aki;
776 }
777 if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
778 throw new OptimizationException(
779 "unable to perform Q.R decomposition on the {0}x{1} jacobian matrix",
780 rows, cols);
781 }
782 if (norm2 > ak2) {
783 nextColumn = i;
784 ak2 = norm2;
785 }
786 }
787 if (ak2 == 0) {
788 rank = k;
789 return;
790 }
791 int pk = permutation[nextColumn];
792 permutation[nextColumn] = permutation[k];
793 permutation[k] = pk;
794
795 // choose alpha such that Hk.u = alpha ek
796 double akk = jacobian[k][pk];
797 double alpha = (akk > 0) ? -Math.sqrt(ak2) : Math.sqrt(ak2);
798 double betak = 1.0 / (ak2 - akk * alpha);
799 beta[pk] = betak;
800
801 // transform the current column
802 diagR[pk] = alpha;
803 jacobian[k][pk] -= alpha;
804
805 // transform the remaining columns
806 for (int dk = cols - 1 - k; dk > 0; --dk) {
807 double gamma = 0;
808 for (int j = k; j < jacobian.length; ++j) {
809 gamma += jacobian[j][pk] * jacobian[j][permutation[k + dk]];
810 }
811 gamma *= betak;
812 for (int j = k; j < jacobian.length; ++j) {
813 jacobian[j][permutation[k + dk]] -= gamma * jacobian[j][pk];
814 }
815 }
816
817 }
818
819 rank = solvedCols;
820
821 }
822
823 /**
824 * Compute the product Qt.y for some Q.R. decomposition.
825 *
826 * @param y vector to multiply (will be overwritten with the result)
827 */
828 private void qTy(double[] y) {
829 for (int k = 0; k < cols; ++k) {
830 int pk = permutation[k];
831 double gamma = 0;
832 for (int i = k; i < rows; ++i) {
833 gamma += jacobian[i][pk] * y[i];
834 }
835 gamma *= beta[pk];
836 for (int i = k; i < rows; ++i) {
837 y[i] -= gamma * jacobian[i][pk];
838 }
839 }
840 }
841
842 }