| /* |
| * 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 |
| * |
| * http://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. |
| */ |
| |
| package org.apache.commons.math.optimization; |
| |
| import org.apache.commons.math.FunctionEvaluationException; |
| import org.apache.commons.math.MathRuntimeException; |
| import org.apache.commons.math.analysis.MultivariateRealFunction; |
| import org.apache.commons.math.analysis.MultivariateVectorialFunction; |
| import org.apache.commons.math.exception.util.LocalizedFormats; |
| import org.apache.commons.math.linear.RealMatrix; |
| |
| /** This class converts {@link MultivariateVectorialFunction vectorial |
| * objective functions} to {@link MultivariateRealFunction scalar objective functions} |
| * when the goal is to minimize them. |
| * <p> |
| * This class is mostly used when the vectorial objective function represents |
| * a theoretical result computed from a point set applied to a model and |
| * the models point must be adjusted to fit the theoretical result to some |
| * reference observations. The observations may be obtained for example from |
| * physical measurements whether the model is built from theoretical |
| * considerations. |
| * </p> |
| * <p> |
| * This class computes a possibly weighted squared sum of the residuals, which is |
| * a scalar value. The residuals are the difference between the theoretical model |
| * (i.e. the output of the vectorial objective function) and the observations. The |
| * class implements the {@link MultivariateRealFunction} interface and can therefore be |
| * minimized by any optimizer supporting scalar objectives functions.This is one way |
| * to perform a least square estimation. There are other ways to do this without using |
| * this converter, as some optimization algorithms directly support vectorial objective |
| * functions. |
| * </p> |
| * <p> |
| * This class support combination of residuals with or without weights and correlations. |
| * </p> |
| * |
| * @see MultivariateRealFunction |
| * @see MultivariateVectorialFunction |
| * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ |
| * @since 2.0 |
| */ |
| |
| public class LeastSquaresConverter implements MultivariateRealFunction { |
| |
| /** Underlying vectorial function. */ |
| private final MultivariateVectorialFunction function; |
| |
| /** Observations to be compared to objective function to compute residuals. */ |
| private final double[] observations; |
| |
| /** Optional weights for the residuals. */ |
| private final double[] weights; |
| |
| /** Optional scaling matrix (weight and correlations) for the residuals. */ |
| private final RealMatrix scale; |
| |
| /** Build a simple converter for uncorrelated residuals with the same weight. |
| * @param function vectorial residuals function to wrap |
| * @param observations observations to be compared to objective function to compute residuals |
| */ |
| public LeastSquaresConverter(final MultivariateVectorialFunction function, |
| final double[] observations) { |
| this.function = function; |
| this.observations = observations.clone(); |
| this.weights = null; |
| this.scale = null; |
| } |
| |
| /** Build a simple converter for uncorrelated residuals with the specific weights. |
| * <p> |
| * The scalar objective function value is computed as: |
| * <pre> |
| * objective = ∑weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup> |
| * </pre> |
| * </p> |
| * <p> |
| * Weights can be used for example to combine residuals with different standard |
| * deviations. As an example, consider a residuals array in which even elements |
| * are angular measurements in degrees with a 0.01° standard deviation and |
| * odd elements are distance measurements in meters with a 15m standard deviation. |
| * In this case, the weights array should be initialized with value |
| * 1.0/(0.01<sup>2</sup>) in the even elements and 1.0/(15.0<sup>2</sup>) in the |
| * odd elements (i.e. reciprocals of variances). |
| * </p> |
| * <p> |
| * The array computed by the objective function, the observations array and the |
| * weights array must have consistent sizes or a {@link FunctionEvaluationException} will be |
| * triggered while computing the scalar objective. |
| * </p> |
| * @param function vectorial residuals function to wrap |
| * @param observations observations to be compared to objective function to compute residuals |
| * @param weights weights to apply to the residuals |
| * @exception IllegalArgumentException if the observations vector and the weights |
| * vector dimensions don't match (objective function dimension is checked only when |
| * the {@link #value(double[])} method is called) |
| */ |
| public LeastSquaresConverter(final MultivariateVectorialFunction function, |
| final double[] observations, final double[] weights) |
| throws IllegalArgumentException { |
| if (observations.length != weights.length) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, |
| observations.length, weights.length); |
| } |
| this.function = function; |
| this.observations = observations.clone(); |
| this.weights = weights.clone(); |
| this.scale = null; |
| } |
| |
| /** Build a simple converter for correlated residuals with the specific weights. |
| * <p> |
| * The scalar objective function value is computed as: |
| * <pre> |
| * objective = y<sup>T</sup>y with y = scale×(observation-objective) |
| * </pre> |
| * </p> |
| * <p> |
| * The array computed by the objective function, the observations array and the |
| * the scaling matrix must have consistent sizes or a {@link FunctionEvaluationException} |
| * will be triggered while computing the scalar objective. |
| * </p> |
| * @param function vectorial residuals function to wrap |
| * @param observations observations to be compared to objective function to compute residuals |
| * @param scale scaling matrix |
| * @exception IllegalArgumentException if the observations vector and the scale |
| * matrix dimensions don't match (objective function dimension is checked only when |
| * the {@link #value(double[])} method is called) |
| */ |
| public LeastSquaresConverter(final MultivariateVectorialFunction function, |
| final double[] observations, final RealMatrix scale) |
| throws IllegalArgumentException { |
| if (observations.length != scale.getColumnDimension()) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, |
| observations.length, scale.getColumnDimension()); |
| } |
| this.function = function; |
| this.observations = observations.clone(); |
| this.weights = null; |
| this.scale = scale.copy(); |
| } |
| |
| /** {@inheritDoc} */ |
| public double value(final double[] point) throws FunctionEvaluationException { |
| |
| // compute residuals |
| final double[] residuals = function.value(point); |
| if (residuals.length != observations.length) { |
| throw new FunctionEvaluationException(point,LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, |
| residuals.length, observations.length); |
| } |
| for (int i = 0; i < residuals.length; ++i) { |
| residuals[i] -= observations[i]; |
| } |
| |
| // compute sum of squares |
| double sumSquares = 0; |
| if (weights != null) { |
| for (int i = 0; i < residuals.length; ++i) { |
| final double ri = residuals[i]; |
| sumSquares += weights[i] * ri * ri; |
| } |
| } else if (scale != null) { |
| for (final double yi : scale.operate(residuals)) { |
| sumSquares += yi * yi; |
| } |
| } else { |
| for (final double ri : residuals) { |
| sumSquares += ri * ri; |
| } |
| } |
| |
| return sumSquares; |
| |
| } |
| |
| } |