Line data Source code
1 : #pragma once
2 :
3 : #include <panoc-alm/inner/detail/panoc-helpers.hpp>
4 : #include <panoc-alm/inner/directions/decl/specialized-lbfgs.hpp>
5 : #include <panoc-alm/inner/directions/lbfgs.hpp>
6 :
7 : #include <cmath>
8 :
9 : namespace pa {
10 :
11 4 : inline void SpecializedLBFGS::initialize(crvec x₀, crvec grad₀) {
12 4 : idx = 0;
13 4 : full = false;
14 4 : x(0) = x₀;
15 4 : g(0) = grad₀;
16 4 : }
17 :
18 : /// Standard L-BFGS update without changing the step size γ.
19 13 : inline bool SpecializedLBFGS::standard_update(crvec xₖ, crvec xₖ₊₁,
20 : crvec pₖ, crvec pₖ₊₁,
21 : crvec gradₖ₊₁) {
22 13 : const auto s = xₖ₊₁ - xₖ;
23 13 : const auto y = pₖ - pₖ₊₁;
24 :
25 13 : real_t yᵀs = y.dot(s);
26 13 : real_t sᵀs = s.squaredNorm();
27 13 : real_t pᵀp = pₖ₊₁.squaredNorm();
28 13 : real_t ρ = 1 / yᵀs;
29 :
30 13 : if (not LBFGS::update_valid(params, yᵀs, sᵀs, pᵀp))
31 0 : return false;
32 :
33 : // Store the new s and y vectors
34 13 : this->s(idx) = s;
35 13 : this->y(idx) = y;
36 13 : this->ρ(idx) = ρ;
37 :
38 : // Store x and the gradient
39 13 : this->x(succ(idx)) = xₖ₊₁;
40 13 : this->g(succ(idx)) = gradₖ₊₁;
41 :
42 : // Increment the index in the circular buffer
43 13 : idx = succ(idx);
44 13 : full |= idx == 0;
45 :
46 13 : return true;
47 13 : }
48 :
49 : /// L-BFGS update when changing the step size γ, recomputing everything.
50 3 : inline bool SpecializedLBFGS::full_update(crvec xₖ, crvec xₖ₊₁,
51 : crvec pₖ_old_γ, crvec pₖ₊₁,
52 : crvec gradₖ₊₁, const Box &C,
53 : real_t γ) {
54 3 : auto &&sₖ = xₖ₊₁ - xₖ;
55 3 : auto &&yₖ = this->w(); // temporary workspace
56 : // Old pₖ is no longer valid, recompute with new γ
57 : (void)pₖ_old_γ;
58 3 : auto &&pₖ = this->p();
59 3 : pₖ = detail::projected_gradient_step(C, γ, x(idx), g(idx));
60 3 : yₖ = pₖ - pₖ₊₁;
61 :
62 3 : assert(x(idx) == xₖ);
63 :
64 3 : real_t yᵀs = yₖ.dot(sₖ);
65 3 : real_t sᵀs = sₖ.squaredNorm();
66 3 : real_t pᵀp = pₖ₊₁.squaredNorm();
67 3 : real_t ρₖ = 1 / yᵀs;
68 :
69 3 : if (not LBFGS::update_valid(params, yᵀs, sᵀs, pᵀp))
70 0 : return false;
71 :
72 : // Recompute all residuals with new γ
73 : // yₖ = pₖ - pₖ₊₁
74 : // pₖ = Π(-γ∇ψ(xₖ), C - xₖ)
75 3 : size_t endidx = full ? idx : pred(0);
76 8 : for (size_t i = pred(idx); i != endidx; i = pred(i)) {
77 5 : this->y(i) = -pₖ /* i+1 */;
78 5 : pₖ = detail::projected_gradient_step(C, γ, this->x(i), this->g(i));
79 5 : this->y(i) += pₖ /* i */;
80 5 : }
81 : // Store the new s and y vectors
82 3 : this->s(idx) = sₖ;
83 3 : this->y(idx) = yₖ;
84 3 : this->ρ(idx) = ρₖ;
85 :
86 : // Store x and the gradient
87 3 : this->x(succ(idx)) = xₖ₊₁;
88 3 : this->g(succ(idx)) = gradₖ₊₁;
89 :
90 : // Increment the index in the circular buffer
91 3 : idx = succ(idx);
92 3 : full |= idx == 0;
93 :
94 3 : return true;
95 3 : }
96 :
97 16 : inline bool SpecializedLBFGS::update(crvec xₖ, crvec xₖ₊₁,
98 : crvec pₖ, crvec pₖ₊₁,
99 : crvec gradₖ₊₁, const Box &C,
100 : real_t γ) {
101 32 : bool ret = (γ == old_γ || old_γ == 0)
102 13 : ? standard_update(xₖ, xₖ₊₁, pₖ, pₖ₊₁, gradₖ₊₁)
103 3 : : full_update(xₖ, xₖ₊₁, pₖ, pₖ₊₁, gradₖ₊₁, C, γ);
104 16 : old_γ = γ;
105 16 : return ret;
106 16 : }
107 :
108 : template <class Vec>
109 : void SpecializedLBFGS::apply(Vec &&q) {
110 : // TODO: dry, reuse standard LBFGS::apply
111 : auto update1 = [&](size_t i) {
112 : α(i) = ρ(i) * (s(i).dot(q));
113 : q -= α(i) * y(i);
114 : };
115 : if (idx)
116 : for (size_t i = idx; i-- > 0;)
117 : update1(i);
118 : if (full)
119 : for (size_t i = history(); i-- > idx;)
120 : update1(i);
121 :
122 : // q = H₀ * q; // TODO: diagonal matrix H₀?
123 :
124 : auto update2 = [&](size_t i) {
125 : real_t β = ρ(i) * (y(i).dot(q));
126 : q += (α(i) - β) * s(i);
127 : };
128 : if (full)
129 : for (size_t i = idx; i < history(); ++i)
130 : update2(i);
131 : for (size_t i = 0; i < idx; ++i)
132 : update2(i);
133 : }
134 :
135 4 : inline void SpecializedLBFGS::resize(size_t n, size_t history) {
136 4 : sto.resize(n + 1, history * 4 + 2);
137 4 : sto.fill(std::numeric_limits<real_t>::quiet_NaN());
138 4 : idx = 0;
139 4 : full = false;
140 4 : }
141 :
142 : inline void SpecializedLBFGS::reset() {
143 : x(0) = x(idx);
144 : g(0) = x(idx);
145 : idx = 0;
146 : full = false;
147 : }
148 :
149 : } // namespace pa
150 :
151 : #include <panoc-alm/inner/directions/decl/panoc-direction-update.hpp>
152 :
153 : namespace pa {
154 :
155 : template <>
156 : struct PANOCDirection<SpecializedLBFGS> {
157 :
158 : static void initialize(SpecializedLBFGS &lbfgs, crvec x₀,
159 : crvec x̂₀, crvec p₀, crvec grad₀) {
160 : lbfgs.initialize(x₀, grad₀);
161 : (void)x̂₀;
162 : (void)p₀;
163 : }
164 :
165 : static bool update(SpecializedLBFGS &lbfgs, crvec xₖ, crvec xₖ₊₁,
166 : crvec pₖ, crvec pₖ₊₁, crvec gradₖ₊₁,
167 : const Box &C, real_t γ) {
168 : return lbfgs.update(xₖ, xₖ₊₁, pₖ, pₖ₊₁, gradₖ₊₁, C, γ);
169 : }
170 :
171 : static bool apply(SpecializedLBFGS &lbfgs, crvec xₖ, crvec x̂ₖ,
172 : crvec pₖ, real_t γ, rvec qₖ) {
173 : (void)xₖ;
174 : (void)x̂ₖ;
175 : (void)γ; // TODO: add this parameter to SLBFGS
176 : qₖ = pₖ;
177 : lbfgs.apply(qₖ);
178 : return true;
179 : }
180 :
181 : static void changed_γ(SpecializedLBFGS &lbfgs, real_t γₖ, real_t old_γₖ) {
182 : (void)lbfgs;
183 : (void)γₖ;
184 : (void)old_γₖ;
185 : }
186 : };
187 :
188 : } // namespace pa
|
