Line data Source code
1 : /*
2 : * File: LBFGSBuffer.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on April 13, 2016, 2:32 PM
6 : *
7 : * Original source code:
8 : * https://github.com/lostella/libLBFGS
9 : *
10 : * ForBES is free software: you can redistribute it and/or modify
11 : * it under the terms of the GNU Lesser General Public License as published by
12 : * the Free Software Foundation, either version 3 of the License, or
13 : * (at your option) any later version.
14 : *
15 : * ForBES is distributed in the hope that it will be useful,
16 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 : * GNU Lesser General Public License for more details.
19 : *
20 : * You should have received a copy of the GNU Lesser General Public License
21 : * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
22 : */
23 :
24 : #include "LBFGSBuffer.h"
25 : #include "Matrix.h"
26 : #include "MatrixFactory.h"
27 : #include <iostream>
28 : #include <string.h>
29 : #include <complex>
30 :
31 : const long LBFGSBuffer::LBFGS_BUFFER_EOB = -1;
32 :
33 5 : LBFGSBuffer::LBFGSBuffer(size_t n, size_t mem) : m_mem(mem) {
34 : // initialize the LBFGS buffer
35 5 : m_current_mem = 0;
36 5 : m_idx = -1;
37 : // allocate memory for matrices
38 5 : m_S = new Matrix(n, mem);
39 5 : m_Y = new Matrix(n, mem);
40 5 : m_Ys = new Matrix(mem, 1);
41 5 : m_alphas = new Matrix(mem, 1);
42 5 : }
43 :
44 10 : LBFGSBuffer::~LBFGSBuffer() {
45 5 : if (m_S != NULL) {
46 5 : delete m_S;
47 : }
48 5 : if (m_Y != NULL) {
49 5 : delete m_Y;
50 : }
51 5 : if (m_Ys != NULL) {
52 5 : delete m_Ys;
53 : }
54 5 : if (m_alphas != NULL) {
55 5 : delete m_alphas;
56 : }
57 5 : m_mem = 0;
58 5 : m_current_mem = 0;
59 5 : m_idx = 0;
60 10 : }
61 :
62 11 : Matrix* LBFGSBuffer::get_S() const {
63 11 : return m_S;
64 : }
65 :
66 10 : Matrix* LBFGSBuffer::get_Y() const {
67 10 : return m_Y;
68 : }
69 :
70 12 : Matrix* LBFGSBuffer::get_Ys() const {
71 12 : return m_Ys;
72 : }
73 :
74 2 : Matrix* LBFGSBuffer::get_alphas() const {
75 2 : return m_alphas;
76 : }
77 :
78 11 : long LBFGSBuffer::cursor() const {
79 11 : return m_idx;
80 : }
81 :
82 32 : int LBFGSBuffer::push(Matrix* s, Matrix* y) {
83 32 : m_idx++;
84 32 : if (m_idx >= m_mem) {
85 9 : m_idx = 0; // start over
86 : }
87 32 : m_current_mem++; // increase current memory
88 32 : if (m_current_mem > m_mem) {
89 17 : m_current_mem = m_mem; // buffer is now full
90 : }
91 :
92 32 : size_t n = m_Y->getNrows();
93 :
94 : // copy data into m_Y
95 32 : double * m_Y_data = m_Y->getData();
96 32 : double * y_data = y->getData();
97 32 : m_Y_data += m_idx*n;
98 32 : memcpy(m_Y_data, y_data, n * sizeof (double));
99 :
100 : // copy data into m_S
101 32 : double * m_S_data = m_S->getData();
102 32 : double * s_data = s->getData();
103 32 : m_S_data += m_idx*n;
104 32 : memcpy(m_S_data, s_data, n * sizeof (double));
105 :
106 32 : Matrix s_mat = MatrixFactory::ShallowVector(s_data, n, 0);
107 64 : Matrix y_mat = MatrixFactory::ShallowVector(y_data, n, 0);
108 64 : Matrix ys = y_mat*s_mat;
109 32 : m_Ys->set(m_idx, 0, ys[0]);
110 64 : return ForBESUtils::STATUS_OK;
111 : }
112 :
113 0 : int LBFGSBuffer::reset() {
114 0 : m_idx = -1;
115 0 : m_current_mem = 0;
116 0 : return ForBESUtils::STATUS_OK;
117 : }
118 :
119 0 : size_t LBFGSBuffer::get_current_mem() const {
120 0 : return m_current_mem;
121 : }
122 :
123 1 : double LBFGSBuffer::hessian_estimate() {
124 : // returns Hk0 = ys_prev / (y_prev'*y_prev)
125 : // or 1.0 is no ys_prev is cached
126 1 : size_t idx_current = cursor();
127 1 : size_t n = m_Y->getNrows();
128 1 : double gamma_k_0 = 1.0;
129 1 : if (m_current_mem > 0) {
130 1 : Matrix curr_y = MatrixFactory::ShallowSubVector(*m_Y, idx_current);
131 1 : double sq_norm_y = 0.0;
132 11 : for (size_t l = 0; l < n; l++) {
133 10 : sq_norm_y += std::pow(curr_y[l], 2);
134 : }
135 1 : gamma_k_0 = (*m_Ys)[idx_current] / sq_norm_y;
136 : }
137 1 : return gamma_k_0;
138 : }
139 :
140 10 : int LBFGSBuffer::update(const Matrix* q, Matrix* r, double& gamma0) {
141 10 : Matrix qq = *q;
142 :
143 : /*
144 : * STEP 1: First loop (compute alpha, update qq)
145 : * [TESTED]
146 : */
147 : long k_minus_j;
148 10 : long j = 0;
149 38 : while ((k_minus_j = get_k_minus_j(j)) != LBFGSBuffer::LBFGS_BUFFER_EOB) {
150 18 : Matrix sq = MatrixFactory::ShallowSubVector(*m_S, k_minus_j) * qq; // <si, q_i>
151 36 : Matrix yi = MatrixFactory::ShallowSubVector(*m_Y, k_minus_j); // yi
152 18 : double alpha_i = sq[0] / (*m_Ys)[k_minus_j];
153 18 : (*m_alphas)[k_minus_j] = alpha_i;
154 18 : Matrix::add(qq, -alpha_i, yi, 1.0);
155 18 : j++;
156 18 : }
157 :
158 : /* Update r */
159 10 : *r = gamma0 * qq;
160 :
161 :
162 :
163 : /*
164 : * STEP 2: Second loop
165 : */
166 10 : j = m_mem;
167 :
168 : // skip empty buffer
169 72 : while (j >= 0 &&
170 52 : (k_minus_j = get_k_minus_j(j)) == LBFGSBuffer::LBFGS_BUFFER_EOB) j--;
171 :
172 :
173 : // iterate backwards
174 38 : while (j >= 0) {
175 18 : k_minus_j = get_k_minus_j(j);
176 18 : Matrix yi = MatrixFactory::ShallowSubVector(*m_Y, k_minus_j);
177 36 : Matrix b = yi * (*r);
178 18 : double beta = b[0] / (*m_Ys)[k_minus_j];
179 36 : Matrix s = MatrixFactory::ShallowSubVector(*m_S, k_minus_j);
180 18 : Matrix::add(*r, (*m_alphas)[k_minus_j] - beta, s, 1.0);
181 18 : j--;
182 18 : }
183 :
184 10 : return ForBESUtils::STATUS_OK;
185 3 : }
186 :
|
