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 */ 017package org.apache.commons.math4.legacy.fitting.leastsquares; 018 019import org.apache.commons.math4.legacy.analysis.MultivariateMatrixFunction; 020import org.apache.commons.math4.legacy.analysis.MultivariateVectorFunction; 021import org.apache.commons.math4.legacy.exception.MathIllegalStateException; 022import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; 023import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem.Evaluation; 024import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix; 025import org.apache.commons.math4.legacy.linear.ArrayRealVector; 026import org.apache.commons.math4.legacy.linear.DiagonalMatrix; 027import org.apache.commons.math4.legacy.linear.EigenDecomposition; 028import org.apache.commons.math4.legacy.linear.RealMatrix; 029import org.apache.commons.math4.legacy.linear.RealVector; 030import org.apache.commons.math4.legacy.optim.AbstractOptimizationProblem; 031import org.apache.commons.math4.legacy.optim.ConvergenceChecker; 032import org.apache.commons.math4.legacy.optim.PointVectorValuePair; 033import org.apache.commons.math4.core.jdkmath.JdkMath; 034import org.apache.commons.math4.legacy.core.IntegerSequence; 035import org.apache.commons.math4.legacy.core.Pair; 036 037/** 038 * A Factory for creating {@link LeastSquaresProblem}s. 039 * 040 * @since 3.3 041 */ 042public final class LeastSquaresFactory { 043 044 /** Prevent instantiation. */ 045 private LeastSquaresFactory() {} 046 047 /** 048 * Create a {@link org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem} 049 * from the given elements. 050 * 051 * @param model the model function. Produces the computed values. 052 * @param observed the observed (target) values 053 * @param start the initial guess. 054 * @param weight the weight matrix 055 * @param checker convergence checker 056 * @param maxEvaluations the maximum number of times to evaluate the model 057 * @param maxIterations the maximum number to times to iterate in the algorithm 058 * @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)} 059 * will defer the evaluation until access to the value is requested. 060 * @param paramValidator Model parameters validator. 061 * @return the specified General Least Squares problem. 062 * 063 * @since 3.4 064 */ 065 public static LeastSquaresProblem create(final MultivariateJacobianFunction model, 066 final RealVector observed, 067 final RealVector start, 068 final RealMatrix weight, 069 final ConvergenceChecker<Evaluation> checker, 070 final int maxEvaluations, 071 final int maxIterations, 072 final boolean lazyEvaluation, 073 final ParameterValidator paramValidator) { 074 final LeastSquaresProblem p = new LocalLeastSquaresProblem(model, 075 observed, 076 start, 077 checker, 078 maxEvaluations, 079 maxIterations, 080 lazyEvaluation, 081 paramValidator); 082 if (weight != null) { 083 return weightMatrix(p, weight); 084 } else { 085 return p; 086 } 087 } 088 089 /** 090 * Create a {@link org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem} 091 * from the given elements. There will be no weights applied (unit weights). 092 * 093 * @param model the model function. Produces the computed values. 094 * @param observed the observed (target) values 095 * @param start the initial guess. 096 * @param checker convergence checker 097 * @param maxEvaluations the maximum number of times to evaluate the model 098 * @param maxIterations the maximum number to times to iterate in the algorithm 099 * @return the specified General Least Squares problem. 100 */ 101 public static LeastSquaresProblem create(final MultivariateJacobianFunction model, 102 final RealVector observed, 103 final RealVector start, 104 final ConvergenceChecker<Evaluation> checker, 105 final int maxEvaluations, 106 final int maxIterations) { 107 return create(model, 108 observed, 109 start, 110 null, 111 checker, 112 maxEvaluations, 113 maxIterations, 114 false, 115 null); 116 } 117 118 /** 119 * Create a {@link org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem} 120 * from the given elements. 121 * 122 * @param model the model function. Produces the computed values. 123 * @param observed the observed (target) values 124 * @param start the initial guess. 125 * @param weight the weight matrix 126 * @param checker convergence checker 127 * @param maxEvaluations the maximum number of times to evaluate the model 128 * @param maxIterations the maximum number to times to iterate in the algorithm 129 * @return the specified General Least Squares problem. 130 */ 131 public static LeastSquaresProblem create(final MultivariateJacobianFunction model, 132 final RealVector observed, 133 final RealVector start, 134 final RealMatrix weight, 135 final ConvergenceChecker<Evaluation> checker, 136 final int maxEvaluations, 137 final int maxIterations) { 138 return weightMatrix(create(model, 139 observed, 140 start, 141 checker, 142 maxEvaluations, 143 maxIterations), 144 weight); 145 } 146 147 /** 148 * Create a {@link org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem} 149 * from the given elements. 150 * <p> 151 * This factory method is provided for continuity with previous interfaces. Newer 152 * applications should use {@link #create(MultivariateJacobianFunction, RealVector, 153 * RealVector, ConvergenceChecker, int, int)}, or {@link #create(MultivariateJacobianFunction, 154 * RealVector, RealVector, RealMatrix, ConvergenceChecker, int, int)}. 155 * 156 * @param model the model function. Produces the computed values. 157 * @param jacobian the jacobian of the model with respect to the parameters 158 * @param observed the observed (target) values 159 * @param start the initial guess. 160 * @param weight the weight matrix 161 * @param checker convergence checker 162 * @param maxEvaluations the maximum number of times to evaluate the model 163 * @param maxIterations the maximum number to times to iterate in the algorithm 164 * @return the specified General Least Squares problem. 165 */ 166 public static LeastSquaresProblem create(final MultivariateVectorFunction model, 167 final MultivariateMatrixFunction jacobian, 168 final double[] observed, 169 final double[] start, 170 final RealMatrix weight, 171 final ConvergenceChecker<Evaluation> checker, 172 final int maxEvaluations, 173 final int maxIterations) { 174 return create(model(model, jacobian), 175 new ArrayRealVector(observed, false), 176 new ArrayRealVector(start, false), 177 weight, 178 checker, 179 maxEvaluations, 180 maxIterations); 181 } 182 183 /** 184 * Apply a dense weight matrix to the {@link LeastSquaresProblem}. 185 * 186 * @param problem the unweighted problem 187 * @param weights the matrix of weights 188 * @return a new {@link LeastSquaresProblem} with the weights applied. The original 189 * {@code problem} is not modified. 190 */ 191 public static LeastSquaresProblem weightMatrix(final LeastSquaresProblem problem, 192 final RealMatrix weights) { 193 final RealMatrix weightSquareRoot = squareRoot(weights); 194 return new LeastSquaresAdapter(problem) { 195 /** {@inheritDoc} */ 196 @Override 197 public Evaluation evaluate(final RealVector point) { 198 return new DenseWeightedEvaluation(super.evaluate(point), weightSquareRoot); 199 } 200 }; 201 } 202 203 /** 204 * Apply a diagonal weight matrix to the {@link LeastSquaresProblem}. 205 * 206 * @param problem the unweighted problem 207 * @param weights the diagonal of the weight matrix 208 * @return a new {@link LeastSquaresProblem} with the weights applied. The original 209 * {@code problem} is not modified. 210 */ 211 public static LeastSquaresProblem weightDiagonal(final LeastSquaresProblem problem, 212 final RealVector weights) { 213 // TODO more efficient implementation 214 return weightMatrix(problem, new DiagonalMatrix(weights.toArray())); 215 } 216 217 /** 218 * Count the evaluations of a particular problem. The {@code counter} will be 219 * incremented every time {@link LeastSquaresProblem#evaluate(RealVector)} is called on 220 * the <em>returned</em> problem. 221 * 222 * @param problem the problem to track. 223 * @param counter the counter to increment. 224 * @return a least squares problem that tracks evaluations 225 */ 226 public static LeastSquaresProblem countEvaluations(final LeastSquaresProblem problem, 227 final IntegerSequence.Incrementor counter) { 228 return new LeastSquaresAdapter(problem) { 229 230 /** {@inheritDoc} */ 231 @Override 232 public Evaluation evaluate(final RealVector point) { 233 counter.increment(); 234 return super.evaluate(point); 235 } 236 237 // Delegate the rest. 238 }; 239 } 240 241 /** 242 * View a convergence checker specified for a {@link PointVectorValuePair} as one 243 * specified for an {@link Evaluation}. 244 * 245 * @param checker the convergence checker to adapt. 246 * @return a convergence checker that delegates to {@code checker}. 247 */ 248 public static ConvergenceChecker<Evaluation> evaluationChecker(final ConvergenceChecker<PointVectorValuePair> checker) { 249 return new ConvergenceChecker<Evaluation>() { 250 /** {@inheritDoc} */ 251 @Override 252 public boolean converged(final int iteration, 253 final Evaluation previous, 254 final Evaluation current) { 255 return checker.converged( 256 iteration, 257 new PointVectorValuePair( 258 previous.getPoint().toArray(), 259 previous.getResiduals().toArray(), 260 false), 261 new PointVectorValuePair( 262 current.getPoint().toArray(), 263 current.getResiduals().toArray(), 264 false) 265 ); 266 } 267 }; 268 } 269 270 /** 271 * Computes the square-root of the weight matrix. 272 * 273 * @param m Symmetric, positive-definite (weight) matrix. 274 * @return the square-root of the weight matrix. 275 */ 276 private static RealMatrix squareRoot(final RealMatrix m) { 277 if (m instanceof DiagonalMatrix) { 278 final int dim = m.getRowDimension(); 279 final RealMatrix sqrtM = new DiagonalMatrix(dim); 280 for (int i = 0; i < dim; i++) { 281 sqrtM.setEntry(i, i, JdkMath.sqrt(m.getEntry(i, i))); 282 } 283 return sqrtM; 284 } else { 285 final EigenDecomposition dec = new EigenDecomposition(m); 286 return dec.getSquareRoot(); 287 } 288 } 289 290 /** 291 * Combine a {@link MultivariateVectorFunction} with a {@link 292 * MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}. 293 * 294 * @param value the vector value function 295 * @param jacobian the Jacobian function 296 * @return a function that computes both at the same time 297 */ 298 public static MultivariateJacobianFunction model(final MultivariateVectorFunction value, 299 final MultivariateMatrixFunction jacobian) { 300 return new LocalValueAndJacobianFunction(value, jacobian); 301 } 302 303 /** 304 * Combine a {@link MultivariateVectorFunction} with a {@link 305 * MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}. 306 */ 307 private static class LocalValueAndJacobianFunction 308 implements ValueAndJacobianFunction { 309 /** Model. */ 310 private final MultivariateVectorFunction value; 311 /** Model's Jacobian. */ 312 private final MultivariateMatrixFunction jacobian; 313 314 /** 315 * @param value Model function. 316 * @param jacobian Model's Jacobian function. 317 */ 318 LocalValueAndJacobianFunction(final MultivariateVectorFunction value, 319 final MultivariateMatrixFunction jacobian) { 320 this.value = value; 321 this.jacobian = jacobian; 322 } 323 324 /** {@inheritDoc} */ 325 @Override 326 public Pair<RealVector, RealMatrix> value(final RealVector point) { 327 //TODO get array from RealVector without copying? 328 final double[] p = point.toArray(); 329 330 // Evaluate. 331 return new Pair<>(computeValue(p), computeJacobian(p)); 332 } 333 334 /** {@inheritDoc} */ 335 @Override 336 public RealVector computeValue(final double[] params) { 337 return new ArrayRealVector(value.value(params), false); 338 } 339 340 /** {@inheritDoc} */ 341 @Override 342 public RealMatrix computeJacobian(final double[] params) { 343 return new Array2DRowRealMatrix(jacobian.value(params), false); 344 } 345 } 346 347 348 /** 349 * A private, "field" immutable (not "real" immutable) implementation of {@link 350 * LeastSquaresProblem}. 351 * @since 3.3 352 */ 353 private static class LocalLeastSquaresProblem 354 extends AbstractOptimizationProblem<Evaluation> 355 implements LeastSquaresProblem { 356 357 /** Target values for the model function at optimum. */ 358 private final RealVector target; 359 /** Model function. */ 360 private final MultivariateJacobianFunction model; 361 /** Initial guess. */ 362 private final RealVector start; 363 /** Whether to use lazy evaluation. */ 364 private final boolean lazyEvaluation; 365 /** Model parameters validator. */ 366 private final ParameterValidator paramValidator; 367 368 /** 369 * Create a {@link LeastSquaresProblem} from the given data. 370 * 371 * @param model the model function 372 * @param target the observed data 373 * @param start the initial guess 374 * @param checker the convergence checker 375 * @param maxEvaluations the allowed evaluations 376 * @param maxIterations the allowed iterations 377 * @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)} 378 * will defer the evaluation until access to the value is requested. 379 * @param paramValidator Model parameters validator. 380 */ 381 LocalLeastSquaresProblem(final MultivariateJacobianFunction model, 382 final RealVector target, 383 final RealVector start, 384 final ConvergenceChecker<Evaluation> checker, 385 final int maxEvaluations, 386 final int maxIterations, 387 final boolean lazyEvaluation, 388 final ParameterValidator paramValidator) { 389 super(maxEvaluations, maxIterations, checker); 390 this.target = target; 391 this.model = model; 392 this.start = start; 393 this.lazyEvaluation = lazyEvaluation; 394 this.paramValidator = paramValidator; 395 396 if (lazyEvaluation && 397 !(model instanceof ValueAndJacobianFunction)) { 398 // Lazy evaluation requires that value and Jacobian 399 // can be computed separately. 400 throw new MathIllegalStateException(LocalizedFormats.INVALID_IMPLEMENTATION, 401 model.getClass().getName()); 402 } 403 } 404 405 /** {@inheritDoc} */ 406 @Override 407 public int getObservationSize() { 408 return target.getDimension(); 409 } 410 411 /** {@inheritDoc} */ 412 @Override 413 public int getParameterSize() { 414 return start.getDimension(); 415 } 416 417 /** {@inheritDoc} */ 418 @Override 419 public RealVector getStart() { 420 return start == null ? null : start.copy(); 421 } 422 423 /** {@inheritDoc} */ 424 @Override 425 public Evaluation evaluate(final RealVector point) { 426 // Copy so optimizer can change point without changing our instance. 427 final RealVector p = paramValidator == null ? 428 point.copy() : 429 paramValidator.validate(point.copy()); 430 431 if (lazyEvaluation) { 432 return new LazyUnweightedEvaluation((ValueAndJacobianFunction) model, 433 target, 434 p); 435 } else { 436 // Evaluate value and jacobian in one function call. 437 final Pair<RealVector, RealMatrix> value = model.value(p); 438 return new UnweightedEvaluation(value.getFirst(), 439 value.getSecond(), 440 target, 441 p); 442 } 443 } 444 445 /** 446 * Container with the model evaluation at a particular point. 447 */ 448 private static final class UnweightedEvaluation extends AbstractEvaluation { 449 /** Point of evaluation. */ 450 private final RealVector point; 451 /** Derivative at point. */ 452 private final RealMatrix jacobian; 453 /** Computed residuals. */ 454 private final RealVector residuals; 455 456 /** 457 * Create an {@link Evaluation} with no weights. 458 * 459 * @param values the computed function values 460 * @param jacobian the computed function Jacobian 461 * @param target the observed values 462 * @param point the abscissa 463 */ 464 private UnweightedEvaluation(final RealVector values, 465 final RealMatrix jacobian, 466 final RealVector target, 467 final RealVector point) { 468 super(target.getDimension()); 469 this.jacobian = jacobian; 470 this.point = point; 471 this.residuals = target.subtract(values); 472 } 473 474 /** {@inheritDoc} */ 475 @Override 476 public RealMatrix getJacobian() { 477 return jacobian; 478 } 479 480 /** {@inheritDoc} */ 481 @Override 482 public RealVector getPoint() { 483 return point; 484 } 485 486 /** {@inheritDoc} */ 487 @Override 488 public RealVector getResiduals() { 489 return residuals; 490 } 491 } 492 493 /** 494 * Container with the model <em>lazy</em> evaluation at a particular point. 495 */ 496 private static final class LazyUnweightedEvaluation extends AbstractEvaluation { 497 /** Point of evaluation. */ 498 private final RealVector point; 499 /** Model and Jacobian functions. */ 500 private final ValueAndJacobianFunction model; 501 /** Target values for the model function at optimum. */ 502 private final RealVector target; 503 504 /** 505 * Create an {@link Evaluation} with no weights. 506 * 507 * @param model the model function 508 * @param target the observed values 509 * @param point the abscissa 510 */ 511 private LazyUnweightedEvaluation(final ValueAndJacobianFunction model, 512 final RealVector target, 513 final RealVector point) { 514 super(target.getDimension()); 515 // Safe to cast as long as we control usage of this class. 516 this.model = model; 517 this.point = point; 518 this.target = target; 519 } 520 521 /** {@inheritDoc} */ 522 @Override 523 public RealMatrix getJacobian() { 524 return model.computeJacobian(point.toArray()); 525 } 526 527 /** {@inheritDoc} */ 528 @Override 529 public RealVector getPoint() { 530 return point; 531 } 532 533 /** {@inheritDoc} */ 534 @Override 535 public RealVector getResiduals() { 536 return target.subtract(model.computeValue(point.toArray())); 537 } 538 } 539 } 540} 541