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    package org.apache.commons.math3.optim.nonlinear.vector.jacobian;
018    
019    import org.apache.commons.math3.exception.DimensionMismatchException;
020    import org.apache.commons.math3.exception.TooManyEvaluationsException;
021    import org.apache.commons.math3.linear.ArrayRealVector;
022    import org.apache.commons.math3.linear.RealMatrix;
023    import org.apache.commons.math3.linear.DecompositionSolver;
024    import org.apache.commons.math3.linear.MatrixUtils;
025    import org.apache.commons.math3.linear.QRDecomposition;
026    import org.apache.commons.math3.linear.EigenDecomposition;
027    import org.apache.commons.math3.optim.OptimizationData;
028    import org.apache.commons.math3.optim.ConvergenceChecker;
029    import org.apache.commons.math3.optim.PointVectorValuePair;
030    import org.apache.commons.math3.optim.nonlinear.vector.Weight;
031    import org.apache.commons.math3.optim.nonlinear.vector.JacobianMultivariateVectorOptimizer;
032    import org.apache.commons.math3.util.FastMath;
033    
034    /**
035     * Base class for implementing least-squares optimizers.
036     * It provides methods for error estimation.
037     *
038     * @version $Id$
039     * @since 3.1
040     */
041    public abstract class AbstractLeastSquaresOptimizer
042        extends JacobianMultivariateVectorOptimizer {
043        /** Square-root of the weight matrix. */
044        private RealMatrix weightMatrixSqrt;
045        /** Cost value (square root of the sum of the residuals). */
046        private double cost;
047    
048        /**
049         * @param checker Convergence checker.
050         */
051        protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
052            super(checker);
053        }
054    
055        /**
056         * Computes the weighted Jacobian matrix.
057         *
058         * @param params Model parameters at which to compute the Jacobian.
059         * @return the weighted Jacobian: W<sup>1/2</sup> J.
060         * @throws DimensionMismatchException if the Jacobian dimension does not
061         * match problem dimension.
062         */
063        protected RealMatrix computeWeightedJacobian(double[] params) {
064            return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params)));
065        }
066    
067        /**
068         * Computes the cost.
069         *
070         * @param residuals Residuals.
071         * @return the cost.
072         * @see #computeResiduals(double[])
073         */
074        protected double computeCost(double[] residuals) {
075            final ArrayRealVector r = new ArrayRealVector(residuals);
076            return FastMath.sqrt(r.dotProduct(getWeight().operate(r)));
077        }
078    
079        /**
080         * Gets the root-mean-square (RMS) value.
081         *
082         * The RMS the root of the arithmetic mean of the square of all weighted
083         * residuals.
084         * This is related to the criterion that is minimized by the optimizer
085         * as follows: If <em>c</em> if the criterion, and <em>n</em> is the
086         * number of measurements, then the RMS is <em>sqrt (c/n)</em>.
087         *
088         * @return the RMS value.
089         */
090        public double getRMS() {
091            return FastMath.sqrt(getChiSquare() / getTargetSize());
092        }
093    
094        /**
095         * Get a Chi-Square-like value assuming the N residuals follow N
096         * distinct normal distributions centered on 0 and whose variances are
097         * the reciprocal of the weights.
098         * @return chi-square value
099         */
100        public double getChiSquare() {
101            return cost * cost;
102        }
103    
104        /**
105         * Gets the square-root of the weight matrix.
106         *
107         * @return the square-root of the weight matrix.
108         */
109        public RealMatrix getWeightSquareRoot() {
110            return weightMatrixSqrt.copy();
111        }
112    
113        /**
114         * Sets the cost.
115         *
116         * @param cost Cost value.
117         */
118        protected void setCost(double cost) {
119            this.cost = cost;
120        }
121    
122        /**
123         * Get the covariance matrix of the optimized parameters.
124         * <br/>
125         * Note that this operation involves the inversion of the
126         * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
127         * Jacobian matrix.
128         * The {@code threshold} parameter is a way for the caller to specify
129         * that the result of this computation should be considered meaningless,
130         * and thus trigger an exception.
131         *
132         * @param params Model parameters.
133         * @param threshold Singularity threshold.
134         * @return the covariance matrix.
135         * @throws org.apache.commons.math3.linear.SingularMatrixException
136         * if the covariance matrix cannot be computed (singular problem).
137         */
138        public double[][] computeCovariances(double[] params,
139                                             double threshold) {
140            // Set up the Jacobian.
141            final RealMatrix j = computeWeightedJacobian(params);
142    
143            // Compute transpose(J)J.
144            final RealMatrix jTj = j.transpose().multiply(j);
145    
146            // Compute the covariances matrix.
147            final DecompositionSolver solver
148                = new QRDecomposition(jTj, threshold).getSolver();
149            return solver.getInverse().getData();
150        }
151    
152        /**
153         * Computes an estimate of the standard deviation of the parameters. The
154         * returned values are the square root of the diagonal coefficients of the
155         * covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]}
156         * is the optimized value of the {@code i}-th parameter, and {@code C} is
157         * the covariance matrix.
158         *
159         * @param params Model parameters.
160         * @param covarianceSingularityThreshold Singularity threshold (see
161         * {@link #computeCovariances(double[],double) computeCovariances}).
162         * @return an estimate of the standard deviation of the optimized parameters
163         * @throws org.apache.commons.math3.linear.SingularMatrixException
164         * if the covariance matrix cannot be computed.
165         */
166        public double[] computeSigma(double[] params,
167                                     double covarianceSingularityThreshold) {
168            final int nC = params.length;
169            final double[] sig = new double[nC];
170            final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
171            for (int i = 0; i < nC; ++i) {
172                sig[i] = FastMath.sqrt(cov[i][i]);
173            }
174            return sig;
175        }
176    
177        /**
178         * {@inheritDoc}
179         *
180         * @param optData Optimization data. The following data will be looked for:
181         * <ul>
182         *  <li>{@link org.apache.commons.math3.optim.MaxEval}</li>
183         *  <li>{@link org.apache.commons.math3.optim.InitialGuess}</li>
184         *  <li>{@link org.apache.commons.math3.optim.SimpleBounds}</li>
185         *  <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Target}</li>
186         *  <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Weight}</li>
187         *  <li>{@link org.apache.commons.math3.optim.nonlinear.vector.ModelFunction}</li>
188         *  <li>{@link org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian}</li>
189         * </ul>
190         * @return {@inheritDoc}
191         * @throws TooManyEvaluationsException if the maximal number of
192         * evaluations is exceeded.
193         * @throws DimensionMismatchException if the initial guess, target, and weight
194         * arguments have inconsistent dimensions.
195         */
196        @Override
197        public PointVectorValuePair optimize(OptimizationData... optData)
198            throws TooManyEvaluationsException {
199            // Retrieve settings.
200            parseOptimizationData(optData);
201            // Set up base class and perform computation.
202            return super.optimize(optData);
203        }
204    
205        /**
206         * Computes the residuals.
207         * The residual is the difference between the observed (target)
208         * values and the model (objective function) value.
209         * There is one residual for each element of the vector-valued
210         * function.
211         *
212         * @param objectiveValue Value of the the objective function. This is
213         * the value returned from a call to
214         * {@link #computeObjectiveValue(double[]) computeObjectiveValue}
215         * (whose array argument contains the model parameters).
216         * @return the residuals.
217         * @throws DimensionMismatchException if {@code params} has a wrong
218         * length.
219         */
220        protected double[] computeResiduals(double[] objectiveValue) {
221            final double[] target = getTarget();
222            if (objectiveValue.length != target.length) {
223                throw new DimensionMismatchException(target.length,
224                                                     objectiveValue.length);
225            }
226    
227            final double[] residuals = new double[target.length];
228            for (int i = 0; i < target.length; i++) {
229                residuals[i] = target[i] - objectiveValue[i];
230            }
231    
232            return residuals;
233        }
234    
235        /**
236         * Scans the list of (required and optional) optimization data that
237         * characterize the problem.
238         * If the weight matrix is specified, the {@link #weightMatrixSqrt}
239         * field is recomputed.
240         *
241         * @param optData Optimization data. The following data will be looked for:
242         * <ul>
243         *  <li>{@link Weight}</li>
244         * </ul>
245         */
246        private void parseOptimizationData(OptimizationData... optData) {
247            // The existing values (as set by the previous call) are reused if
248            // not provided in the argument list.
249            for (OptimizationData data : optData) {
250                if (data instanceof Weight) {
251                    weightMatrixSqrt = squareRoot(((Weight) data).getWeight());
252                    // If more data must be parsed, this statement _must_ be
253                    // changed to "continue".
254                    break;
255                }
256            }
257        }
258    
259        /**
260         * Computes the square-root of the weight matrix.
261         *
262         * @param m Symmetric, positive-definite (weight) matrix.
263         * @return the square-root of the weight matrix.
264         */
265        private RealMatrix squareRoot(RealMatrix m) {
266            final EigenDecomposition dec = new EigenDecomposition(m);
267            return dec.getSquareRoot();
268        }
269    }