Java example source code file (AbstractLeastSquaresOptimizer.java)

This example Java source code file (AbstractLeastSquaresOptimizer.java) is included in the alvinalexander.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" ^TM.

Learn more about this Java project at its project page.

Java - Java tags/keywords

abstractleastsquaresoptimizer, arrayrealvector, deprecated, derivativestructure, diagonalmatrix, differentiablemultivariatevectoroptimizer, eigendecomposition, initialguess, multivariatedifferentiablevectorfunction, override, pointvectorvaluepair, realmatrix, target, weight

The AbstractLeastSquaresOptimizer.java Java example source code

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math3.optimization.general;

import org.apache.commons.math3.analysis.DifferentiableMultivariateVectorFunction;
import org.apache.commons.math3.analysis.FunctionUtils;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.analysis.differentiation.MultivariateDifferentiableVectorFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.optimization.OptimizationData;
import org.apache.commons.math3.optimization.InitialGuess;
import org.apache.commons.math3.optimization.Target;
import org.apache.commons.math3.optimization.Weight;
import org.apache.commons.math3.optimization.ConvergenceChecker;
import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer;
import org.apache.commons.math3.optimization.PointVectorValuePair;
import org.apache.commons.math3.optimization.direct.BaseAbstractMultivariateVectorOptimizer;
import org.apache.commons.math3.util.FastMath;

/**
 * Base class for implementing least squares optimizers.
 * It handles the boilerplate methods associated to thresholds settings,
 * Jacobian and error estimation.
 * <br/>
 * This class constructs the Jacobian matrix of the function argument in method
 * {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,
 * org.apache.commons.math3.analysis.MultivariateVectorFunction,OptimizationData[])
 * optimize} and assumes that the rows of that matrix iterate on the model
 * functions while the columns iterate on the parameters; thus, the numbers
 * of rows is equal to the dimension of the
 * {@link org.apache.commons.math3.optimization.Target Target} while
 * the number of columns is equal to the dimension of the
 * {@link org.apache.commons.math3.optimization.InitialGuess InitialGuess}.
 *
 * @deprecated As of 3.1 (to be removed in 4.0).
 * @since 1.2
 */
