Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include <panoc-alm/inner/detail/limited-memory-qr.hpp>
       4             : 
       5             : namespace pa {
       6             : 
       7             : /**
       8             :  * @brief   Solve one step of Anderson acceleration to find a fixed point of a 
       9             :  *          function g(x):
      10             :  * 
      11             :  * @f$ g(x^\star) - x^\star = 0 @f$
      12             :  * 
      13             :  * Updates the QR factorization of @f$ \mathcal{R}_k = QR @f$, solves the least 
      14             :  * squares problem to find @f$ \gamma_\text{LS} @f$, computes the next 
      15             :  * iterate @f$ x_{k+1} @f$, and stores the current function value @f$ g_k @f$
      16             :  * in the matrix @f$ G @f$, which is used as a circular buffer.
      17             :  * @f[ \begin{aligned}
      18             :  * \mathcal{R}_k &= \begin{pmatrix} \Delta r_{k-m_k} & \dots & \Delta r_{k-1} \end{pmatrix} \in \mathbb{R}^{n\times m_k} \\
      19             :  * \Delta r_i &= r_{i+1} - r_i \\
      20             :  * r_i &= g_i - x_i \\
      21             :  * g_i &= g(x_i) \\
      22             :  * \DeclareMathOperator*{\argmin}{argmin}
      23             :  * \gamma_\text{LS} &= \argmin_\gamma \left\| \mathcal{R}_k \gamma - r_k \right\|_2 \\
      24             :  * \alpha_i &= \begin{cases} \gamma_\text{LS}[0] & i = 0 \\
      25             :  *                           \gamma_\text{LS}[i] - \gamma_\text{LS}[i-1] & 0 < i < m_k \\
      26             :  *                           1 - \gamma_\text{LS}[m_k - 1] & i = m_k \end{cases} \\
      27             :  * x_{k+1} &= \sum_{i=0}^{m_k} \alpha_i\,g_i
      28             :  * \end{aligned} @f]
      29             :  */
      30           4 : inline void minimize_update_anderson(LimitedMemoryQR &qr, rmat G, crvec rₖ,
      31             :                                      crvec rₖ₋₁, crvec gₖ, rvec γ_LS,
      32             :                                      rvec xₖ_aa) {
      33             :     // Update QR factorization for Anderson acceleration
      34           4 :     if (qr.num_columns() == qr.m()) // if the history buffer is full
      35           3 :         qr.remove_column();
      36           4 :     qr.add_column(rₖ - rₖ₋₁);
      37             : 
      38             :     // Solve least squares problem Anderson acceleration
      39             :     // γ = argmin ‖ ΔR γ - rₖ ‖²
      40           4 :     qr.solve_col(rₖ, γ_LS);
      41             : 
      42             :     // Iterate over columns of G, whose indices match the indices of the matrix
      43             :     // R in the QR factorization, stored as a circular buffer.
      44           4 :     auto g_it  = qr.ring_iter().begin();
      45           4 :     auto g_end = qr.ring_iter().end();
      46           4 :     assert(g_it != g_end);
      47             : 
      48             :     // Compute Anderson acceleration next iterate yₑₓₜ = ∑ₙ₌₀ αₙ gₙ
      49             :     // α₀ = γ₀             if n = 0
      50             :     // αₙ = γₙ - γₙ₋₁      if 0 < n < mₖ
      51             :     // αₘ = 1 - γₘ₋₁       if n = mₖ
      52           4 :     real_t α = γ_LS(0);
      53           4 :     xₖ_aa    = α * G.col((*g_it).circular);
      54          12 :     while (++g_it != g_end) {
      55          32 :         auto [i, g_idx] = *g_it; // [zero based index, circular index]
      56           8 :         α               = γ_LS(i) - γ_LS(i - 1);
      57           8 :         xₖ_aa += α * G.col(g_idx);
      58           8 :     }
      59           4 :     α = 1 - γ_LS(qr.num_columns() - 1);
      60           4 :     xₖ_aa += α * gₖ;
      61             : 
      62             :     // Add the new column to G
      63           4 :     G.col(qr.ring_tail()) = gₖ; // TODO: avoid copy, make G an array of vectors
      64           4 : }
      65             : 
      66             : } // namespace pa