vl_bfgs.py 11.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
Theo Steininger's avatar
Theo Steininger committed
13
14
15
16
17
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
theos's avatar
theos committed
18

Martin Reinecke's avatar
Martin Reinecke committed
19
20
21
from __future__ import division
from builtins import range
from builtins import object
22
23
import numpy as np

24
from .descent_minimizer import DescentMinimizer
25
from .line_searching import LineSearchStrongWolfe
theos's avatar
theos committed
26
27


28
class VL_BFGS(DescentMinimizer):
29
30
    def __init__(self, line_searcher=LineSearchStrongWolfe(), callback=None,
                 convergence_tolerance=1E-4, convergence_level=3,
31
                 iteration_limit=None, max_history_length=5):
32
33
34
35
36
37
38
39
40
41

        super(VL_BFGS, self).__init__(
                                line_searcher=line_searcher,
                                callback=callback,
                                convergence_tolerance=convergence_tolerance,
                                convergence_level=convergence_level,
                                iteration_limit=iteration_limit)

        self.max_history_length = max_history_length

42
    def __call__(self, energy):
43
        self._information_store = None
44
        return super(VL_BFGS, self).__call__(energy)
45

46
47
48
49
50
51
52
53
54
    def get_descend_direction(self, energy):
        """Implementation of the Vector-free L-BFGS minimization scheme.

        Find the descent direction by using the inverse Hessian.
        Instead of storing the whole matrix, it stores only the last few
        updates, which are used to do operations requiring the inverse
        Hessian product. The updates are represented in a new basis to optimize
        the algorithm.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
55
56
        Parameters
        ----------
57
58
59
60
        energy : Energy
            An instance of the Energy class which shall be minized. The
            position of `energy` is used as the starting point of minization.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
61
62
63
        Returns
        -------
        descend_direction : Field
64
65
66
67
68
69
70
            Returns the descent direction.

        References
        ----------
        W. Chen, Z. Wang, J. Zhou, "Large-scale L-BFGS using MapReduce", 2014,
        Microsoft

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
71
        """
72

73
74
        x = energy.position
        gradient = energy.gradient
75
76
77
78
79
80
81
82
83
84
85
86
        # initialize the information store if it doesn't already exist
        try:
            self._information_store.add_new_point(x, gradient)
        except AttributeError:
            self._information_store = InformationStore(self.max_history_length,
                                                       x0=x,
                                                       gradient=gradient)

        b = self._information_store.b
        delta = self._information_store.delta

        descend_direction = delta[0] * b[0]
Martin Reinecke's avatar
Martin Reinecke committed
87
        for i in range(1, len(delta)):
88
89
90
            descend_direction += delta[i] * b[i]

        return descend_direction
theos's avatar
theos committed
91
92
93


class InformationStore(object):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
94
    """Class for storing a list of past updates.
95

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
96
97
98
99
100
101
102
103
    Parameters
    ----------
    max_history_length : integer
        Maximum number of stored past updates.
    x0 : Field
        Initial position in variable space.
    gradient : Field
        Gradient at position x0.
104

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
    Attributes
    ----------
    max_history_length : integer
        Maximum number of stored past updates.
    s : List
        List of past position differences, which are Fields.
    y : List
        List of past gradient differences, which are Fields.
    last_x : Field
        Initial position in variable space.
    last_gradient : Field
        Gradient at initial position.
    k : integer
        Number of currently stored past updates.
    _ss_store : dictionary
        Dictionary of scalar products between different elements of s.
    _sy_store : dictionary
        Dictionary of scalar products between elements of s and y.
    _yy_store : dictionary
        Dictionary of scalar products between different elements of y.
125

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
126
    """
127
128
129
130
    def __init__(self, max_history_length, x0, gradient):
        self.max_history_length = max_history_length
        self.s = LimitedList(max_history_length)
        self.y = LimitedList(max_history_length)
131
132
        self.last_x = x0.copy()
        self.last_gradient = gradient.copy()
theos's avatar
theos committed
133
        self.k = 0
134
135
136
137
138
139
140

        self._ss_store = {}
        self._sy_store = {}
        self._yy_store = {}

    @property
    def history_length(self):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
141
        """Returns the number of currently stored updates.
142

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
143
        """
