View Javadoc
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    *      https://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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.optim.nonlinear.vector.leastsquares;
23  
24  import org.hipparchus.analysis.MultivariateMatrixFunction;
25  import org.hipparchus.analysis.MultivariateVectorFunction;
26  import org.hipparchus.exception.MathIllegalStateException;
27  import org.hipparchus.linear.Array2DRowRealMatrix;
28  import org.hipparchus.linear.ArrayRealVector;
29  import org.hipparchus.linear.DiagonalMatrix;
30  import org.hipparchus.linear.EigenDecompositionSymmetric;
31  import org.hipparchus.linear.RealMatrix;
32  import org.hipparchus.linear.RealVector;
33  import org.hipparchus.optim.AbstractOptimizationProblem;
34  import org.hipparchus.optim.ConvergenceChecker;
35  import org.hipparchus.optim.LocalizedOptimFormats;
36  import org.hipparchus.optim.PointVectorValuePair;
37  import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
38  import org.hipparchus.util.FastMath;
39  import org.hipparchus.util.Incrementor;
40  import org.hipparchus.util.Pair;
41  
42  /**
43   * A Factory for creating {@link LeastSquaresProblem}s.
44   *
45   */
46  public class LeastSquaresFactory {
47  
48      /** Prevent instantiation. */
49      private LeastSquaresFactory() {}
50  
51      /**
52       * Create a {@link org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem}
53       * from the given elements. There will be no weights applied (unit weights).
54       *
55       * @param model          the model function. Produces the computed values.
56       * @param observed       the observed (target) values
57       * @param start          the initial guess.
58       * @param weight         the weight matrix
59       * @param checker        convergence checker
60       * @param maxEvaluations the maximum number of times to evaluate the model
61       * @param maxIterations  the maximum number to times to iterate in the algorithm
62       * @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)}
63       * will defer the evaluation until access to the value is requested.
64       * @param paramValidator Model parameters validator.
65       * @return the specified General Least Squares problem.
66       *
67       */
68      public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
69                                               final RealVector observed,
70                                               final RealVector start,
71                                               final RealMatrix weight,
72                                               final ConvergenceChecker<Evaluation> checker,
73                                               final int maxEvaluations,
74                                               final int maxIterations,
75                                               final boolean lazyEvaluation,
76                                               final ParameterValidator paramValidator) {
77          final LeastSquaresProblem p = new LocalLeastSquaresProblem(model,
78                                                                     observed,
79                                                                     start,
80                                                                     checker,
81                                                                     maxEvaluations,
82                                                                     maxIterations,
83                                                                     lazyEvaluation,
84                                                                     paramValidator);
85          if (weight != null) {
86              return weightMatrix(p, weight);
87          } else {
88              return p;
89          }
90      }
91  
92      /**
93       * Create a {@link org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem}
94       * from the given elements. There will be no weights applied (unit weights).
95       *
96       * @param model          the model function. Produces the computed values.
97       * @param observed       the observed (target) values
98       * @param start          the initial guess.
99       * @param checker        convergence checker
100      * @param maxEvaluations the maximum number of times to evaluate the model
101      * @param maxIterations  the maximum number to times to iterate in the algorithm
102      * @return the specified General Least Squares problem.
103      */
104     public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
105                                              final RealVector observed,
106                                              final RealVector start,
107                                              final ConvergenceChecker<Evaluation> checker,
108                                              final int maxEvaluations,
109                                              final int maxIterations) {
110         return create(model,
111                       observed,
112                       start,
113                       null,
114                       checker,
115                       maxEvaluations,
116                       maxIterations,
117                       false,
118                       null);
119     }
120 
121     /**
122      * Create a {@link org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem}
123      * from the given elements.
124      *
125      * @param model          the model function. Produces the computed values.
126      * @param observed       the observed (target) values
127      * @param start          the initial guess.
128      * @param weight         the weight matrix
129      * @param checker        convergence checker
130      * @param maxEvaluations the maximum number of times to evaluate the model
131      * @param maxIterations  the maximum number to times to iterate in the algorithm
132      * @return the specified General Least Squares problem.
133      */
134     public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
135                                              final RealVector observed,
136                                              final RealVector start,
137                                              final RealMatrix weight,
138                                              final ConvergenceChecker<Evaluation> checker,
139                                              final int maxEvaluations,
140                                              final int maxIterations) {
141         return weightMatrix(create(model,
142                                    observed,
143                                    start,
144                                    checker,
145                                    maxEvaluations,
146                                    maxIterations),
147                             weight);
148     }
149 
150     /**
151      * Create a {@link org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem}
152      * from the given elements.
153      * <p>
154      * This factory method is provided for continuity with previous interfaces. Newer
155      * applications should use {@link #create(MultivariateJacobianFunction, RealVector,
156      * RealVector, ConvergenceChecker, int, int)}, or {@link #create(MultivariateJacobianFunction,
157      * RealVector, RealVector, RealMatrix, ConvergenceChecker, int, int)}.
158      *
159      * @param model          the model function. Produces the computed values.
160      * @param jacobian       the jacobian of the model with respect to the parameters
161      * @param observed       the observed (target) values
162      * @param start          the initial guess.
163      * @param weight         the weight matrix
164      * @param checker        convergence checker
165      * @param maxEvaluations the maximum number of times to evaluate the model
166      * @param maxIterations  the maximum number to times to iterate in the algorithm
167      * @return the specified General Least Squares problem.
168      */
169     public static LeastSquaresProblem create(final MultivariateVectorFunction model,
170                                              final MultivariateMatrixFunction jacobian,
171                                              final double[] observed,
172                                              final double[] start,
173                                              final RealMatrix weight,
174                                              final ConvergenceChecker<Evaluation> checker,
175                                              final int maxEvaluations,
176                                              final int maxIterations) {
177         return create(model(model, jacobian),
178                       new ArrayRealVector(observed, false),
179                       new ArrayRealVector(start, false),
180                       weight,
181                       checker,
182                       maxEvaluations,
183                       maxIterations);
184     }
185 
186     /**
187      * Apply a dense weight matrix to the {@link LeastSquaresProblem}.
188      *
189      * @param problem the unweighted problem
190      * @param weights the matrix of weights
191      * @return a new {@link LeastSquaresProblem} with the weights applied. The original
192      *         {@code problem} is not modified.
193      */
194     public static LeastSquaresProblem weightMatrix(final LeastSquaresProblem problem,
195                                                    final RealMatrix weights) {
196         final RealMatrix weightSquareRoot = squareRoot(weights);
197         return new LeastSquaresAdapter(problem) {
198             /** {@inheritDoc} */
199             @Override
200             public Evaluation evaluate(final RealVector point) {
201                 return new DenseWeightedEvaluation(super.evaluate(point), weightSquareRoot);
202             }
203         };
204     }
205 
206     /**
207      * Apply a diagonal weight matrix to the {@link LeastSquaresProblem}.
208      *
209      * @param problem the unweighted problem
210      * @param weights the diagonal of the weight matrix
211      * @return a new {@link LeastSquaresProblem} with the weights applied. The original
212      *         {@code problem} is not modified.
213      */
214     public static LeastSquaresProblem weightDiagonal(final LeastSquaresProblem problem,
215                                                      final RealVector weights) {
216         // TODO more efficient implementation
217         return weightMatrix(problem, new DiagonalMatrix(weights.toArray()));
218     }
219 
220     /**
221      * Count the evaluations of a particular problem. The {@code counter} will be
222      * incremented every time {@link LeastSquaresProblem#evaluate(RealVector)} is called on
223      * the <em>returned</em> problem.
224      *
225      * @param problem the problem to track.
226      * @param counter the counter to increment.
227      * @return a least squares problem that tracks evaluations
228      */
229     public static LeastSquaresProblem countEvaluations(final LeastSquaresProblem problem,
230                                                        final Incrementor counter) {
231         return new LeastSquaresAdapter(problem) {
232 
233             /** {@inheritDoc} */
234             @Override
235             public Evaluation evaluate(final RealVector point) {
236                 counter.increment();
237                 return super.evaluate(point);
238             }
239 
240             // Delegate the rest.
241         };
242     }
243 
244     /**
245      * View a convergence checker specified for a {@link PointVectorValuePair} as one
246      * specified for an {@link Evaluation}.
247      *
248      * @param checker the convergence checker to adapt.
249      * @return a convergence checker that delegates to {@code checker}.
250      */
251     public static ConvergenceChecker<Evaluation> evaluationChecker(final ConvergenceChecker<PointVectorValuePair> checker) {
252         return new ConvergenceChecker<Evaluation>() {
253             /** {@inheritDoc} */
254             @Override
255             public boolean converged(final int iteration,
256                                      final Evaluation previous,
257                                      final Evaluation current) {
258                 return checker.converged(
259                         iteration,
260                         new PointVectorValuePair(
261                                 previous.getPoint().toArray(),
262                                 previous.getResiduals().toArray(),
263                                 false),
264                         new PointVectorValuePair(
265                                 current.getPoint().toArray(),
266                                 current.getResiduals().toArray(),
267                                 false)
268                 );
269             }
270         };
271     }
272 
273     /**
274      * Computes the square-root of the weight matrix.
275      *
276      * @param m Symmetric, positive-definite (weight) matrix.
277      * @return the square-root of the weight matrix.
278      */
279     private static RealMatrix squareRoot(final RealMatrix m) {
280         if (m instanceof DiagonalMatrix) {
281             final int dim = m.getRowDimension();
282             final RealMatrix sqrtM = new DiagonalMatrix(dim);
283             for (int i = 0; i < dim; i++) {
284                 sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
285             }
286             return sqrtM;
287         } else {
288             final EigenDecompositionSymmetric dec = new EigenDecompositionSymmetric(m);
289             return dec.getSquareRoot();
290         }
291     }
292 
293     /**
294      * Combine a {@link MultivariateVectorFunction} with a {@link
295      * MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}.
296      *
297      * @param value    the vector value function
298      * @param jacobian the Jacobian function
299      * @return a function that computes both at the same time
300      */
301     public static MultivariateJacobianFunction model(final MultivariateVectorFunction value,
302                                                      final MultivariateMatrixFunction jacobian) {
303         return new LocalValueAndJacobianFunction(value, jacobian);
304     }
305 
306     /**
307      * Combine a {@link MultivariateVectorFunction} with a {@link
308      * MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}.
309      */
310     private static class LocalValueAndJacobianFunction
311         implements ValueAndJacobianFunction {
312         /** Model. */
313         private final MultivariateVectorFunction value;
314         /** Model's Jacobian. */
315         private final MultivariateMatrixFunction jacobian;
316 
317         /**
318          * @param value Model function.
319          * @param jacobian Model's Jacobian function.
320          */
321         LocalValueAndJacobianFunction(final MultivariateVectorFunction value,
322                                       final MultivariateMatrixFunction jacobian) {
323             this.value = value;
324             this.jacobian = jacobian;
325         }
326 
327         /** {@inheritDoc} */
328         @Override
329         public Pair<RealVector, RealMatrix> value(final RealVector point) {
330             //TODO get array from RealVector without copying?
331             final double[] p = point.toArray();
332 
333             // Evaluate.
334             return new Pair<RealVector, RealMatrix>(computeValue(p),
335                                                     computeJacobian(p));
336         }
337 
338         /** {@inheritDoc} */
339         @Override
340         public RealVector computeValue(final double[] params) {
341             return new ArrayRealVector(value.value(params), false);
342         }
343 
344         /** {@inheritDoc} */
345         @Override
346         public RealMatrix computeJacobian(final double[] params) {
347             return new Array2DRowRealMatrix(jacobian.value(params), false);
348         }
349     }
350 
351 
352     /**
353      * A private, "field" immutable (not "real" immutable) implementation of {@link
354      * LeastSquaresProblem}.
355      */
356     private static class LocalLeastSquaresProblem
357             extends AbstractOptimizationProblem<Evaluation>
358             implements LeastSquaresProblem {
359 
360         /** Target values for the model function at optimum. */
361         private final RealVector target;
362         /** Model function. */
363         private final MultivariateJacobianFunction model;
364         /** Initial guess. */
365         private final RealVector start;
366         /** Whether to use lazy evaluation. */
367         private final boolean lazyEvaluation;
368         /** Model parameters validator. */
369         private final ParameterValidator paramValidator;
370 
371         /**
372          * Create a {@link LeastSquaresProblem} from the given data.
373          *
374          * @param model          the model function
375          * @param target         the observed data
376          * @param start          the initial guess
377          * @param checker        the convergence checker
378          * @param maxEvaluations the allowed evaluations
379          * @param maxIterations  the allowed iterations
380          * @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)}
381          * will defer the evaluation until access to the value is requested.
382          * @param paramValidator Model parameters validator.
383          */
384         LocalLeastSquaresProblem(final MultivariateJacobianFunction model,
385                                  final RealVector target,
386                                  final RealVector start,
387                                  final ConvergenceChecker<Evaluation> checker,
388                                  final int maxEvaluations,
389                                  final int maxIterations,
390                                  final boolean lazyEvaluation,
391                                  final ParameterValidator paramValidator) {
392             super(maxEvaluations, maxIterations, checker);
393             this.target = target;
394             this.model = model;
395             this.start = start;
396             this.lazyEvaluation = lazyEvaluation;
397             this.paramValidator = paramValidator;
398 
399             if (lazyEvaluation &&
400                 !(model instanceof ValueAndJacobianFunction)) {
401                 // Lazy evaluation requires that value and Jacobian
402                 // can be computed separately.
403                 throw new MathIllegalStateException(LocalizedOptimFormats.INVALID_IMPLEMENTATION,
404                                                     model.getClass().getName());
405             }
406         }
407 
408         /** {@inheritDoc} */
409         @Override
410         public int getObservationSize() {
411             return target.getDimension();
412         }
413 
414         /** {@inheritDoc} */
415         @Override
416         public int getParameterSize() {
417             return start.getDimension();
418         }
419 
420         /** {@inheritDoc} */
421         @Override
422         public RealVector getStart() {
423             return start == null ? null : start.copy();
424         }
425 
426         /** {@inheritDoc} */
427         @Override
428         public Evaluation evaluate(final RealVector point) {
429             // Copy so optimizer can change point without changing our instance.
430             final RealVector p = paramValidator == null ?
431                 point.copy() :
432                 paramValidator.validate(point.copy());
433 
434             if (lazyEvaluation) {
435                 return new LazyUnweightedEvaluation((ValueAndJacobianFunction) model,
436                                                     target,
437                                                     p);
438             } else {
439                 // Evaluate value and jacobian in one function call.
440                 final Pair<RealVector, RealMatrix> value = model.value(p);
441                 return new UnweightedEvaluation(value.getFirst(),
442                                                 value.getSecond(),
443                                                 target,
444                                                 p);
445             }
446         }
447 
448         /**
449          * Container with the model evaluation at a particular point.
450          */
451         private static class UnweightedEvaluation extends AbstractEvaluation {
452             /** Point of evaluation. */
453             private final RealVector point;
454             /** Derivative at point. */
455             private final RealMatrix jacobian;
456             /** Computed residuals. */
457             private final RealVector residuals;
458 
459             /**
460              * Create an {@link Evaluation} with no weights.
461              *
462              * @param values   the computed function values
463              * @param jacobian the computed function Jacobian
464              * @param target   the observed values
465              * @param point    the abscissa
466              */
467             private UnweightedEvaluation(final RealVector values,
468                                          final RealMatrix jacobian,
469                                          final RealVector target,
470                                          final RealVector point) {
471                 super(target.getDimension());
472                 this.jacobian = jacobian;
473                 this.point = point;
474                 this.residuals = target.subtract(values);
475             }
476 
477             /** {@inheritDoc} */
478             @Override
479             public RealMatrix getJacobian() {
480                 return jacobian;
481             }
482 
483             /** {@inheritDoc} */
484             @Override
485             public RealVector getPoint() {
486                 return point;
487             }
488 
489             /** {@inheritDoc} */
490             @Override
491             public RealVector getResiduals() {
492                 return residuals;
493             }
494         }
495 
496         /**
497          * Container with the model <em>lazy</em> evaluation at a particular point.
498          */
499         private static class LazyUnweightedEvaluation extends AbstractEvaluation {
500             /** Point of evaluation. */
501             private final RealVector point;
502             /** Model and Jacobian functions. */
503             private final ValueAndJacobianFunction model;
504             /** Target values for the model function at optimum. */
505             private final RealVector target;
506 
507             /**
508              * Create an {@link Evaluation} with no weights.
509              *
510              * @param model  the model function
511              * @param target the observed values
512              * @param point  the abscissa
513              */
514             private LazyUnweightedEvaluation(final ValueAndJacobianFunction model,
515                                              final RealVector target,
516                                              final RealVector point) {
517                 super(target.getDimension());
518                 // Safe to cast as long as we control usage of this class.
519                 this.model = model;
520                 this.point = point;
521                 this.target = target;
522             }
523 
524             /** {@inheritDoc} */
525             @Override
526             public RealMatrix getJacobian() {
527                 return model.computeJacobian(point.toArray());
528             }
529 
530             /** {@inheritDoc} */
531             @Override
532             public RealVector getPoint() {
533                 return point;
534             }
535 
536             /** {@inheritDoc} */
537             @Override
538             public RealVector getResiduals() {
539                 return target.subtract(model.computeValue(point.toArray()));
540             }
541         }
542     }
543 }
544