Gradientenverfahren

Beim Gradientenverfahren handelt es sich um ein Verfahren zur Minimierung einer Funktion f vom Typ Rn -> R.

Gradientenverfahren beginnen bei einem Startwert als Näherung und verbessern diesen schrittweise. Für eine Näherung x wird der Gradient grad(x) von f in x bestimmt. Das klassische Gradientenverfahren vermutet das Minimum der Funktion in Richtung d des stärksten Abstiegs, d.h. d = -grad(x). Als neue Näherung an das Minimum wird x+ = x + a * d gesetzt. Dabei ist a ist die Schrittweite, die z.B. mit der Armijoregel bestimmt werden kann.

Als Pseudo Code sieht das Verfahren wie folgt aus:

procedure minimize(nDimensions:vectorIndex;
                   const param;
                   f:n1function; (* Funktion von param und vector *)
                   var x:vector;
                   var y:extended;
                   nMaxSteps:longint;
                   var nSteps:longint;
                   epsilon:extended;
                   var fail:boolean);

  procedure vectorCopy(const v:vector;
                       var w:vector);
  var i:vectorIndex;
  begin
    for i:=1 to nDimensions do
      w[i]:=v[i];
  end;

  procedure gradient(const x0:vector;
                     const y0:extended;
                     var g:vector);
  const nMaxPoints=2;
  var nPoints:integer;
      dxiList,dyList:array [1..2] of extended;

    procedure savePoint(const dx,dy:extended);
    begin
      if dy <> 0 then begin
        inc(nPoints);
        dxiList[nPoints]:=dx;
        dyList[nPoints]:=dy;
      end;
    end;

  const minh=1e-15;
        maxh=1e-7;
  var i:vectorIndex;
      x:vector;
      h,xi,y:extended;
  begin
    vectorCopy(x0, x);
    for i:=1 to nDimensions do begin
      h:=minh;
      xi:=x[i];
      nPoints:=0;
      (* Suche nMaxPoints-viele Stellen in der Nähe
         von xi, für die gilt (y - y0) / dx <> 0. *)
      repeat
        x[i]:=xi + h;
        f(param, x, y, fail);
        if fail then begin
          x[i]:=xi - h;
          f(param, x, y, fail);
          if fail then exit;
          savePoint(-h, y - y0);
        end else begin
          savePoint(h, y - y0);
        end;
        h:=h * 2;
      until (nPoints >= nMaxPoints)
         or (h > maxh);
      (* Vektor x restaurieren *)
      x[i]:=xi;
      (* Die (nPoints+1)-vielen Stellen werden als
         Stützpunkte eines Polynoms verwendet.
         g[i] ist die Steigung des Polynoms für
         x=0. *)
      case nPoints of
        0: g[i]:=0;
        (* Beachten: (0,0) ist immer ein Punkt des Polynoms *)
        1: g[i]:=dyList[1] / dxiList[1];
        2: g[i]:=(dyList[2] * dxiList[1] / dxiList[2]
                 -dyList[1] * dxiList[2] / dxiList[1])
                 / (dxiList[1] - dxiList[2]);
      end;
    end;
  end;

  function armijo(const alpha:extended;
                  const x0,d:vector):extended;
  var i:vectorIndex;
      x:vector;
      y:extended;
  begin
    for i:=1 to nDimensions do begin
      x[i]:=x0[i] - alpha * d[i];
    end;
    f(param, x, y, fail);
    if fail then exit;
    armijo:=y;
  end;

const sigma=0.33;
      p=2;
var nRestSteps:longint;
    y_old,r,beta,a:extended;
    x_old,g:vector;
    i:vectorIndex;
begin
  nRestSteps:=nMaxSteps;
  f(param, x, y, fail);
  if fail then exit;
  repeat
    y_old:=y;
    vectorCopy(x, x_old);
    (* Gradient bestimmen *)
    gradient(x, y, g);
    if fail then exit;
    (* Schrittweite bestimmen *)
    begin
      (* Richtungsableitung in Richtung des Gradienten bestimmen *)
      r:=0;
      for i:=1 to nDimensions do begin
        r:=r + sqr(g[i]);
      end;
      (* Anfangsschrittweite *)
      beta:=1;
      if r > 1e-30 then begin
        (* beta * r darf nicht zu klein sein *)
        a:=armijo(beta, x, g);
        while (not fail)
          and ((y - a) / (beta * r) >= sigma) do begin
          beta:=beta * p;
          a:=armijo(beta, x, g);
        end;
        (* beta halbieren bis Bedingung erfüllt oder beta * r zu klein ist *)
        while (beta > 1e-30)
          and (((y - a) / (beta * r) < sigma)
               or fail) do begin
          beta:=beta / p;
          a:=armijo(beta, x, g);
        end;
      end;
    end;
    (* Schritt durchführen *)
    for i:=1 to nDimensions do begin
      x[i]:=x[i] - beta * g[i];
    end;
    f(param, x,y, fail);
    if fail then exit;
    dec(nRestSteps);
  until (nRestSteps <= 0)
     or (y < epsilon)
     or (y >= y_old);
  if y > y_old then begin
    vectorCopy(x_old, x);
    y:=y_old;
  end;
  nSteps:=nMaxSteps - nRestSteps;
end;