NumLib::AndersonAcceleration Class Reference

Detailed Description

Anderson acceleration of the damped Picard fixpoint iteration.

Owns a sliding window of the last depth damped steps and the Gram matrix of that window, and mixes them into an accelerated iterate. History vectors are taken from the global vector provider on demand and returned in the destructor (RAII), so no manual release is required at the call site.

The window and its Gram matrix are only ever mutated together; keeping them in one object makes "history and Gram stay in sync" a class invariant instead of a convention spread across the solver loop.

A depth below min_mixing_depth admits no mixing (the sum-to-one constraint forces unit weight on the single stored step), so such an instance is an inert no-op: it takes no vector from the global vector provider and leaves the iterate untouched. This is the plain-Picard case.

Definition at line 30 of file AndersonAcceleration.h.

#include <AndersonAcceleration.h>

Classes
struct	HistoryEntry

Public Member Functions
	AndersonAcceleration (int depth)
	~AndersonAcceleration ()
	AndersonAcceleration (AndersonAcceleration const &)=delete
AndersonAcceleration &	operator= (AndersonAcceleration const &)=delete
void	accelerate (GlobalVector const &x_old, GlobalVector &x_new)
void	dropLastStep ()

Static Public Attributes
static constexpr int	min_mixing_depth = 2

Static Private Member Functions
static void	releaseHistoryEntry (HistoryEntry const &entry)
	Returns entry's vectors to the global vector provider.

Private Attributes
int const	_depth
std::vector< HistoryEntry >	_history
	Circular buffer of history entries, oldest first (size <= _depth).
Eigen::MatrixXd	_gram

Constructor & Destructor Documentation

AndersonAcceleration() [1/2]

NumLib::AndersonAcceleration::AndersonAcceleration ( int depth )

explicit

Parameters

depth

number of previous iterates retained for mixing; a value below min_mixing_depth makes the instance an inert no-op (plain Picard).

Definition at line 150 of file AndersonAcceleration.cpp.

    : _depth(depth), _gram(depth, depth)
{
    if (_depth >= min_mixing_depth)
    {
        _history.reserve(_depth);
    }
}

References _depth, _gram, _history, and min_mixing_depth.

Referenced by AndersonAcceleration(), and operator=().

~AndersonAcceleration()

NumLib::AndersonAcceleration::~AndersonAcceleration

(

)

Definition at line 159 of file AndersonAcceleration.cpp.

{
    for (auto const& entry : _history)
    {
        releaseHistoryEntry(entry);
    }
}

References _history, and releaseHistoryEntry().

AndersonAcceleration() [2/2]

NumLib::AndersonAcceleration::AndersonAcceleration ( AndersonAcceleration const & )

delete

References AndersonAcceleration().

Member Function Documentation

accelerate()

void NumLib::AndersonAcceleration::accelerate	(	GlobalVector const &	x_old,
		GlobalVector &	x_new )

Records the newest damped step \( x_{\rm old} \to x_{\rm new} \) and overwrites x_new in place with the Anderson-mixed iterate.

While fewer than two steps are stored, or when the mixture is rejected as untrustworthy (see detail::computeAndersonWeights), x_new is left unchanged. For a no-op instance (depth < min_mixing_depth) this does nothing.

Parameters

x_old	the iterate entering this step.
x_new	in: the damped Picard step \( f = \beta(g(x_{\rm old}) - x_{\rm old}) \) added to x_old; out: the mixed iterate.

Definition at line 173 of file AndersonAcceleration.cpp.

{
    namespace LinAlg = MathLib::LinAlg;
 
    // A depth below min_mixing_depth admits no mixing (plain Picard); nothing
    // is stored.
    if (_depth < min_mixing_depth)
    {
        return;
    }
 
    // Additionally mixes the last _depth damped steps
    //   f_i = beta*(g(x_i) - x_i) (i.e. x_new - x_old computed after the beta
    //   relaxation already applied by the caller) to find the optimal theta
    //   minimising ||sum theta_i f_i|| s.t. sum theta_i = 1, then sets
    //   x_new = sum theta_i*(x_i + f_i).
 
    // Whether the circular buffer is full and the oldest entry is about to be
    // evicted (needed for the incremental Gram update).
    bool const rotated = static_cast<int>(_history.size()) == _depth;
    if (!rotated)
    {
        // The id out-params are unused: the provider allocates a fresh vector
        // on every call and never re-fetches by id.
        std::size_t x_id = 0u;
        std::size_t f_id = 0u;
        _history.push_back(
            {&NumLib::GlobalVectorProvider::provider.getVector(x_id),
             &NumLib::GlobalVectorProvider::provider.getVector(f_id)});
    }
    else
    {
        // Recycle the oldest entry as the newest one.
        std::rotate(_history.begin(), _history.begin() + 1, _history.end());
    }
 
    auto const& newest = _history.back();
 
    // x = x_old, f = x_new - x_old
    LinAlg::copy(x_old, *newest.x);
    LinAlg::copy(x_new, *newest.f);
    LinAlg::axpy(*newest.f, -1.0, x_old);
 
    // Actual window size, <= _depth while the buffer fills.
    int const history_size = static_cast<int>(_history.size());
 
    // Incrementally maintain the (history_size x history_size) Gram matrix
    // G = F^T F whose columns are the stored damped steps f_0 ...
    // f_{history_size-1}. All steps but the newest are unchanged from the
    // previous iteration, so only the last row/column is recomputed -
    // history_size dot products instead of a full
    // history_size*(history_size+1)/2 rebuild. On a rotate the oldest entry
    // (index 0) was evicted, so the cached block is first shifted up-left by
    // one.
    if (rotated)
    {
        _gram.topLeftCorner(history_size - 1, history_size - 1) =
            _gram.block(1, 1, history_size - 1, history_size - 1).eval();
    }
    for (int i = 0; i < history_size; ++i)
    {
        double const d = LinAlg::dot(*_history[i].f, *newest.f);
        _gram(i, history_size - 1) = d;
        _gram(history_size - 1, i) = d;
    }
 
    // A single stored step needs no mixing: the sum-to-one constraint forces
    // theta = (1), which just reproduces the damped step already held in x_new.
    if (history_size < min_mixing_depth)
    {
        return;
    }
 
    // Solve G theta = e (least-squares) with the constraint sum theta_i = 1 via
    // a simple Lagrange formulation:
    //
    //   [ G  1 ] [ theta  ] = [ 0 ]
    //   [ 1  0 ] [ lambda ]   [ 1 ]
    //
    // The beta factor scales G by beta^2 and cancels in theta, so the weights
    // are identical to the undamped case. The Anderson update is then:
    //   x_anderson = sum_i theta_i * (x_i + f_i)
    //              = sum_i theta_i * (x_i + beta*(g(x_i)-x_i))
    //              = sum_i theta_i * ((1-beta)*x_i + beta*g(x_i))
    Eigen::MatrixXd const G = _gram.topLeftCorner(history_size, history_size);
 
    Eigen::VectorXd const theta = detail::computeAndersonWeights(G);
 
    // Accumulate the Anderson mixed iterate directly into x_new. Its previous
    // value is no longer needed: the newest step was already extracted from it
    // above, and it is not aliased by any history entry (those are independent
    // copies).
    x_new.setZero();
    for (int i = 0; i < history_size; ++i)
    {
        // x_new += theta_i * (x_i + f_i)
        LinAlg::axpy(x_new, theta(i), *_history[i].x);
        LinAlg::axpy(x_new, theta(i), *_history[i].f);
    }
 
    DBUG("Picard/Anderson: history size {:d}, theta=[{:.4g}]", history_size,
         fmt::join(theta.data(), theta.data() + history_size, ", "));
}

