vl_bfgs.py 11.5 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.
Theo Steininger's avatar
Theo Steininger committed
18

19
20
import numpy as np

21
from .descent_minimizer import DescentMinimizer
22
from .line_searching import LineSearchStrongWolfe
Theo Steininger's avatar
Theo Steininger committed
23
24


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

        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

39
    def __call__(self, energy):
40
        self._information_store = None
41
        return super(VL_BFGS, self).__call__(energy)
42

43
44
45
46
47
48
49
50
51
    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
52
53
        Parameters
        ----------
54
55
56
57
        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
58
59
60
        Returns
        -------
        descend_direction : Field
61
62
63
64
65
66
67
            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
68
        """
69

70
71
        x = energy.position
        gradient = energy.gradient
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
        # 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]
        for i in xrange(1, len(delta)):
            descend_direction += delta[i] * b[i]

        return descend_direction
Theo Steininger's avatar
Theo Steininger committed
88
89
90


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

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
    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.
122

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
123
    """
124
125
126
127
    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)
128
129
        self.last_x = x0.copy()
        self.last_gradient = gradient.copy()
Theo Steininger's avatar
Theo Steininger committed
130
        self.k = 0
131
132
133
134
135
136
137

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

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
140
        """
141
142
143
144
        return min(self.k, self.max_history_length)

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
147
148
149
150
        Returns
        -------
        result : List
            List of new basis vectors.
151

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
152
        """
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
        result = []
        m = self.history_length
        k = self.k

        s = self.s
        for i in xrange(m):
            result.append(s[k-m+i])

        y = self.y
        for i in xrange(m):
            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
171
        """Generates the (2m+1) * (2m+1) scalar matrix.
172

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
176
177
        Returns
        -------
178
        result : numpy.ndarray
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
179
            Scalar matrix.
180

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
181
        """
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
        m = self.history_length
        k = self.k
        result = np.empty((2*m+1, 2*m+1), dtype=np.float)

        for i in xrange(m):
            for j in xrange(m):
                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
Theo Steininger's avatar
Theo Steininger committed
207
208

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

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
217
        """
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
        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)

        for j in xrange(m-1, -1, -1):
            delta_b_b = sum([delta[l] * b_dot_b[l, j] for l in xrange(2*m+1)])
            alpha[j] = delta_b_b/b_dot_b[j, m+j]
            delta[m+j] -= alpha[j]

        for i in xrange(2*m+1):
            delta[i] *= b_dot_b[m-1, 2*m-1]/b_dot_b[2*m-1, 2*m-1]

        for j in xrange(m-1, -1, -1):
            delta_b_b = sum([delta[l]*b_dot_b[m+j, l] for l in xrange(2*m+1)])
            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
242
        """Updates the dictionary _ss_store with a new scalar product.
243
244
245

        Returns the scalar product of s_i and s_j.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
246
247
248
249
250
251
        Parameters
        ----------
        i : integer
            s index.
        j : integer
            s index.
252

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

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

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

        Returns the scalar product of s_i and y_j.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
269
270
271
272
273
274
        Parameters
        ----------
        i : integer
            s index.
        j : integer
            y index.
275

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

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

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

        Returns the scalar product of y_i and y_j.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
292
293
294
295
296
297
        Parameters
        ----------
        i : integer
            y index.
        j : integer
            y index.
298

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

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

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
313
314
315
316
        Returns
        -------
        scalar product : float
            Scalar product.
317

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
318
        """
Martin Reinecke's avatar
Martin Reinecke committed
319
        return self.s[i].vdot(self.last_gradient)
320
321

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
324
325
326
327
        Returns
        -------
        scalar product : float
            Scalar product.
328

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
329
        """
Martin Reinecke's avatar
Martin Reinecke committed
330
        return self.y[i].vdot(self.last_gradient)
331
332

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
335
336
337
338
        Returns
        -------
        scalar product : float
            Scalar product.
339

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
340
        """
Martin Reinecke's avatar
Martin Reinecke committed
341
        return self.last_gradient.vdot(self.last_gradient)
Theo Steininger's avatar
Theo Steininger committed
342
343

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

346
347
348
        Calculates the new position and gradient differences and adds them to
        the respective list.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
349
        """
Theo Steininger's avatar
Theo Steininger committed
350
351
352
353
354
355
356
357
        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)

358
359
        self.last_x = x.copy()
        self.last_gradient = gradient.copy()
Theo Steininger's avatar
Theo Steininger committed
360
361
362


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

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
370
371
372
373
374
375
376
377
378
    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.
379

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
380
    """
Theo Steininger's avatar
Theo Steininger committed
381
382
383
384
385
386
    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
387
        """Returns the element with index [index-offset].
388

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
389
390
391
392
        Parameters
        ----------
        index : integer
            Index of the selected element.
393

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
394
395
396
397
        Returns
        -------
            selected element
        """
Theo Steininger's avatar
Theo Steininger committed
398
399
400
        return self._storage[index-self._offset]

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

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
406
407
408
409
        Parameters
        ----------
        value : anything
            New element in the list.
410

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
411
        """
Theo Steininger's avatar
Theo Steininger committed
412
413
414
415
        if len(self._storage) == self.history_length:
            self._storage.pop(0)
            self._offset += 1
        self._storage.append(value)