Как найти минимум функции python

Using, for instance, scipy‘s fmin (which contains an implementation of the Nelder-Mead algorithm), you can try this:

import numpy as np
from scipy.optimize import fmin
import math

def f(x):
    exp = (math.pow(x[0], 2) + math.pow(x[1], 2)) * -1
    return math.exp(exp) * math.cos(x[0] * x[1]) * math.sin(x[0] * x[1])

fmin(f,np.array([0,0]))

which yields the following output:

Optimization terminated successfully.
     Current function value: -0.161198
     Iterations: 60
     Function evaluations: 113 
array([ 0.62665701, -0.62663095])

Please keep in mind that:

1) with scipy you need to convert your function into a function accepting an array (I showed how to do it in the example above);

2) fmin uses, like most of its pairs, an iterative algorithm, therefore you must provide a starting point (in my example, I provided (0,0)). You can provide different starting points to obtain different minima/maxima.

Время на прочтение
5 мин

Количество просмотров 56K

Метод Нелдера — Мида — метод оптимизации (поиска минимума) функции от нескольких переменных. Простой и в тоже время эффективный метод, позволяющий оптимизировать функции без использования градиентов. Метод надежен и, как правило, показывает хорошие результаты, хотя и отсутствует теория сходимости. Может использоваться в функции optimize из модуля scipy.optimize популярной библиотеки для языка python, которая используется для математических расчетов.

Алгоритм заключается в формировании симплекса (simplex) и последующего его деформирования в направлении минимума, посредством трех операций:

1) Отражение (reflection);
2) Растяжение (expansion);
3) Сжатие (contract);

Симплекс представляет из себя геометрическую фигуру, являющуюся n — мерным обобщением треугольника. Для одномерного пространства — это отрезок, для двумерного — треугольник. Таким образом n — мерный симплекс имеет n + 1 вершину.

Алгоритм

1) Пусть

функция, которую необходимо оптимизировать. На первом шаге выбираем три случайные точки (об этом чуть позже) и формируем симплекс (треугольник). Вычисляем значение функции в каждой точке:

,

,

.

Сортируем точки по значениям функции

в этих точках, таким образом получаем двойное неравенство:

Мы ищем минимум функции, а следовательно, на данном шаге лучшей будет та точка, в которой значение функции минимально. Для удобства переобозначим точки следующим образом:

b =

, g =

, w =

, где best, good, worst — соответственно.

2) На следующем шаге находим середину отрезка, точками которого являются g и b. Т.к. координаты середины отрезка равны полусумме координат его концов, получаем:

В более общем виде можно записать так:

3) Применяем операцию отражения:
Находим точку

, следующим образом:

Т.е. фактически отражаем точку w относительно mid. В качестве коэффициента берут как правило 1. Проверяем нашу точку: если

, то это хорошая точка. А теперь попробуем расстояние увеличить в 2 раза, вдруг нам повезет и мы найдем точку еще лучше.

4) Применяем операцию растяжения:
Находим точку

следующим образом:

В качестве γ принимаем γ = 2, т.е. расстояние увеличиваем в 2 раза.

Проверяем точку

:

Если

, то нам повезло и мы нашли точку лучше, чем есть на данный момент, если бы этого не произошло, мы бы остановились на точке

.

Далее заменяем точку w на

, в итоге получаем:

5) Если же нам совсем не повезло и мы не нашли хороших точек, пробуем операцию сжатия.
Как следует из названия операции мы будем уменьшать наш отрезок и искать хорошие точки внутри треугольника.

Пробуем найти хорошую точку

:

Коэффициент β принимаем равным 0.5, т.е. точка

на середине отрезка wmid.

Существует еще одна операция — shrink (сокращение). В данном случае, мы переопределяем весь симплекс. Оставляем только «лучшую» точку, остальные определяем следующим образом:

Коэффициент δ берут равным 0.5.

По существу передвигаем точки по направлению к текущей «лучшей» точке. Преобразование выглядит следующим образом:

Необходимо отметить, что данная операция дорого обходится, поскольку необходимо заменять точки в симплексе. К счастью было установлено, при проведении большого количества экспериментов, что shrink — трансформация редко случается на практике.

Алгоритм заканчивается, когда:

1) Было выполнено необходимое количество итераций.
2) Площадь симплекса достигла определенной величины.
3) Текущее лучшее решение достигло необходимой точности.

Как и в большинстве эвристических методов, не существует идеального способа выбора инициализирующих точек. Как уже было сказано, можно брать случайные точки, находящиеся недалеко друг от друга для формирования симплекса; но есть решение и получше, которое используется в реализации алгоритма в MATHLAB:

