alpaqa 1.0.0a10
Nonconvex constrained optimization
Loading...
Searching...
No Matches
anderson-helpers.hpp
Go to the documentation of this file.
1#pragma once
2
5#include <alpaqa/export.hpp>
6
7namespace alpaqa {
8
9/**
10 * @brief Solve one step of Anderson acceleration to find a fixed point of a
11 * function g(x):
12 *
13 * @f$ g(x^\star) - x^\star = 0 @f$
14 *
15 * Updates the QR factorization of @f$ \mathcal{R}_k = QR @f$, solves the least
16 * squares problem to find @f$ \gamma_\text{LS} @f$, computes the next
17 * iterate @f$ x_{k+1} @f$, and stores the current function value @f$ g_k @f$
18 * in the matrix @f$ G @f$, which is used as a circular buffer.
19 * @f[ \begin{aligned}
20 * \def\gammaLS{\gamma_\text{LS}}
21 * m_k &= \min \{k, m\} \\
22 * g_i &= g(x_i) \\
23 * r_i &= r(x_i) g_i - x_i \\
24 * \Delta r_i &= r_i - r_{i-1} \\
25 * \mathcal{R}_k &= \begin{pmatrix} \Delta r_{k-m_k+1} & \dots & \Delta r_k \end{pmatrix} \in \R^{n\times m_k}
26 * \\
27 * \gammaLS &= \argmin_{\gamma \in \R^{m_k}}\; \norm{\mathcal{R}_k \gamma - r_k}_2 \\
28 * \alpha_i &= \begin{cases} \gammaLS[0] & i = 0 \\
29 * \gammaLS[i] - \gammaLS[i-1] & 0 \lt i \lt m_k \\
30 * 1 - \gammaLS[m_k - 1] & i = m_k \end{cases} \\
31 * \tilde G_k &= \begin{pmatrix} g_{k - m_k} & \dots & g_{k-1} \end{pmatrix} \\
32 * G_k &= \begin{pmatrix} g_{k - m_k} & \dots & g_{k} \end{pmatrix} \\
33 * &= \begin{pmatrix} \tilde G_k & g_{k} \end{pmatrix} \\
34 * x_{k+1} &= \sum_{i=0}^{m_k} \alpha_i\,g_{k - m_k + i} \\
35 * &= G_k \alpha \\
36 * \end{aligned} @f]
37 */
38template <Config Conf = DefaultConfig>
40 /// [inout] QR factorization of @f$ \mathcal{R}_k @f$
42 /// [inout] Matrix of previous function values @f$ \tilde G_k @f$
43 /// (stored as ring buffer with the same indices as `qr`)
44 rmat<Conf> G̃,
45 /// [in] Current residual @f$ r_k @f$
46 crvec<Conf> rₖ,
47 /// [in] Previous residual @f$ r_{k-1} @f$
48 crvec<Conf> rₗₐₛₜ,
49 /// [in] Current function value @f$ g_k @f$
50 crvec<Conf> gₖ,
51 /// [in] Minimum divisor when solving close to singular systems,
52 /// scaled by the maximum eigenvalue of R
53 real_t<Conf> min_div_fac,
54 /// [out] Solution to the least squares system
55 rvec<Conf> γ_LS,
56 /// [out] Next Anderson iterate
57 rvec<Conf> xₖ_aa) {
58
59 // Update QR factorization for Anderson acceleration
60 if (qr.num_columns() == qr.m()) // if the history buffer is full
61 qr.remove_column();
62 qr.add_column(rₖ - rₗₐₛₜ);
63
64 // Solve least squares problem Anderson acceleration
65 // γ = argmin ‖ ΔR γ - rₖ ‖²
66 qr.solve_col(rₖ, γ_LS, qr.get_max_eig() * min_div_fac);
67
68 // Iterate over columns of G, whose indices match the indices of the matrix
69 // R in the QR factorization, stored as a circular buffer.
70 auto g_it = qr.ring_iter().begin();
71 auto g_end = qr.ring_iter().end();
72 assert(g_it != g_end);
73
74 // Compute Anderson acceleration next iterate yₑₓₜ = ∑ₙ₌₀ αₙ gₙ
75 // α₀ = γ₀ if n = 0
76 // αₙ = γₙ - γₙ₋₁ if 0 < n < mₖ
77 // αₘ = 1 - γₘ₋₁ if n = mₖ
78 auto α = γ_LS(0);
79 xₖ_aa = α * G̃.col((*g_it).circular);
80 while (++g_it != g_end) {
81 auto [i, g_idx] = *g_it; // [zero based index, circular index]
82 α = γ_LS(i) - γ_LS(i - 1);
83 xₖ_aa += α * G̃.col(g_idx);
84 }
85 α = 1 - γ_LS(qr.num_columns() - 1);
86 xₖ_aa += α * gₖ;
87
88 // Add the new column to G
89 G̃.col(qr.ring_tail()) = gₖ; // TODO: avoid copy, make G an array of vectors
90}
91
92} // namespace alpaqa
Incremental QR factorization using modified Gram-Schmidt with reorthogonalization.
CircularRange< index_t > ring_iter() const
Get iterators in the circular buffer.
length_t num_columns() const
Get the number of columns that are currently stored.
void remove_column()
Remove the leftmost column.
void add_column(const VecV &v)
Add the given column to the right.
void solve_col(const VecB &b, VecX &x, real_t tol=0) const
Solve the least squares problem Ax = b.
index_t ring_tail() const
Get the tail index of the circular buffer (points to one past the most recent element).
real_t get_max_eig() const
Get the maximum eigenvalue of R.
typename Conf::rmat rmat
Definition: config.hpp:72
typename Conf::real_t real_t
Definition: config.hpp:63
typename Conf::rvec rvec
Definition: config.hpp:67
typename Conf::crvec crvec
Definition: config.hpp:68
void minimize_update_anderson(LimitedMemoryQR< Conf > &qr, rmat< Conf > G̃, crvec< Conf > rₖ, crvec< Conf > rₗₐₛₜ, crvec< Conf > gₖ, real_t< Conf > min_div_fac, rvec< Conf > γ_LS, rvec< Conf > xₖ_aa)
Solve one step of Anderson acceleration to find a fixed point of a function g(x):