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

[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, спасибо.