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