Помогите в градиентные методы оптимизации

[B@R5uk]

мат-ламер в сообщении #1682501 писал(а):

Но тут вопрос - а как выбрать этот множитель? Если нет никакой информации насчёт функции, то таки надо делать хотя бы неполную оптимизацию вдоль выбранного направления.

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

$$y=\frac{1}{2}x^TAx,\qquad\qquad A^T=A$$

$$g=\frac{dy}{dx}=Ax$$

$$x_{k+1}=x_k+\alpha g_k$$

Если потребовать минимум функции вдоль направления спуска (стратегия наискорейшего спуска, Steepest Descent) и проварьировать по параметру, то придём к формулам:

$$\alpha=\frac{g_k^Tg_k}{g_k^TAg_k}=\frac{x_k^TA^2x_k}{x_k^TA^3x_k}$$

$$x_{k+1}=x_k-\frac{g_k^Tg_k}{g_k^TAg_k}g_k=x_k-\frac{x_k^TA^2x_k}{x_k^TA^3x_k}Ax_k$$

Если же потребовать минимум модуля (или, что более удобно квадрата модуля) градиента (стратегия минимального градиента, Minimum Gradient), то получится так:

$$\alpha=\frac{g_k^TAg_k}{g_k^TA^2g_k}=\frac{x_k^TA^3x_k}{x_k^TA^4x_k}$$

$$x_{k+1}=x_k-\frac{g_k^TAg_k}{g_k^TA^2g_k}g_k=x_k-\frac{x_k^TA^3x_k}{x_k^TA^4x_k}Ax_k$$

Очевидно, оператор шага S в обоих стратегиях не является линейным (является дробно-рациональным?):

$$S_{SD}\bigl(x\bigr)=x-\frac{x^TA^2x}{x^TA^3x}Ax,\qquad\qquad S_{MG}\bigl(x\bigr)=x-\frac{x^TA^3x}{x^TA^4x}Ax$$

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

Сразу говорю, на картинке нарисована неправда: проекции окружности после первой итерации перевёрнуты по вертикали. Я это сделал, чтобы продемонстрировать, что существует выделенное направление, векторы вдоль которого при применении оператора шага S, во-первых укорачиваются медленно, а во-вторых не меняют своего направления (с учётом смены знака по координате y). Траектории, для них — это (жёлтая) линия, вторая против часовой стрелки от горизонтальной белой на левом графике и первая, розовая — на втором. Левый график соответствует стратегии максимального спуска, правый — минимального градиента. Синими окружностями помечены точки, в которые отображается исходная единичная окружность стратегией с фиксированным шагом. Получаются окружности, а не эллипсы, потому что шаг выбран оптимальным, а это даёт одинаковые собственные числа для оператора S, соответствующие минимальному и максимальному собственным числам матрицы A квадратичной формы. Поскольку задача двумерная, других чисел быть не может. Поскольку я перевернул результат первой итерации по вертикали, траектории, отмеченные разноцветными тонкими линиями, перестали быть ортогональными линиям уровня целевой функции, если первую итерацию перевернуть назад, то будет всё нормально. Видно, что методы наискорейшего спуска и минимального градиента в самом худшем случае чуть-чуть хуже метода с фиксированным градиентом с оптимально подобранным шагом. По крайней мере в двумерии. Так же, если очень присмотреться, то можно заметить, что векторы в процессе итераций прижимаются к наихудшему направлению спуска. Ну и очевидно, что собственные вектора оператора A схлопываются в ноль за одну итерацию.

У меня есть подозрение, что в трёхмерии и пространствах большей размерности может получиться что-то более интересное. Надо будет попробоать.

Программа:

Используется синтаксис Matlab M

%   File: study_GD_2.m
%   dxdy.ru/topic160189.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
%format long
format short

num_iter = 3;

%   Diagonal bilinear form
A = diag ([1 3]);

func_objective = @(x, A) 0.5 * x' * A * x;
func_gradient  = @(x, A)            A * x;

%   Starting circle and space for data
ppp = 2 * pi * (0 : 0.001 : 1);
num_points = numel (ppp);
xxx = zeros (num_iter, num_points);
yyy = zeros (num_iter, num_points);
xxx (1, :) = cos (ppp);
yyy (1, :) = sin (ppp);
xxx2 = xxx;
yyy2 = yyy;
xxx3 = xxx;
yyy3 = yyy;

%   Gradient descent: steepest descent strategy
for k = 2 : num_iter
    for l = 1 : num_points
        x = [xxx(k - 1, l); yyy(k - 1, l)];
        g = func_gradient (x, A);
        xx = x - g * g' * g / (g' * A * g);
        xxx (k, l) = xx (1);
        yyy (k, l) = xx (2);
    end
