Оптимизация функции многих переменных
   
Метод Ньютона Метод Нелдера-Мида
   


Метод Ньютона

Задача: Методом Ньютона вычислить минимум функции двух переменных.

Для примера, возьмем функцию Розенброка:
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 (по Нелдеру) !!!