144
145
146
147
        return min(self.k, self.max_history_length)

    @property
    def b(self):
148
149
        """Combines s, y and gradient to form the new base vectors b.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
150
151
152
153
        Returns
        -------
        result : List
            List of new basis vectors.
154

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
155
        """
156
157
158
159
160
        result = []
        m = self.history_length
        k = self.k

        s = self.s
Martin Reinecke's avatar
Martin Reinecke committed
161
        for i in range(m):
162
163
164
            result.append(s[k-m+i])

        y = self.y
Martin Reinecke's avatar
Martin Reinecke committed
165
        for i in range(m):
166
167
168
169
170
171
172
173
            result.append(y[k-m+i])

        result.append(self.last_gradient)

        return result

    @property
    def b_dot_b(self):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
174
        """Generates the (2m+1) * (2m+1) scalar matrix.
175

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
176
        The i,j-th element of the matrix is a scalar product between the i-th
177
178
        and j-th base vector.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
179
180
        Returns
        -------
181
        result : numpy.ndarray
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
182
            Scalar matrix.
183

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
184
        """
185
186
187
188
        m = self.history_length
        k = self.k
        result = np.empty((2*m+1, 2*m+1), dtype=np.float)

Martin Reinecke's avatar
Martin Reinecke committed
189
190
        for i in range(m):
            for j in range(m):
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
                result[i, j] = self.ss_store(k-m+i, k-m+j)

                sy_ij = self.sy_store(k-m+i, k-m+j)
                result[i, m+j] = sy_ij
                result[m+j, i] = sy_ij

                result[m+i, m+j] = self.yy_store(k-m+i, k-m+j)

            sgrad_i = self.sgrad_store(k-m+i)
            result[2*m, i] = sgrad_i
            result[i, 2*m] = sgrad_i

            ygrad_i = self.ygrad_store(k-m+i)
            result[2*m, m+i] = ygrad_i
            result[m+i, 2*m] = ygrad_i

        result[2*m, 2*m] = self.gradgrad_store()

        return result
theos's avatar
theos committed
210
211

    @property
212
    def delta(self):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
213
        """Calculates the new scalar coefficients (deltas).
214

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
215
216
217
218
        Returns
        -------
        delta : List
            List of the new scalar coefficients (deltas).
219

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
220
        """
221
222
223
224
225
226
227
228
        m = self.history_length
        b_dot_b = self.b_dot_b

        delta = np.zeros(2*m+1, dtype=np.float)
        delta[2*m] = -1

        alpha = np.empty(m, dtype=np.float)

Martin Reinecke's avatar
Martin Reinecke committed
229
230
        for j in range(m-1, -1, -1):
            delta_b_b = sum([delta[l] * b_dot_b[l, j] for l in range(2*m+1)])
231
232
233
            alpha[j] = delta_b_b/b_dot_b[j, m+j]
            delta[m+j] -= alpha[j]

Martin Reinecke's avatar
Martin Reinecke committed
234
        for i in range(2*m+1):
235
236
            delta[i] *= b_dot_b[m-1, 2*m-1]/b_dot_b[2*m-1, 2*m-1]

Martin Reinecke's avatar
Martin Reinecke committed
237
238
        for j in range(m-1, -1, -1):
            delta_b_b = sum([delta[l]*b_dot_b[m+j, l] for l in range(2*m+1)])
239
240
241
242
243
244
            beta = delta_b_b/b_dot_b[j, m+j]
            delta[j] += (alpha[j] - beta)

        return delta

    def ss_store(self, i, j):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
245
        """Updates the dictionary _ss_store with a new scalar product.
246
247
248

        Returns the scalar product of s_i and s_j.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
249
250
251
252
253
254
        Parameters
        ----------
        i : integer
            s index.
        j : integer
            s index.
255

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
256
257
258
259
        Returns
        -------
        _ss_store[key] : float
            Scalar product of s_i and s_j.
260

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
261
        """
262
263
        key = tuple(sorted((i, j)))
        if key not in self._ss_store:
Martin Reinecke's avatar
Martin Reinecke committed
264
            self._ss_store[key] = self.s[i].vdot(self.s[j])
265
266
267
        return self._ss_store[key]

    def sy_store(self, i, j):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