Выбор первой точки

поручаем пользователю, если он имеет некоторое представление о возможном хорошем решении, в противном случае выбирается случайным образом. Остальные точки выбираются исходя из

, на небольшом расстоянии вдоль направления каждого измерения:

где

— единичный вектор.

определяется таким образом:

= 0.05, если коэффициент при

в определении

не нулевой.

= 0.00025, если коэффициент при

в определении нулевой.

Пример:

Найти экстремум следующей функции:

В качестве начальных возьмем точки:

Вычислим значение функции в каждой точке:

Переобозначим точки следующим образом:

Находим середину отрезка bg:

Находим точку

(операция отражения):

если α=1, тогда:

Проверяем точку

:

, т.к.

пробуем увеличить отрезок (операция растяжения).

если γ = 2, тогда:

Проверяем значение функции в точке

:

Оказалось, что точка

«лучше» точки b. Следовательно мы получаем новые вершины:

И алгоритм начинается сначала.

Таблица значений для 10 итераций:

Best	Good	Worst

Аналитически находим экстремум функции, он достигается в точке

.
После 10 итераций мы получаем достаточно точное приближение:

Еще о методе:

Алгоритм Нелдера — Мида в основном используется для выбора параметра в машинном обучении. В сущности, симплекс-метод используется для оптимизации параметров модели. Это связано с тем, что данный метод оптимизирует целевую функцию довольно быстро и эффективно (особенно там, где не используется shrink — модификация).

С другой стороны, в силу отсутствия теории сходимости, на практике метод может приводить к неверному ответу даже на гладких (непрерывно дифференцируемых) функциях. Также возможна ситуация, когда рабочий симплекс находится далеко от оптимальной точки, а алгоритм производит большое число итераций, при этом мало изменяя значения функции. Эвристический метод решения этой проблемы заключается в запуске алгоритма несколько раз и ограничении числа итераций.

Реализация на языке программирования python:

Создаем вспомогательный класс Vector и перегружаем операторы для возможности производить с векторами базовые операции. Я намерено не использовал вспомогательные библиотеки для реализации алгоритма, т.к. в таком случае зачастую снижается восприятие.

#!/usr/bin/python
# -*- coding: utf-8 -*-
class Vector(object):
    def __init__(self, x, y):
        """ Create a vector, example: v = Vector(1,2) """
        self.x = x
        self.y = y

    def __repr__(self):
        return "({0}, {1})".format(self.x, self.y)

    def __add__(self, other):
        x = self.x + other.x
        y = self.y + other.y
        return Vector(x, y)

    def __sub__(self, other):
        x = self.x - other.x
        y = self.y - other.y
        return Vector(x, y)

    def __rmul__(self, other):
        x = self.x * other
        y = self.y * other
        return Vector(x, y)

    def __truediv__(self, other):
        x = self.x / other
        y = self.y / other
        return Vector(x, y)

    def c(self):
        return (self.x, self.y)
        
# objective function
def f(point):
    x, y = point
    return x**2 + x*y + y**2 - 6*x - 9*y

def nelder_mead(alpha=1, beta=0.5, gamma=2, maxiter=10):
    
    # initialization
    v1 = Vector(0, 0)
    v2 = Vector(1.0, 0)
    v3 = Vector(0, 1)

    for i in range(maxiter):
        adict = {v1:f(v1.c()), v2:f(v2.c()), v3:f(v3.c())}
        points = sorted(adict.items(), key=lambda x: x[1])
        
        b = points[0][0]
        g = points[1][0]
        w = points[2][0]
        
        
        mid = (g + b)/2

        # reflection
        xr = mid + alpha * (mid - w)
        if f(xr.c()) < f(g.c()):
            w = xr
        else:
            if f(xr.c()) < f(w.c()):
                w = xr
            c = (w + mid)/2
            if f(c.c()) < f(w.c()):
                w = c
        if f(xr.c()) < f(b.c()):

            # expansion
            xe = mid + gamma * (xr - mid)
            if f(xe.c()) < f(xr.c()):
                w = xe
            else:
                w = xr
        if f(xr.c()) > f(g.c()):
            
            # contraction
            xc = mid + beta * (w - mid)
            if f(xc.c()) < f(w.c()):
                w = xc

        # update points
        v1 = w
        v2 = g
        v3 = b
    return b

print("Result of Nelder-Mead algorithm: ")
xk = nelder_mead()
print("Best poits is: %s"%(xk))

Спасибо за чтение статьи. Надеюсь она была Вам полезна и Вы узнали много нового.
С вами был FUNNYDMAN. Удачной оптимизации!)

