anderson.jl

module Anderson

export anderson

function anderson(g,x0,m::Int,niter::Int,eps::Real,warming=0)
    @assert length(size(x0)) == 1 #1D array input
    N = size(x0,1)
    T = eltype(x0)
    #xs: ring buffer storing the iterates, from newest to oldest
    xs = zeros(T,N,m+1)
    gs = zeros(T,N,m+1)
    xs[:,1] = x0
    errs = zeros(niter)
    for n = 1:niter
        gs[:,1] = g(xs[:,1])
        err = norm(gs[:,1] - xs[:,1],Inf)
        errs[n] = err
        # println("$n $err")
        if(err < eps)
            break
        end
        m_eff = min(n-1,m)
        new_x = gs[:,1]
        if m_eff > 0 && n > warming
            mat = (gs[:,2:m_eff+1] - xs[:,2:m_eff+1]) .- (gs[:,1] - xs[:,1])
            alphas = -mat \ (gs[:,1] - xs[:,1])
            # alphas = -(mat'*mat) \ (mat'* (gs[:,1] - xs[:,1]))
            for i = 1:m_eff
                new_x .+= alphas[i].*(gs[:,i+1] .- gs[:,1])
            end
        end

        xs = circshift(xs,(0,1))
        gs = circshift(gs,(0,1))
        xs[:,1] = new_x
    end
    gs[:,1],errs
end
end


module Anderson_test
using Anderson
function test_anderson()
    N = 100
    A = randn(N,N)
    A = A'*A
    b = randn(N)
    x0 = randn(N)
    g(x) = x - .00001(A*x - b)
    x = Anderson.anderson(g,x0,200,200,1e-10)
    println(norm(A*x-b)/norm(A*x0-b))
end
end