end

%   Gradient descent: minimum gradient strategy
for k = 2 : num_iter
    for l = 1 : num_points
        x = [xxx2(k - 1, l); yyy2(k - 1, l)];
        g = func_gradient (x, A);
        ag = A * g;
        xx = x - g * g' * ag / (ag' * ag);
        xxx2 (k, l) = xx (1);
        yyy2 (k, l) = xx (2);
    end
end

%   Gradient descent: fixed step strategy (best step selected)
step = 2 / trace (A);
for k = 2 : num_iter
    for l = 1 : num_points
        x = [xxx3(k - 1, l); yyy3(k - 1, l)];
        g = func_gradient (x, A);
        xx = x - step * g;
        xxx3 (k, l) = xx (1);
        yyy3 (k, l) = xx (2);
    end
end

%   Objective function levels values
xx = -1 : 0.01 : 1;
zz = zeros (size (xx, 2));
for k = 1 : numel (xx)
    for l = 1 : numel (xx)
        zz (l, k) = func_objective ([xx(k); xx(l)], A);
    end
end

cl = 'wmycg';

%   SD graph
subplot (121)
contourf (xx, xx, zz, -1 : 0.1 : 2)
hold on
for k = 1 : num_iter
    plot (xxx3 (k, :), yyy3 (k, :), 'LineWidth', 3)
end
for k = 1 : num_iter
    %   Odd circle image is flipped vertically
    plot (xxx (k, :), (-1) ^ (k - 1) * yyy (k, :), 'r', 'LineWidth', 3)
end
for k = 0 : 9
    n = 1 + k * 26;
    for l = 1 : 2
        %   Odd circle image is flipped vertically
        plot ([xxx(l, n) xxx(l + 1, n)], (-1) ^ l * [-yyy(l, n) yyy(l + 1, n)], cl (mod (k, numel (cl)) + 1))
    end
end
n = 251;
l = 1;
%   Vertical eigenvector trajectory
plot ([xxx(l, n) xxx(l + 1, n)], (-1) ^ l * [-yyy(l, n) yyy(l + 1, n)], 'w')
hold off
axis equal
xlim ([-1 1])
ylim ([-1 1])

%   MG graph
subplot (122)
contourf (xx, xx, zz, -1 : 0.1 : 2)
hold on
for k = 1 : num_iter
    plot (xxx3 (k, :), yyy3 (k, :), 'LineWidth', 3)
end
for k = 1 : num_iter
    %   Odd circle image is flipped vertically
    plot (xxx2 (k, :), (-1) ^ (k - 1) * yyy2 (k, :), 'r', 'LineWidth', 3)
end
for k = 0 : 7
    n = 1 + k * 32;
    for l = 1 : 2
        %   Odd circle image is flipped vertically
        plot ([xxx2(l, n) xxx2(l + 1, n)], (-1) ^ l * [-yyy2(l, n) yyy2(l + 1, n)], cl (mod (k, numel (cl)) + 1))
    end
end
n = 251;
l = 1;
%   Vertical eigenvector trajectory
plot ([xxx2(l, n) xxx2(l + 1, n)], (-1) ^ l * [-yyy2(l, n) yyy2(l + 1, n)], 'w')
hold off
axis equal
xlim ([-1 1])
ylim ([-1 1])

%   Possible colour maps:
%   hsv gray hot cool bone copper pink flag
colormap (pink)
 

В программе несколькими постами ранее, я облажался и вместо zz (l, k) = func_objective ... написал zz (k, l) = func_objective ... Из-за этого на картинке в том же посте (и при работе самой программы) линии уровня целевой функции отображены неправильно: они отражены вдоль диагонали первого квадранта. Другими словами, иксы и игреки перепутаны местами. Надеюсь, в этой программе нет опечаток.

[B@R5uk]

Осваиваю тут потихоньку Jorge Nocedal, Stephen J. Wright - Numerical Optimization. zgemm, спасибо за наводку в соседней теме, отличная книжка. В частности прочитал последним, что такое условия Вольфе, откуда у них ноги растут и почему они лучше критерия Голдстейна, хотя оба преследуют одну и ту же цель: обеспечить глобальную сходимость алгоритмов (то бишь к локальному экстремуму, а не просто к стационарной точке). В книге это шикарно объясняется с наглядными картинками. Недостаток критерия Голдстейна в том, что он может исключить из потенциального шага алгоритма все локальные минимумы. Тем не менее, его часто используют с Ньютоновскими методами, но не используют с квази-Ньютоновскими. Такого рода практические рекомендации тоже достоинство книжки.

