1/*2* Licensed to the Apache Software Foundation (ASF) under one or more3* contributor license agreements. See the NOTICE file distributed with4* this work for additional information regarding copyright ownership.5* The ASF licenses this file to You under the Apache License, Version 2.06* (the "License"); you may not use this file except in compliance with7* the License. You may obtain a copy of the License at8*9* http://www.apache.org/licenses/LICENSE-2.010*11* Unless required by applicable law or agreed to in writing, software12* 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 and15* limitations under the License.16*/17 18packageorg.apache.commons.math3.optim.nonlinear.scalar; 19 20importorg.apache.commons.math3.analysis.MultivariateFunction; 21importorg.apache.commons.math3.analysis.MultivariateVectorFunction; 22importorg.apache.commons.math3.exception.DimensionMismatchException; 23importorg.apache.commons.math3.linear.RealMatrix; 24 25/**26* This class converts27* {@link MultivariateVectorFunction vectorial objective functions} to28* {@link MultivariateFunction scalar objective functions}29* when the goal is to minimize them.30* <br/>31* This class is mostly used when the vectorial objective function represents32* a theoretical result computed from a point set applied to a model and33* the models point must be adjusted to fit the theoretical result to some34* reference observations. The observations may be obtained for example from35* physical measurements whether the model is built from theoretical36* considerations.37* <br/>38* This class computes a possibly weighted squared sum of the residuals, which is39* a scalar value. The residuals are the difference between the theoretical model40* (i.e. the output of the vectorial objective function) and the observations. The41* class implements the {@link MultivariateFunction} interface and can therefore be42* minimized by any optimizer supporting scalar objectives functions.This is one way43* to perform a least square estimation. There are other ways to do this without using44* this converter, as some optimization algorithms directly support vectorial objective45* functions.46* <br/>47* This class support combination of residuals with or without weights and correlations.48*49* @see MultivariateFunction50* @see MultivariateVectorFunction51* @since 2.052*/53 54publicclassLeastSquaresConverterimplementsMultivariateFunction { 55/** Underlying vectorial function. */56privatefinalMultivariateVectorFunction function; 57/** Observations to be compared to objective function to compute residuals. */58privatefinaldouble[] observations; 59/** Optional weights for the residuals. */60privatefinaldouble[] weights; 61/** Optional scaling matrix (weight and correlations) for the residuals. */62privatefinalRealMatrix scale; 63 64/**65* Builds a simple converter for uncorrelated residuals with identical66* weights.67*68* @param function vectorial residuals function to wrap69* @param observations observations to be compared to objective function to compute residuals70*/71publicLeastSquaresConverter(finalMultivariateVectorFunction function, 72finaldouble[] observations) { 73this.function = function; 74this.observations = observations.clone(); 75this.weights =null; 76this.scale =null; 77 } 78 79/**80* Builds a simple converter for uncorrelated residuals with the81* specified weights.82* <p>83* The scalar objective function value is computed as:84* <pre>85* objective = ∑weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup>86* </pre>87* </p>88* <p>89* Weights can be used for example to combine residuals with different standard90* deviations. As an example, consider a residuals array in which even elements91* are angular measurements in degrees with a 0.01° standard deviation and92* odd elements are distance measurements in meters with a 15m standard deviation.93* In this case, the weights array should be initialized with value94* 1.0/(0.01<sup>2</sup>) in the even elements and 1.0/(15.0<sup>2</sup>) in the95* odd elements (i.e. reciprocals of variances).96* </p>97* <p>98* The array computed by the objective function, the observations array and the99* weights array must have consistent sizes or a {@link DimensionMismatchException}100* will be triggered while computing the scalar objective.101* </p>102*103* @param function vectorial residuals function to wrap104* @param observations observations to be compared to objective function to compute residuals105* @param weights weights to apply to the residuals106* @throws DimensionMismatchException if the observations vector and the weights107* vector dimensions do not match (objective function dimension is checked only when108* the {@link #value(double[])} method is called)109*/110publicLeastSquaresConverter(finalMultivariateVectorFunction function, 111finaldouble[] observations, 112finaldouble[] weights) { 113if(observations.length != weights.length) { 114thrownewDimensionMismatchException(observations.length, weights.length); 115 } 116this.function = function; 117this.observations = observations.clone(); 118this.weights = weights.clone(); 119this.scale =null; 120 } 121 122/**123* Builds a simple converter for correlated residuals with the124* specified weights.125* <p>126* The scalar objective function value is computed as:127* <pre>128* objective = y<sup>T</sup>y with y = scale×(observation-objective)129* </pre>130* </p>131* <p>132* The array computed by the objective function, the observations array and the133* the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}134* will be triggered while computing the scalar objective.135* </p>136*137* @param function vectorial residuals function to wrap138* @param observations observations to be compared to objective function to compute residuals139* @param scale scaling matrix140* @throws DimensionMismatchException if the observations vector and the scale141* matrix dimensions do not match (objective function dimension is checked only when142* the {@link #value(double[])} method is called)143*/144publicLeastSquaresConverter(finalMultivariateVectorFunction function, 145finaldouble[] observations, 146finalRealMatrix scale) { 147if(observations.length != scale.getColumnDimension()) { 148thrownewDimensionMismatchException(observations.length, scale.getColumnDimension()); 149 } 150this.function = function; 151this.observations = observations.clone(); 152this.weights =null; 153this.scale = scale.copy(); 154 } 155 156/** {@inheritDoc} */157publicdoublevalue(finaldouble[] point) { 158// compute residuals159finaldouble[] residuals = function.value(point); 160if(residuals.length != observations.length) { 161thrownewDimensionMismatchException(residuals.length, observations.length); 162 } 163for(inti = 0; i < residuals.length; ++i) { 164 residuals[i] -= observations[i]; 165 } 166 167// compute sum of squares168doublesumSquares = 0; 169if(weights !=null) { 170for(inti = 0; i < residuals.length; ++i) { 171finaldoubleri = residuals[i]; 172 sumSquares += weights[i] * ri * ri; 173 } 174 }elseif(scale !=null) { 175for(finaldoubleyi : scale.operate(residuals)) { 176 sumSquares += yi * yi; 177 } 178 }else{ 179for(finaldoubleri : residuals) { 180 sumSquares += ri * ri; 181 } 182 } 183 184returnsumSquares; 185 } 186 }