Java tutorial
/** * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies * * Please see distribution for license. */ package com.opengamma.analytics.math.minimization; import org.apache.commons.lang.NotImplementedException; import com.opengamma.analytics.math.MathException; import com.opengamma.analytics.math.function.Function1D; import com.opengamma.analytics.math.matrix.DoubleMatrix1D; import com.opengamma.analytics.math.matrix.DoubleMatrix2D; import com.opengamma.analytics.math.matrix.DoubleMatrixUtils; import com.opengamma.analytics.math.matrix.MatrixAlgebra; import com.opengamma.analytics.math.matrix.MatrixAlgebraFactory; import com.opengamma.util.ArgumentChecker; /** * */ public class QuasiNewtonVectorMinimizer implements MinimizerWithGradient<Function1D<DoubleMatrix1D, Double>, Function1D<DoubleMatrix1D, DoubleMatrix1D>, DoubleMatrix1D> { private static final int RESET_FREQ = 200; private static final double ALPHA = 1e-4; private static final double BETA = 1.5; private static final double EPS = 1e-8; private static final int DEF_MAX_STEPS = 200; private static final MatrixAlgebra MA = MatrixAlgebraFactory.getMatrixAlgebra("OG"); private static final QuasiNewtonInverseHessianUpdate DEF_UPDATER = new BroydenFletcherGoldfarbShannoInverseHessianUpdate(); private final double _absoluteTol, _relativeTol; private final int _maxSteps; private final QuasiNewtonInverseHessianUpdate _hessainUpdater; public QuasiNewtonVectorMinimizer() { this(EPS, EPS, DEF_MAX_STEPS); } public QuasiNewtonVectorMinimizer(final double absTolerance, final double relTolerance, final int maxInterations) { this(absTolerance, relTolerance, maxInterations, DEF_UPDATER); } public QuasiNewtonVectorMinimizer(final double absoluteTol, final double relativeTol, final int maxInterations, final QuasiNewtonInverseHessianUpdate hessianUpdater) { ArgumentChecker.notNull(hessianUpdater, "null updater"); ArgumentChecker.notNegative(absoluteTol, "absolute tolerance"); ArgumentChecker.notNegative(relativeTol, "relative tolerance"); ArgumentChecker.notNegative(maxInterations, "maxSteps"); _absoluteTol = absoluteTol; _relativeTol = relativeTol; _maxSteps = maxInterations; _hessainUpdater = hessianUpdater; } /** * Disabled because not working properly (see JIRA issue) * @param function The function * @param startPosition The start position * @return The minimum */ @Override public DoubleMatrix1D minimize(final Function1D<DoubleMatrix1D, Double> function, final DoubleMatrix1D startPosition) { throw new NotImplementedException("Please supply gradient function or use ConjugateGradient"); } @Override public DoubleMatrix1D minimize(final Function1D<DoubleMatrix1D, Double> function, final Function1D<DoubleMatrix1D, DoubleMatrix1D> grad, final DoubleMatrix1D startPosition) { final DataBundle data = new DataBundle(); final double y = function.evaluate(startPosition); data.setX(startPosition); data.setG0(y * y); data.setGrad(grad.evaluate(startPosition)); data.setInverseHessianEsimate(getInitializedMatrix(startPosition)); if (!getNextPosition(function, grad, data)) { throw new MathException("Cannot work with this starting position. Please choose another point"); } int count = 0; int resetCount = 1; while (!isConverged(data)) { if ((resetCount) % RESET_FREQ == 0) { data.setInverseHessianEsimate(getInitializedMatrix(startPosition)); resetCount = 1; } else { _hessainUpdater.update(data); } if (!getNextPosition(function, grad, data)) { data.setInverseHessianEsimate(getInitializedMatrix(startPosition)); resetCount = 1; if (!getNextPosition(function, grad, data)) { throw new MathException("Failed to converge in backtracking"); } } count++; resetCount++; if (count > _maxSteps) { throw new MathException("Failed to converge after " + _maxSteps + " iterations. Final point reached: " + data.getX().toString()); } } return data.getX(); } private DoubleMatrix2D getInitializedMatrix(final DoubleMatrix1D startPosition) { return DoubleMatrixUtils.getIdentityMatrix2D(startPosition.getNumberOfElements()); } private DoubleMatrix1D getDirection(final DataBundle data) { return (DoubleMatrix1D) MA.multiply(data.getInverseHessianEsimate(), MA.scale(data.getGrad(), -1.0)); } private boolean getNextPosition(final Function1D<DoubleMatrix1D, Double> function, final Function1D<DoubleMatrix1D, DoubleMatrix1D> grad, final DataBundle data) { final DoubleMatrix1D p = getDirection(data); if (data.getLambda0() < 1.0) { data.setLambda0(1.0); } else { data.setLambda0(data.getLambda0() * BETA); } updatePosition(p, function, data); final double g1 = data.getG1(); // the function is invalid at the new position, try to recover if (Double.isInfinite(g1) || Double.isNaN(g1)) { bisectBacktrack(p, function, data); } if (data.getG1() > data.getG0() / (1 + ALPHA * data.getLambda0())) { quadraticBacktrack(p, function, data); int count = 0; while (data.getG1() > data.getG0() / (1 + ALPHA * data.getLambda0())) { if (count > 5) { return false; } cubicBacktrack(p, function, data); count++; } } final DoubleMatrix1D deltaX = data.getDeltaX(); data.setX((DoubleMatrix1D) MA.add(data.getX(), deltaX)); data.setG0(data.getG1()); final DoubleMatrix1D gradNew = grad.evaluate(data.getX()); data.setDeltaGrad((DoubleMatrix1D) MA.subtract(gradNew, data.getGrad())); data.setGrad(gradNew); return true; } protected void updatePosition(final DoubleMatrix1D p, final Function1D<DoubleMatrix1D, Double> function, final DataBundle data) { final double lambda0 = data.getLambda0(); final DoubleMatrix1D deltaX = (DoubleMatrix1D) MA.scale(p, lambda0); final DoubleMatrix1D xNew = (DoubleMatrix1D) MA.add(data.getX(), deltaX); data.setDeltaX(deltaX); data.setG2(data.getG1()); final double y = function.evaluate(xNew); data.setG1(y * y); } private void bisectBacktrack(final DoubleMatrix1D p, final Function1D<DoubleMatrix1D, Double> function, final DataBundle data) { do { data.setLambda0(data.getLambda0() * 0.1); updatePosition(p, function, data); } while (Double.isNaN(data.getG1()) || Double.isInfinite(data.getG1()) || Double.isNaN(data.getG2()) || Double.isInfinite(data.getG2())); } private void quadraticBacktrack(final DoubleMatrix1D p, final Function1D<DoubleMatrix1D, Double> function, final DataBundle data) { final double lambda0 = data.getLambda0(); final double g0 = data.getG0(); final double lambda = Math.max(0.01 * lambda0, g0 * lambda0 * lambda0 / (data.getG1() + g0 * (2 * lambda0 - 1))); data.swapLambdaAndReplace(lambda); updatePosition(p, function, data); } private void cubicBacktrack(final DoubleMatrix1D p, final Function1D<DoubleMatrix1D, Double> function, final DataBundle data) { double temp1, temp2, temp3, temp4, temp5; final double lambda0 = data.getLambda0(); final double lambda1 = data.getLambda1(); final double g0 = data.getG0(); temp1 = 1.0 / lambda0 / lambda0; temp2 = 1.0 / lambda1 / lambda1; temp3 = data.getG1() + g0 * (2 * lambda0 - 1.0); temp4 = data.getG2() + g0 * (2 * lambda1 - 1.0); temp5 = 1.0 / (lambda0 - lambda1); final double a = temp5 * (temp1 * temp3 - temp2 * temp4); final double b = temp5 * (-lambda1 * temp1 * temp3 + lambda0 * temp2 * temp4); double lambda = (-b + Math.sqrt(b * b + 6 * a * g0)) / 3 / a; lambda = Math.min(Math.max(lambda, 0.01 * lambda0), 0.75 * lambda1); // make sure new lambda is between 1% & 75% of old value data.swapLambdaAndReplace(lambda); updatePosition(p, function, data); } private boolean isConverged(final DataBundle data) { final DoubleMatrix1D deltaX = data.getDeltaX(); final DoubleMatrix1D x = data.getX(); final int n = deltaX.getNumberOfElements(); double diff, scale; for (int i = 0; i < n; i++) { diff = Math.abs(deltaX.getEntry(i)); scale = Math.abs(x.getEntry(i)); if (diff > _absoluteTol + scale * _relativeTol) { return false; } } return (MA.getNorm2(data.getGrad()) < _absoluteTol); } /** * Data bundle for intermediate data */ public static class DataBundle { private double _g0; private double _g1; private double _g2; private double _lambda0; private double _lambda1; private DoubleMatrix1D _deltaGrad; private DoubleMatrix1D _grad; private DoubleMatrix1D _deltaX; private DoubleMatrix1D _x; private DoubleMatrix2D _h; public double getG0() { return _g0; } public double getG1() { return _g1; } public double getG2() { return _g2; } public double getLambda0() { return _lambda0; } public double getLambda1() { return _lambda1; } public DoubleMatrix1D getDeltaGrad() { return _deltaGrad; } public DoubleMatrix1D getGrad() { return _grad; } public DoubleMatrix1D getDeltaX() { return _deltaX; } public DoubleMatrix1D getX() { return _x; } public void setG0(final double g0) { _g0 = g0; } public void setG1(final double g1) { _g1 = g1; } public void setG2(final double g2) { _g2 = g2; } public void setLambda0(final double lambda0) { _lambda0 = lambda0; } public void setDeltaGrad(final DoubleMatrix1D deltaGrad) { _deltaGrad = deltaGrad; } public void setGrad(final DoubleMatrix1D grad) { _grad = grad; } public void setDeltaX(final DoubleMatrix1D deltaX) { _deltaX = deltaX; } public void setX(final DoubleMatrix1D x) { _x = x; } /** * Inverse Hessian matrix * @return The inverse Hessian Matrix */ public DoubleMatrix2D getInverseHessianEsimate() { return _h; } public void setInverseHessianEsimate(final DoubleMatrix2D estimate) { _h = estimate; } public void swapLambdaAndReplace(final double lambda0) { _lambda1 = _lambda0; _lambda0 = lambda0; } } }