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