anderson-helpers.hpp

Go to the documentation of this file.

#pragma once
 
#include <panoc-alm/inner/detail/limited-memory-qr.hpp>
 
namespace pa {
 
inline void minimize_update_anderson(LimitedMemoryQR &qr, rmat G, crvec rₖ,
                                     crvec rₖ₋₁, crvec gₖ, rvec γ_LS,
                                     rvec xₖ_aa) {
    // Update QR factorization for Anderson acceleration
    if (qr.num_columns() == qr.m()) // if the history buffer is full
        qr.remove_column();
    qr.add_column(rₖ - rₖ₋₁);
 
    // Solve least squares problem Anderson acceleration
    // γ = argmin ‖ ΔR γ - rₖ ‖²
    qr.solve_col(rₖ, γ_LS);
 
    // Iterate over columns of G, whose indices match the indices of the matrix
    // R in the QR factorization, stored as a circular buffer.
    auto g_it  = qr.ring_iter().begin();
    auto g_end = qr.ring_iter().end();
    assert(g_it != g_end);
 
    // Compute Anderson acceleration next iterate yₑₓₜ = ∑ₙ₌₀ αₙ gₙ
    // α₀ = γ₀             if n = 0
    // αₙ = γₙ - γₙ₋₁      if 0 < n < mₖ
    // αₘ = 1 - γₘ₋₁       if n = mₖ
    real_t α = γ_LS(0);
    xₖ_aa    = α * G.col((*g_it).circular);
    while (++g_it != g_end) {
        auto [i, g_idx] = *g_it; // [zero based index, circular index]
        α               = γ_LS(i) - γ_LS(i - 1);
        xₖ_aa += α * G.col(g_idx);
    }
    α = 1 - γ_LS(qr.num_columns() - 1);
    xₖ_aa += α * gₖ;
 
    // Add the new column to G
    G.col(qr.ring_tail()) = gₖ; // TODO: avoid copy, make G an array of vectors
}
 
} // namespace pa