1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.optimization.general; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ConvergenceException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.FunctionEvaluationException; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.UnivariateRealFunction; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.solvers.BrentSolver; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.solvers.UnivariateRealSolver; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.optimization.GoalType; 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.optimization.OptimizationException; 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.optimization.RealPointValuePair; 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Non-linear conjugate gradient optimizer. 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This class supports both the Fletcher-Reeves and the Polak-Ribière 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * update formulas for the conjugate search directions. It also supports 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * optional preconditioning. 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.0 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class NonLinearConjugateGradientOptimizer 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond extends AbstractScalarDifferentiableOptimizer { 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Update formula for the beta parameter. */ 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final ConjugateGradientFormula updateFormula; 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Preconditioner (may be null). */ 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private Preconditioner preconditioner; 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** solver to use in the line search (may be null). */ 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private UnivariateRealSolver solver; 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Initial step used to bracket the optimum in line search. */ 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double initialStep; 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor with default settings. 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The convergence check is set to a {@link 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * org.apache.commons.math.optimization.SimpleVectorialValueChecker} 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * and the maximal number of iterations is set to 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link AbstractScalarDifferentiableOptimizer#DEFAULT_MAX_ITERATIONS}. 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param updateFormula formula to use for updating the β parameter, 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ConjugateGradientFormula#POLAK_RIBIERE} 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula) { 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.updateFormula = updateFormula; 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond preconditioner = null; 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond solver = null; 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond initialStep = 1.0; 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Set the preconditioner. 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param preconditioner preconditioner to use for next optimization, 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * may be null to remove an already registered preconditioner 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setPreconditioner(final Preconditioner preconditioner) { 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.preconditioner = preconditioner; 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Set the solver to use during line search. 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param lineSearchSolver solver to use during line search, may be null 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to remove an already registered solver and fall back to the 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * default {@link BrentSolver Brent solver}. 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setLineSearchSolver(final UnivariateRealSolver lineSearchSolver) { 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.solver = lineSearchSolver; 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Set the initial step used to bracket the optimum in line search. 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The initial step is a factor with respect to the search direction, 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * which itself is roughly related to the gradient of the function 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param initialStep initial step used to bracket the optimum in line search, 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * if a non-positive value is used, the initial step is reset to its 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * default value of 1.0 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setInitialStep(final double initialStep) { 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (initialStep <= 0) { 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.initialStep = 1.0; 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.initialStep = initialStep; 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected RealPointValuePair doOptimize() 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws FunctionEvaluationException, OptimizationException, IllegalArgumentException { 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // initialization 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (preconditioner == null) { 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond preconditioner = new IdentityPreconditioner(); 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (solver == null) { 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond solver = new BrentSolver(); 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = point.length; 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[] r = computeObjectiveGradient(point); 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (goal == GoalType.MINIMIZE) { 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond r[i] = -r[i]; 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // initial search direction 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[] steepestDescent = preconditioner.precondition(point, r); 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[] searchDirection = steepestDescent.clone(); 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double delta = 0; 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond delta += r[i] * searchDirection[i]; 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealPointValuePair current = null; 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond while (true) { 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double objective = computeObjectiveValue(point); 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealPointValuePair previous = current; 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond current = new RealPointValuePair(point, objective); 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (previous != null) { 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (checker.converged(getIterations(), previous, current)) { 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // we have found an optimum 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return current; 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond incrementIterationsCounter(); 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double dTd = 0; 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (final double di : searchDirection) { 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond dTd += di * di; 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // find the optimal step in the search direction 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final UnivariateRealFunction lsf = new LineSearchFunction(searchDirection); 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double step = solver.solve(lsf, 0, findUpperBound(lsf, 0, initialStep)); 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // validate new point 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < point.length; ++i) { 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond point[i] += step * searchDirection[i]; 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond r = computeObjectiveGradient(point); 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (goal == GoalType.MINIMIZE) { 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond r[i] = -r[i]; 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute beta 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double deltaOld = delta; 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] newSteepestDescent = preconditioner.precondition(point, r); 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond delta = 0; 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond delta += r[i] * newSteepestDescent[i]; 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double beta; 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (updateFormula == ConjugateGradientFormula.FLETCHER_REEVES) { 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond beta = delta / deltaOld; 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double deltaMid = 0; 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < r.length; ++i) { 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond deltaMid += r[i] * steepestDescent[i]; 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond beta = (delta - deltaMid) / deltaOld; 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond steepestDescent = newSteepestDescent; 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute conjugate search direction 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if ((getIterations() % n == 0) || (beta < 0)) { 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // break conjugation: reset search direction 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond searchDirection = steepestDescent.clone(); 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute new conjugate search direction 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond searchDirection[i] = steepestDescent[i] + beta * searchDirection[i]; 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (ConvergenceException ce) { 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OptimizationException(ce); 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Find the upper bound b ensuring bracketing of a root between a and b 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f function whose root must be bracketed 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param a lower bound of the interval 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param h initial step to try 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return b such that f(a) and f(b) have opposite signs 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception FunctionEvaluationException if the function cannot be computed 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception OptimizationException if no bracket can be found 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double findUpperBound(final UnivariateRealFunction f, 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double a, final double h) 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws FunctionEvaluationException, OptimizationException { 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yA = f.value(a); 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double yB = yA; 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (double step = h; step < Double.MAX_VALUE; step *= FastMath.max(2, yA / yB)) { 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double b = a + step; 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yB = f.value(b); 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (yA * yB <= 0) { 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return b; 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OptimizationException(LocalizedFormats.UNABLE_TO_BRACKET_OPTIMUM_IN_LINE_SEARCH); 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Default identity preconditioner. */ 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static class IdentityPreconditioner implements Preconditioner { 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] precondition(double[] variables, double[] r) { 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return r.clone(); 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Internal class for line search. 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The function represented by this class is the dot product of 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the objective function gradient and the search direction. Its 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * value is zero when the gradient is orthogonal to the search 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * direction, i.e. when the objective function value is a local 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * extremum along the search direction. 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private class LineSearchFunction implements UnivariateRealFunction { 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Search direction. */ 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] searchDirection; 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor. 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param searchDirection search direction 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public LineSearchFunction(final double[] searchDirection) { 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.searchDirection = searchDirection; 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(double x) throws FunctionEvaluationException { 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // current point in the search direction 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] shiftedPoint = point.clone(); 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < shiftedPoint.length; ++i) { 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond shiftedPoint[i] += x * searchDirection[i]; 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // gradient of the objective function 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] gradient; 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond gradient = computeObjectiveGradient(shiftedPoint); 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // dot product with the search direction 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double dotProduct = 0; 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < gradient.length; ++i) { 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond dotProduct += gradient[i] * searchDirection[i]; 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return dotProduct; 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 295