diff --git a/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java b/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java index 72d932f6b..7f0b3421c 100644 --- a/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java +++ b/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java @@ -26,7 +26,6 @@ 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.Array2DRowRealMatrix; import org.apache.commons.math3.linear.DecompositionSolver; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.QRDecomposition; @@ -61,7 +60,11 @@ import org.apache.commons.math3.util.FastMath; public abstract class AbstractLeastSquaresOptimizer extends BaseAbstractMultivariateVectorOptimizer implements DifferentiableMultivariateVectorOptimizer { - /** Singularity threshold (cf. {@link #getCovariances(double)}). */ + /** + * 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. @@ -69,19 +72,41 @@ public abstract class AbstractLeastSquaresOptimizer * {@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 #computeJacobian(double[])} instead. */ + @Deprecated protected double[][] weightedResidualJacobian; - /** Number of columns of the jacobian matrix. */ + /** Number of columns of the jacobian matrix. + * @deprecated As of 3.1. + */ + @Deprecated protected int cols; - /** Number of rows of the jacobian matrix. */ + /** Number of rows of the jacobian matrix. + * @deprecated As of 3.1. + */ + @Deprecated protected int rows; - /** Current point. */ + /** Current point. + * @deprecated As of 3.1. + */ + @Deprecated protected double[] point; - /** Current objective function value. */ + /** Current objective function value. + * @deprecated As of 3.1. + */ + @Deprecated protected double[] objective; - /** Weighted residuals */ + /** Weighted residuals + * @deprecated As of 3.1. + */ + @Deprecated protected double[] weightedResiduals; - /** Cost value (square root of the sum of the residuals). */ + /** 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; @@ -118,23 +143,39 @@ public abstract class AbstractLeastSquaresOptimizer * * @throws DimensionMismatchException if the Jacobian dimension does not * match problem dimension. + * @deprecated As of 3.1. Please use {@link #computeJacobian(double[])} + * instead. */ + @Deprecated protected void updateJacobian() { + computeJacobian(point); + } + + /** + * Computes the Jacobian matrix. + * + * @param params Model parameters at which to compute the Jacobian. + * @return the weighted Jacobian: -(W1/2 J). + * @throws DimensionMismatchException if the Jacobian dimension does not + * match problem dimension. + * @since 3.1 + */ + protected RealMatrix computeJacobian(double[] params) { ++jacobianEvaluations; - DerivativeStructure[] dsPoint = new DerivativeStructure[point.length]; - for (int i = 0; i < point.length; ++i) { - dsPoint[i] = new DerivativeStructure(point.length, 1, i, point[i]); - } - DerivativeStructure[] dsValue = jF.value(dsPoint); - if (dsValue.length != rows) { - throw new DimensionMismatchException(dsValue.length, rows); + 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; - final int nC = point.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[point.length]; + int[] orders = new int[nC]; for (int j = 0; j < nC; ++j) { orders[j] = 1; jacobianData[i][j] = dsValue[i].getPartialDerivative(orders); @@ -142,8 +183,15 @@ public abstract class AbstractLeastSquaresOptimizer } } - weightedResidualJacobian - = weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(jacobianData)).scalarMultiply(-1).getData(); + // XXX What is the purpose of the multiplication by -1? + final RealMatrix weightedJacobian + = weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(jacobianData)).scalarMultiply(-1); + + // XXX For backwards-compatibility (field "weightedResidualJacobian" + // must be removed in 4.0). + weightedResidualJacobian = weightedJacobian.getData(); + + return weightedJacobian; } /** @@ -152,18 +200,36 @@ public abstract class AbstractLeastSquaresOptimizer * 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[])} and {@link #computeCost(double[])} + * instead. */ + @Deprecated protected void updateResidualsAndCost() { - final double[] res = computeResidual(point); - final ArrayRealVector residuals = new ArrayRealVector(res); - final RealMatrix weight = getWeight(); + objective = computeObjectiveValue(point); + final double[] res = computeResiduals(objective); // Compute cost. - cost = FastMath.sqrt(residuals.dotProduct(weight.operate(residuals))); - // Compute weighted residuals. + cost = computeCost(res); + + // Compute weighted residuals. XXX To be moved to "LevenbergMarquardtOptimizer". + 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 @@ -188,15 +254,27 @@ public abstract class AbstractLeastSquaresOptimizer return cost * cost; } + /** + * Sets the cost. + * + * @param cost Cost value. + * @since 3.1 + */ + public 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); } @@ -215,14 +293,37 @@ public abstract class AbstractLeastSquaresOptimizer * @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) { - // Set up the jacobian. - updateJacobian(); + return computeCovariances(point, threshold); + } + + /** + * Get the covariance matrix of the optimized parameters. + *
+ * Note that this operation involves the inversion of the + * 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). + */ + public double[][] computeCovariances(double[] params, + double threshold) { + // Set up the Jacobian. + final RealMatrix j = computeJacobian(params); // Compute transpose(J)J. - final RealMatrix wrj = new Array2DRowRealMatrix(weightedResidualJacobian); - final RealMatrix jTj = wrj.transpose().multiply(wrj); + final RealMatrix jTj = j.transpose().multiply(j); // Compute the covariances matrix. final DecompositionSolver solver @@ -255,9 +356,9 @@ public abstract class AbstractLeastSquaresOptimizer * @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 #getSigma()} should be used - * instead. It should be emphasized that {@link #guessParametersErrors()} and - * {@link #getSigma()} are not strictly equivalent. + * @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 not strictly equivalent. */ @Deprecated public double[] guessParametersErrors() { @@ -267,7 +368,7 @@ public abstract class AbstractLeastSquaresOptimizer } double[] errors = new double[cols]; final double c = FastMath.sqrt(getChiSquare() / (rows - cols)); - double[][] covar = getCovariances(); + double[][] covar = computeCovariances(point, 1e-14); for (int i = 0; i < errors.length; ++i) { errors[i] = FastMath.sqrt(covar[i][i]) * c; } @@ -275,22 +376,25 @@ public abstract class AbstractLeastSquaresOptimizer } /** - *

- * Returns an estimate of the standard deviation of the parameters. The + * 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. */ - public double[] getSigma() { - final double[] sig = new double[cols]; - final double[][] cov = getCovariances(); - for (int i = 0; i < sig.length; ++i) { + 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; @@ -409,6 +513,36 @@ public abstract class AbstractLeastSquaresOptimizer 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. + */ + 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. * @@ -419,34 +553,4 @@ public abstract class AbstractLeastSquaresOptimizer final EigenDecomposition dec = new EigenDecomposition(m); return dec.getSquareRoot(); } - - /** - * Computes the residuals. - * The residual is the difference between the observed (target) - * values and the model (objective function) value, for the given - * parameters. - * There is one residual for each element of the vector-valued - * function. - * - * @param params Parameters of the model. - * @return the residuals. - * @throws DimensionMismatchException if {@code params} has a wrong - * length. - */ - private double[] computeResidual(double[] params) { - if (params.length != getStartPoint().length) { - throw new DimensionMismatchException(params.length, - getStartPoint().length); - } - - objective = computeObjectiveValue(params); - final double[] target = getTarget(); - - final double[] residuals = new double[target.length]; - for (int i = 0; i < target.length; i++) { - residuals[i] = target[i] - objective[i]; - } - - return residuals; - } } diff --git a/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java b/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java index 685757048..9a5abf43e 100644 --- a/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java +++ b/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java @@ -364,14 +364,11 @@ public abstract class AbstractLeastSquaresOptimizerAbstractTest { Assert.assertEquals(69.96016176931406, circle.getRadius(center), 1.0e-6); Assert.assertEquals(96.07590211815305, center.getX(), 1.0e-6); Assert.assertEquals(48.13516790438953, center.getY(), 1.0e-6); - double[][] cov = optimizer.getCovariances(); + double[][] cov = optimizer.computeCovariances(optimum.getPoint(), 1e-14); Assert.assertEquals(1.839, cov[0][0], 0.001); Assert.assertEquals(0.731, cov[0][1], 0.001); Assert.assertEquals(cov[0][1], cov[1][0], 1.0e-14); Assert.assertEquals(0.786, cov[1][1], 0.001); - double[] errors = optimizer.guessParametersErrors(); - Assert.assertEquals(1.384, errors[0], 0.001); - Assert.assertEquals(0.905, errors[1], 0.001); // add perfect measurements and check errors are reduced double r = circle.getRadius(center); @@ -382,15 +379,12 @@ public abstract class AbstractLeastSquaresOptimizerAbstractTest { Arrays.fill(target, 0.0); double[] weights = new double[circle.getN()]; Arrays.fill(weights, 2.0); - optimizer.optimize(100, circle, target, weights, new double[] { 98.680, 47.345 }); - cov = optimizer.getCovariances(); + optimum = optimizer.optimize(100, circle, target, weights, new double[] { 98.680, 47.345 }); + cov = optimizer.computeCovariances(optimum.getPoint(), 1e-14); Assert.assertEquals(0.0016, cov[0][0], 0.001); Assert.assertEquals(3.2e-7, cov[0][1], 1.0e-9); Assert.assertEquals(cov[0][1], cov[1][0], 1.0e-14); Assert.assertEquals(0.0016, cov[1][1], 0.001); - errors = optimizer.guessParametersErrors(); - Assert.assertEquals(0.004, errors[0], 0.001); - Assert.assertEquals(0.004, errors[1], 0.001); } @Test @@ -481,16 +475,11 @@ public abstract class AbstractLeastSquaresOptimizerAbstractTest { optimum = optimizer.optimize(100, problem, data[1], w, initial); final double[] actual = optimum.getPoint(); - final double[] actualSig = optimizer.guessParametersErrors(); for (int i = 0; i < actual.length; i++) { double expected = dataset.getParameter(i); double delta = FastMath.abs(errParams * expected); Assert.assertEquals(dataset.getName() + ", param #" + i, expected, actual[i], delta); - expected = dataset.getParameterStandardDeviation(i); - delta = FastMath.abs(errParamsSd * expected); - Assert.assertEquals(dataset.getName() + ", sd of param #" + i, - expected, actualSig[i], delta); } } diff --git a/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerTest.java b/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerTest.java index 8c95ba621..64ace2f72 100644 --- a/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerTest.java +++ b/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerTest.java @@ -29,9 +29,10 @@ public class AbstractLeastSquaresOptimizerTest { @Override protected PointVectorValuePair doOptimize() { - updateResidualsAndCost(); - updateJacobian(); - return null; + final double[] params = getStartPoint(); + final double[] res = computeResiduals(computeObjectiveValue(params)); + setCost(computeCost(res)); + return new PointVectorValuePair(params, null); } }; } @@ -75,7 +76,7 @@ public class AbstractLeastSquaresOptimizerTest { } @Test - public void testGetSigma() throws IOException { + public void testComputeSigma() throws IOException { final StatisticalReferenceDataset dataset; dataset = StatisticalReferenceDatasetFactory.createKirby2(); final AbstractLeastSquaresOptimizer optimizer; @@ -85,14 +86,14 @@ public class AbstractLeastSquaresOptimizerTest { final double[] w = new double[y.length]; Arrays.fill(w, 1.0); - final int dof = y.length-a.length; - optimizer.optimize(1, dataset.getLeastSquaresProblem(), y, w, a); - final double[] sig = optimizer.getSigma(); + final int dof = y.length - a.length; + final PointVectorValuePair optimum = optimizer.optimize(1, dataset.getLeastSquaresProblem(), y, w, a); + final double[] sig = optimizer.computeSigma(optimum.getPoint(), 1e-14); final double[] expected = dataset.getParametersStandardDeviations(); for (int i = 0; i < sig.length; i++) { - final double actual = FastMath.sqrt(optimizer.getChiSquare()/dof)*sig[i]; + final double actual = FastMath.sqrt(optimizer.getChiSquare() / dof) * sig[i]; Assert.assertEquals(dataset.getName() + ", parameter #" + i, - expected[i], actual, 1.3e-8 * expected[i]); + expected[i], actual, 1e-7 * expected[i]); } } } diff --git a/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerTestValidation.java b/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerTestValidation.java index 532f443f6..34dc3ddd8 100644 --- a/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerTestValidation.java +++ b/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerTestValidation.java @@ -115,9 +115,9 @@ public class AbstractLeastSquaresOptimizerTestValidation { // Estimation of the standard deviation (diagonal elements of the // covariance matrix). - optim.optimize(Integer.MAX_VALUE, + final PointVectorValuePair optimum = optim.optimize(Integer.MAX_VALUE, problem, problem.target(), problem.weight(), init); - final double[] sigma = optim.getSigma(); + final double[] sigma = optim.computeSigma(optimum.getPoint(), 1e-14); // Accumulate statistics. for (int i = 0; i < numParams; i++) { @@ -220,7 +220,7 @@ public class AbstractLeastSquaresOptimizerTestValidation { // Get chi-square of the best parameters set for the given set of // observations. final double bestChi2N = getChi2N(optim, problem, regress); - final double[] sigma = optim.getSigma(); + final double[] sigma = optim.computeSigma(regress, 1e-14); // Monte-Carlo (generates a grid of parameters). final int mcRepeat = MONTE_CARLO_RUNS; @@ -312,10 +312,9 @@ class DummyOptimizer extends AbstractLeastSquaresOptimizer { */ @Override public PointVectorValuePair doOptimize() { - // In order to be able to access the chi-square. - updateResidualsAndCost(); - - // Dummy value. - return null; + final double[] params = getStartPoint(); + final double[] res = computeResiduals(computeObjectiveValue(params)); + setCost(computeCost(res)); + return new PointVectorValuePair(params, null); } }