    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    }