Question: Programming Anderson Acceleration in Maple

Hello! 

I am currently programming the Anderson Acceleration for fixed point iteration in Maple. The algorithm comes from Walker et al. 2011 (ANDERSON ACCELERATION FOR FIXED-POINT ITERATIONS) (if you are interested in this problem please read the short paper). The code that Walker supplies runs fine in Matlab, with qrdelete as a built-in function. However in Maple I have decided to skip operations on QR, and instead opted to create a new QR every time I increase or decrease the amount of residuals df. However, here comes the kicker, somehow Maple decides to turn a vector or matrix into a procedure when I Concatenate, or DeleteColumn. I could really use a working Anderson Acceleration code for my research (my research is not based on AA or root solvers in general, but spectral methods). I will paste the entire code here. This is my attempt at getting Walker's original Matlab code to work in Maple.

I could use some pointers and tips. Can you program this in a more efficient way? I would be happy to learn. *Notice that the code works for a host of different equations, but not all. Feel free to let me know if this question or inquiry is inappropriate and I will of course delete the post.

restart:
Digits:=12:
Runtime:=time():
with(LinearAlgebra):
with(plots):
with(orthopoly):
with(ArrayTools):
    

phix := Vector(2):
X    := Vector(2):
phix[1] := cos(x[2]):
phix[2] := 3*cos(x[1]):
X[1] := 0.0:
X[2] := 0.0:

#Code AndersonAcceleration.
AndersonAcceleration:=proc(N,phi,X0)
global x, xS, here1, here2;
local mMax, itmax, atol, rtol, droptol, beta, AAstart, res_hist, df, DGg, gamma;
local DG, mAA, iter, gval,fval,res_norm, tol, f_old, g_old, y, i, k, Q, R, QRg, dfT, DGT;

(*
%------------------------------------------------------------------------
% Function xS = AndersonAcceleration(N,phi,x0).
%
% Fixed-point iteration with or without Anderson acceleration.
% 'phi' is the fixed point iteration map and 'xS' is the 
% solution, so that xS = phi(xS).
%
% Input arguments:
%    X0 = Initial value solution. Note: In this function this variable 
%        must be a column vector.
%       1. 'mMax' = maximum number of stored residuals (non-negative integer).
%       NOTE: 'mMax' = 0 => no acceleration. default=1000
%       2. 'itmax' = maximum allowable number of iterations. default=1000
%       3. 'atol' = absolute error tolerance. default=1.0e-6
%       4. 'rtol' = relative error tolerance. default=1.0e-3
%       5. 'droptol' = tolerance for dropping stored residual vectors to 
%       improve conditioning: If 'droptol' > 0, drop residuals if the
%       condition number exceeds droptol; if droptol <= 0,
%       do not drop residuals.
%       6. 'beta' = damping factor: If 'beta' > 0 (and beta ~= 1), then 
%       the step is damped by beta; otherwise, the step is not damped.
%       NOTE: 'beta' can be a function handle; form beta(iter), where iter 
%       is the iteration number and 0 < beta(iter) <= 1.
%       7. 'AAstart' = acceleration delay factor: If 'AAstart' > 0, start 
%       acceleration when iter = AAstart.
%
% Output:
% xS = Solution vector.
%
% The notation used is that of H.F. Walker: Anderson Acceleration:
% Algorithms and implementation
%------------------------------------------------------------------------
*)

mMax    := 1000:
itmax   := 1000:
atol    := 1.0e-8:
rtol    := 1.0e-12:
droptol := 1.0e4:
beta    := 1.0:
AAstart := 0:

# Initialize storage arrays and number of stored residuals.
DG := Matrix():
df := Matrix():
DGg := Vector(N);
QRg := Vector(N);
mAA := 0:

for iter from 0 to itmax do

   x:=X0:
   gval := Vector(phi):
   fval := gval - X0:
   res_norm := norm(fval,2):
   print(res_norm);    
   # Set the residual tolerance on the initial iteration.
   if iter = 0 then
      tol := max(atol,rtol*res_norm):
   fi:
    
   # Convergence test, if converged the loop stops.
   if res_norm <= tol then
      print(res_norm);
      break;   # Breaks for-loop
   fi:
    
   # If resnorm is larger than 1e8 at iter > 5, problem stops
   if res_norm >1e8 and iter > 5 then
      print(res_norm);
      break; # Breaks for-loop, diverged
   fi:

   # Fixed point iteration without acceleration, if mMax == 0.
   if mMax = 0 or iter < AAstart then
      # We update E <- g(E) to obtain the next approximate solution.
      for i from 1 to N do
         X0[i] := gval[i]:
      od:
   else
      # With Anderson acceleration.
      # Update the df vector and the DG array.
      if iter > AAstart then
         if mAA < mMax or Size(df,2) = 1 then
            print(whattype(df));
            df := Concatenate(2,df,fval-f_old):
            DG := Concatenate(2,DG,gval-g_old):
         else 
            df := Concatenate(2,df[..,-1],fval-f_old):
            DG := Concatenate(2,DG[..,-1],gval-g_old):   
         fi:
         mAA := mAA + 1:
      fi:   # iter
      print(iter,mAA);
      print(whattype(df));
      # We define the old g and f values for the next iteration
      f_old := fval;
      g_old := gval;
      
      if mAA = 0 then
         # Initialization
         # If mAA == 0, update X <- g(X) to obtain themah next approximate
         # solution. No least-squares problem is solved for iter = 0
         for i from 1 to N do
            X0[i] := gval[i]:
         od:
      else
         if mAA > 1 then
            Q,R := QRDecomposition(df,datatype=float);
            while ConditionNumber(R) > droptol do
                if mAA = 2 then
                   print('here1'):
                   df := convert(DeleteColumn(df,1),Vector);
                   DG := convert(DeleteColumn(DG,1),Vector);
                else
                   df := DeleteColumn(df,1);
                   DG := DeleteColumn(DG,1);
                fi:
                Q,R := QRDecomposition(df,datatype=float);
                mAA := mAA - 1;
                print(Q,R,mAA);
            od:
          
            if Size(df,2) > 1 then
               print(Q,R);
               gamma := LeastSquares([Q,R],fval);
               print(gamma);
            else
               R := norm(df,2);
               Q := MTM[mldivide](R,df);
               gamma := MTM[mldivide](R,Transpose(Q).fval);
               print(gamma);
            fi:
         else
            R := norm(df,2);
            Q := MTM[mldivide](R,df);
            gamma := MTM[mldivide](R,Transpose(Q).fval);
         fi:
 
         if Size(gamma,1) > 1 then
            DGg:=DG.gamma:
         else
            DGg:=DG*gamma;
         fi:

         # Update the approximate solution.
         for i from 1 to N do
            X0[i] := gval[i] - DGg[i];
         od:
         print('Sol',X0);
         
         # Damping for non-zero beta
         if beta > 0 and beta <> 1 then
            if mAA = 1 then
               QRg := Q*R*gamma;
            else
               QRg := df.gamma;
            fi:
            for i from 1 to N do
               X0[i] := X0[i] - (1-beta)*(fval[i] - QRg[i]);
            od:
         fi:# isa(beta ...
         print(iter,mAA);
      fi: # mAA = 0
   fi:# mMax == 0

od:
xS := Vector(N);
for i from 1 to N do xS[i]:=X0[i]: od:
return xS
end:

AndersonAcceleration(2,phix,X):