Метод Ньютона
Задача: Методом Ньютона вычислить минимум функции двух переменных.
Для примера, возьмем функцию Розенброка:
f(x, y) = 100(y - x2)2 + (1 - x)2
Очевидно, что минимум достигается в точке (1, 1).
Алгоритм: Пусть (x0, y0) – начальное приближение для точки минимума.
Очередное приближение вычисляется по формуле:
|xk+1| |xk| |h1,1 h1,2|-1 |ƒ(xk)|
| | = | | - | | | |
|yk+1| |yk| |h2,1 h2,2| |ƒ(yk)|
где (ƒx, ƒy)' – градиент функции ƒ(x, y), а hi, j - элементы матрицы Гессе (матрицы вторых производных):
h1, 1 = ƒxx
h1, 2 = h2, 1 = ƒxy
h2, 2 = ƒyy
В покомпонентном виде приведенная выше формула имеет вид:
xk+1 = xk - (h2,2 * ƒxk - h1,2 * ƒyk) / det(h)
yk+1 = yk - (- h2,1 * ƒxk + h1,1 * ƒyk) / det(h)
где det(h) - определитель матрицы h.
Cчет ведется итерациями до тех пор, пока два последовательных приближения не будут отличаться по третьей норме больше, чем на ε.
Программа, вычисляющая минимум функции методом Ньютона, приведена ниже:
function _f(x, y: double): double;
begin
_f := (100.0*sqr(y-x*x) + sqr(1.0-x));
end;
function _Fx(x, y: double): double;
begin
_Fx := (-400.0*(y-x*x)*x-2.0* (1.0-x));
end;
function _Fy(x, y: double): double;
begin
_Fy := (200.0*(y-x*x));
end;
function _H11(x, y: double): double;
begin
_H11 := (-400.0*(y-x*x)+ 800.0*x*x+2.0);
end;
function _H12(x, y: double): double;
begin
_H12 := (-400.0*x);
end;
function _H22(x, y: double): double;
begin
_H22 := (200.0);
end;
const
eps = 0.001;
var
x0, y0, f: double;
fx, fy, h11, h12, h22: double;
det, dx, dy: double;
iter: integer;
begin
iter := 0;
x0 := -1.2; y0 := 1;
f := _f(x0, y0);
writeln(' iter=', iter:3, ' x=', x0:9:4,
' y=', y0:9:4, ' f(x,y)=', f:9:4);
repeat
inc(iter);
h11 := _H11(x0,y0);
h12 := _H12(x0,y0);
h22 := _H22(x0,y0);
fx := _Fx(x0,y0); fy := _Fy(x0,y0);
det := h11*h22-h12*h12;
dx := (h22*fx-h12*fy)/det;
dy := (h11*fy-h12*fx)/det;
x0 := x0 - dx; y0 := y0 - dy; f := _f(x0,y0);
writeln(' iter=', iter:3, ' x=', x0:9:4,
' y=', y0:9:4, ' f(x,y)=', f:9:4);
until (sqrt(sqr(dx)+sqr(dy)) < eps);
end.Результат работы программы:
iter= 0 x= -1.2000 y= 1.0000 f(x,y)= 24.2000
iter= 1 x= -1.1753 y= 1.3807 f(x,y)= 4.7319
iter= 2 x= 0.7631 y= -3.1750 f(x,y)=1411.8452
iter= 3 x= 0.7634 y= 0.5828 f(x,y)= 0.0560
iter= 4 x= 1.0000 y= 0.9440 f(x,y)= 0.3132
iter= 5 x= 1.0000 y= 1.0000 f(x,y)= 0.0000
iter= 6 x= 1.0000 y= 1.0000 f(x,y)= 0.0000
Как видим, результат найден всего за 6 итераций.
Но у этого метода есть явный недостаток: приходится вручную вычислять k первых производных, и k*(k+1)/2 вторых
производных, что при размерности пространства, скажем, равном 7, приведет к вычислению 7 производных первого порядка и 28 - второго;
шансы совершить ошибку просто огромны...
Метод Нелдера-Мида
Подобного недостатка лишен другой метод оптимизации - метод Нелдера-Мида.
Этот метод является одним из наиболее быстрых и наиболее надежных не градиентных методов многомерной оптимизации.
Идея этого метода состоит в следующем: в пространстве оптимизируемых параметров генерируется случайная точка. Затем строится
n-мерный симплекс с центром в этой точке, и длиной стороны L. Далее в каждой из вершин симплекса вычисляется
значение оценки. Выбирается вершина с наибольшей оценкой. Вычисляется центр тяжести остальных n вершин. Проводится
оптимизация шага в направлении от наихудшей вершины к центру тяжести остальных вершин.
Эта процедура повторяется до тех пор, пока не окажется, что оптимизация не изменяет положения вершины. После этого выбирается
вершина с наилучшей оценкой и вокруг нее снова строится симплекс с меньшими размерами (например L/2).
Процедура продолжается до тех пор, пока размер симплекса, который необходимо построить, не окажется меньше требуемой точности.
Программа, реализующая поиск минимума функции Розенброка данным методом:
const
eps = 0.000001; { Точность определения минимума }
alfa = 1.0; { Коэффициент отражения }
beta = 0.5; { Коэффициент сжатия }
gamma = 2.0; { Коэффициент растяжения }
t = 0.2;
n_max = 6; { Макс. количество переменных }
var
x, y: array[0 .. n_max] of double;
fh, fl, f4, f5, f6: double;
h, l, it: integer;
{ *** Оптимизируемая функция *** }
function f(k: integer): double;
var r1, r2: double;
begin
r1 := y[k]-x[k]*x[k]; r2 := 1-x[k];
f := 100*r1*r1+r2*r2;
end;
function maxf: double;
var f1, f2, r: double;
begin
f1 := f(1); f2 := f(2); r := f(0);
h := 0;
if r < f1 then begin
r := f1; h := 1;
end;
if r < f2 then begin
r := f2; h := 2;
end;
maxf := r;
end;
function minf: double;
var f1, f2, r: double;
begin
f1 := f(1); f2 := f(2); r := f(0);
l := 0;
if f1 < r then begin
r := f1; l := 1;
end;
if f2 < r then begin
r := f2; l := 2;
end;
minf := r;
end;
label TheEnd;
var
i, flag: integer;
r, r1, x0, y0, x1, y1: double;
begin
x0 := -1.2; y0 := 1;
x[0] := x0-0.5*t; y[0] := y0-t*sqrt(3)/6;
x[1] := x0; y[1] := y0; r := f(1); y[1] := y0+t*sqrt(3)/3;
x[2] := x0+0.5*t; y[2] := y[0];
it := 0;
writeln(' it=', it:5, ' x=', x0:8:4, ' y=', y0:8:4, ' f=', r:8:4);
repeat
fh := maxf; fl := minf;
x[3] := 0.5*(x[0]+x[1]+x[2]-x[h]);
y[3] := 0.5*(y[0]+y[1]+y[2]-y[h]);
x[4] := (1+alfa)*x[3]-alfa*x[h];
y[4] := (1+alfa)*y[3]-alfa*y[h];
f4 := f(4);
if (f4 < fl) then begin
x[5]:=(1-gamma)*x[3]+gamma*x[4];
y[5]:=(1-gamma)*y[3]+gamma*y[4];
f5:=f(5);
if (f5 < fl) then begin x[h]:=x[5]; y[h]:=y[5]; end
else begin x[h]:=x[4]; y[h]:=y[4]; end;
goto TheEnd;
end;
flag := 0;
for i := 0 to pred(3) do begin
if ((i <> h) and (f4 > f(i))) then inc(flag);
end;
if (flag = 2) then begin
x[6]:=beta*x[h]+(1-beta)*x[3]; y[6]:=beta*y[h]+(1-beta)*y[3];
x[h]:=x[6]; y[h]:=y[6];
goto TheEnd;
end;
if (f(4) < fh) then begin
for i := 0 to pred(3) do begin
x[i]:=0.5*(x[i]+x[l]); y[i]:=0.5*(y[i]+y[l]);
end;
end
else begin
x[h]:=x[4]; y[h]:=y[4];
end;
TheEnd:;
r := 0;
for i := 0 to pred(3) do begin
r1 := f(i)-f(3); r := r + r1*r1;
end;
r:=sqrt(r/3);
inc(it);
r1:=(f(0)+f(1)+f(2))/3;
x0:=(x[0]+x[1]+x[2])/3;
y0:=(y[0]+y[1]+y[2])/3;
if it mod 20 = 0 then
writeln(' it=', it:5, ' x=', x0:8:4, ' y=', y0:8:4, ' f=', r:8:4);
until (r < eps);
writeln(' it=', it:5, ' x=', x0:8:4, ' y=', y0:8:4, ' f=', r:8:4);
end.Запустим эту программу и посмотрим на возвращаемые результаты:
it= 0 x= -1.2000 y= 1.0000 f= 24.2000
it= 20 x= -0.5415 y= 0.2818 f= 0.0164
it= 40 x= -0.2912 y= 0.1231 f= 0.0291
it= 60 x= 0.0123 y= -0.0106 f= 0.0060
it= 80 x= 0.1233 y= -0.0075 f= 0.0039
it= 100 x= 0.4136 y= 0.1820 f= 0.1305
it= 120 x= 0.6089 y= 0.3575 f= 0.0019
it= 140 x= 0.6456 y= 0.4065 f= 0.0095
it= 160 x= 0.8838 y= 0.7762 f= 0.0017
it= 179 x= 1.0003 y= 1.0006 f= 0.0000
Как видим, с одной стороны метод Нелдера-Мида имеет преимущество перед методом Ньютона
(для того, чтобы иайти минимум, скажем, вот такой функции:
f(x, y) = (1/2)x2 - 4xy + 9y2 + 3x - 14y + (1/2)
, надо всего лишь заменить оптимизируемую
функцию:
function f(k: integer): double;
begin
f := (1/2)*sqr(x[k]) - 4*x[k]*y[k] + 9*sqr(y[k]) + 3*x[k] - 14*y[k] + (1/2);
end;
, и программа выдает верный результат без необходимости высчитывать производные первого и второго порядков:
it= 0 x= -1.2000 y= 1.0000 f= -2.5800
it= 20 x= -0.3740 y= 0.6605 f= 0.0017
it= 40 x= 0.8577 y= 0.9712 f= 0.0044
it= 53 x= 1.0012 y= 1.0005 f= 0.0000
), но посмотрите, насколько отличается число итераций, необходимое для получения
результата в этих двух методах! 6 итераций (по Ньютону) против
179 (по Нелдеру) !!!