2.7. Mathematical optimization: finding minima of functions¶

Authors: Gaël Varoquaux

Mathematical optimization deals with the
problem of finding numerically minimums (or maximums or zeros) of
a function. In this context, the function is called cost function, or
objective function, or energy.

Here, we are interested in using scipy.optimize for black-box
optimization: we do not rely on the mathematical expression of the
function that we are optimizing. Note that this expression can often be
used for more efficient, non black-box, optimization.

2.7.1. Knowing your problem¶

Not all optimization problems are equal. Knowing your problem enables you
to choose the right tool.

Dimensionality of the problem

The scale of an optimization problem is pretty much set by the
dimensionality of the problem, i.e. the number of scalar variables
on which the search is performed.

2.7.1.1. Convex versus non-convex optimization¶

A convex function: f is above all its tangents. equivalently, for two point A, B, f(C) lies below the segment [f(A), f(B])], if A < C < B	A non-convex function

Optimizing convex functions is easy. Optimizing non-convex functions can
be very hard.

Note

It can be proven that for a convex function a local minimum is
also a global minimum. Then, in some sense, the minimum is unique.

2.7.1.2. Smooth and non-smooth problems¶

A smooth function: The gradient is defined everywhere, and is a continuous function	A non-smooth function

Optimizing smooth functions is easier
(true in the context of black-box optimization, otherwise
Linear Programming
is an example of methods which deal very efficiently with
piece-wise linear functions).

2.7.1.3. Noisy versus exact cost functions¶

Noisy (blue) and non-noisy (green) functions

Noisy gradients

Many optimization methods rely on gradients of the objective function.
If the gradient function is not given, they are computed numerically,
which induces errors. In such situation, even if the objective
function is not noisy, a gradient-based optimization may be a noisy
optimization.

2.7.1.4. Constraints¶

Optimizations under constraints Here:

2.7.2. A review of the different optimizers¶

2.7.2.1. Getting started: 1D optimization¶

Let’s get started by finding the minimum of the scalar function
. scipy.optimize.minimize_scalar() uses
Brent’s method to find the minimum of a function:

>>> from scipy import optimize
>>> def f(x):
...     return -np.exp(-(x - 0.7)**2)
>>> result = optimize.minimize_scalar(f)
>>> result.success # check if solver was successful
True
>>> x_min = result.x
>>> x_min 
0.699999999...
>>> x_min - 0.7 
-2.16...e-10

Brent’s method on a quadratic function: it converges in 3 iterations, as the quadratic approximation is then exact.
Brent’s method on a non-convex function: note that the fact that the optimizer avoided the local minimum is a matter of luck.

Note

You can use different solvers using the parameter method.

2.7.2.2. Gradient based methods¶

Some intuitions about gradient descent¶

Here we focus on intuitions, not code. Code will follow.

Gradient descent
basically consists in taking small steps in the direction of the
gradient, that is the direction of the steepest descent.

Fixed step gradient descent¶

A well-conditioned quadratic function.
An ill-conditioned quadratic function. The core problem of gradient-methods on ill-conditioned problems is that the gradient tends not to point in the direction of the minimum.

We can see that very anisotropic (ill-conditioned) functions are harder
to optimize.

Take home message: conditioning number and preconditioning

If you know natural scaling for your variables, prescale them so that
they behave similarly. This is related to preconditioning.

Also, it clearly can be advantageous to take bigger steps. This
is done in gradient descent code using a
line search.

Adaptive step gradient descent¶

A well-conditioned quadratic function.
An ill-conditioned quadratic function.
An ill-conditioned non-quadratic function.
An ill-conditioned very non-quadratic function.

The more a function looks like a quadratic function (elliptic
iso-curves), the easier it is to optimize.

Conjugate gradient descent¶

The gradient descent algorithms above are toys not to be used on real
problems.

As can be seen from the above experiments, one of the problems of the
simple gradient descent algorithms, is that it tends to oscillate across
a valley, each time following the direction of the gradient, that makes
it cross the valley. The conjugate gradient solves this problem by adding
a friction term: each step depends on the two last values of the
gradient and sharp turns are reduced.

Conjugate gradient descent¶

An ill-conditioned non-quadratic function.
An ill-conditioned very non-quadratic function.

scipy provides scipy.optimize.minimize() to find the minimum of scalar
functions of one or more variables. The simple conjugate gradient method can
be used by setting the parameter method to CG

>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> optimize.minimize(f, [2, -1], method="CG")    
     fun: 1.6...e-11
     jac: array([-6.15...e-06,   2.53...e-07])
 message: ...'Optimization terminated successfully.'
    nfev: 108
     nit: 13
    njev: 27
  status: 0
 success: True
       x: array([0.99999...,  0.99998...])

