Локальная минимизация по нескольким переменным#
Поднимаем ставки. Пусть теперь \(f:\mathbb{R}^m \rightarrow \mathbb{R}\) - дифференциируемая функция нескольких переменных. Также хотим найти локальный минимум.
В основе всего и вся будет лежать градиент (хотя есть и неградиентные методы, см. симплекс-метод).
Градиентный спуск#
Начнём с самого простого и очевидного - обычный градиентный спуск.
Его алгоритм уже описывался, но мы любезно напомним. В общем случае, он выглядит вот так:
Алгоритм градиентного спуска
Задают начальное приближение и точность расчёта \(\vec{x}_{0}, \varepsilon\)
Рассчитывают \(\vec{x}_{n+1}=\vec{x}_n-\alpha \cdot \nabla f\left(\vec{x}_n\right)\), где \(\alpha\) - шаг градиентного спуска или learning rate .
Проверяют условия остановки:
Если \(\left|\vec{x}_{n+1}-\vec{x}_n\right|>\varepsilon,\;\; \left|f\left(\vec{x}_{n+1}\right)-f\left(\vec{x}_n\right)\right|>\varepsilon\) или \(\left\|\nabla f\left(\vec{x}_{n+1}\right)\right\|>\varepsilon\) (выбирают одно из условий), то \(n=n+1\) и переход к шагу 2.
Иначе \(\vec{x}=\vec{x}_n\) и остановка.
Про него уже много было сказано, так что погнали сразу модифицировать.
Наискорейший спуск#
Самая очевидная модификация - сделать шаг градиентного спуска зависимым от текущей точки \(\alpha \rightarrow \alpha_n\), чтобы «не топтаться на месте». С другой стороны, слишком далеко тоже шагнуть не хочется.
Это наводит нас на мысль, что \(\alpha_n\) можно выбирать из условия минимизации функции одной переменной (sic!):
Эту самую минимизацию под одной переменной можно проводить методами из раздела выше, которые, как правило, сходятся очень быстро.
Итого, алгоритм метода наискорейшего спуска:
Алгоритм наискорейшего спуска
Задают начальное приближение и точность расчёта \(\vec{x}_{0}, \varepsilon\)
Рассчитывают \(\vec{x}_{n+1}=\vec{x}_n-\alpha_n \cdot \nabla f\left(\vec{x}_n\right)\), где \(\alpha_n\) считают из однопараметрической минимизации: $\( \alpha_n \stackrel{\text { def }}{=}\underset{\alpha}{\operatorname{argmin}} \; f\left(\vec{x}_n-\alpha \cdot \nabla f\left(\vec{x}_n\right)\right) \)$
Проверяют условия остановки:
Если \(\left|\vec{x}_{n+1}-\vec{x}_n\right|>\varepsilon,\;\; \left|f\left(\vec{x}_{n+1}\right)-f\left(\vec{x}_n\right)\right|>\varepsilon\) или \(\left\|\nabla f\left(\vec{x}_{n+1}\right)\right\|>\varepsilon\) (выбирают одно из условий), то \(n=n+1\) и переход к шагу 2.
Иначе \(\vec{x}=\vec{x}_n\) и остановка.
Работает стабильнее и быстрее обычного град. спуска. Но есть методы и покруче.
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
def f(x):
return 5*x[0]**2 + x[1]**2 # Пример функции (квадратичная)
def grad_f(x):
return np.array([10*x[0], 2*x[1]]) # Градиент функции
def rapid_descent(starting_point, epsilon):
x = starting_point
iterations = 0
steps = [x] # Список для хранения шагов алгоритма
while True:
phi = lambda a : f(x - a * grad_f(x))
alpha = minimize_scalar(phi).x
x_new = x - alpha * grad_f(x)
iterations += 1
steps.append(x_new)
if np.linalg.norm(x_new - x) < epsilon:
x = x_new
break
x = x_new
return x, iterations, np.array(steps)
# Генерация сетки для построения линий уровня
x = np.linspace(-2, 2, 400)
y = np.linspace(-2, 2, 400)
X, Y = np.meshgrid(x, y)
Z = 5*X**2 + Y**2
# Вызов функции градиентного спуска
starting_point = np.array([2, 2])
epsilon = 0.001
min_point, iterations, steps = rapid_descent(starting_point, epsilon)
# Построение графика
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, Z, levels=10, cmap='viridis')
plt.colorbar()
# Построение линий уровня
contour = plt.contour(X, Y, Z, levels=10, colors='black')
# Подписи к линиям уровня
plt.clabel(contour, colors='k', fmt='%2.1f', fontsize=12)
# Построение шагов алгоритма
for i in range(len(steps)-1):
plt.arrow(steps[i][0], steps[i][1], steps[i+1][0] - steps[i][0], steps[i+1][1] - steps[i][1],
shape='full', lw=1, color='red', length_includes_head=True, head_width=0.05, zorder=4)
# Подпись с количеством итераций
plt.text(-1.5, 1.8, f'Iterations: {iterations}', color='white', fontsize=12)
# Настройка осей и отображение
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Rapid Descent Visualization')
plt.show()
Сопряженные градиенты#
Наискорейший спуск работает хорошо, но можно ещё лучше. В наискорейшем спуске между соседними итерациями мы двигались в ортогональных направлениях - так получалось из минимизации функции вдоль текущего градиента. Но почему стоит двигаться именно ортогонально?
Если мы находимся в окрестности локального минимума, то вся функция приближается некоторой положительно определенной квадратичной формой, которую можно сопаставить с симметричной матрицей \(G\). Оказывается, можно ориентироваться на ортогональность градиентов именно в терминах этой матрицы (т.н. сопряженность):
Алгоритм, использующий этот подход называется методом сопряженных градиентов. Как правило, работает гораздо лучше наискорейшего. Собсна, алгоритм:
Алгоритм метода сопряженных градиентов
Пусть \(\vec{x}_0\) - начальная точка, \(\vec{r}_0=-\nabla f(\vec{x}_0\) - антиградиент.
Положим \(\vec{S}_0=\vec{r}_0\) и найдём минимум вдоль направления \(\vec{S}_0\). Обозначим точку минимума \(\vec{x}_1\).
Пусть на некотором шаге мы находимся в точке \(\vec{x}_k\), и \(\vec{r}_k\) - направление антиградиента. Положим \(\vec{S}_k=\vec{r}_k+\beta_k \vec{S}_{k-1}\), где \(\beta_k\) выбирают
либо \(\frac{\left(\vec{r}_k, \vec{r}_k\right)}{\left(\vec{r}_{k-1}, \vec{r}_{k-1}\right)}\) (стандартный алгоритм - Флетчера-Ривса),
либо \(\max \left(0, \frac{\left(\vec{r}_k, \vec{r}_k-\vec{r}_{k-1}\right)}{\left(\vec{r}_{k-1}, \vec{r}_{k-1}\right)}\right)\) (алгоритм Полака-Рибьера, иногда циклится).
После чего найдём минимум в направлении \(\overrightarrow{S_k}\) и обозначим точку минимума \(\vec{x}_{k+1}\).
Если в вычисленном направлении выполнилось условие остановки, то нужно забыть предыдущее направление, положив \(\beta_k=0\) и повторив шаг 2. Если оно выполнилось второй раз подряд, то остановка алгоритма.
Как правило, работает во много раз лучше метода наискорейшего спуска, но необходимо хорошее начальное приближение.
Примечание. Вообще, есть много вариантов выбора сопряженных градиентов. Вот некоторые из них:
Формула |
Авторы |
Год |
|---|---|---|
\(\beta_k^{HS} = \cfrac{\langle r_{k}, r_k - r_{k-1} \rangle}{\langle S_k, r_k - r_{k-1} \rangle}\) |
Хестенс, Штифель |
1952 г. |
\(\beta_k^{FR} = \cfrac{|r_k|^2}{|r_{k-1}|^2}\) |
Флетчер, Ривз |
1964 г. |
\(\beta_k^{PR} = \cfrac{\langle r_{k}, r_k - r_{k-1} \rangle}{|r_{k-1}|^2}\) |
Полак, Рибьер |
1969 г. |
\(\beta_k^{CD} = \cfrac{|r_k|^2}{\langle -S_k, r_{k-1} \rangle}\) |
Флетчер |
1987 г. |
\(\beta_k^{LS} = \cfrac{\langle r_k, r_k - r_{k-1} \rangle}{\langle -S_k, r_{k-1} \rangle}\) |
Лиу, Стори |
1991 г. |
\(\beta_k^{DY} = \cfrac{|r_k|^2}{\langle S_k, r_k - r_{k-1} \rangle}\) |
Дай, Юань |
1999 г. |
\(\beta_k^{HZ} = \left\langle r_k - r_{k-1} - 2 S_k \cfrac{|r_k - r_{k-1}|^2}{\left\langle S_k, r_k - r_{k-1} \right\rangle}, \cfrac{r_k}{\left\langle S_k, r_k - r_{k-1} \right\rangle} \right\rangle\) |
Хагер, Жанг |
2005 г. |
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return np.sin(x[0]*x[1])
def grad_f(x):
return np.array([np.cos(x[0]*x[1])*x[1], np.cos(x[0]*x[1])*x[0]])
def conjugate_gradients(starting_point, epsilon):
x = starting_point
iterations = 0
steps = [x]
flag = 0
r = -grad_f(x)
S = r
while True:
phi = lambda a: f(x + a * S)
alpha = minimize_scalar(phi).x
x_new = x + alpha * S
r_new = -grad_f(x_new)
beta = np.dot(r_new, r_new) / np.dot(r, r)
if flag:
beta = 0
S = r_new + beta * S
iterations += 1
steps.append(x_new)
if (np.linalg.norm(x_new - x) < epsilon):
if flag:
x = x_new
break
else:
flag = 1
x = x_new
r = r_new
continue
flag = 0
x = x_new
r = r_new
return x, iterations, np.array(steps)
# Генерация сетки для построения линий уровня
x = np.linspace(-2, 2, 400)
y = np.linspace(-2, 2, 400)
X, Y = np.meshgrid(x, y)
Z = np.sin(X*Y)
# Вызов функции градиентного спуска
starting_point = np.array([0.5, -0.1])
epsilon = 0.001
min_point, iterations, steps = conjugate_gradients(starting_point, epsilon)
# Построение графика
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, Z, levels=10, cmap='viridis')
plt.colorbar()
# Построение линий уровня
contour = plt.contour(X, Y, Z, levels=10, colors='black')
# Подписи к линиям уровня
plt.clabel(contour, colors='k', fmt='%2.1f', fontsize=12)
# Построение шагов алгоритма
for i in range(len(steps)-1):
plt.arrow(steps[i][0], steps[i][1], steps[i+1][0] - steps[i][0], steps[i+1][1] - steps[i][1],
shape='full', lw=1, color='red', length_includes_head=True, head_width=0.05, zorder=4)
# Подпись с количеством итераций
plt.text(-1.5, 1.8, f'Iterations: {iterations}', color='white', fontsize=12)
# Настройка осей и отображение
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Conjugate Gradients Visualization')
print("Сошлись за 2 итерации")
plt.show()
Сошлись за 2 итерации
