C++ API Reference¶
-
namespace quala¶
Typedefs
-
using real_t = double¶
Default floating point type.
-
using realmat = Eigen::Matrix<real_t, Eigen::Dynamic, Eigen::Dynamic>¶
Default type for floating point matrices.
-
using index_t = Eigen::Index¶
Default type for vector indices.
Functions
-
inline void minimize_update_anderson(LimitedMemoryQR &qr, rmat G, crvec rₖ, crvec rₖ₋₁, crvec gₖ, rvec γ_LS, rvec xₖ_aa)¶
Solve one step of Anderson acceleration to find a fixed point of a function g(x):
\( g(x^\star) - x^\star = 0 \)
Updates the QR factorization of \( \mathcal{R}_k = QR \), solves the least squares problem to find \( \gamma_\text{LS} \), computes the next iterate \( x_{k+1} \), and stores the current function value \( g_k \) in the matrix \( G \), which is used as a circular buffer.
\[\begin{split} \begin{aligned} \mathcal{R}_k &= \begin{pmatrix} \Delta r_{k-m_k} & \dots & \Delta r_{k-1} \end{pmatrix} \in \mathbb{R}^{n\times m_k} \\ \Delta r_i &= r_{i+1} - r_i \\ r_i &= g_i - x_i \\ g_i &= g(x_i) \\ \DeclareMathOperator*{\argmin}{argmin} \gamma_\text{LS} &= \argmin_\gamma \left\| \mathcal{R}_k \gamma - r_k \right\|_2 \\ \alpha_i &= \begin{cases} \gamma_\text{LS}[0] & i = 0 \\ \gamma_\text{LS}[i] - \gamma_\text{LS}[i-1] & 0 < i < m_k \\ 1 - \gamma_\text{LS}[m_k - 1] & i = m_k \end{cases} \\ x_{k+1} &= \sum_{i=0}^{m_k} \alpha_i\,g_i \end{aligned} \end{split}\]
Variables
-
class AndersonAccel¶
- #include <quala/anderson-acceleration.hpp>
Anderson Acceleration.
- Todo:
Condition estimation of the QR factorization.
Public Types
-
using Params = AndersonAccelParams¶
Public Functions
-
inline AndersonAccel(Params params, length_t n)¶
- Parameters
params – Parameters.
n – Problem dimension (size of the vectors).
-
inline void resize(length_t n)¶
Change the problem dimension.
Flushes the history.
- Parameters
n – Problem dimension (size of the vectors).
-
inline void initialize(crvec g₀, vec r₀)¶
Call this function on the first iteration to initialize the accelerator.
-
inline void compute(crvec gₖ, crvec rₖ, rvec xₖ_aa)¶
Compute the accelerated iterate \( x^k_\text{AA} \), given the function value at the current iterate \( g^k = g(x^k) \) and the corresponding residual \( r^k = g^k - x^k \).
-
inline void compute(crvec gₖ, vec &&rₖ, rvec xₖ_aa)¶
Compute the accelerated iterate \( x^k_\text{AA} \), given the function value at the current iterate \( g^k = g(x^k) \) and the corresponding residual \( r^k = g^k - x^k \).
-
inline void reset()¶
Reset the accelerator (but keep the last function value and residual, so calling initialize is not necessary).
-
struct AndersonAccelParams¶
- #include <quala/anderson-acceleration.hpp>
Parameters for the AndersonAccel class.
-
class BroydenGood¶
- #include <quala/broyden-good.hpp>
Broyden’s “Good” method for solving systems of nonlinear equations \( F(x) = 0 \).
Where \( B_k \) approximates the Jacobian of \( F(x_k) \) and \( H_k \triangleq B_k^{-1} \).\[\begin{split} \begin{aligned} B_{k+1} &= B_k + \frac{\left(y_k - B_k s_k\right) s_k^\top} {s_k^\top s_k} \\ H_{k+1} &= H_k + \frac{\left(s_k - H_k y_k\right) \left(s_k^\top H_k\right)} {s_k^\top H_k y_k} \\ s_k &\triangleq x_{k+1} - x_k \\ y_k &\triangleq F(x_{k+1}) - F(x_k) \\ \end{aligned} \end{split}\]- Todo:
Damping.
Public Types
-
using Params = BroydenGoodParams¶
Public Functions
-
template<class VecS, class VecY>
bool update_sy(const anymat<VecS> &s, const anymat<VecY> &y, bool forced = false)¶ Update the inverse Jacobian approximation using the new vectors sₖ = xₖ₊₁ - xₖ and yₖ = pₖ₊₁ - pₖ.
-
inline bool update(crvec xₖ, crvec xₖ₊₁, crvec pₖ, crvec pₖ₊₁, bool forced = false)¶
Update the inverse Jacobian approximation using the new vectors xₖ₊₁ and pₖ₊₁.
-
inline bool apply(rvec q, real_t γ)¶
Apply the inverse Jacobian approximation to the given vector q, i.e.
\( q \leftarrow H_k q \). Initial inverse Hessian approximation is set to \( H_0 = \gamma I \). The result is scaled by a factor
γ. Ifγis negative, the result is scaled by \( \frac{s^\top y}{y^\top y} \).
-
inline void reset()¶
Throw away the approximation and all previous vectors s and y.
-
inline void resize(length_t n)¶
Re-allocate storage for a problem with a different size.
Causes a reset.
-
inline index_t succ(index_t i) const¶
Get the next index in the circular buffer of previous s and y vectors.
-
inline index_t pred(index_t i) const¶
Get the previous index in the circular buffer of s and y vectors.
-
inline length_t current_history() const¶
Get the number of previous s and y vectors currently stored in the buffer.
- inline auto s̃ (index_t i)
- inline auto s̃ (index_t i) const
-
struct BroydenGoodParams¶
- #include <quala/broyden-good.hpp>
Parameters for the BroydenGood class.
Public Members
-
bool force_pos_def = false¶
If set to true, the inverse Jacobian estimate should remain definite.
-
bool restarted = true¶
If set to true, the buffer is cleared after
memoryiterations.If set to false, a circular buffer is used, replacing the oldest pair of vectors, as a limited-memory type algorithm. However, since each vector s̃ depends on the previous vectors, this algorithm is not theoretically correct, because the s̃ vectors are not updated when the old history is overwritten.
-
bool force_pos_def = false¶
-
struct BroydenStorage¶
- #include <quala/broyden-good.hpp>
Layout:
┌───── 2 m + 1 ─────┐ ┌ ┌───┬───┬───┬───┬───┐ │ │ │ │ │ │ │ n │ │ s │ s̃ │ s │ s̃ │ w │ │ │ │ │ │ │ │ └ └───┴───┴───┴───┴───┘
Public Types
-
template<class IndexT = size_t>
struct CircularIndexIterator¶ - #include <quala/util/ringbuffer.hpp>
Public Types
-
using Indices = CircularIndices<Index>¶
-
using reference = value_type¶
-
using difference_type = std::ptrdiff_t¶
-
using pointer = void¶
-
using iterator_category = std::input_iterator_tag¶
Public Functions
-
inline CircularIndexIterator()¶
-
inline CircularIndexIterator &operator++()¶
-
inline CircularIndexIterator &operator--()¶
-
inline CircularIndexIterator operator++(int)¶
-
inline CircularIndexIterator operator--(int)¶
-
template<class IndexT>
bool operator==(CircularIndexIterator<IndexT> a, CircularIndexIterator<IndexT> b)¶ Note
Only valid for two indices in the same range.
-
template<class IndexT>
bool operator!=(CircularIndexIterator<IndexT> a, CircularIndexIterator<IndexT> b)¶ Note
Only valid for two indices in the same range.
-
using Indices = CircularIndices<Index>¶
-
template<class IndexT = size_t>
struct CircularIndices¶ - #include <quala/util/ringbuffer.hpp>
-
template<class IndexT>
bool operator==(CircularIndices<IndexT> a, CircularIndices<IndexT> b)¶ Note
Only valid for two indices in the same range.
-
template<class IndexT>
bool operator!=(CircularIndices<IndexT> a, CircularIndices<IndexT> b)¶ Note
Only valid for two indices in the same range.
-
template<class IndexT>
-
template<class IndexT>
class CircularRange¶ - #include <quala/util/ringbuffer.hpp>
Public Types
-
using Indices = CircularIndices<Index>¶
-
using const_iterator = CircularIndexIterator<Index>¶
-
using iterator = const_iterator¶
-
using const_reverse_iterator = ReverseCircularIndexIterator<Index>¶
-
using reverse_iterator = const_reverse_iterator¶
Public Functions
-
inline const_iterator cbegin() const¶
-
inline const_iterator cend() const¶
-
inline reverse_iterator rbegin() const¶
-
inline reverse_iterator rend() const¶
-
inline const_reverse_iterator crbegin() const¶
-
inline const_reverse_iterator crend() const¶
-
using Indices = CircularIndices<Index>¶
-
class LBFGS¶
- #include <quala/decl/lbfgs.hpp>
Limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm.
Public Types
-
enum class Sign¶
The sign of the vectors \( p \) passed to the LBFGS::update method.
Values:
-
enumerator Positive¶
\( p \sim \nabla \psi(x) \)
-
enumerator Negative¶
\( p \sim -\nabla \psi(x) \)
-
enumerator Positive¶
-
using Params = LBFGSParams¶
Public Functions
-
template<class VecS, class VecY>
inline bool update_sy(const anymat<VecS> &s, const anymat<VecY> &y, real_t pₖ₊₁ᵀpₖ₊₁, bool forced = false)¶ Update the inverse Hessian approximation using the new vectors sₖ = xₖ₊₁ - xₖ and yₖ = pₖ₊₁ - pₖ.
-
inline bool update(crvec xₖ, crvec xₖ₊₁, crvec pₖ, crvec pₖ₊₁, Sign sign = Sign::Positive, bool forced = false)¶
Update the inverse Hessian approximation using the new vectors xₖ₊₁ and pₖ₊₁.
-
inline bool apply(rvec q, real_t γ = -1)¶
Apply the inverse Hessian approximation to the given vector q.
Initial inverse Hessian approximation is set to \( H_0 = \gamma I \). If
γis negative, \( H_0 = \frac{s^\top y}{y^\top y} I \).
-
template<class IndexVec>
bool apply(rvec q, real_t γ, const IndexVec &J)¶ Apply the inverse Hessian approximation to the given vector q, applying only the columns and rows of the Hessian in the index set J.
-
inline void reset()¶
Throw away the approximation and all previous vectors s and y.
-
inline void resize(length_t n)¶
Re-allocate storage for a problem with a different size.
Causes a reset.
-
inline index_t succ(index_t i) const¶
Get the next index in the circular buffer of previous s and y vectors.
-
inline index_t pred(index_t i) const¶
Get the previous index in the circular buffer of s and y vectors.
-
inline length_t current_history() const¶
Get the number of previous s and y vectors currently stored in the buffer.
Public Static Functions
-
static inline bool update_valid(const LBFGSParams ¶ms, real_t yᵀs, real_t sᵀs, real_t pᵀp)¶
Check if the new vectors s and y allow for a valid BFGS update that preserves the positive definiteness of the Hessian approximation.
-
enum class Sign¶
-
struct LBFGSParams¶
- #include <quala/decl/lbfgs.hpp>
Parameters for the LBFGS class.
Public Members
-
CBFGSParams cbfgs¶
Parameters in the cautious BFGS update condition.
\[ \frac{y^\top s}{s^\top s} \ge \epsilon \| g \|^\alpha \]
-
bool force_pos_def = true¶
If set to true, the inverse Hessian estimate should remain definite, i.e.
a check is performed that rejects the update if \( y^\top s \le \text{min_div_fac} \cdot s^\top s \). If set to false, just try to prevent a singular Hessian by rejecting the update if \( \left| y^\top s \right| \le \text{min_div_fac} \cdot s^\top s \).
-
CBFGSParams cbfgs¶
-
struct LBFGSStorage¶
- #include <quala/decl/lbfgs.hpp>
Layout:
┌───── 2 m ─────┐ ┌ ┌───┬───┬───┬───┐ │ │ │ │ │ │ │ │ s │ y │ s │ y │ n+1 │ │ │ │ │ │ │ ├───┼───┼───┼───┤ │ │ ρ │ α │ ρ │ α │ └ └───┴───┴───┴───┘Public Types
-
class LimitedMemoryQR¶
- #include <quala/detail/limited-memory-qr.hpp>
Incremental QR factorization using modified Gram-Schmidt with reorthogonalization.
Computes A = QR while allowing efficient removal of the first column of A or adding new columns at the end of A.
Public Types
Public Functions
-
LimitedMemoryQR() = default¶
-
inline LimitedMemoryQR(length_t n, length_t m)¶
The maximum dimensions of Q are n×m and the maximum dimensions of R are m×m.
- Parameters
n – The size of the vectors, the number of rows of A.
m – The maximum number of columns of A.
-
inline void remove_column()¶
Remove the leftmost column.
-
template<class VecB, class VecX>
inline void solve_col(const VecB &b, VecX &x) const¶ Solve the least squares problem Ax = b.
-
template<class MatB, class MatX>
inline void solve(const MatB &B, MatX &X) const¶ Solve the least squares problem AX = B.
-
template<class Derived>
inline solve_ret_t<Derived> solve(const Eigen::DenseBase<Derived> &B)¶ Solve the least squares problem AX = B.
-
inline const mat &get_raw_R() const¶
Get the full, raw storage for the upper triangular matrix R.
The columns of this matrix are permuted because it’s stored as a circular buffer for efficiently appending columns to the end and popping columns from the front.
-
inline mat get_full_R() const¶
Get the full storage for the upper triangular matrix R but with the columns in the correct order.
Note
Meant for tests only, creates a permuted copy.
-
inline mat get_R() const¶
Get the matrix R such that Q times R is the original matrix.
Note
Meant for tests only, creates a permuted copy.
-
inline mat get_Q() const¶
Get the matrix Q such that Q times R is the original matrix.
Note
Meant for tests only, creates a copy.
-
inline unsigned long get_reorth_count() const¶
Get the number of MGS reorthogonalizations.
-
inline void clear_reorth_count()¶
Reset the number of MGS reorthogonalizations.
-
inline void reset()¶
Reset all indices, clearing the Q and R matrices.
-
inline void resize(length_t n, length_t m)¶
Re-allocate storage for a problem with a different size.
Causes a reset.
-
inline index_t ring_head() const¶
Get the head index of the circular buffer (points to the oldest element).
-
inline index_t ring_tail() const¶
Get the tail index of the circular buffer (points to one past the most recent element).
-
inline CircularRange<index_t> ring_iter() const¶
Get iterators in the circular buffer.
-
inline ReverseCircularRange<index_t> ring_reverse_iter() const¶
Get reverse iterators in the circular buffer.
-
LimitedMemoryQR() = default¶
-
template<class T>
class MaxHistory¶ - #include <quala/util/ringbuffer.hpp>
-
template<class IndexT = size_t>
struct ReverseCircularIndexIterator¶ - #include <quala/util/ringbuffer.hpp>
Public Types
-
using ForwardIterator = CircularIndexIterator<IndexT>¶
-
using Index = typename ForwardIterator::Index¶
-
using Indices = typename ForwardIterator::Indices¶
-
using reference = value_type¶
-
using difference_type = std::ptrdiff_t¶
-
using pointer = void¶
-
using iterator_category = std::input_iterator_tag¶
Public Functions
-
inline ReverseCircularIndexIterator()¶
-
inline ReverseCircularIndexIterator(ForwardIterator forwardit)¶
-
inline ReverseCircularIndexIterator &operator++()¶
-
inline ReverseCircularIndexIterator &operator--()¶
-
inline ReverseCircularIndexIterator operator++(int)¶
-
inline ReverseCircularIndexIterator operator--(int)¶
Public Members
-
ForwardIterator forwardit¶
-
template<class IndexT>
bool operator==(ReverseCircularIndexIterator<IndexT> a, ReverseCircularIndexIterator<IndexT> b)¶ Note
Only valid for two indices in the same range.
-
template<class IndexT>
bool operator!=(ReverseCircularIndexIterator<IndexT> a, ReverseCircularIndexIterator<IndexT> b)¶ Note
Only valid for two indices in the same range.
-
using ForwardIterator = CircularIndexIterator<IndexT>¶
-
template<class IndexT>
class ReverseCircularRange¶ - #include <quala/util/ringbuffer.hpp>
Public Types
-
using ForwardRange = CircularRange<IndexT>¶
-
using Index = typename ForwardRange::Index¶
-
using Indices = typename ForwardRange::Indices¶
-
using const_iterator = typename ForwardRange::const_reverse_iterator¶
-
using iterator = typename ForwardRange::reverse_iterator¶
-
using const_reverse_iterator = typename ForwardRange::const_iterator¶
-
using reverse_iterator = typename ForwardRange::iterator¶
Public Functions
-
inline ReverseCircularRange(const ForwardRange &forwardrange)¶
-
inline const_iterator cbegin() const¶
-
inline const_iterator cend() const¶
-
inline reverse_iterator rbegin() const¶
-
inline reverse_iterator rend() const¶
-
inline const_reverse_iterator crbegin() const¶
-
inline const_reverse_iterator crend() const¶
-
using ForwardRange = CircularRange<IndexT>¶
-
class vec_allocator¶
- #include <quala/util/alloc.hpp>
Public Functions
-
inline vec_allocator(size_t num_vec, Eigen::Index n)¶
-
vec_allocator(const vec_allocator&) = delete¶
-
vec_allocator(vec_allocator&&) = delete¶
-
vec_allocator &operator=(const vec_allocator&) = delete¶
-
vec_allocator &operator=(vec_allocator&&) = delete¶
-
inline auto alloc()¶
-
inline alloc_raii_wrapper alloc_raii()¶
-
inline size_t size() const¶
-
inline size_t used_space() const¶
-
inline Eigen::Index vector_size() const¶
-
inline size_t highwater() const¶
-
struct alloc_raii_wrapper¶
- #include <quala/util/alloc.hpp>
Public Functions
-
inline alloc_raii_wrapper(real_t *dptr, Eigen::Index n, vec_allocator *alloc)¶
-
inline alloc_raii_wrapper(mvec &&v, vec_allocator *alloc)¶
-
inline ~alloc_raii_wrapper()¶
-
alloc_raii_wrapper(const alloc_raii_wrapper&) = delete¶
-
inline alloc_raii_wrapper(alloc_raii_wrapper &&o) noexcept¶
-
alloc_raii_wrapper &operator=(const alloc_raii_wrapper&) = delete¶
-
inline alloc_raii_wrapper &operator=(alloc_raii_wrapper &&o) noexcept¶
-
inline alloc_raii_wrapper &operator=(crvec v)¶
-
inline alloc_raii_wrapper(real_t *dptr, Eigen::Index n, vec_allocator *alloc)¶
-
inline vec_allocator(size_t num_vec, Eigen::Index n)¶
-
namespace vec_util¶
Functions
-
template<class V, class M>
auto norm_squared_weighted(V &&v, M &&Σ)¶ Get the Σ norm squared of a given vector, with Σ a diagonal matrix.
- Returns
\( \langle v, \Sigma v \rangle \)
-
template<class V, class M>
-
using real_t = double¶