How to solve set of nonlinear fixed-point equations using Anderson acceleration in NLsolve.jl?

I have a set of two nonlinear fixed-point equations like

x = g1(x, y)
y = g2(x, y)

where x and y are n-dimensional arrays (n=1,2,3), length(x) and length(y) can be very large, e.g. 128x128x128 for 3D problem. How can I solve it using anderson algorithm in the NLsolve.jl package which seems to be designed for x = g(x) only?

Note that your problem has precisely the form z = g(z) if you were to consider z := (x, y) and g := (g1, g2).

The only thing you have to do is to concatenate both variables and functions. Assuming x and y have the same size and dimensions, one way to do this would be
z = cat(x, y; dims=ndims(x)+1) # concatenate along a new dimension
And then a similar procedure for your functions.

But personally I am not a fan of this because explicit concatenation with these cat functions is not very efficient, especially for big arrays like yours. It also only works properly if x and y have the same size and typeof. On top of that, you’d need to throw @views everywhere to avoid further allocations.

So something much nicer would be to use an array of arrays such as
z = [x, y]
g(z) = [g1(z[1], z[2]), g2(z[1], z[2])]

But nlsolve can’t really operate on such containers.

However, it can work on slightly more refined containers such as the ones in RecursiveArraysTools, which are designed specifically for these tasks.

using RecursiveArrayTools
z = VectorOfArray([x, y])
g(z) = VectorOfArray([g1(z[1], z[2]), g2(z[1], z[2])])
nlsolve(g, z; method=:anderson)

Great! I will try them out and mark this post as the solution if everything goes well. Thanks!

By studying the source code of NLsolve, I am confused on which function should I serve the nlsolve.
It seems that nlsolve expect f(x) = 0 to solve. Thus, should I do

f(z) = g(z) - z
nlsolve(f, z; method=:anderson)


Another issue I have encountered is the way NLsolve defines its AndersonCache object. The line 24 in NLsolve.jl/src/solvers/anderson.jl is

mmax = min(length(x), m)

which is too restrictive for VectorOfArray objects. Here x is a VectorOfArray object, thus length(x) is the number of AbstractArrays in the object (it is 2 in this example). However, we expect mmax can be much larger because the length of each AbstractArray in the object can be very large, e.g. 128x128x128 for 3D problems.

At present, I wrote a constructor for AndersonCache by myself to allow any m.

A fixed-point equation f(z) = z can always be mapped to g(z) = 0 by precisely that recipe!
Note that NLsolve.jl has a wrapper for this, where instead of calling nlsolve you call fixedpoint and this is taken care for you
But personally I would keep using the standard nlsolve function call: being more explicit helps in the long run.

Also note that m is NOT related with the size of your solution z (otherwise that min operation could have you end up with an underdetermined system!). Here you can find explained what m is for Anderson acceleration - Wikipedia

So if my fixed point is x=g(x), I should input f(x) = g(x) - x to nlsolve, right?

I know exactly what m does in Anderson mixing. The issue here is the Base.length defined by VectorOfArray only returns the number of vectors it contains. Thus If I use the NLsove provided version of AndersonCache constructor, m is restricted to be less than 2 (it always choose the minmum of m given and length(VectorOfArray)). I have written my own version of Anderson mixing long ago in C++ and I can use m as large as 100 to accelerate the convergence.

I marked this post as a solution as it provides a very good approach to represent a set of fixed-point equations involving many high-dimensional data. However, this approach cannot directly apply to NLsolve.anderson which deals with VectorOfArray objects in a way that it only allows <= length(VectorOfArray) rather than the length of multi-dimensional data.

I am ending up rewriting my own version of Anderson mixing which closely follows Walker’s pseudocode. (ie. its way to treat β relaxation). As Q*R*γ results in a flatten vector, we have to figure out a way to copy this long vector back into VectorOfArray{MyType} object. At present, I overload Base.copyto! to accomplish the task. But it is not universal.

Thanks for your input @FMeirinhos !

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

function anderson_(df::Union{NonDifferentiable, OnceDifferentiable},
    @. 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

        # 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])
        res /= length(initial_x)

        if iter % 100 == 0
            @show iter, res

        # check convergence
        converged = (res < ftol) ? true : false
        if any(isnan, cache.x) || any(isnan, fx)
        # converged = x_converged || f_converged
        if converged
            @show iter, res

        # 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 -
                # circularly shift differences of g
                Δg1 = cache.Δgs[1]
                for i in 1:(m-1)
                    cache.Δgs[i] = cache.Δgs[i+1]
                cache.Δgs[m] = Δg1
                @. cache.Δgs[m] = cache.g -

        @. cache.fxold = fx
        @. = cache.g

        if m_eff == 0 || m == 0
            if 0 < β ≤ 1.0
                @. cache.x += β * fx
                @. cache.x = cache.g
            if m_eff == 1
                NLsolve.qradd!(cache.Q, cache.R, vec(Δf), m_eff)
                if m_eff > m
                    # delete left-most column of QR decomposition
                    NLsolve.qrdelete!(cache.Q, cache.R, m)
                    m_eff -= 1
                NLsolve.qradd!(cache.Q, cache.R, vec(Δf), m_eff)

            # 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))

            # 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]
            @. cache.x = cache.g
            if 0 < β < 1
                copyto!(Δf, Q * R * γs)
                @. cache.x -= (1-β) * (fx - Δf)
        iter += 1

    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)