268
        """Updates the dictionary _sy_store with a new scalar product.
269
270
271

        Returns the scalar product of s_i and y_j.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
272
273
274
275
276
277
        Parameters
        ----------
        i : integer
            s index.
        j : integer
            y index.
278

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
279
280
281
282
        Returns
        -------
        _sy_store[key] : float
            Scalar product of s_i and y_j.
283

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
284
        """
285
286
        key = (i, j)
        if key not in self._sy_store:
Martin Reinecke's avatar
Martin Reinecke committed
287
            self._sy_store[key] = self.s[i].vdot(self.y[j])
288
289
290
        return self._sy_store[key]

    def yy_store(self, i, j):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
291
        """Updates the dictionary _yy_store with a new scalar product.
292
293
294

        Returns the scalar product of y_i and y_j.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
295
296
297
298
299
300
        Parameters
        ----------
        i : integer
            y index.
        j : integer
            y index.
301

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
302
303
304
305
        Returns
        ------
        _yy_store[key] : float
            Scalar product of y_i and y_j.
306

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
307
        """
308
309
        key = tuple(sorted((i, j)))
        if key not in self._yy_store:
Martin Reinecke's avatar
Martin Reinecke committed
310
            self._yy_store[key] = self.y[i].vdot(self.y[j])
311
312
313
        return self._yy_store[key]

    def sgrad_store(self, i):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
314
        """Returns scalar product between s_i and gradient on initial position.
315

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
316
317
318
319
        Returns
        -------
        scalar product : float
            Scalar product.
320

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
321
        """
Martin Reinecke's avatar
Martin Reinecke committed
322
        return self.s[i].vdot(self.last_gradient)
323
324

    def ygrad_store(self, i):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
325
        """Returns scalar product between y_i and gradient on initial position.
326

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
327
328
329
330
        Returns
        -------
        scalar product : float
            Scalar product.
331

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
332
        """
Martin Reinecke's avatar
Martin Reinecke committed
333
        return self.y[i].vdot(self.last_gradient)
334
335

    def gradgrad_store(self):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
336
        """Returns scalar product of gradient on initial position with itself.
337

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
338
339
340
341
        Returns
        -------
        scalar product : float
            Scalar product.
342

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
343
        """
Martin Reinecke's avatar
Martin Reinecke committed
344
        return self.last_gradient.vdot(self.last_gradient)
theos's avatar
theos committed
345
346

    def add_new_point(self, x, gradient):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
347
348
        """Updates the s list and y list.

349
350
351
        Calculates the new position and gradient differences and adds them to
        the respective list.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
352
        """
theos's avatar
theos committed
353
354
355
356
357
358
359
360
        self.k += 1

        new_s = x - self.last_x
        self.s.add(new_s)

        new_y = gradient - self.last_gradient
        self.y.add(new_y)

361
362
        self.last_x = x.copy()
        self.last_gradient = gradient.copy()
theos's avatar
theos committed
363
364
365


class LimitedList(object):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
366
    """Class for creating a list of limited length.
367

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
368
369
370
371
    Parameters
    ----------
    history_length : integer
        Maximum number of stored past updates.
372

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
373
374
375
376
377
378
379
380
381
    Attributes
    ----------
    history_length : integer
        Maximum number of stored past updates.
    _offset : integer
        Offset to correct the indices which are bigger than maximum history.
        length.
    _storage : list
        List where input values are stored.
382

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
383
    """
theos's avatar
theos committed
384
385
386
387
388
389
    def __init__(self, history_length):
        self.history_length = int(history_length)
        self._offset = 0
        self._storage = []

    def __getitem__(self, index):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
390
        """Returns the element with index [index-offset].
391

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
392
393
394
395
        Parameters
        ----------
        index : integer
            Index of the selected element.
396

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
397
398
399
400
        Returns
        -------
            selected element
        """
theos's avatar
theos committed
401
402
403
        return self._storage[index-self._offset]

    def add(self, value):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
404
        """Adds a new element to the list.
405

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
406
407
        If the list is of length maximum history then it removes the first
        element first.
408

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
409
410
411
412
        Parameters
        ----------
        value : anything
            New element in the list.
413

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
414
        """
theos's avatar
theos committed
415
416
417
418
        if len(self._storage) == self.history_length:
            self._storage.pop(0)
            self._offset += 1
        self._storage.append(value)