And here is my implementation (the intial_x has to be `AbstractVectorOfArray`

):

```
function anderson_(df::Union{NonDifferentiable, OnceDifferentiable},
initial_x::AbstractArray,
xtol::Real,
ftol::Real,
iterations::Integer,
store_trace::Bool,
show_trace::Bool,
extended_trace::Bool,
β::Real,
aa_start::Integer,
droptol::Real,
cache::NLsolve.AndersonCache)
@. cache.x = initial_x
tr = NLsolve.SolverTrace()
tracing = store_trace || show_trace || extended_trace
aa_iteration = cache.γs !== nothing
x_converged, f_converged, converged = false, false, false
m = aa_iteration ? length(cache.γs) : 0
aa_iteration && (m_eff = 0)
# Preallocate Δf = f_{k+1}(x) - f_{k}(x)
Δf = deepcopy(cache.fxold)
iter = 0
while iter < iterations
# evaluate function: f(x) = g(x) - x
value!!(df, cache.x)
# Give a name to value(df). It will not allocate.
fx = value(df)
# check that all values are finite
NLsolve.check_isfinite(fx)
# Compute g(x) = x + f(x)
@. cache.g = cache.x + fx
res = zero(eltype(initial_x))
for i in 1:length(initial_x)
res += vecnorm1(fx[i]) / vecnorm1(cache.x[i])
end
res /= length(initial_x)
if iter % 100 == 0
@show iter, res
end
# check convergence
converged = (res < ftol) ? true : false
if any(isnan, cache.x) || any(isnan, fx)
break
end
# converged = x_converged || f_converged
if converged
@show iter, res
break
end
# perform Anderson acceleration
if iter > aa_start && m > 0
# increase size of history
m_eff += 1
@. Δf = fx - cache.fxold
if m_eff < (m+1)
@. cache.Δgs[m_eff] = cache.g - cache.gold
else
# circularly shift differences of g
Δg1 = cache.Δgs[1]
for i in 1:(m-1)
cache.Δgs[i] = cache.Δgs[i+1]
end
cache.Δgs[m] = Δg1
@. cache.Δgs[m] = cache.g - cache.gold
end
end
@. cache.fxold = fx
@. cache.gold = cache.g
if m_eff == 0 || m == 0
if 0 < β ≤ 1.0
@. cache.x += β * fx
else
@. cache.x = cache.g
end
else
if m_eff == 1
NLsolve.qradd!(cache.Q, cache.R, vec(Δf), m_eff)
else
if m_eff > m
# delete left-most column of QR decomposition
NLsolve.qrdelete!(cache.Q, cache.R, m)
m_eff -= 1
end
NLsolve.qradd!(cache.Q, cache.R, vec(Δf), m_eff)
end
# define current Q and R matrices
Q, R = view(cache.Q, :, 1:m_eff), UpperTriangular(view(cache.R, 1:m_eff, 1:m_eff))
# check condition
if droptol > 0
while cond(R) > droptol && m_eff > 1
NLsolve.qrdelete!(cache.Q, cache.R, m_eff)
m_eff -= 1
Q, R = view(cache.Q, :, 1:m_eff), UpperTriangular(view(cache.R, 1:m_eff, 1:m_eff))
end
end
# solve least squares problem
γs = view(cache.γs, 1:m_eff)
ldiv!(R, mul!(γs, Q', vec(fx)))
# update next iterate
for i in 1:m_eff
@. cache.g -= cache.γs[i] * cache.Δgs[i]
end
@. cache.x = cache.g
if 0 < β < 1
copyto!(Δf, Q * R * γs)
@. cache.x -= (1-β) * (fx - Δf)
end
end
iter += 1
end
return NLsolve.SolverResults("Anderson m=$m β=$β aa_start=$aa_start droptol=$droptol",
initial_x, deepcopy(cache.x), norm(value(df), Inf),
iter, x_converged, xtol, f_converged, ftol, tr,
first(df.f_calls), 0)
end
```