suanPan
LBFGS.hpp
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
41#ifndef LBFGS_HPP
42#define LBFGS_HPP
43
44#include <suanPan.h>
45#include <deque>
46
47using std::deque;
48using std::vector;
49
50template<typename T> concept Differentiable = requires(T t, const vec& x) {
51 t.evaluate_residual(x); t.evaluate_jacobian(x);
52};
53
54class LBFGS final {
55 deque<vec> hist_ninja, hist_residual;
56 deque<double> hist_factor;
57 vector<double> alpha;
58
59 const unsigned max_hist;
60 const unsigned max_iteration;
61 const double abs_tol;
62 const double rel_tol;
63
64public:
65 explicit LBFGS(const unsigned MH = 10, const unsigned MI = 500, const double AT = 1E-12, const double RT = 1E-12)
66 : max_hist(MH)
67 , max_iteration(MI)
68 , abs_tol(AT)
69 , rel_tol(RT) {}
70
71 template<Differentiable F> int optimize(F& func, vec& x) {
72 // clear container
73 hist_ninja.clear();
74 hist_residual.clear();
75 hist_factor.clear();
76
77 vec ninja;
78 const auto ref_magnitude = norm(x);
79
80 // iteration counter
81 auto counter = 0u;
82 while(true) {
83 const auto residual = func.evaluate_residual(x);
84
85 if(0 == counter) ninja = solve(func.evaluate_jacobian(x), residual);
86 else {
87 // clear temporary factor container
88 alpha.clear();
89 alpha.reserve(hist_ninja.size());
90 // commit current residual
91 hist_residual.back() += residual;
92 // commit current factor after obtaining ninja and residual
93 hist_factor.emplace_back(dot(hist_ninja.back(), hist_residual.back()));
94 // copy current residual to ninja
95 ninja = residual;
96 // perform two-step recursive loop
97 // right side loop
98 for(auto J = static_cast<int>(hist_factor.size()) - 1; J >= 0; --J) {
99 // compute and commit alpha
100 alpha.emplace_back(dot(hist_ninja[J], ninja) / hist_factor[J]);
101 // update ninja
102 ninja -= alpha.back() * hist_residual[J];
103 }
104 // apply the Hessian from the factorization in the first iteration
105 ninja *= dot(hist_residual.back(), hist_ninja.back()) / dot(hist_residual.back(), hist_residual.back());
106 // left side loop
107 for(size_t I = 0, J = hist_factor.size() - 1; I < hist_factor.size(); ++I, --J) ninja += (alpha[J] - dot(hist_residual[I], ninja) / hist_factor[I]) * hist_ninja[I];
108 }
109
110 // commit current displacement increment
111 hist_ninja.emplace_back(-ninja);
112 hist_residual.emplace_back(-residual);
113
114 const auto error = norm(ninja);
115 const auto ref_error = error / ref_magnitude;
116 suanpan_debug("Local iteration error: {:.5E}.\n", ref_error);
117 if(error <= abs_tol && ref_error <= rel_tol) return SUANPAN_SUCCESS;
118 if(++counter > max_iteration) return SUANPAN_FAIL;
119
120 // check if the maximum record number is hit (LBFGS)
121 if(counter > max_hist) {
122 hist_ninja.pop_front();
123 hist_residual.pop_front();
124 hist_factor.pop_front();
125 }
126
127 x -= ninja;
128 }
129 }
130};
131
132#endif
133
The LBFGS class defines a solver using LBFGS iteration method.
Definition: LBFGS.hpp:54
int optimize(F &func, vec &x)
Definition: LBFGS.hpp:71
LBFGS(const unsigned MH=10, const unsigned MI=500, const double AT=1E-12, const double RT=1E-12)
Definition: LBFGS.hpp:65
Definition: LBFGS.hpp:50
std::vector< T > vector
Definition: container.h:53
void error(const std::string_view file_name, const int line, const std::string_view format_str, const T &... args)
Definition: suanPan.h:234
double norm(const vec &)
Definition: tensor.cpp:370
#define suanpan_debug(...)
Definition: suanPan.h:307
constexpr auto SUANPAN_SUCCESS
Definition: suanPan.h:172
constexpr auto SUANPAN_FAIL
Definition: suanPan.h:173