# Implement gradient method in functional programming style

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 = 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[#[], #[], #[], #[], #[], #[]] &, {0, \[Alpha],
x, g, NormG, v}, #[] <= iter && #[]^2 >= eps &];

Prepend[
"Object Value"}] // TableForm
]  • 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}. Oct 25, 2019 at 11:57
• Problem has been solved. Somthing wrong with ctria in NestWhileList. It should be #[] not #[]^2. Oct 26, 2019 at 5:05
• Please post the arguments used in the example in copyable code. Tx. Oct 26, 2019 at 14:39

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

Expressions like #[] 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},
{step + 1,                   (*update Iteration counts*)
stepfactor,                 (*α step size*)
xnext,                      (*update x*)
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*)
1/2 x0.A.x0 + b.x0}, (*object function value*)
test]
];


Example:

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

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

ClearAll[gmQuadraticF, igmQuadraticF];
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},
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;
{step + 1,                   (*update Iteration counts*)
stepfactor,                 (*α step size*)
xnext,                      (*update x*)
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*)
`