LeastSquaresFactory.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
/*
* This is not the original file distributed by the Apache Software Foundation
* It has been modified by the Hipparchus project
*/
package org.hipparchus.optim.nonlinear.vector.leastsquares;
import org.hipparchus.analysis.MultivariateMatrixFunction;
import org.hipparchus.analysis.MultivariateVectorFunction;
import org.hipparchus.exception.MathIllegalStateException;
import org.hipparchus.linear.Array2DRowRealMatrix;
import org.hipparchus.linear.ArrayRealVector;
import org.hipparchus.linear.DiagonalMatrix;
import org.hipparchus.linear.EigenDecompositionSymmetric;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.linear.RealVector;
import org.hipparchus.optim.AbstractOptimizationProblem;
import org.hipparchus.optim.ConvergenceChecker;
import org.hipparchus.optim.LocalizedOptimFormats;
import org.hipparchus.optim.PointVectorValuePair;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Incrementor;
import org.hipparchus.util.Pair;
/**
* A Factory for creating {@link LeastSquaresProblem}s.
*
*/
public class LeastSquaresFactory {
/** Prevent instantiation. */
private LeastSquaresFactory() {}
/**
* Create a {@link org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem}
* from the given elements. There will be no weights applied (unit weights).
*
* @param model the model function. Produces the computed values.
* @param observed the observed (target) values
* @param start the initial guess.
* @param weight the weight matrix
* @param checker convergence checker
* @param maxEvaluations the maximum number of times to evaluate the model
* @param maxIterations the maximum number to times to iterate in the algorithm
* @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)}
* will defer the evaluation until access to the value is requested.
* @param paramValidator Model parameters validator.
* @return the specified General Least Squares problem.
*
*/
public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
final RealVector observed,
final RealVector start,
final RealMatrix weight,
final ConvergenceChecker<Evaluation> checker,
final int maxEvaluations,
final int maxIterations,
final boolean lazyEvaluation,
final ParameterValidator paramValidator) {
final LeastSquaresProblem p = new LocalLeastSquaresProblem(model,
observed,
start,
checker,
maxEvaluations,
maxIterations,
lazyEvaluation,
paramValidator);
if (weight != null) {
return weightMatrix(p, weight);
} else {
return p;
}
}
/**
* Create a {@link org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem}
* from the given elements. There will be no weights applied (unit weights).
*
* @param model the model function. Produces the computed values.
* @param observed the observed (target) values
* @param start the initial guess.
* @param checker convergence checker
* @param maxEvaluations the maximum number of times to evaluate the model
* @param maxIterations the maximum number to times to iterate in the algorithm
* @return the specified General Least Squares problem.
*/
public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
final RealVector observed,
final RealVector start,
final ConvergenceChecker<Evaluation> checker,
final int maxEvaluations,
final int maxIterations) {
return create(model,
observed,
start,
null,
checker,
maxEvaluations,
maxIterations,
false,
null);
}
/**
* Create a {@link org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem}
* from the given elements.
*
* @param model the model function. Produces the computed values.
* @param observed the observed (target) values
* @param start the initial guess.
* @param weight the weight matrix
* @param checker convergence checker
* @param maxEvaluations the maximum number of times to evaluate the model
* @param maxIterations the maximum number to times to iterate in the algorithm
* @return the specified General Least Squares problem.
*/
public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
final RealVector observed,
final RealVector start,
final RealMatrix weight,
final ConvergenceChecker<Evaluation> checker,
final int maxEvaluations,
final int maxIterations) {
return weightMatrix(create(model,
observed,
start,
checker,
maxEvaluations,
maxIterations),
weight);
}
/**
* Create a {@link org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem}
* from the given elements.
* <p>
* This factory method is provided for continuity with previous interfaces. Newer
* applications should use {@link #create(MultivariateJacobianFunction, RealVector,
* RealVector, ConvergenceChecker, int, int)}, or {@link #create(MultivariateJacobianFunction,
* RealVector, RealVector, RealMatrix, ConvergenceChecker, int, int)}.
*
* @param model the model function. Produces the computed values.
* @param jacobian the jacobian of the model with respect to the parameters
* @param observed the observed (target) values
* @param start the initial guess.
* @param weight the weight matrix
* @param checker convergence checker
* @param maxEvaluations the maximum number of times to evaluate the model
* @param maxIterations the maximum number to times to iterate in the algorithm
* @return the specified General Least Squares problem.
*/
public static LeastSquaresProblem create(final MultivariateVectorFunction model,
final MultivariateMatrixFunction jacobian,
final double[] observed,
final double[] start,
final RealMatrix weight,
final ConvergenceChecker<Evaluation> checker,
final int maxEvaluations,
final int maxIterations) {
return create(model(model, jacobian),
new ArrayRealVector(observed, false),
new ArrayRealVector(start, false),
weight,
checker,
maxEvaluations,
maxIterations);
}
/**
* Apply a dense weight matrix to the {@link LeastSquaresProblem}.
*
* @param problem the unweighted problem
* @param weights the matrix of weights
* @return a new {@link LeastSquaresProblem} with the weights applied. The original
* {@code problem} is not modified.
*/
public static LeastSquaresProblem weightMatrix(final LeastSquaresProblem problem,
final RealMatrix weights) {
final RealMatrix weightSquareRoot = squareRoot(weights);
return new LeastSquaresAdapter(problem) {
/** {@inheritDoc} */
@Override
public Evaluation evaluate(final RealVector point) {
return new DenseWeightedEvaluation(super.evaluate(point), weightSquareRoot);
}
};
}
/**
* Apply a diagonal weight matrix to the {@link LeastSquaresProblem}.
*
* @param problem the unweighted problem
* @param weights the diagonal of the weight matrix
* @return a new {@link LeastSquaresProblem} with the weights applied. The original
* {@code problem} is not modified.
*/
public static LeastSquaresProblem weightDiagonal(final LeastSquaresProblem problem,
final RealVector weights) {
// TODO more efficient implementation
return weightMatrix(problem, new DiagonalMatrix(weights.toArray()));
}
/**
* Count the evaluations of a particular problem. The {@code counter} will be
* incremented every time {@link LeastSquaresProblem#evaluate(RealVector)} is called on
* the <em>returned</em> problem.
*
* @param problem the problem to track.
* @param counter the counter to increment.
* @return a least squares problem that tracks evaluations
*/
public static LeastSquaresProblem countEvaluations(final LeastSquaresProblem problem,
final Incrementor counter) {
return new LeastSquaresAdapter(problem) {
/** {@inheritDoc} */
@Override
public Evaluation evaluate(final RealVector point) {
counter.increment();
return super.evaluate(point);
}
// Delegate the rest.
};
}
/**
* View a convergence checker specified for a {@link PointVectorValuePair} as one
* specified for an {@link Evaluation}.
*
* @param checker the convergence checker to adapt.
* @return a convergence checker that delegates to {@code checker}.
*/
public static ConvergenceChecker<Evaluation> evaluationChecker(final ConvergenceChecker<PointVectorValuePair> checker) {
return new ConvergenceChecker<Evaluation>() {
/** {@inheritDoc} */
@Override
public boolean converged(final int iteration,
final Evaluation previous,
final Evaluation current) {
return checker.converged(
iteration,
new PointVectorValuePair(
previous.getPoint().toArray(),
previous.getResiduals().toArray(),
false),
new PointVectorValuePair(
current.getPoint().toArray(),
current.getResiduals().toArray(),
false)
);
}
};
}
/**
* Computes the square-root of the weight matrix.
*
* @param m Symmetric, positive-definite (weight) matrix.
* @return the square-root of the weight matrix.
*/
private static RealMatrix squareRoot(final RealMatrix m) {
if (m instanceof DiagonalMatrix) {
final int dim = m.getRowDimension();
final RealMatrix sqrtM = new DiagonalMatrix(dim);
for (int i = 0; i < dim; i++) {
sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
}
return sqrtM;
} else {
final EigenDecompositionSymmetric dec = new EigenDecompositionSymmetric(m);
return dec.getSquareRoot();
}
}
/**
* Combine a {@link MultivariateVectorFunction} with a {@link
* MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}.
*
* @param value the vector value function
* @param jacobian the Jacobian function
* @return a function that computes both at the same time
*/
public static MultivariateJacobianFunction model(final MultivariateVectorFunction value,
final MultivariateMatrixFunction jacobian) {
return new LocalValueAndJacobianFunction(value, jacobian);
}
/**
* Combine a {@link MultivariateVectorFunction} with a {@link
* MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}.
*/
private static class LocalValueAndJacobianFunction
implements ValueAndJacobianFunction {
/** Model. */
private final MultivariateVectorFunction value;
/** Model's Jacobian. */
private final MultivariateMatrixFunction jacobian;
/**
* @param value Model function.
* @param jacobian Model's Jacobian function.
*/
LocalValueAndJacobianFunction(final MultivariateVectorFunction value,
final MultivariateMatrixFunction jacobian) {
this.value = value;
this.jacobian = jacobian;
}
/** {@inheritDoc} */
@Override
public Pair<RealVector, RealMatrix> value(final RealVector point) {
//TODO get array from RealVector without copying?
final double[] p = point.toArray();
// Evaluate.
return new Pair<RealVector, RealMatrix>(computeValue(p),
computeJacobian(p));
}
/** {@inheritDoc} */
@Override
public RealVector computeValue(final double[] params) {
return new ArrayRealVector(value.value(params), false);
}
/** {@inheritDoc} */
@Override
public RealMatrix computeJacobian(final double[] params) {
return new Array2DRowRealMatrix(jacobian.value(params), false);
}
}
/**
* A private, "field" immutable (not "real" immutable) implementation of {@link
* LeastSquaresProblem}.
*/
private static class LocalLeastSquaresProblem
extends AbstractOptimizationProblem<Evaluation>
implements LeastSquaresProblem {
/** Target values for the model function at optimum. */
private final RealVector target;
/** Model function. */
private final MultivariateJacobianFunction model;
/** Initial guess. */
private final RealVector start;
/** Whether to use lazy evaluation. */
private final boolean lazyEvaluation;
/** Model parameters validator. */
private final ParameterValidator paramValidator;
/**
* Create a {@link LeastSquaresProblem} from the given data.
*
* @param model the model function
* @param target the observed data
* @param start the initial guess
* @param checker the convergence checker
* @param maxEvaluations the allowed evaluations
* @param maxIterations the allowed iterations
* @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)}
* will defer the evaluation until access to the value is requested.
* @param paramValidator Model parameters validator.
*/
LocalLeastSquaresProblem(final MultivariateJacobianFunction model,
final RealVector target,
final RealVector start,
final ConvergenceChecker<Evaluation> checker,
final int maxEvaluations,
final int maxIterations,
final boolean lazyEvaluation,
final ParameterValidator paramValidator) {
super(maxEvaluations, maxIterations, checker);
this.target = target;
this.model = model;
this.start = start;
this.lazyEvaluation = lazyEvaluation;
this.paramValidator = paramValidator;
if (lazyEvaluation &&
!(model instanceof ValueAndJacobianFunction)) {
// Lazy evaluation requires that value and Jacobian
// can be computed separately.
throw new MathIllegalStateException(LocalizedOptimFormats.INVALID_IMPLEMENTATION,
model.getClass().getName());
}
}
/** {@inheritDoc} */
@Override
public int getObservationSize() {
return target.getDimension();
}
/** {@inheritDoc} */
@Override
public int getParameterSize() {
return start.getDimension();
}
/** {@inheritDoc} */
@Override
public RealVector getStart() {
return start == null ? null : start.copy();
}
/** {@inheritDoc} */
@Override
public Evaluation evaluate(final RealVector point) {
// Copy so optimizer can change point without changing our instance.
final RealVector p = paramValidator == null ?
point.copy() :
paramValidator.validate(point.copy());
if (lazyEvaluation) {
return new LazyUnweightedEvaluation((ValueAndJacobianFunction) model,
target,
p);
} else {
// Evaluate value and jacobian in one function call.
final Pair<RealVector, RealMatrix> value = model.value(p);
return new UnweightedEvaluation(value.getFirst(),
value.getSecond(),
target,
p);
}
}
/**
* Container with the model evaluation at a particular point.
*/
private static class UnweightedEvaluation extends AbstractEvaluation {
/** Point of evaluation. */
private final RealVector point;
/** Derivative at point. */
private final RealMatrix jacobian;
/** Computed residuals. */
private final RealVector residuals;
/**
* Create an {@link Evaluation} with no weights.
*
* @param values the computed function values
* @param jacobian the computed function Jacobian
* @param target the observed values
* @param point the abscissa
*/
private UnweightedEvaluation(final RealVector values,
final RealMatrix jacobian,
final RealVector target,
final RealVector point) {
super(target.getDimension());
this.jacobian = jacobian;
this.point = point;
this.residuals = target.subtract(values);
}
/** {@inheritDoc} */
@Override
public RealMatrix getJacobian() {
return jacobian;
}
/** {@inheritDoc} */
@Override
public RealVector getPoint() {
return point;
}
/** {@inheritDoc} */
@Override
public RealVector getResiduals() {
return residuals;
}
}
/**
* Container with the model <em>lazy</em> evaluation at a particular point.
*/
private static class LazyUnweightedEvaluation extends AbstractEvaluation {
/** Point of evaluation. */
private final RealVector point;
/** Model and Jacobian functions. */
private final ValueAndJacobianFunction model;
/** Target values for the model function at optimum. */
private final RealVector target;
/**
* Create an {@link Evaluation} with no weights.
*
* @param model the model function
* @param target the observed values
* @param point the abscissa
*/
private LazyUnweightedEvaluation(final ValueAndJacobianFunction model,
final RealVector target,
final RealVector point) {
super(target.getDimension());
// Safe to cast as long as we control usage of this class.
this.model = model;
this.point = point;
this.target = target;
}
/** {@inheritDoc} */
@Override
public RealMatrix getJacobian() {
return model.computeJacobian(point.toArray());
}
/** {@inheritDoc} */
@Override
public RealVector getPoint() {
return point;
}
/** {@inheritDoc} */
@Override
public RealVector getResiduals() {
return target.subtract(model.computeValue(point.toArray()));
}
}
}
}