AbstractLeastSquaresOptimizer.java

  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. package org.apache.commons.math3.optim.nonlinear.vector.jacobian;

  19. import org.apache.commons.math3.exception.DimensionMismatchException;
  20. import org.apache.commons.math3.exception.TooManyEvaluationsException;
  21. import org.apache.commons.math3.linear.ArrayRealVector;
  22. import org.apache.commons.math3.linear.RealMatrix;
  23. import org.apache.commons.math3.linear.DiagonalMatrix;
  24. import org.apache.commons.math3.linear.DecompositionSolver;
  25. import org.apache.commons.math3.linear.MatrixUtils;
  26. import org.apache.commons.math3.linear.QRDecomposition;
  27. import org.apache.commons.math3.linear.EigenDecomposition;
  28. import org.apache.commons.math3.optim.OptimizationData;
  29. import org.apache.commons.math3.optim.ConvergenceChecker;
  30. import org.apache.commons.math3.optim.PointVectorValuePair;
  31. import org.apache.commons.math3.optim.nonlinear.vector.Weight;
  32. import org.apache.commons.math3.optim.nonlinear.vector.JacobianMultivariateVectorOptimizer;
  33. import org.apache.commons.math3.util.FastMath;

  35. /**
  36.  * Base class for implementing least-squares optimizers.
  37.  * It provides methods for error estimation.
  38.  *
  39.  * @since 3.1
  40.  * @deprecated All classes and interfaces in this package are deprecated.
  41.  * The optimizers that were provided here were moved to the
  42.  * {@link org.apache.commons.math3.fitting.leastsquares} package
  43.  * (cf. MATH-1008).
  44.  */
  45. @Deprecated
  46. public abstract class AbstractLeastSquaresOptimizer
  47.     extends JacobianMultivariateVectorOptimizer {
  48.     /** Square-root of the weight matrix. */
  49.     private RealMatrix weightMatrixSqrt;
  50.     /** Cost value (square root of the sum of the residuals). */
  51.     private double cost;

  53.     /**
  54.      * @param checker Convergence checker.
  55.      */
  56.     protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
  57.         super(checker);
  58.     }

  60.     /**
  61.      * Computes the weighted Jacobian matrix.
  62.      *
  63.      * @param params Model parameters at which to compute the Jacobian.
  64.      * @return the weighted Jacobian: W<sup>1/2</sup> J.
  65.      * @throws DimensionMismatchException if the Jacobian dimension does not
  66.      * match problem dimension.
  67.      */
  68.     protected RealMatrix computeWeightedJacobian(double[] params) {
  69.         return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params)));
  70.     }

  72.     /**
  73.      * Computes the cost.
  74.      *
  75.      * @param residuals Residuals.
  76.      * @return the cost.
  77.      * @see #computeResiduals(double[])
  78.      */
  79.     protected double computeCost(double[] residuals) {
  80.         final ArrayRealVector r = new ArrayRealVector(residuals);
  81.         return FastMath.sqrt(r.dotProduct(getWeight().operate(r)));
  82.     }

  84.     /**
  85.      * Gets the root-mean-square (RMS) value.
  86.      *
  87.      * The RMS the root of the arithmetic mean of the square of all weighted
  88.      * residuals.
  89.      * This is related to the criterion that is minimized by the optimizer
  90.      * as follows: If <em>c</em> if the criterion, and <em>n</em> is the
  91.      * number of measurements, then the RMS is <em>sqrt (c/n)</em>.
  92.      *
  93.      * @return the RMS value.
  94.      */
  95.     public double getRMS() {
  96.         return FastMath.sqrt(getChiSquare() / getTargetSize());
  97.     }

  99.     /**
  100.      * Get a Chi-Square-like value assuming the N residuals follow N
  101.      * distinct normal distributions centered on 0 and whose variances are
  102.      * the reciprocal of the weights.
  103.      * @return chi-square value
  104.      */
  105.     public double getChiSquare() {
  106.         return cost * cost;
  107.     }

  109.     /**
  110.      * Gets the square-root of the weight matrix.
  111.      *
  112.      * @return the square-root of the weight matrix.
  113.      */
  114.     public RealMatrix getWeightSquareRoot() {
  115.         return weightMatrixSqrt.copy();
  116.     }

  118.     /**
  119.      * Sets the cost.
  120.      *
  121.      * @param cost Cost value.
  122.      */
  123.     protected void setCost(double cost) {
  124.         this.cost = cost;
  125.     }

  127.     /**
  128.      * Get the covariance matrix of the optimized parameters.
  129.      * <br/>
  130.      * Note that this operation involves the inversion of the
  131.      * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
  132.      * Jacobian matrix.
  133.      * The {@code threshold} parameter is a way for the caller to specify
  134.      * that the result of this computation should be considered meaningless,
  135.      * and thus trigger an exception.
  136.      *
  137.      * @param params Model parameters.
  138.      * @param threshold Singularity threshold.
  139.      * @return the covariance matrix.
  140.      * @throws org.apache.commons.math3.linear.SingularMatrixException
  141.      * if the covariance matrix cannot be computed (singular problem).
  142.      */
  143.     public double[][] computeCovariances(double[] params,
  144.                                          double threshold) {
  145.         // Set up the Jacobian.
  146.         final RealMatrix j = computeWeightedJacobian(params);

  148.         // Compute transpose(J)J.
  149.         final RealMatrix jTj = j.transpose().multiply(j);

  151.         // Compute the covariances matrix.
  152.         final DecompositionSolver solver
  153.             = new QRDecomposition(jTj, threshold).getSolver();
  154.         return solver.getInverse().getData();
  155.     }

  157.     /**
  158.      * Computes an estimate of the standard deviation of the parameters. The
  159.      * returned values are the square root of the diagonal coefficients of the
  160.      * covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]}
  161.      * is the optimized value of the {@code i}-th parameter, and {@code C} is
  162.      * the covariance matrix.
  163.      *
  164.      * @param params Model parameters.
  165.      * @param covarianceSingularityThreshold Singularity threshold (see
  166.      * {@link #computeCovariances(double[],double) computeCovariances}).
  167.      * @return an estimate of the standard deviation of the optimized parameters
  168.      * @throws org.apache.commons.math3.linear.SingularMatrixException
  169.      * if the covariance matrix cannot be computed.
  170.      */
  171.     public double[] computeSigma(double[] params,
  172.                                  double covarianceSingularityThreshold) {
  173.         final int nC = params.length;
  174.         final double[] sig = new double[nC];
  175.         final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
  176.         for (int i = 0; i < nC; ++i) {
  177.             sig[i] = FastMath.sqrt(cov[i][i]);
  178.         }
  179.         return sig;
  180.     }

  182.     /**
  183.      * {@inheritDoc}
  184.      *
  185.      * @param optData Optimization data. In addition to those documented in
  186.      * {@link JacobianMultivariateVectorOptimizer#parseOptimizationData(OptimizationData[])
  187.      * JacobianMultivariateVectorOptimizer}, this method will register the following data:
  188.      * <ul>
  189.      *  <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Weight}</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.         // Set up base class and perform computation.
  201.         return super.optimize(optData);
  202.     }

  204.     /**
  205.      * Computes the residuals.
  206.      * The residual is the difference between the observed (target)
  207.      * values and the model (objective function) value.
  208.      * There is one residual for each element of the vector-valued
  209.      * function.
  210.      *
  211.      * @param objectiveValue Value of the the objective function. This is
  212.      * the value returned from a call to
  213.      * {@link #computeObjectiveValue(double[]) computeObjectiveValue}
  214.      * (whose array argument contains the model parameters).
  215.      * @return the residuals.
  216.      * @throws DimensionMismatchException if {@code params} has a wrong
  217.      * length.
  218.      */
  219.     protected double[] computeResiduals(double[] objectiveValue) {
  220.         final double[] target = getTarget();
  221.         if (objectiveValue.length != target.length) {
  222.             throw new DimensionMismatchException(target.length,
  223.                                                  objectiveValue.length);
  224.         }

  226.         final double[] residuals = new double[target.length];
  227.         for (int i = 0; i < target.length; i++) {
  228.             residuals[i] = target[i] - objectiveValue[i];
  229.         }

  231.         return residuals;
  232.     }

  234.     /**
  235.      * Scans the list of (required and optional) optimization data that
  236.      * characterize the problem.
  237.      * If the weight matrix is specified, the {@link #weightMatrixSqrt}
  238.      * field is recomputed.
  239.      *
  240.      * @param optData Optimization data. The following data will be looked for:
  241.      * <ul>
  242.      *  <li>{@link Weight}</li>
  243.      * </ul>
  244.      */
  245.     @Override
  246.     protected void parseOptimizationData(OptimizationData... optData) {
  247.         // Allow base class to register its own data.
  248.         super.parseOptimizationData(optData);

  250.         // The existing values (as set by the previous call) are reused if
  251.         // not provided in the argument list.
  252.         for (OptimizationData data : optData) {
  253.             if (data instanceof Weight) {
  254.                 weightMatrixSqrt = squareRoot(((Weight) data).getWeight());
  255.                 // If more data must be parsed, this statement _must_ be
  256.                 // changed to "continue".
  257.                 break;
  258.             }
  259.         }
  260.     }

  262.     /**
  263.      * Computes the square-root of the weight matrix.
  264.      *
  265.      * @param m Symmetric, positive-definite (weight) matrix.
  266.      * @return the square-root of the weight matrix.
  267.      */
  268.     private RealMatrix squareRoot(RealMatrix m) {
  269.         if (m instanceof DiagonalMatrix) {
  270.             final int dim = m.getRowDimension();
  271.             final RealMatrix sqrtM = new DiagonalMatrix(dim);
  272.             for (int i = 0; i < dim; i++) {
  273.                 sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
  274.             }
  275.             return sqrtM;
  276.         } else {
  277.             final EigenDecomposition dec = new EigenDecomposition(m);
  278.             return dec.getSquareRoot();
  279.         }
  280.     }
  281. }
Created with JaCoCo 0.7.2.201409121644