Gradient methods need the Jacobian (gradient) of the function. They can compute it
numerically, but will perform better if you can pass them the gradient:

>>> def jacobian(x):
...     return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.minimize(f, [2, 1], method="CG", jac=jacobian)    
     fun: 2.957...e-14
     jac: array([ 7.1825...e-07,  -2.9903...e-07])
 message: 'Optimization terminated successfully.'
    nfev: 16
     nit: 8
    njev: 16
  status: 0
 success: True
       x: array([1.0000...,  1.0000...])

Note that the function has only been evaluated 27 times, compared to 108
without the gradient.

2.7.2.3. Newton and quasi-newton methods¶

Newton methods: using the Hessian (2nd differential)¶

Newton methods use a
local quadratic approximation to compute the jump direction. For this
purpose, they rely on the 2 first derivative of the function: the
gradient and the Hessian.

An ill-conditioned quadratic function: Note that, as the quadratic approximation is exact, the Newton method is blazing fast
An ill-conditioned non-quadratic function: Here we are optimizing a Gaussian, which is always below its quadratic approximation. As a result, the Newton method overshoots and leads to oscillations.
An ill-conditioned very non-quadratic function:

In scipy, you can use the Newton method by setting method to Newton-CG in
scipy.optimize.minimize(). Here, CG refers to the fact that an internal
inversion of the Hessian is performed by conjugate gradient

>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> def jacobian(x):
...     return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.minimize(f, [2,-1], method="Newton-CG", jac=jacobian)    
     fun: 1.5...e-15
     jac: array([  1.0575...e-07,  -7.4832...e-08])
 message: ...'Optimization terminated successfully.'
    nfev: 11
    nhev: 0
     nit: 10
    njev: 52
  status: 0
 success: True
       x: array([0.99999...,  0.99999...])

Note that compared to a conjugate gradient (above), Newton’s method has
required less function evaluations, but more gradient evaluations, as it
uses it to approximate the Hessian. Let’s compute the Hessian and pass it
to the algorithm:

>>> def hessian(x): # Computed with sympy
...     return np.array(((1 - 4*x[1] + 12*x[0]**2, -4*x[0]), (-4*x[0], 2)))
>>> optimize.minimize(f, [2,-1], method="Newton-CG", jac=jacobian, hess=hessian)    
     fun: 1.6277...e-15
     jac: array([  1.1104...e-07,  -7.7809...e-08])
 message: ...'Optimization terminated successfully.'
    nfev: 11
    nhev: 10
     nit: 10
    njev: 20
  status: 0
 success: True
       x: array([0.99999...,  0.99999...])

Note

At very high-dimension, the inversion of the Hessian can be costly
and unstable (large scale > 250).

Note

Newton optimizers should not to be confused with Newton’s root finding
method, based on the same principles, scipy.optimize.newton().

Quasi-Newton methods: approximating the Hessian on the fly¶

BFGS: BFGS (Broyden-Fletcher-Goldfarb-Shanno algorithm) refines at
each step an approximation of the Hessian.

2.7.3. Full code examples¶

2.7.4. Examples for the mathematical optimization chapter¶

Gallery generated by Sphinx-Gallery

An ill-conditioned quadratic function: On a exactly quadratic function, BFGS is not as fast as Newton’s method, but still very fast.
An ill-conditioned non-quadratic function: Here BFGS does better than Newton, as its empirical estimate of the curvature is better than that given by the Hessian.
An ill-conditioned very non-quadratic function:

>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> def jacobian(x):
...     return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.minimize(f, [2, -1], method="BFGS", jac=jacobian)    
      fun: 2.6306...e-16
 hess_inv: array([[0.99986...,  2.0000...],
       [2.0000...,  4.498...]])
      jac: array([  6.7089...e-08,  -3.2222...e-08])
  message: ...'Optimization terminated successfully.'
     nfev: 10
      nit: 8
     njev: 10
   status: 0
  success: True
        x: array([1.        ,  0.99999...])

L-BFGS: Limited-memory BFGS Sits between BFGS and conjugate gradient:
in very high dimensions (> 250) the Hessian matrix is too costly to
compute and invert. L-BFGS keeps a low-rank version. In addition, box bounds
are also supported by L-BFGS-B:

