From 6dc6cf5b2fba50eb602283101de4910787ad5579 Mon Sep 17 00:00:00 2001 From: adrian Date: Tue, 1 May 2018 14:06:50 -0400 Subject: [PATCH] MATH-1459: Create a way to automatically calculate a Jacobian matrix using a differentiator. --- ...rentiatorMultivariateJacobianFunction.java | 103 ++++++++++++ ...torVectorMultivariateJacobianFunction.java | 106 ++++++++++++ .../leastsquares/BevingtonProblem.java | 63 +++++++ ...ectorMultivariateJacobianFunctionTest.java | 154 ++++++++++++++++++ .../LevenbergMarquardtOptimizerTest.java | 64 -------- 5 files changed, 426 insertions(+), 64 deletions(-) create mode 100644 src/main/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorMultivariateJacobianFunction.java create mode 100644 src/main/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunction.java create mode 100644 src/test/java/org/apache/commons/math4/fitting/leastsquares/BevingtonProblem.java create mode 100644 src/test/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunctionTest.java diff --git a/src/main/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorMultivariateJacobianFunction.java b/src/main/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorMultivariateJacobianFunction.java new file mode 100644 index 000000000..313571d38 --- /dev/null +++ b/src/main/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorMultivariateJacobianFunction.java @@ -0,0 +1,103 @@ +/* + * 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.math4.fitting.leastsquares; + +import org.apache.commons.math4.analysis.MultivariateFunction; +import org.apache.commons.math4.analysis.UnivariateFunction; +import org.apache.commons.math4.analysis.differentiation.DerivativeStructure; +import org.apache.commons.math4.analysis.differentiation.UnivariateFunctionDifferentiator; +import org.apache.commons.math4.linear.Array2DRowRealMatrix; +import org.apache.commons.math4.linear.ArrayRealVector; +import org.apache.commons.math4.linear.RealMatrix; +import org.apache.commons.math4.linear.RealVector; +import org.apache.commons.math4.util.Pair; + +/** + * A MultivariateJacobianFunction (a thing that requires a derivative) + * combined with the thing that can find derivatives. + * + * Can be used with a LeastSquaresProblem, a LeastSquaresFactory, or a LeastSquaresBuilder. + * + * This version that works with MultivariateFunction + * @see DifferentiatorVectorMultivariateJacobianFunction for version that works with MultivariateVectorFunction + */ +public class DifferentiatorMultivariateJacobianFunction implements MultivariateJacobianFunction { + /** + * The input function to find a jacobian for. + */ + private final MultivariateFunction function; + /** + * The differentiator to use to find the jacobian. + */ + private final UnivariateFunctionDifferentiator differentiator; + + /** + * Build the jacobian function using a differentiator. + * + * @param function the function to turn into a jacobian + * @param differentiator the differentiator to find the derivative + * + * This version that works with MultivariateFunction + * @see DifferentiatorVectorMultivariateJacobianFunction for version that works with MultivariateVectorFunction + */ + public DifferentiatorMultivariateJacobianFunction(MultivariateFunction function, UnivariateFunctionDifferentiator differentiator) { + this.function = function; + this.differentiator = differentiator; + } + + /** + * @inheritDoc + */ + @Override + public Pair value(RealVector point) { + ArrayRealVector value = new ArrayRealVector(1); + value.setEntry(0, function.value(point.toArray())); + RealMatrix jacobian = new Array2DRowRealMatrix(1, point.getDimension()); + + for(int column = 0; column < point.getDimension(); column++) { + final int columnFinal = column; + double originalPoint = point.getEntry(column); + double partialDerivative = getPartialDerivative(testPoint -> { + + point.setEntry(columnFinal, testPoint); + + double testPointOutput = function.value(point.toArray()); + + point.setEntry(columnFinal, originalPoint); //set it back + + return testPointOutput; + }, originalPoint); + + jacobian.setEntry(0, column, partialDerivative); + } + + return new Pair<>(value, jacobian); + } + + /** + * Returns first order derivative for the function passed in using a differentiator + * @param univariateFunction the function to differentiate + * @param atParameterValue the point at which to differentiate it at + * @return the slope at that point + */ + private double getPartialDerivative(UnivariateFunction univariateFunction, double atParameterValue) { + return differentiator + .differentiate(univariateFunction) + .value(new DerivativeStructure(1, 1, 0, atParameterValue)) + .getPartialDerivative(1); + } +} diff --git a/src/main/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunction.java b/src/main/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunction.java new file mode 100644 index 000000000..2f8295aa8 --- /dev/null +++ b/src/main/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunction.java @@ -0,0 +1,106 @@ +/* + * 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.math4.fitting.leastsquares; + +import org.apache.commons.math4.analysis.MultivariateVectorFunction; +import org.apache.commons.math4.analysis.UnivariateVectorFunction; +import org.apache.commons.math4.analysis.differentiation.DerivativeStructure; +import org.apache.commons.math4.analysis.differentiation.UnivariateVectorFunctionDifferentiator; +import org.apache.commons.math4.linear.Array2DRowRealMatrix; +import org.apache.commons.math4.linear.ArrayRealVector; +import org.apache.commons.math4.linear.RealMatrix; +import org.apache.commons.math4.linear.RealVector; +import org.apache.commons.math4.util.Pair; + +/** + * A MultivariateJacobianFunction (a thing that requires a derivative) + * combined with the thing that can find derivatives. + * + * Can be used with a LeastSquaresProblem, a LeastSquaresFactory, or a LeastSquaresBuilder. + * + * This version that works with MultivariateVectorFunction + * @see DifferentiatorMultivariateJacobianFunction for version that works with MultivariateFunction + */ +public class DifferentiatorVectorMultivariateJacobianFunction implements MultivariateJacobianFunction { + /** + * The input function to find a jacobian for. + */ + private final MultivariateVectorFunction function; + /** + * The differentiator to use to find the jacobian. + */ + private final UnivariateVectorFunctionDifferentiator differentiator; + + /** + * Build the jacobian function using a differentiator. + * + * @param function the function to turn into a jacobian + * @param differentiator the differentiator to find the derivative + * + * This version that works with MultivariateFunction + * @see DifferentiatorVectorMultivariateJacobianFunction for version that works with MultivariateVectorFunction + */ + public DifferentiatorVectorMultivariateJacobianFunction(MultivariateVectorFunction function, UnivariateVectorFunctionDifferentiator differentiator) { + this.function = function; + this.differentiator = differentiator; + } + + /** + * @inheritDoc + */ + @Override + public Pair value(RealVector point) { + RealVector value = new ArrayRealVector(function.value(point.toArray())); + RealMatrix jacobian = new Array2DRowRealMatrix(value.getDimension(), point.getDimension()); + + for(int column = 0; column < point.getDimension(); column++) { + final int columnFinal = column; + double originalPoint = point.getEntry(column); + double[] partialDerivatives = getPartialDerivative(testPoint -> { + + point.setEntry(columnFinal, testPoint); + + double[] testPointValue = function.value(point.toArray()); + + point.setEntry(columnFinal, originalPoint); //set it back + + return testPointValue; + }, originalPoint); + + jacobian.setColumn(column, partialDerivatives); + } + + return new Pair<>(value, jacobian); + } + + /** + * Returns first order derivative for the function passed in using a differentiator + * @param univariateVectorFunction the function to differentiate + * @param atParameterValue the point at which to differentiate it at + * @return the slopes at that point + */ + private double[] getPartialDerivative(UnivariateVectorFunction univariateVectorFunction, double atParameterValue) { + DerivativeStructure[] derivatives = differentiator + .differentiate(univariateVectorFunction) + .value(new DerivativeStructure(1, 1, 0, atParameterValue)); + double[] derivativesOut = new double[derivatives.length]; + for(int index=0;index time; + private List count; + + public BevingtonProblem() { + time = new ArrayList<>(); + count = new ArrayList<>(); + } + + public void addPoint(double t, double c) { + time.add(t); + count.add(c); + } + + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + @Override + public double[] value(double[] params) { + double[] values = new double[time.size()]; + for (int i = 0; i < values.length; ++i) { + final double t = time.get(i); + values[i] = params[0] + + params[1] * FastMath.exp(-t / params[3]) + + params[2] * FastMath.exp(-t / params[4]); + } + return values; + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + @Override + public double[][] value(double[] params) { + double[][] jacobian = new double[time.size()][5]; + + for (int i = 0; i < jacobian.length; ++i) { + final double t = time.get(i); + jacobian[i][0] = 1; + + final double p3 = params[3]; + final double p4 = params[4]; + final double tOp3 = t / p3; + final double tOp4 = t / p4; + jacobian[i][1] = FastMath.exp(-tOp3); + jacobian[i][2] = FastMath.exp(-tOp4); + jacobian[i][3] = params[1] * FastMath.exp(-tOp3) * tOp3 / p3; + jacobian[i][4] = params[2] * FastMath.exp(-tOp4) * tOp4 / p4; + } + return jacobian; + } + }; + } +} diff --git a/src/test/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunctionTest.java b/src/test/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunctionTest.java new file mode 100644 index 000000000..3321ec6db --- /dev/null +++ b/src/test/java/org/apache/commons/math4/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunctionTest.java @@ -0,0 +1,154 @@ +/* + * 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.math4.fitting.leastsquares; + +import org.apache.commons.math4.analysis.differentiation.FiniteDifferencesDifferentiator; +import org.apache.commons.math4.analysis.differentiation.UnivariateVectorFunctionDifferentiator; +import org.apache.commons.math4.linear.DiagonalMatrix; +import org.apache.commons.math4.linear.RealMatrix; +import org.apache.commons.math4.linear.RealVector; +import org.apache.commons.math4.optim.SimpleVectorValueChecker; +import org.apache.commons.math4.util.FastMath; +import org.junit.Assert; +import org.junit.Test; + +public class DifferentiatorVectorMultivariateJacobianFunctionTest { + + private static final int POINTS = 20; + private static final double STEP_SIZE = 0.2; + + private final UnivariateVectorFunctionDifferentiator differentiator = new FiniteDifferencesDifferentiator(POINTS, STEP_SIZE); + private final LeastSquaresOptimizer optimizer = this.getOptimizer(); + + public LeastSquaresBuilder base() { + return new LeastSquaresBuilder() + .checkerPair(new SimpleVectorValueChecker(1e-6, 1e-6)) + .maxEvaluations(100) + .maxIterations(getMaxIterations()); + } + + public LeastSquaresBuilder builder(BevingtonProblem problem, boolean automatic){ + if(automatic) { + return base() + .model(new DifferentiatorVectorMultivariateJacobianFunction(problem.getModelFunction(), differentiator)); + } + else { + return base() + .model(problem.getModelFunction(), problem.getModelFunctionJacobian()); + } + } + + public int getMaxIterations() { + return 25; + } + + public LeastSquaresOptimizer getOptimizer() { + return new LevenbergMarquardtOptimizer(); + } + + /** + * Non-linear test case: fitting of decay curve (from Chapter 8 of + * Bevington's textbook, "Data reduction and analysis for the physical sciences"). + */ + @Test + public void testBevington() { + + // the analytical optimum to compare to + final LeastSquaresOptimizer.Optimum analyticalOptimum = findBevington(false); + final RealVector analyticalSolution = analyticalOptimum.getPoint(); + final RealMatrix analyticalCovarianceMatrix = analyticalOptimum.getCovariances(1e-14); + + // the automatic DifferentiatorVectorMultivariateJacobianFunction optimum + final LeastSquaresOptimizer.Optimum automaticOptimum = findBevington(true); + final RealVector automaticSolution = automaticOptimum.getPoint(); + final RealMatrix automaticCovarianceMatrix = automaticOptimum.getCovariances(1e-14); + + final int numParams = analyticalOptimum.getPoint().getDimension(); + + // Check that the automatic solution is within the reference error range. + for (int i = 0; i < numParams; i++) { + final double error = FastMath.sqrt(analyticalCovarianceMatrix.getEntry(i, i)); + Assert.assertEquals("Parameter " + i, analyticalSolution.getEntry(i), automaticSolution.getEntry(i), error); + } + + // Check that each entry of the computed covariance matrix is within 1% + // of the reference analytical matrix entry. + for (int i = 0; i < numParams; i++) { + for (int j = 0; j < numParams; j++) { + Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]", + analyticalCovarianceMatrix.getEntry(i, j), + automaticCovarianceMatrix.getEntry(i, j), + FastMath.abs(0.01 * analyticalCovarianceMatrix.getEntry(i, j))); + } + } + + // Check various measures of goodness-of-fit. + final double tol = 1e-40; + Assert.assertEquals(analyticalOptimum.getChiSquare(), automaticOptimum.getChiSquare(), tol); + Assert.assertEquals(analyticalOptimum.getCost(), automaticOptimum.getCost(), tol); + Assert.assertEquals(analyticalOptimum.getRMS(), automaticOptimum.getRMS(), tol); + Assert.assertEquals(analyticalOptimum.getReducedChiSquare(automaticOptimum.getPoint().getDimension()), automaticOptimum.getReducedChiSquare(automaticOptimum.getPoint().getDimension()), tol); + } + + /** + * Build the problem and return the optimum, doesn't actually test the results. + * + * Pass in if you want to test using analytical derivatives, + * or the automatic {@link DifferentiatorVectorMultivariateJacobianFunction} + * + * @param automatic automatic {@link DifferentiatorVectorMultivariateJacobianFunction}, as opposed to analytical + * @return the optimum for this test + */ + private LeastSquaresOptimizer.Optimum findBevington(boolean automatic) { + final double[][] dataPoints = { + // column 1 = times + { 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, + 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, + 315, 330, 345, 360, 375, 390, 405, 420, 435, 450, + 465, 480, 495, 510, 525, 540, 555, 570, 585, 600, + 615, 630, 645, 660, 675, 690, 705, 720, 735, 750, + 765, 780, 795, 810, 825, 840, 855, 870, 885, }, + // column 2 = measured counts + { 775, 479, 380, 302, 185, 157, 137, 119, 110, 89, + 74, 61, 66, 68, 48, 54, 51, 46, 55, 29, + 28, 37, 49, 26, 35, 29, 31, 24, 25, 35, + 24, 30, 26, 28, 21, 18, 20, 27, 17, 17, + 14, 17, 24, 11, 22, 17, 12, 10, 13, 16, + 9, 9, 14, 21, 17, 13, 12, 18, 10, }, + }; + final double[] start = {10, 900, 80, 27, 225}; + + final BevingtonProblem problem = new BevingtonProblem(); + + final int len = dataPoints[0].length; + final double[] weights = new double[len]; + for (int i = 0; i < len; i++) { + problem.addPoint(dataPoints[0][i], + dataPoints[1][i]); + + weights[i] = 1 / dataPoints[1][i]; + } + + return optimizer.optimize( + builder(problem, automatic) + .target(dataPoints[1]) + .weight(new DiagonalMatrix(weights)) + .start(start) + .build() + ); + } +} diff --git a/src/test/java/org/apache/commons/math4/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java b/src/test/java/org/apache/commons/math4/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java index d0ee97f3e..0be1116f4 100644 --- a/src/test/java/org/apache/commons/math4/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java +++ b/src/test/java/org/apache/commons/math4/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java @@ -17,15 +17,8 @@ package org.apache.commons.math4.fitting.leastsquares; -import org.apache.commons.math4.analysis.MultivariateMatrixFunction; -import org.apache.commons.math4.analysis.MultivariateVectorFunction; import org.apache.commons.math4.exception.DimensionMismatchException; import org.apache.commons.math4.exception.TooManyEvaluationsException; -import org.apache.commons.math4.fitting.leastsquares.LeastSquaresBuilder; -import org.apache.commons.math4.fitting.leastsquares.LeastSquaresOptimizer; -import org.apache.commons.math4.fitting.leastsquares.LeastSquaresProblem; -import org.apache.commons.math4.fitting.leastsquares.LevenbergMarquardtOptimizer; -import org.apache.commons.math4.fitting.leastsquares.ParameterValidator; import org.apache.commons.math4.fitting.leastsquares.LeastSquaresOptimizer.Optimum; import org.apache.commons.math4.fitting.leastsquares.LeastSquaresProblem.Evaluation; import org.apache.commons.math4.geometry.euclidean.twod.Cartesian2D; @@ -39,9 +32,6 @@ import org.apache.commons.numbers.core.Precision; import org.junit.Assert; import org.junit.Test; -import java.util.ArrayList; -import java.util.List; - import static org.hamcrest.CoreMatchers.is; /** @@ -357,58 +347,4 @@ public class LevenbergMarquardtOptimizerTest Assert.assertThat(optimum.getEvaluations(), is(2)); } - private static class BevingtonProblem { - private List time; - private List count; - - public BevingtonProblem() { - time = new ArrayList<>(); - count = new ArrayList<>(); - } - - public void addPoint(double t, double c) { - time.add(t); - count.add(c); - } - - public MultivariateVectorFunction getModelFunction() { - return new MultivariateVectorFunction() { - @Override - public double[] value(double[] params) { - double[] values = new double[time.size()]; - for (int i = 0; i < values.length; ++i) { - final double t = time.get(i); - values[i] = params[0] + - params[1] * FastMath.exp(-t / params[3]) + - params[2] * FastMath.exp(-t / params[4]); - } - return values; - } - }; - } - - public MultivariateMatrixFunction getModelFunctionJacobian() { - return new MultivariateMatrixFunction() { - @Override - public double[][] value(double[] params) { - double[][] jacobian = new double[time.size()][5]; - - for (int i = 0; i < jacobian.length; ++i) { - final double t = time.get(i); - jacobian[i][0] = 1; - - final double p3 = params[3]; - final double p4 = params[4]; - final double tOp3 = t / p3; - final double tOp4 = t / p4; - jacobian[i][1] = FastMath.exp(-tOp3); - jacobian[i][2] = FastMath.exp(-tOp4); - jacobian[i][3] = params[1] * FastMath.exp(-tOp3) * tOp3 / p3; - jacobian[i][4] = params[2] * FastMath.exp(-tOp4) * tOp4 / p4; - } - return jacobian; - } - }; - } - } }