Implement gradient method in functional programming style

Question

I implemented the gm for quadratic function in Mathematica, but it looks like more in an OOP style, as the pic shows.

How to rewrite it in Mathematica style? I tried with following code:

gmQuadraticF[A_, b_, x0_, eps_, iter_] :=Module[{x, temp, path},(*g[k] = x[k]+b*)
temp = RecurrenceTable[
       {x[k + 1] == x[k] - N[Norm[A.x[k] + b]^2/(A.x[k] + b).A.(A.x[k] + b)]*(A.x[k] + b), 
       x[0] = x0 }, x, {k, 0, iter_}];(* with other codes ommitted*)];

I came across two problems. The first is flow control with recurrencetable. I want to stop the iteration when Norm[A.x[k]+b] is less than eps (tolerence rate). I tried to calculate the whole list first than apply TakeWhile to select but didn't work Second, the "b" caused the recuurencetable return a numerical exception. Any hint? Or any function should I apply instead?

-------------------------Update-----------------------------------------------

Thanks to Thies Heidecke. I've partially get it done. Two implementations return same result in first several steps. However, the implementation in functional way always stops earlier than OOP style. How to solve that?

gmQuadraticF[A_, b_, x0_, eps_, iter_] :=
 Module[
  {x, g, \[Alpha], NormG, f, v, path}, 
  x = x0;
  g = A.x + b;
  v = 1/2 x0.A.x0 + b.x0;
  \[Alpha] = N[ (Norm[g]^2)/g.A.g];
  NormG = Norm[g];
  f[step_, \[Alpha]_, x_, g_, sqaureNormG_, v_] := {
    step + 1,(*update Iteration counts*)

    N[Norm[g]^2/g.A.g],(*\[Alpha] step size*)

    x - N[Norm[g]^2/g.A.g]*g,(*update x*)

    A.(x - N[Norm[g]^2/g.A.g]*g) + b, (*update gradient*)

    Norm[A.(x - N[Norm[g]^2/g.A.g]*g) + b],(*update gradient norm*)

      1/2 (x - N[Norm[g]^2/g.A.g]*g).A.(x - N[Norm[g]^2/g.A.g]*g) + 
     b.(x - N[Norm[g]^2/g.A.g]*g)(*object function value*)};

  path = NestWhileList[
    f[#[[1]], #[[2]], #[[3]], #[[4]], #[[5]], #[[6]]] &, {0, \[Alpha],
      x, g, NormG, v}, #[[1]] <= iter && #[[5]]^2 >= eps &];

  Prepend[
    path, {"Iteration", "Step Size", "X", "Gradient", "Gradient Norm",
      "Object Value"}] // TableForm
  ]

Answer 1

Here's a slight refactoring of the OP's solution.

Expressions like #[[5]] are hard to read and therefore hard to debug. One can use a function with symbolic names for arguments to make the code easier to understand.

I think it's better to return raw output data of a function and leave it formatting etc. to processing the returned data. So I moved TableForm outside the function gmQuadraticF.

gmQuadraticF[A_, b_, x0_, eps_, iter_] := Module[{g, NormG, f, test},
   f[{step_, α_, x_, g_, squareNormG_, v_}] :=
    With[{stepfactor = Norm[g]^2/g.A.g},
     With[{ng = stepfactor*g},
      With[{xnext = x - ng},
       With[{grad = A.xnext + b},
        {step + 1,                   (*update Iteration counts*)
         stepfactor,                 (*α step size*)
         xnext,                      (*update x*)
         grad,                       (*update gradient*)
         Norm[grad],                 (*update gradient norm*)
         1/2 xnext.A.xnext + b.xnext (*object function value*)}
        ]]]];
   test[{step_, α_, x_, g_, squareNormG_, v_}] := 
    step <= iter && squareNormG >= eps;
   g = A.x0 + b;
   NormG = Norm[g];
   NestWhileList[
    f,
    {0,                   (* initial step number *)
     (NormG^2)/g.A.g,     (*α step size*)
     x0,                  (*initial x*)
     g,                   (*intial gradient*)
     NormG,               (*intial gradient norm*)
     1/2 x0.A.x0 + b.x0}, (*object function value*)
    test]
   ];
gmQuadraticF[TableHeadings] = {"Iteration", "Step Size", "X", 
   "Gradient", "Gradient Norm", "Object Value"};

Example:

SeedRandom[0];
dim = 2;
aa = RandomReal[{-1, 1}, {dim, dim}];
aa = aa\[Transpose].aa;
bb = RandomReal[1, dim];
x0 = Developer`ToPackedArray@{13., -31.};
eps = 1.*^-12;

TableForm[
 gmQuadraticF[aa, bb, x0, eps, 200],
 TableHeadings -> {None, gmQuadraticF[TableHeadings]}
 ]

Here's a Mathematica-like version with precision control:

ClearAll[gmQuadraticF, igmQuadraticF];
Options@gmQuadraticF = {WorkingPrecision -> Automatic};
gmQuadraticF[A_, b_, x0_, eps_, iter_Integer, OptionsPattern[]] := 
  Module[{wp},
   wp = OptionValue[WorkingPrecision];
   wp = wp /. 
     Automatic -> (Precision[{A, b, x0}] /. 
        Infinity -> MachinePrecision);
   Block[{$MaxPrecision = wp, $MinPrecision = wp},
    igmQuadraticF[
     SetPrecision[A, wp],
     SetPrecision[b, wp],
     SetPrecision[x0, wp],
     eps, iter]
    ]
   ];
igmQuadraticF[A_, b_, x0_, eps_, iter_] := 
  Module[{g, NormG, f, test, stepfactor, ng, xnext, grad},
   f[{step_, α_, x_, g_, sqaureNormG_, v_}] := (
     stepfactor = Norm[g]^2/g.A.g;
     ng = stepfactor*g;
     xnext = x - ng;
     grad = A.xnext + b;
     {step + 1,                   (*update Iteration counts*)
      stepfactor,                 (*α step size*)
      xnext,                      (*update x*)
      grad,                       (*update gradient*)
      Norm[grad],                 (*update gradient norm*)
      1/2 xnext.A.xnext + b.xnext (*object function value*)});
   test[{step_, α_, x_, g_, squareNormG_, v_}] := 
    step <= iter && squareNormG >= eps;
   g = A.x0 + b;
   NormG = Norm[g];
   NestWhileList[
    f,
    {0,                   (* initial step number *)
     (NormG^2)/g.A.g,     (*α step size*)
     x0,                  (*initial x*)
     g,                   (*intial gradient*)
     NormG,               (*intial gradient norm*)
     1/2 x0.A.x0 + b.x0}, (*object function value*)
    test]
   ];