>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> def jacobian(x):
...     return np.array((-2*.5*(1 - x[0]) - 4*x[0]*(x[1] - x[0]**2), 2*(x[1] - x[0]**2)))
>>> optimize.minimize(f, [2, 2], method="L-BFGS-B", jac=jacobian)    
      fun: 1.4417...e-15
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
      jac: array([  1.0233...e-07,  -2.5929...e-08])
  message: ...'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 17
      nit: 16
   status: 0
  success: True
        x: array([1.0000...,  1.0000...])

2.7.4.12. Gradient-less methods¶

A shooting method: the Powell algorithm¶

Almost a gradient approach

An ill-conditioned quadratic function: Powell’s method isn’t too sensitive to local ill-conditionning in low dimensions
An ill-conditioned very non-quadratic function:

Simplex method: the Nelder-Mead¶

The Nelder-Mead algorithms is a generalization of dichotomy approaches to
high-dimensional spaces. The algorithm works by refining a simplex, the generalization of intervals
and triangles to high-dimensional spaces, to bracket the minimum.

Strong points: it is robust to noise, as it does not rely on
computing gradients. Thus it can work on functions that are not locally
smooth such as experimental data points, as long as they display a
large-scale bell-shape behavior. However it is slower than gradient-based
methods on smooth, non-noisy functions.

An ill-conditioned non-quadratic function:
An ill-conditioned very non-quadratic function:

Using the Nelder-Mead solver in scipy.optimize.minimize():

>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> optimize.minimize(f, [2, -1], method="Nelder-Mead") 
 final_simplex: (array([[1.0000...,  1.0000...],
       [0.99998... ,  0.99996... ],
       [1.0000...,  1.0000... ]]), array([1.1152...e-10,   1.5367...e-10,   4.9883...e-10]))
           fun: 1.1152...e-10
       message: ...'Optimization terminated successfully.'
          nfev: 111
           nit: 58
        status: 0
       success: True
             x: array([1.0000...,  1.0000...])

2.7.4.13. Global optimizers¶

If your problem does not admit a unique local minimum (which can be hard
to test unless the function is convex), and you do not have prior
information to initialize the optimization close to the solution, you
may need a global optimizer.

Brute force: a grid search¶

scipy.optimize.brute() evaluates the function on a given grid of
parameters and returns the parameters corresponding to the minimum
value. The parameters are specified with ranges given to
numpy.mgrid. By default, 20 steps are taken in each direction:

>>> def f(x):   # The rosenbrock function
...     return .5*(1 - x[0])**2 + (x[1] - x[0]**2)**2
>>> optimize.brute(f, ((-1, 2), (-1, 2))) 
array([1.0000...,  1.0000...])

>>> optimize.check_grad(f, jacobian, [2, -1])
2.384185791015625e-07

np.random.seed(0)
K = np.random.normal(size=(100, 100))

def f(x):
    return np.sum((np.dot(K, x - 1))**2) + np.sum(x**2)**2

>>> def f(x):
...     return np.arctan(x) - np.arctan(np.linspace(0, 1, len(x)))

>>> x0 = np.zeros(10)
>>> optimize.leastsq(f, x0)
(array([0.        ,  0.11111111,  0.22222222,  0.33333333,  0.44444444,
        0.55555556,  0.66666667,  0.77777778,  0.88888889,  1.        ]), 2)

>>> def g(x):
...     return np.sum(f(x)**2)
>>> optimize.minimize(g, x0, method="BFGS")   
  fun: 2.6940...e-11
 hess_inv: array([[...
 ...
           ...]])
      jac: array([...
 ...
           ...])
  message: ...'Optimization terminated successfully.'
     nfev: 144
      nit: 11
     njev: 12
   status: 0
  success: True
        x: array([-7.3...e-09,   1.1111...e-01,   2.2222...e-01, 3.3333...e-01,
        4.4444...e-01,   5.5555...e-01, 6.6666...e-01,   7.7777...e-01,
        8.8889...e-01, 1.0000...e+00])

>>> def f(t, omega, phi):
...     return np.cos(omega * t + phi)

>>> x = np.linspace(0, 3, 50)
>>> y = f(x, 1.5, 1) + .1*np.random.normal(size=50)

>>> optimize.curve_fit(f, x, y)   
(array([1.5185...,  0.92665...]), array([[ 0.00037..., -0.00056...],
       [-0.0005...,  0.00123...]]))

>>> def f(x):
...    return np.sqrt((x[0] - 3)**2 + (x[1] - 2)**2)
>>> optimize.minimize(f, np.array([0, 0]), bounds=((-1.5, 1.5), (-1.5, 1.5)))   
      fun: 1.5811...
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
      jac: array([-0.94868..., -0.31622...])
  message: ...'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 9
      nit: 2
   status: 0
  success: True
        x: array([1.5,  1.5])

2.7.8. Full code examples¶