@Deprecated
public abstract class AbstractLeastSquaresOptimizer
    extends BaseAbstractMultivariateVectorOptimizer<DifferentiableMultivariateVectorFunction>
    implements DifferentiableMultivariateVectorOptimizer {
    /**
     * Singularity threshold (cf. {@link #getCovariances(double)}).
     * @deprecated As of 3.1.
     */
    @Deprecated
    private static final double DEFAULT_SINGULARITY_THRESHOLD = 1e-14;
    /**
     * Jacobian matrix of the weighted residuals.
     * This matrix is in canonical form just after the calls to
     * {@link #updateJacobian()}, but may be modified by the solver
     * in the derived class (the {@link LevenbergMarquardtOptimizer
     * Levenberg-Marquardt optimizer} does this).
     * @deprecated As of 3.1. To be removed in 4.0. Please use
     * {@link #computeWeightedJacobian(double[])} instead.
     */
    @Deprecated
    protected double[][] weightedResidualJacobian;
    /** Number of columns of the jacobian matrix.
     * @deprecated As of 3.1.
     */
    @Deprecated
    protected int cols;
    /** Number of rows of the jacobian matrix.
     * @deprecated As of 3.1.
     */
    @Deprecated
    protected int rows;
    /** Current point.
     * @deprecated As of 3.1.
     */
    @Deprecated
    protected double[] point;
    /** Current objective function value.
     * @deprecated As of 3.1.
     */
    @Deprecated
    protected double[] objective;
    /** Weighted residuals
     * @deprecated As of 3.1.
     */
    @Deprecated
    protected double[] weightedResiduals;
    /** Cost value (square root of the sum of the residuals).
     * @deprecated As of 3.1. Field to become "private" in 4.0.
     * Please use {@link #setCost(double)}.
     */
    @Deprecated
    protected double cost;
    /** Objective function derivatives. */
    private MultivariateDifferentiableVectorFunction jF;
    /** Number of evaluations of the Jacobian. */
    private int jacobianEvaluations;
    /** Square-root of the weight matrix. */
    private RealMatrix weightMatrixSqrt;

    /**
     * Simple constructor with default settings.
     * The convergence check is set to a {@link
     * org.apache.commons.math3.optimization.SimpleVectorValueChecker}.
     * @deprecated See {@link org.apache.commons.math3.optimization.SimpleValueChecker#SimpleValueChecker()}
     */
    @Deprecated
    protected AbstractLeastSquaresOptimizer() {}

    /**
     * @param checker Convergence checker.
     */
    protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
        super(checker);
    }

    /**
     * @return the number of evaluations of the Jacobian function.
     */
    public int getJacobianEvaluations() {
        return jacobianEvaluations;
    }

    /**
     * Update the jacobian matrix.
     *
     * @throws DimensionMismatchException if the Jacobian dimension does not
     * match problem dimension.
     * @deprecated As of 3.1. Please use {@link #computeWeightedJacobian(double[])}
     * instead.
     */
    @Deprecated
    protected void updateJacobian() {
        final RealMatrix weightedJacobian = computeWeightedJacobian(point);
        weightedResidualJacobian = weightedJacobian.scalarMultiply(-1).getData();
    }

    /**
     * Computes the Jacobian matrix.
     *
     * @param params Model parameters at which to compute the Jacobian.
     * @return the weighted Jacobian: W<sup>1/2 J.
     * @throws DimensionMismatchException if the Jacobian dimension does not
     * match problem dimension.
     * @since 3.1
     */
    protected RealMatrix computeWeightedJacobian(double[] params) {
        ++jacobianEvaluations;

        final DerivativeStructure[] dsPoint = new DerivativeStructure[params.length];
        final int nC = params.length;
        for (int i = 0; i < nC; ++i) {
            dsPoint[i] = new DerivativeStructure(nC, 1, i, params[i]);
        }
        final DerivativeStructure[] dsValue = jF.value(dsPoint);
        final int nR = getTarget().length;
        if (dsValue.length != nR) {
            throw new DimensionMismatchException(dsValue.length, nR);
        }
        final double[][] jacobianData = new double[nR][nC];
        for (int i = 0; i < nR; ++i) {
            int[] orders = new int[nC];
            for (int j = 0; j < nC; ++j) {
                orders[j] = 1;
                jacobianData[i][j] = dsValue[i].getPartialDerivative(orders);
                orders[j] = 0;
            }
        }

        return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(jacobianData));
    }

    /**
     * Update the residuals array and cost function value.
     * @throws DimensionMismatchException if the dimension does not match the
     * problem dimension.
     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
     * if the maximal number of evaluations is exceeded.
     * @deprecated As of 3.1. Please use {@link #computeResiduals(double[])},
     * {@link #computeObjectiveValue(double[])}, {@link #computeCost(double[])}
     * and {@link #setCost(double)} instead.
     */
    @Deprecated
    protected void updateResidualsAndCost() {
        objective = computeObjectiveValue(point);
        final double[] res = computeResiduals(objective);

        // Compute cost.
        cost = computeCost(res);

        // Compute weighted residuals.
        final ArrayRealVector residuals = new ArrayRealVector(res);
        weightedResiduals = weightMatrixSqrt.operate(residuals).toArray();
    }

    /**
     * Computes the cost.
     *
     * @param residuals Residuals.
     * @return the cost.
     * @see #computeResiduals(double[])
     * @since 3.1
     */
    protected double computeCost(double[] residuals) {
        final ArrayRealVector r = new ArrayRealVector(residuals);
        return FastMath.sqrt(r.dotProduct(getWeight().operate(r)));
    }

    /**
     * Get the Root Mean Square value.
     * Get the Root Mean Square value, i.e. the root of the arithmetic
     * mean of the square of all weighted residuals. This is related to the
     * criterion that is minimized by the optimizer as follows: if
     * <em>c if the criterion, and n is the number of
     * measurements, then the RMS is <em>sqrt (c/n).
     *
     * @return RMS value
     */
    public double getRMS() {
        return FastMath.sqrt(getChiSquare() / rows);
    }

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

    /**
     * Gets the square-root of the weight matrix.
     *
     * @return the square-root of the weight matrix.
     * @since 3.1
     */
    public RealMatrix getWeightSquareRoot() {
        return weightMatrixSqrt.copy();
    }

    /**
     * Sets the cost.
     *
     * @param cost Cost value.
     * @since 3.1
     */
    protected void setCost(double cost) {
        this.cost = cost;
    }

    /**
     * Get the covariance matrix of the optimized parameters.
     *
     * @return the covariance matrix.
     * @throws org.apache.commons.math3.linear.SingularMatrixException
     * if the covariance matrix cannot be computed (singular problem).
     * @see #getCovariances(double)
     * @deprecated As of 3.1. Please use {@link #computeCovariances(double[],double)}
     * instead.
     */
    @Deprecated
    public double[][] getCovariances() {
        return getCovariances(DEFAULT_SINGULARITY_THRESHOLD);
    }

    /**
     * Get the covariance matrix of the optimized parameters.
     * <br/>
     * Note that this operation involves the inversion of the
     * <code>JTJ matrix, where {@code J} is the
     * Jacobian matrix.
     * The {@code threshold} parameter is a way for the caller to specify
     * that the result of this computation should be considered meaningless,
     * and thus trigger an exception.
     *
     * @param threshold Singularity threshold.
     * @return the covariance matrix.
     * @throws org.apache.commons.math3.linear.SingularMatrixException
     * if the covariance matrix cannot be computed (singular problem).
     * @deprecated As of 3.1. Please use {@link #computeCovariances(double[],double)}
     * instead.
     */
    @Deprecated
    public double[][] getCovariances(double threshold) {
        return computeCovariances(point, threshold);
    }

    /**
     * Get the covariance matrix of the optimized parameters.
     * <br/>
     * Note that this operation involves the inversion of the
     * <code>JTJ matrix, where {@code J} is the
     * Jacobian matrix.
     * The {@code threshold} parameter is a way for the caller to specify
     * that the result of this computation should be considered meaningless,
     * and thus trigger an exception.
     *
     * @param params Model parameters.
     * @param threshold Singularity threshold.
     * @return the covariance matrix.
     * @throws org.apache.commons.math3.linear.SingularMatrixException
     * if the covariance matrix cannot be computed (singular problem).
     * @since 3.1
     */
    public double[][] computeCovariances(double[] params,
                                         double threshold) {
        // Set up the Jacobian.
        final RealMatrix j = computeWeightedJacobian(params);

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

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

    /**
     * <p>
     * Returns an estimate of the standard deviation of each parameter. The
     * returned values are the so-called (asymptotic) standard errors on the
     * parameters, defined as {@code sd(a[i]) = sqrt(S / (n - m) * C[i][i])},
     * where {@code a[i]} is the optimized value of the {@code i}-th parameter,
     * {@code S} is the minimized value of the sum of squares objective function
     * (as returned by {@link #getChiSquare()}), {@code n} is the number of
     * observations, {@code m} is the number of parameters and {@code C} is the
     * covariance matrix.
     * </p>
     * <p>
     * See also
     * <a href="http://en.wikipedia.org/wiki/Least_squares">Wikipedia,
     * or
     * <a href="http://mathworld.wolfram.com/LeastSquaresFitting.html">MathWorld,
     * equations (34) and (35) for a particular case.
     * </p>
     *
     * @return an estimate of the standard deviation of the optimized parameters
     * @throws org.apache.commons.math3.linear.SingularMatrixException
     * if the covariance matrix cannot be computed.
     * @throws NumberIsTooSmallException if the number of degrees of freedom is not
     * positive, i.e. the number of measurements is less or equal to the number of
     * parameters.
     * @deprecated as of version 3.1, {@link #computeSigma(double[],double)} should be used
     * instead. It should be emphasized that {@code guessParametersErrors} and
     * {@code computeSigma} are <em>not strictly equivalent.
     */
    @Deprecated
    public double[] guessParametersErrors() {
        if (rows <= cols) {
            throw new NumberIsTooSmallException(LocalizedFormats.NO_DEGREES_OF_FREEDOM,
                                                rows, cols, false);
        }
        double[] errors = new double[cols];
        final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
        double[][] covar = computeCovariances(point, 1e-14);
        for (int i = 0; i < errors.length; ++i) {
            errors[i] = FastMath.sqrt(covar[i][i]) * c;
        }
        return errors;
    }

    /**
     * Computes an estimate of the standard deviation of the parameters. The
     * returned values are the square root of the diagonal coefficients of the
     * covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]}
     * is the optimized value of the {@code i}-th parameter, and {@code C} is
     * the covariance matrix.
     *
     * @param params Model parameters.
     * @param covarianceSingularityThreshold Singularity threshold (see
     * {@link #computeCovariances(double[],double) computeCovariances}).
     * @return an estimate of the standard deviation of the optimized parameters
     * @throws org.apache.commons.math3.linear.SingularMatrixException
     * if the covariance matrix cannot be computed.
     * @since 3.1
     */
    public double[] computeSigma(double[] params,
                                 double covarianceSingularityThreshold) {
        final int nC = params.length;
        final double[] sig = new double[nC];
        final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
        for (int i = 0; i < nC; ++i) {
            sig[i] = FastMath.sqrt(cov[i][i]);
        }
        return sig;
    }

    /** {@inheritDoc}
     * @deprecated As of 3.1. Please use
     * {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,
     * org.apache.commons.math3.analysis.MultivariateVectorFunction,OptimizationData[])
     * optimize(int,MultivariateDifferentiableVectorFunction,OptimizationData...)}
     * instead.
     */
    @Override
    @Deprecated
    public PointVectorValuePair optimize(int maxEval,
                                         final DifferentiableMultivariateVectorFunction f,
                                         final double[] target, final double[] weights,
                                         final double[] startPoint) {
        return optimizeInternal(maxEval,
                                FunctionUtils.toMultivariateDifferentiableVectorFunction(f),
                                new Target(target),
                                new Weight(weights),
                                new InitialGuess(startPoint));
    }

    /**
     * Optimize an objective function.
     * Optimization is considered to be a weighted least-squares minimization.
     * The cost function to be minimized is
     * <code>?weighti(objectivei - targeti)2
     *
     * @param f Objective function.
     * @param target Target value for the objective functions at optimum.
     * @param weights Weights for the least squares cost computation.
     * @param startPoint Start point for optimization.
     * @return the point/value pair giving the optimal value for objective
     * function.
     * @param maxEval Maximum number of function evaluations.
     * @throws org.apache.commons.math3.exception.DimensionMismatchException
     * if the start point dimension is wrong.
     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
     * if the maximal number of evaluations is exceeded.
     * @throws org.apache.commons.math3.exception.NullArgumentException if
     * any argument is {@code null}.
     * @deprecated As of 3.1. Please use
     * {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,
     * org.apache.commons.math3.analysis.MultivariateVectorFunction,OptimizationData[])
     * optimize(int,MultivariateDifferentiableVectorFunction,OptimizationData...)}
     * instead.
     */
    @Deprecated
    public PointVectorValuePair optimize(final int maxEval,
                                         final MultivariateDifferentiableVectorFunction f,
                                         final double[] target, final double[] weights,
                                         final double[] startPoint) {
        return optimizeInternal(maxEval, f,
                                new Target(target),
                                new Weight(weights),
                                new InitialGuess(startPoint));
    }

    /**
     * Optimize an objective function.
     * Optimization is considered to be a weighted least-squares minimization.
     * The cost function to be minimized is
     * <code>?weighti(objectivei - targeti)2
     *
     * @param maxEval Allowed number of evaluations of the objective function.
     * @param f Objective function.
     * @param optData Optimization data. The following data will be looked for:
     * <ul>
     *  <li>{@link Target}
     *  <li>{@link Weight}
     *  <li>{@link InitialGuess}
     * </ul>
     * @return the point/value pair giving the optimal value of the objective
     * function.
     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException if
     * the maximal number of evaluations is exceeded.
     * @throws DimensionMismatchException if the target, and weight arguments
     * have inconsistent dimensions.
     * @see BaseAbstractMultivariateVectorOptimizer#optimizeInternal(int,
     * org.apache.commons.math3.analysis.MultivariateVectorFunction,OptimizationData[])
     * @since 3.1
     * @deprecated As of 3.1. Override is necessary only until this class's generic
     * argument is changed to {@code MultivariateDifferentiableVectorFunction}.
     */
    @Deprecated
    protected PointVectorValuePair optimizeInternal(final int maxEval,
                                                    final MultivariateDifferentiableVectorFunction f,
                                                    OptimizationData... optData) {
        // XXX Conversion will be removed when the generic argument of the
        // base class becomes "MultivariateDifferentiableVectorFunction".
        return super.optimizeInternal(maxEval, FunctionUtils.toDifferentiableMultivariateVectorFunction(f), optData);
    }

    /** {@inheritDoc} */
    @Override
    protected void setUp() {
        super.setUp();

        // Reset counter.
        jacobianEvaluations = 0;

        // Square-root of the weight matrix.
        weightMatrixSqrt = squareRoot(getWeight());

        // Store least squares problem characteristics.
        // XXX The conversion won't be necessary when the generic argument of
        // the base class becomes "MultivariateDifferentiableVectorFunction".
        // XXX "jF" is not strictly necessary anymore but is currently more
        // efficient than converting the value returned from "getObjectiveFunction()"
        // every time it is used.
        jF = FunctionUtils.toMultivariateDifferentiableVectorFunction((DifferentiableMultivariateVectorFunction) getObjectiveFunction());

        // Arrays shared with "private" and "protected" methods.
        point = getStartPoint();
        rows = getTarget().length;
        cols = point.length;
    }

    /**
     * Computes the residuals.
     * The residual is the difference between the observed (target)
     * values and the model (objective function) value.
     * There is one residual for each element of the vector-valued
     * function.
     *
     * @param objectiveValue Value of the the objective function. This is
     * the value returned from a call to
     * {@link #computeObjectiveValue(double[]) computeObjectiveValue}
     * (whose array argument contains the model parameters).
     * @return the residuals.
     * @throws DimensionMismatchException if {@code params} has a wrong
     * length.
     * @since 3.1
     */
    protected double[] computeResiduals(double[] objectiveValue) {
        final double[] target = getTarget();
        if (objectiveValue.length != target.length) {
            throw new DimensionMismatchException(target.length,
                                                 objectiveValue.length);
        }

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

        return residuals;
    }

    /**
     * Computes the square-root of the weight matrix.
     *
     * @param m Symmetric, positive-definite (weight) matrix.
     * @return the square-root of the weight matrix.
     */
    private RealMatrix squareRoot(RealMatrix m) {
        if (m instanceof DiagonalMatrix) {
            final int dim = m.getRowDimension();
            final RealMatrix sqrtM = new DiagonalMatrix(dim);
            for (int i = 0; i < dim; i++) {
               sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
            }
            return sqrtM;
        } else {
            final EigenDecomposition dec = new EigenDecomposition(m);
            return dec.getSquareRoot();
        }
    }
}