001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.optim.nonlinear.scalar.noderiv;
018
019import java.util.Arrays;
020import org.apache.commons.math4.exception.MathUnsupportedOperationException;
021import org.apache.commons.math4.exception.NotStrictlyPositiveException;
022import org.apache.commons.math4.exception.NumberIsTooSmallException;
023import org.apache.commons.math4.exception.util.LocalizedFormats;
024import org.apache.commons.math4.optim.ConvergenceChecker;
025import org.apache.commons.math4.optim.PointValuePair;
026import org.apache.commons.math4.optim.nonlinear.scalar.GoalType;
027import org.apache.commons.math4.optim.nonlinear.scalar.LineSearch;
028import org.apache.commons.math4.optim.nonlinear.scalar.MultivariateOptimizer;
029import org.apache.commons.math4.optim.univariate.UnivariatePointValuePair;
030import org.apache.commons.math4.util.FastMath;
031
032/**
033 * Powell's algorithm.
034 * This code is translated and adapted from the Python version of this
035 * algorithm (as implemented in module {@code optimize.py} v0.5 of
036 * <em>SciPy</em>).
037 * <br>
038 * The default stopping criterion is based on the differences of the
039 * function value between two successive iterations. It is however possible
040 * to define a custom convergence checker that might terminate the algorithm
041 * earlier.
042 * <br>
043 * Line search is performed by the {@link LineSearch} class.
044 * <br>
045 * Constraints are not supported: the call to
046 * {@link #optimize(org.apache.commons.math4.optim.OptimizationData...)} will throw
047 * {@link MathUnsupportedOperationException} if bounds are passed to it.
048 * In order to impose simple constraints, the objective function must be
049 * wrapped in an adapter like
050 * {@link org.apache.commons.math4.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
051 * MultivariateFunctionMappingAdapter} or
052 * {@link org.apache.commons.math4.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
053 * MultivariateFunctionPenaltyAdapter}.
054 *
055 * @since 2.2
056 */
057public class PowellOptimizer
058    extends MultivariateOptimizer {
059    /**
060     * Minimum relative tolerance.
061     */
062    private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
063    /**
064     * Relative threshold.
065     */
066    private final double relativeThreshold;
067    /**
068     * Absolute threshold.
069     */
070    private final double absoluteThreshold;
071    /**
072     * Line search.
073     */
074    private final LineSearch line;
075
076    /**
077     * This constructor allows to specify a user-defined convergence checker,
078     * in addition to the parameters that control the default convergence
079     * checking procedure.
080     * <br>
081     * The internal line search tolerances are set to the square-root of their
082     * corresponding value in the multivariate optimizer.
083     *
084     * @param rel Relative threshold.
085     * @param abs Absolute threshold.
086     * @param checker Convergence checker.
087     * @throws NotStrictlyPositiveException if {@code abs <= 0}.
088     * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
089     */
090    public PowellOptimizer(double rel,
091                           double abs,
092                           ConvergenceChecker<PointValuePair> checker) {
093        this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
094    }
095
096    /**
097     * This constructor allows to specify a user-defined convergence checker,
098     * in addition to the parameters that control the default convergence
099     * checking procedure and the line search tolerances.
100     *
101     * @param rel Relative threshold for this optimizer.
102     * @param abs Absolute threshold for this optimizer.
103     * @param lineRel Relative threshold for the internal line search optimizer.
104     * @param lineAbs Absolute threshold for the internal line search optimizer.
105     * @param checker Convergence checker.
106     * @throws NotStrictlyPositiveException if {@code abs <= 0}.
107     * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
108     */
109    public PowellOptimizer(double rel,
110                           double abs,
111                           double lineRel,
112                           double lineAbs,
113                           ConvergenceChecker<PointValuePair> checker) {
114        super(checker);
115
116        if (rel < MIN_RELATIVE_TOLERANCE) {
117            throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
118        }
119        if (abs <= 0) {
120            throw new NotStrictlyPositiveException(abs);
121        }
122        relativeThreshold = rel;
123        absoluteThreshold = abs;
124
125        // Create the line search optimizer.
126        line = new LineSearch(this,
127                              lineRel,
128                              lineAbs,
129                              1d);
130    }
131
132    /**
133     * The parameters control the default convergence checking procedure.
134     * <br>
135     * The internal line search tolerances are set to the square-root of their
136     * corresponding value in the multivariate optimizer.
137     *
138     * @param rel Relative threshold.
139     * @param abs Absolute threshold.
140     * @throws NotStrictlyPositiveException if {@code abs <= 0}.
141     * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
142     */
143    public PowellOptimizer(double rel,
144                           double abs) {
145        this(rel, abs, null);
146    }
147
148    /**
149     * Builds an instance with the default convergence checking procedure.
150     *
151     * @param rel Relative threshold.
152     * @param abs Absolute threshold.
153     * @param lineRel Relative threshold for the internal line search optimizer.
154     * @param lineAbs Absolute threshold for the internal line search optimizer.
155     * @throws NotStrictlyPositiveException if {@code abs <= 0}.
156     * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
157     */
158    public PowellOptimizer(double rel,
159                           double abs,
160                           double lineRel,
161                           double lineAbs) {
162        this(rel, abs, lineRel, lineAbs, null);
163    }
164
165    /** {@inheritDoc} */
166    @Override
167    protected PointValuePair doOptimize() {
168        checkParameters();
169
170        final GoalType goal = getGoalType();
171        final double[] guess = getStartPoint();
172        final int n = guess.length;
173
174        final double[][] direc = new double[n][n];
175        for (int i = 0; i < n; i++) {
176            direc[i][i] = 1;
177        }
178
179        final ConvergenceChecker<PointValuePair> checker
180            = getConvergenceChecker();
181
182        double[] x = guess;
183        double fVal = computeObjectiveValue(x);
184        double[] x1 = x.clone();
185        while (true) {
186            incrementIterationCount();
187
188            double fX = fVal;
189            double fX2 = 0;
190            double delta = 0;
191            int bigInd = 0;
192            double alphaMin = 0;
193
194            for (int i = 0; i < n; i++) {
195                final double[] d = Arrays.copyOf(direc[i], direc[i].length);
196
197                fX2 = fVal;
198
199                final UnivariatePointValuePair optimum = line.search(x, d);
200                fVal = optimum.getValue();
201                alphaMin = optimum.getPoint();
202                final double[][] result = newPointAndDirection(x, d, alphaMin);
203                x = result[0];
204
205                if ((fX2 - fVal) > delta) {
206                    delta = fX2 - fVal;
207                    bigInd = i;
208                }
209            }
210
211            // Default convergence check.
212            boolean stop = 2 * (fX - fVal) <=
213                (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
214                 absoluteThreshold);
215
216            final PointValuePair previous = new PointValuePair(x1, fX);
217            final PointValuePair current = new PointValuePair(x, fVal);
218            if (!stop && checker != null) { // User-defined stopping criteria.
219                stop = checker.converged(getIterations(), previous, current);
220            }
221            if (stop) {
222                if (goal == GoalType.MINIMIZE) {
223                    return (fVal < fX) ? current : previous;
224                } else {
225                    return (fVal > fX) ? current : previous;
226                }
227            }
228
229            final double[] d = new double[n];
230            final double[] x2 = new double[n];
231            for (int i = 0; i < n; i++) {
232                d[i] = x[i] - x1[i];
233                x2[i] = 2 * x[i] - x1[i];
234            }
235
236            x1 = x.clone();
237            fX2 = computeObjectiveValue(x2);
238
239            if (fX > fX2) {
240                double t = 2 * (fX + fX2 - 2 * fVal);
241                double temp = fX - fVal - delta;
242                t *= temp * temp;
243                temp = fX - fX2;
244                t -= delta * temp * temp;
245
246                if (t < 0.0) {
247                    final UnivariatePointValuePair optimum = line.search(x, d);
248                    fVal = optimum.getValue();
249                    alphaMin = optimum.getPoint();
250                    final double[][] result = newPointAndDirection(x, d, alphaMin);
251                    x = result[0];
252
253                    final int lastInd = n - 1;
254                    direc[bigInd] = direc[lastInd];
255                    direc[lastInd] = result[1];
256                }
257            }
258        }
259    }
260
261    /**
262     * Compute a new point (in the original space) and a new direction
263     * vector, resulting from the line search.
264     *
265     * @param p Point used in the line search.
266     * @param d Direction used in the line search.
267     * @param optimum Optimum found by the line search.
268     * @return a 2-element array containing the new point (at index 0) and
269     * the new direction (at index 1).
270     */
271    private double[][] newPointAndDirection(double[] p,
272                                            double[] d,
273                                            double optimum) {
274        final int n = p.length;
275        final double[] nP = new double[n];
276        final double[] nD = new double[n];
277        for (int i = 0; i < n; i++) {
278            nD[i] = d[i] * optimum;
279            nP[i] = p[i] + nD[i];
280        }
281
282        final double[][] result = new double[2][];
283        result[0] = nP;
284        result[1] = nD;
285
286        return result;
287    }
288
289    /**
290     * @throws MathUnsupportedOperationException if bounds were passed to the
291     * {@link #optimize(org.apache.commons.math4.optim.OptimizationData[]) optimize} method.
292     */
293    private void checkParameters() {
294        if (getLowerBound() != null ||
295            getUpperBound() != null) {
296            throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
297        }
298    }
299}