001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.optimization; 019 020 import org.apache.commons.math.FunctionEvaluationException; 021 import org.apache.commons.math.MathRuntimeException; 022 import org.apache.commons.math.analysis.MultivariateRealFunction; 023 import org.apache.commons.math.analysis.MultivariateVectorialFunction; 024 import org.apache.commons.math.linear.RealMatrix; 025 026 /** This class converts {@link MultivariateVectorialFunction vectorial 027 * objective functions} to {@link MultivariateRealFunction scalar objective functions} 028 * when the goal is to minimize them. 029 * <p> 030 * This class is mostly used when the vectorial objective function represents 031 * a theoretical result computed from a point set applied to a model and 032 * the models point must be adjusted to fit the theoretical result to some 033 * reference observations. The observations may be obtained for example from 034 * physical measurements whether the model is built from theoretical 035 * considerations. 036 * </p> 037 * <p> 038 * This class computes a possibly weighted squared sum of the residuals, which is 039 * a scalar value. The residuals are the difference between the theoretical model 040 * (i.e. the output of the vectorial objective function) and the observations. The 041 * class implements the {@link MultivariateRealFunction} interface and can therefore be 042 * minimized by any optimizer supporting scalar objectives functions.This is one way 043 * to perform a least square estimation. There are other ways to do this without using 044 * this converter, as some optimization algorithms directly support vectorial objective 045 * functions. 046 * </p> 047 * <p> 048 * This class support combination of residuals with or without weights and correlations. 049 * </p> 050 * 051 * @see MultivariateRealFunction 052 * @see MultivariateVectorialFunction 053 * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $ 054 * @since 2.0 055 */ 056 057 public class LeastSquaresConverter implements MultivariateRealFunction { 058 059 /** Underlying vectorial function. */ 060 private final MultivariateVectorialFunction function; 061 062 /** Observations to be compared to objective function to compute residuals. */ 063 private final double[] observations; 064 065 /** Optional weights for the residuals. */ 066 private final double[] weights; 067 068 /** Optional scaling matrix (weight and correlations) for the residuals. */ 069 private final RealMatrix scale; 070 071 /** Build a simple converter for uncorrelated residuals with the same weight. 072 * @param function vectorial residuals function to wrap 073 * @param observations observations to be compared to objective function to compute residuals 074 */ 075 public LeastSquaresConverter(final MultivariateVectorialFunction function, 076 final double[] observations) { 077 this.function = function; 078 this.observations = observations.clone(); 079 this.weights = null; 080 this.scale = null; 081 } 082 083 /** Build a simple converter for uncorrelated residuals with the specific weights. 084 * <p> 085 * The scalar objective function value is computed as: 086 * <pre> 087 * objective = ∑weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup> 088 * </pre> 089 * </p> 090 * <p> 091 * Weights can be used for example to combine residuals with different standard 092 * deviations. As an example, consider a residuals array in which even elements 093 * are angular measurements in degrees with a 0.01° standard deviation and 094 * odd elements are distance measurements in meters with a 15m standard deviation. 095 * In this case, the weights array should be initialized with value 096 * 1.0/(0.01<sup>2</sup>) in the even elements and 1.0/(15.0<sup>2</sup>) in the 097 * odd elements (i.e. reciprocals of variances). 098 * </p> 099 * <p> 100 * The array computed by the objective function, the observations array and the 101 * weights array must have consistent sizes or a {@link FunctionEvaluationException} will be 102 * triggered while computing the scalar objective. 103 * </p> 104 * @param function vectorial residuals function to wrap 105 * @param observations observations to be compared to objective function to compute residuals 106 * @param weights weights to apply to the residuals 107 * @exception IllegalArgumentException if the observations vector and the weights 108 * vector dimensions don't match (objective function dimension is checked only when 109 * the {@link #value(double[])} method is called) 110 */ 111 public LeastSquaresConverter(final MultivariateVectorialFunction function, 112 final double[] observations, final double[] weights) 113 throws IllegalArgumentException { 114 if (observations.length != weights.length) { 115 throw MathRuntimeException.createIllegalArgumentException( 116 "dimension mismatch {0} != {1}", 117 observations.length, weights.length); 118 } 119 this.function = function; 120 this.observations = observations.clone(); 121 this.weights = weights.clone(); 122 this.scale = null; 123 } 124 125 /** Build a simple converter for correlated residuals with the specific weights. 126 * <p> 127 * The scalar objective function value is computed as: 128 * <pre> 129 * objective = y<sup>T</sup>y with y = scale×(observation-objective) 130 * </pre> 131 * </p> 132 * <p> 133 * The array computed by the objective function, the observations array and the 134 * the scaling matrix must have consistent sizes or a {@link FunctionEvaluationException} 135 * will be triggered while computing the scalar objective. 136 * </p> 137 * @param function vectorial residuals function to wrap 138 * @param observations observations to be compared to objective function to compute residuals 139 * @param scale scaling matrix 140 * @exception IllegalArgumentException if the observations vector and the scale 141 * matrix dimensions don't match (objective function dimension is checked only when 142 * the {@link #value(double[])} method is called) 143 */ 144 public LeastSquaresConverter(final MultivariateVectorialFunction function, 145 final double[] observations, final RealMatrix scale) 146 throws IllegalArgumentException { 147 if (observations.length != scale.getColumnDimension()) { 148 throw MathRuntimeException.createIllegalArgumentException( 149 "dimension mismatch {0} != {1}", 150 observations.length, scale.getColumnDimension()); 151 } 152 this.function = function; 153 this.observations = observations.clone(); 154 this.weights = null; 155 this.scale = scale.copy(); 156 } 157 158 /** {@inheritDoc} */ 159 public double value(final double[] point) throws FunctionEvaluationException { 160 161 // compute residuals 162 final double[] residuals = function.value(point); 163 if (residuals.length != observations.length) { 164 throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}", 165 residuals.length, observations.length); 166 } 167 for (int i = 0; i < residuals.length; ++i) { 168 residuals[i] -= observations[i]; 169 } 170 171 // compute sum of squares 172 double sumSquares = 0; 173 if (weights != null) { 174 for (int i = 0; i < residuals.length; ++i) { 175 final double ri = residuals[i]; 176 sumSquares += weights[i] * ri * ri; 177 } 178 } else if (scale != null) { 179 for (final double yi : scale.operate(residuals)) { 180 sumSquares += yi * yi; 181 } 182 } else { 183 for (final double ri : residuals) { 184 sumSquares += ri * ri; 185 } 186 } 187 188 return sumSquares; 189 190 } 191 192 }