1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.optimization; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.FunctionEvaluationException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.MathRuntimeException; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.MultivariateRealFunction; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.MultivariateVectorialFunction; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.RealMatrix; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** This class converts {@link MultivariateVectorialFunction vectorial 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * objective functions} to {@link MultivariateRealFunction scalar objective functions} 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * when the goal is to minimize them. 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This class is mostly used when the vectorial objective function represents 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a theoretical result computed from a point set applied to a model and 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the models point must be adjusted to fit the theoretical result to some 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * reference observations. The observations may be obtained for example from 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * physical measurements whether the model is built from theoretical 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * considerations. 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This class computes a possibly weighted squared sum of the residuals, which is 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a scalar value. The residuals are the difference between the theoretical model 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (i.e. the output of the vectorial objective function) and the observations. The 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * class implements the {@link MultivariateRealFunction} interface and can therefore be 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * minimized by any optimizer supporting scalar objectives functions.This is one way 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to perform a least square estimation. There are other ways to do this without using 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this converter, as some optimization algorithms directly support vectorial objective 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * functions. 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This class support combination of residuals with or without weights and correlations. 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see MultivariateRealFunction 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see MultivariateVectorialFunction 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.0 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class LeastSquaresConverter implements MultivariateRealFunction { 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Underlying vectorial function. */ 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final MultivariateVectorialFunction function; 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Observations to be compared to objective function to compute residuals. */ 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] observations; 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Optional weights for the residuals. */ 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] weights; 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Optional scaling matrix (weight and correlations) for the residuals. */ 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final RealMatrix scale; 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Build a simple converter for uncorrelated residuals with the same weight. 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param function vectorial residuals function to wrap 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param observations observations to be compared to objective function to compute residuals 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public LeastSquaresConverter(final MultivariateVectorialFunction function, 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] observations) { 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.function = function; 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.observations = observations.clone(); 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.weights = null; 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.scale = null; 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Build a simple converter for uncorrelated residuals with the specific weights. 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The scalar objective function value is computed as: 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * objective = ∑weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup> 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Weights can be used for example to combine residuals with different standard 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * deviations. As an example, consider a residuals array in which even elements 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * are angular measurements in degrees with a 0.01° standard deviation and 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * odd elements are distance measurements in meters with a 15m standard deviation. 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * In this case, the weights array should be initialized with value 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 1.0/(0.01<sup>2</sup>) in the even elements and 1.0/(15.0<sup>2</sup>) in the 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * odd elements (i.e. reciprocals of variances). 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The array computed by the objective function, the observations array and the 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * weights array must have consistent sizes or a {@link FunctionEvaluationException} will be 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * triggered while computing the scalar objective. 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param function vectorial residuals function to wrap 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param observations observations to be compared to objective function to compute residuals 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param weights weights to apply to the residuals 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception IllegalArgumentException if the observations vector and the weights 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * vector dimensions don't match (objective function dimension is checked only when 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the {@link #value(double[])} method is called) 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public LeastSquaresConverter(final MultivariateVectorialFunction function, 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] observations, final double[] weights) 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException { 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (observations.length != weights.length) { 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond observations.length, weights.length); 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.function = function; 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.observations = observations.clone(); 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.weights = weights.clone(); 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.scale = null; 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Build a simple converter for correlated residuals with the specific weights. 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The scalar objective function value is computed as: 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * objective = y<sup>T</sup>y with y = scale×(observation-objective) 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The array computed by the objective function, the observations array and the 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the scaling matrix must have consistent sizes or a {@link FunctionEvaluationException} 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * will be triggered while computing the scalar objective. 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param function vectorial residuals function to wrap 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param observations observations to be compared to objective function to compute residuals 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param scale scaling matrix 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception IllegalArgumentException if the observations vector and the scale 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * matrix dimensions don't match (objective function dimension is checked only when 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the {@link #value(double[])} method is called) 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public LeastSquaresConverter(final MultivariateVectorialFunction function, 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] observations, final RealMatrix scale) 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException { 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (observations.length != scale.getColumnDimension()) { 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond observations.length, scale.getColumnDimension()); 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.function = function; 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.observations = observations.clone(); 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.weights = null; 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.scale = scale.copy(); 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(final double[] point) throws FunctionEvaluationException { 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute residuals 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] residuals = function.value(point); 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (residuals.length != observations.length) { 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new FunctionEvaluationException(point,LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond residuals.length, observations.length); 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < residuals.length; ++i) { 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond residuals[i] -= observations[i]; 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute sum of squares 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sumSquares = 0; 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (weights != null) { 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < residuals.length; ++i) { 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double ri = residuals[i]; 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sumSquares += weights[i] * ri * ri; 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else if (scale != null) { 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (final double yi : scale.operate(residuals)) { 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sumSquares += yi * yi; 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (final double ri : residuals) { 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sumSquares += ri * ri; 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return sumSquares; 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 194