Реализация градиентного спуска в MatLab
Функция, для которой осуществляется поиск минимума (самая простая):
,
гладкая (по крайней мере ), т.е. всюду дважды непрерывно дифференцируема.
Постановка задачи:
.
Исходный код:
clear all; clc;
% Значения коэффициентов
b0 = 0.8;
b1 = 1.3;
b11 = -6.1;
b12 = -0.6;
b2 = 1.7;
b22 = -1.7;
g = 0.05; % постоянная шага
d = 0.01; % дельта
% Начальная точка
x1 = -9;
x2 = 8;
k = 1; % Счетчик шагов
kmax = 100; % Предельное число шагов,
% задается для предотвращения зацикливания
% Массивы для хранения промежуточных координат
x1trace = [x1];
x2trace = [x2];
i = 2;
while k < kmax
% Спуск по обеим координатам сразу
gr1 = b1 + b12*x2 + 2*b11*x1;
gr2 = b2 + b12*x1 + 2*b22*x2;
x1 = x1 + g*gr1;
x2 = x2 + g*gr2;
% Сохранение координат
x1trace(i) = x1;
x2trace(i) = x2;
i = i + 1;
% Проверка условия останова
if sqrt(gr1^2 + gr2^2) <= d;
break; % Выход из цикла в случае выполнения условия
end
k = k + 1;
end
% Построение графика
x = -10:0.1:10;
y = -10:0.1:10;
[X, Y] = meshgrid(x, y);
Z = b0 + b1*X + b2*Y + b12*X.*Y + b11*X.^2 + b22*Y.^2;
[C, h] = contour(X, Y, Z);
clabel(C, h); % Отображение меток на линиях уровня
hold on;
|
|
plot(x1trace, x2trace, '-x');
% Вывод начальной точки на график
text(x1trace(1) + 0.2, x2trace(1) + 0.5, 'M0');
% Вывод решения на график
text(x1 + 2, x2,...
strvcat(['x1 = ' (num2str(x1))],...
['x2 = ' (num2str(x2))],...
['k = ' (num2str(k))]));
Должен получиться вот такой график (приблизительно)
Задание: запустить код в программе и получить график. График прислать мне немедленно!
Более сложная функция:
,
Значения коэффициентов:
gr1 = 2*x + 4*cos(x + 5*y); % производная по х
gr2 = 6*y + 20*cos(x + 5*y); % производная по y
Исходный код:
clear all; clc;
% Значения коэффициентов
g = 0.05; % постоянная шага
d = 0.01; % дельта
% Начальная точка
x1 = -8;
x2 = 8;
k = 1; % Счетчик шагов
kmax = 100; % Предельное число шагов,
% задается для предотвращения зацикливания
% Массивы для хранения промежуточных координат
x1trace = [x1];
x2trace = [x2];
i = 2;
while k < kmax
% Спуск по обеим координатам сразу
gr1 = 2*x1 + 4*cos(x1 + 5*x2);
gr2 = 6*x2 + 20*cos(x1 + 5*x2);
x1 = x1 - g*gr1;
x2 = x2 - g*gr2;
% Сохранение координат
x1trace(i) = x1;
x2trace(i) = x2;
i = i + 1;
% Проверка условия останова
if sqrt(gr1^2 + gr2^2) <= d;
break; % Выход из цикла в случае выполнения условия
end
k = k + 1;
end
% Построение графика
x = -10:0.1:10;
y = -10:0.1:10;
[X, Y] = meshgrid(x, y);
Z = 4*sin(X + 5*X) + X.^2 + 3*Y.^2;
[C, h] = contour(X, Y, Z);
clabel(C, h); % Отображение меток на линиях уровня
hold on;
plot(x1trace, x2trace, '-x');
% Вывод начальной точки на график
text(x1trace(1) + 0.2, x2trace(1) + 0.5, 'M0');
% Вывод решения на график
text(x1 + 2, x2,...
strvcat(['x1 = ' (num2str(x1))],...
['x2 = ' (num2str(x2))],...
['k = ' (num2str(k))]));
Задание: запустить код в программе и получить график. График прислать мне немедленно!
Градиентные методы хороши для гладких функций. Если же функция имеет локальные максимумы (в данном случае минимумы), то они могут проваливаться в эти минимумы. Подбирая шаг можно добиться того, чтобы направление спуска "выпрыгивало" из этих ям, но все равно, в данном случае, чтобы алгоритм сошелся, нужно сильно задрать вверх порог останова, иначе градиент будет кружить бесконечно вокруг одного из максимумов, к которому он придет.
|
|
Исходный код:
clear all; clc;
% f(x,y)=x^2+3y^2+4sin(x+5y)
% Значения коэффициентов
b0 = 1;
b1 = 3;
b2 = 4;
b3 = 5;
g = 0.04; % постоянная шага
d = 3; % дельта
% Начальная точка
x1 = -9;
x2 = 8;
k = 1; % Счетчик шагов
kmax = 100; % Предельное число шагов,
% задается для предотвращения зацикливания
% Массивы для хранения промежуточных координат
x1trace = [x1];
x2trace = [x2];
i = 2;
while k < kmax
% Спуск по обеим координатам сразу
% gr1 = b1 + b12*x2 + 2*b11*x1;
% gr2 = b2 + b12*x1 + 2*b22*x2;
gr1 = 2*b0*x1 + b2*cos(x1+b3*x2);
gr2 = 2*b1*x2 + b2*b3*cos(x1+b3*x2);
x1 = x1 - g*gr1;
x2 = x2 - g*gr2;
% Сохранение координат
x1trace(i) = x1;
x2trace(i) = x2;
i = i + 1;
% Проверка условия останова
sqrt(gr1^2 + gr2^2)
if sqrt(gr1^2 + gr2^2) <= d;
break; % Выход из цикла в случае выполнения условия
end
k = k + 1;
end
% Построение графика
x = -10:0.1:10;
y = -10:0.1:10;
[X, Y] = meshgrid(x, y);
% Z = b0 + b1*X + b2*Y + b12*X.*Y + b11*X.^2 + b22*Y.^2;
Z = b0*X.^2 + b1*Y.^2+ b2*sin(X+b3*X);
[C, h] = contour(X, Y, Z);
clabel(C, h); % Отображение меток на линиях уровня
hold on;
plot(x1trace, x2trace, '-x');
% Вывод начальной точки на график
text(x1trace(1) + 0.2, x2trace(1) + 0.5, 'M0');
% Вывод решения на график
text(x1 + 2, x2,...
strvcat(['x1 = ' (num2str(x1))],...
['x2 = ' (num2str(x2))],...
['k = ' (num2str(k))]));
Тоже прогнать программу, прислать график.