Нахождение локальных максимумов и минимумов пользовательских функций - PullRequest
Новая
       6

Нахождение локальных максимумов и минимумов пользовательских функций

0 голосов
/ 29 апреля 2018

Что я хочу

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

Моя функция выглядит так: 

import numpy as np

def func(x,y):
  return (np.cos(x*10))**2 + (np.sin(y*10))**2

Методы

Вот методы, которые я собираюсь использовать:

  1. На самом деле я уже делаю что-то подобное в Mathematica. Я дифференцирую функцию один раз, а затем дважды. Я смотрю на точки, где первая производная равна 0, вычисляю их значения и их местоположение. Затем я беру вторую производную в этих местах и ​​проверяю, являются ли они минимумами или максимумами.

  2. Мне также интересно, можно ли просто сделать двумерный массив значений функции по x и y и найти максимальное и минимальное значения этого массива. Но для этого нужно знать, как точно определить сетку x и y, чтобы надежно фиксировать поведение функции

Для последнего случая я уже нашел несколько способов, таких как этот .

Я просто хотел знать, какой метод имеет больше смысла с точки зрения эффективности, скорости, точности или даже элегантности в Python?

1 Ответ

0 голосов
/ 30 апреля 2018

найти список стационарных точек, их значений и местоположений, а также того, являются ли они минимумами или максимумами.

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

Символическое решение с SymPy

Для простых функций, таких как ваш пример, SymPy будет работать нормально. Вот полный пример нахождения стационарных точек и их классификации по собственным значениям гессиана. 

import sympy as sym
x, y = sym.symbols("x y")
f = sym.cos(x*10)**2 + sym.sin(y*10)**2
gradient = sym.derive_by_array(f, (x, y))
hessian = sym.Matrix(2, 2, sym.derive_by_array(gradient, (x, y)))

Пока что гессиан является символической матрицей 2 на 2: [[200*sin(10*x)**2 - 200*cos(10*x)**2, 0], [0, -200*sin(10*y)**2 + 200*cos(10*y)**2]]. Затем мы находим стационарные точки, приравнивая gradient к нулю, и вставляем их в гессиан по одному. 

stationary_points = sym.solve(gradient, (x, y))
for p in stationary_points:
    value = f.subs({x: p[0], y: p[1]})
    hess = hessian.subs({x: p[0], y: p[1]})
    eigenvals = hess.eigenvals()
    if all(ev > 0 for ev in eigenvals):
        print("Local minimum at {} with value {}".format(p, value))
    elif all(ev < 0 for ev in eigenvals):
        print("Local maximum at {} with value {}".format(p, value))
    elif any(ev > 0 for ev in eigenvals) and any(ev < 0 for ev in eigenvals):
        print("Saddle point at {} with value {}".format(p, value))
    else:
        print("Could not classify the stationary point at {} with value {}".format(p, value))

Последнее предложение необходимо, потому что, когда гессиан является только полу определенным, мы не можем сказать, какой тип стационарной точки это то, что (x**2 + y**4 и x**2 - y**4 имеют один и тот же гессиан в (0, 0) ) но другое поведение). Выход: 

Saddle point at (0, 0) with value 1
Local maximum at (0, pi/20) with value 2
Saddle point at (0, pi/10) with value 1
Local maximum at (0, 3*pi/20) with value 2
Local minimum at (pi/20, 0) with value 0
Saddle point at (pi/20, pi/20) with value 1
Local minimum at (pi/20, pi/10) with value 0
Saddle point at (pi/20, 3*pi/20) with value 1
Saddle point at (pi/10, 0) with value 1
Local maximum at (pi/10, pi/20) with value 2
Saddle point at (pi/10, pi/10) with value 1
Local maximum at (pi/10, 3*pi/20) with value 2
Local minimum at (3*pi/20, 0) with value 0
Saddle point at (3*pi/20, pi/20) with value 1
Local minimum at (3*pi/20, pi/10) with value 0
Saddle point at (3*pi/20, 3*pi/20) with value 1

Очевидно, solve не нашел всех решений (их бесконечно много). Рассмотрим solve vs solveset , но в любом случае обрабатывать бесконечно много решений сложно.

Числовая оптимизация с SciPy

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

  1. Каждый прогон найдет только один минимум.
  2. Заменив f на -f, вы также можете найти максимум.
  3. Изменение начальной точки поиска (аргумент x0 minimize) может привести к другому максимуму или минимуму. Тем не менее, вы никогда не узнаете, есть ли другие экстремумы, которые вы еще не видели.
  4. Ничего из этого не найдет седловые точки.

Смешанная стратегия

Используя lambdify можно превратить символическое выражение в функцию Python, которую можно передавать в числовые решатели SciPy. 

from scipy.optimize import fsolve
grad = sym.lambdify((x, y), gradient)
fsolve(lambda v: grad(v[0], v[1]), (1, 2))

Возвращает некоторую стационарную точку, [0.9424778 , 2.04203522] в этом примере. Какой это пункт, зависит от первоначального предположения, которое было (1, 2). Обычно (но не всегда) вы получите решение, близкое к первоначальному предположению.

Это имеет преимущество перед подходом прямой минимизации в том, что седловые точки также могут быть обнаружены. Тем не менее, будет сложно найти все решения, поскольку каждый запуск fsolve вызывает только одно.

...