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&eacute;rique matricielle
037     * appliqu&eacute;e &agrave; l'art de l'ing&eacute;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&eacute;rique matricielle appliqu&eacute;e &agrave;
723       * l'art de l'ing&eacute;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    }