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