В книге так же делается упор на разделение одной итерации любого алгоритма на два ключевых этапа:
1) выбор направления спуска;
2) выбор шага вдоль направления.
Несмотря на контринтуитивность, направление наискорейшего спуска, даваемое антиградиентом, далеко не всегда является оптимальным в плане скорости сходимости. Это особенно актуально, когда матрица Гёссе функции имеет слишком высокое число обусловленности. В связи с этим встаёт вопрос: а как это направление можно выбрать более оптимально? Это порождает серию различных алгоритмов, которые различаются на этом этапе итерации. С выбором величины шага тоже не всё так просто, в связи с чем условия Вольфе и еже с ними. Забавно, что совершенно разные алгоритмы могут подходить к этому второму этапу совершенно одинаково, а разные реализации одного и того же алгоритма совершенно по-разному.

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

[zgemm]

B@R5uk в сообщении #1682817 писал(а):

Это особенно актуально, когда матрица Гёссе функции имеет слишком высокое число обусловленности.

Бороться с высоким числом обусловленности и публиковать на этом уйму статей - это целая эпоха в конце 90-х начале 2000х. В большинстве случаев все заканчивалось тем, что матрицу сдвигали вверх на единичную. Сейчас, если говорить от BFGS, то принято решать сразу с обратной, а если норма обратной матрицы сильно растет, то или скалировать, или ограничивать норму факторов.

У меня в моем солвере сделано так:
1. если новый добавляемый член в обратную больше на константу большую единице по норме всего остального, то это плохо, надо делать рестарт.
2. если больше примерно 1/10, надо реортагонализовать и попытаться объединить,
3. иначе - работать по классике.

B@R5uk в сообщении #1682817 писал(а):

1) выбор направления спуска;
2) выбор шага вдоль направления.
... В связи с этим встаёт вопрос: а как это направление можно выбрать более оптимально?

Никогда не пользуйте только одно направление. У Вас всегда должна быть возможность переключения на что-то другое. Я как-то выше уже писал. У меня в солвере сделано так
1. Ньютон (если у пользователя есть callback),
2. Квазиньютон,
3. обычный градиент, ортогональный предыдущему шагу,
4. градиент.

Второй момент - одномерная минимизация.
Квазиньютон очень любит, чтобы эта минимизация была хорошей, а вот сам дает только хорошее направление, но длина обычно бывает очень не правильной.

Выбор метода направления и выбор способа и числа шагов в одномерной минимизации - через минмаличтическую ANN, которая учится на предыдущих итерациях и, если пользователь предоставил, с предыдущих вызов.

[B@R5uk]

zgemm, спасибо.

[B@R5uk]

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

У меня такой вопрос (чтобы не создавать отдельную тему): существует ли численный аналог метода множителей Лагранжа в духе метода Ньютона? А то вдруг я тут изобретаю себе велосипед, когда уже всё придумано. Просто, когда ищется условный экстремум задачи

$$\begin{cases}\widehat{x}=\underset{x}{\arg\min}f(x)\\g(x)=0\end{cases}$$

$$f:\mathbb{R}^n\to\mathbb{R},\qquad g:\mathbb{R}^n\to\mathbb{R}^m\qquad n>m$$

единственный подход, который гарантированно работает из того, что я знаю — это параметризовать каким-то образом пространство

$$g(x)=0$$

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

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

Есть ещё метод внутренней точки, но он больше заточен под поиск экстремума в выпуклой области, заданной неравенствами. Наверное, его можно прикрутить и к задаче выше путём хитрого выбора сходящейся последовательности штрафных функций, но это в целом будет сомнительно и некошерно. Подскажите, пожалуйста, что ещё в этой области разработано и используется?

[dsge]

Если ограничения линейные, то весьма эффективен метод Ньютона с некоторой добавкой, обеспечивающей выполнение ограничений.

Если ограничения нелинейные, то можно применить также Sequential Quadratic Programming (SQP). На каждом шаге целевая функция и ограничения раскладываются в ряды Тэйлора и решается задача квадратичного программирования.

[B@R5uk]

dsge, спасибо за наводку! У Носедаля как раз про SQP есть целая глава, почитаю. Нужно ли перед SQP освоить QP или это две отдельные вещи?

[dsge]

B@R5uk в сообщении #1688752 писал(а):

Нужно ли перед SQP освоить QP или это две отдельные вещи?

QP - это составная часть алгоритма SQP.