References _depth, _gram, _history, MathLib::LinAlg::axpy(), NumLib::detail::computeAndersonWeights(), MathLib::LinAlg::copy(), DBUG(), MathLib::LinAlg::dot(), min_mixing_depth, NumLib::GlobalVectorProvider::provider, and MathLib::EigenVector::setZero().

Referenced by NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::solve().

dropLastStep()

void NumLib::AndersonAcceleration::dropLastStep

(

)

Discards the step recorded by the most recent accelerate() call, used when the current iteration is repeated. The Gram shift performed while recording is deliberately not undone (see implementation).

Definition at line 278 of file AndersonAcceleration.cpp.

{
    // Drop the just-added (now stale) history entry, since the iteration is
    // being repeated. In the full-buffer case this is the recycled slot;
    // releasing it by reference is safe because the provider tracks vectors by
    // pointer, not by the (unused) id.
    //
    // The rotation and the Gram shift performed by accelerate() are not undone,
    // and need not be: dropping the newest entry leaves the buffer holding the
    // remaining entries in order, and the shifted top-left block of the Gram
    // matrix is exactly their Gram matrix. The oldest entry stays evicted,
    // which merely shortens the sliding window by one.
    if (!_history.empty())
    {
        releaseHistoryEntry(_history.back());
        _history.pop_back();
    }
}

References _history, and releaseHistoryEntry().

Referenced by NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::solve().

operator=()

AndersonAcceleration & NumLib::AndersonAcceleration::operator= ( AndersonAcceleration const & )

delete

References AndersonAcceleration().

releaseHistoryEntry()

void NumLib::AndersonAcceleration::releaseHistoryEntry ( HistoryEntry const & entry )

staticprivate

Returns entry's vectors to the global vector provider.

Definition at line 167 of file AndersonAcceleration.cpp.

{
    NumLib::GlobalVectorProvider::provider.releaseVector(*entry.x);
    NumLib::GlobalVectorProvider::provider.releaseVector(*entry.f);
}

References NumLib::AndersonAcceleration::HistoryEntry::f, NumLib::GlobalVectorProvider::provider, and NumLib::AndersonAcceleration::HistoryEntry::x.

Referenced by ~AndersonAcceleration(), and dropLastStep().

Member Data Documentation

_depth

int const NumLib::AndersonAcceleration::_depth

private

Maximum window size; mixing is active only for values >= min_mixing_depth.

Definition at line 81 of file AndersonAcceleration.h.

Referenced by AndersonAcceleration(), and accelerate().

_gram

Eigen::MatrixXd NumLib::AndersonAcceleration::_gram

private

Gram matrix G = F^T F of the stored steps, maintained incrementally across iterations (only the newest step's row/column is recomputed). Sized once to the maximum window (_depth x _depth) so the incremental update never reallocates; only the leading history_size x history_size block is live while the window fills.

Definition at line 91 of file AndersonAcceleration.h.

Referenced by AndersonAcceleration(), and accelerate().

_history

std::vector<HistoryEntry> NumLib::AndersonAcceleration::_history

private

Circular buffer of history entries, oldest first (size <= _depth).

Definition at line 84 of file AndersonAcceleration.h.

Referenced by AndersonAcceleration(), ~AndersonAcceleration(), accelerate(), and dropLastStep().

min_mixing_depth

int NumLib::AndersonAcceleration::min_mixing_depth = 2

staticconstexpr

Smallest depth that admits mixing: below it the sum-to-one constraint forces unit weight on the single stored step, i.e. plain Picard.

Definition at line 35 of file AndersonAcceleration.h.

Referenced by AndersonAcceleration(), accelerate(), NumLib::createNonlinearSolver(), and NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::solve().

---

The documentation for this class was generated from the following files:

• AndersonAcceleration.h
• AndersonAcceleration.cpp