3
$\begingroup$

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

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
  ]

enter image description here

enter image description here

$\endgroup$
3
  • 2
    $\begingroup$ Have a look at NestList/NestWhileList. You can encode all the variables that change with each step as a list, e.g. {x,g,steps}. $\endgroup$ Commented Oct 25, 2019 at 11:57
  • 1
    $\begingroup$ Problem has been solved. Somthing wrong with ctria in NestWhileList. It should be #[[5]] not #[[5]]^2. $\endgroup$
    – Hot Pizza
    Commented Oct 26, 2019 at 5:05
  • $\begingroup$ Please post the arguments used in the example in copyable code. Tx. $\endgroup$
    – Michael E2
    Commented Oct 26, 2019 at 14:39

1 Answer 1

5
$\begingroup$

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]}
 ]

Mathematica graphics

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]
   ];
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.