vl_bfgs.py 9.16 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
    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
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
111
        Latest position in variable space.
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
112
    last_gradient : Field
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
113
        Gradient at latest position.
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
114
    k : integer
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
115
        Number of updates that have taken place
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
116
117
118
119
120
121
    _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

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
132
133
134
135
        hl = max_history_length
        self._ss_store = np.empty((hl,hl),dtype=np.float64)
        self._sy_store = np.empty((hl,hl),dtype=np.float64)
        self._yy_store = np.empty((hl,hl),dtype=np.float64)
136
137
138

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

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

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

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
153
        """
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
        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
172
        """Generates the (2m+1) * (2m+1) scalar matrix.
173

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

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

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

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
187
        # update the stores
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
188
189
190
        if (m>0):
            k1 = (k-1)%m

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
191
        for i in xrange(m):
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
192
193
194
195
196
197
            kmi = (k-m+i)%m
            self._ss_store[kmi,k1] = self._ss_store[k1,kmi] \
                = self.s[k-m+i].vdot(self.s[k-1])
            self._yy_store[kmi,k1] = self._yy_store[k1,kmi] \
                = self.y[k-m+i].vdot(self.y[k-1])
            self._sy_store[kmi,k1] = self.s[k-m+i].vdot(self.y[k-1])
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
198
        for j in xrange(m-1):
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
199
            self._sy_store[k1,(k-m+j)%m] = self.s[k-1].vdot(self.y[k-m+j])
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
200

201
        for i in xrange(m):
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
202
            kmi = (k-m+i)%m
203
            for j in xrange(m):
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
204
205
206
                kmj = (k-m+j)%m
                result[i, j] = self._ss_store[kmi, kmj]
                result[i, m+j] = result[m+j, i] = self._sy_store[kmi, kmj]
Martin Reinecke's avatar
fix  
Martin Reinecke committed
207
                result[m+i, m+j] = self._yy_store[kmi, kmj]
208

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
209
            sgrad_i = self.s[k-m+i].vdot(self.last_gradient)
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
210
            result[2*m, i] = result[i, 2*m] = sgrad_i
211

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
212
            ygrad_i = self.y[k-m+i].vdot(self.last_gradient)
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
213
            result[2*m, m+i] = result[m+i, 2*m] = ygrad_i
214

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
215
        result[2*m, 2*m] = self.last_gradient.norm()
216
217

        return result
Theo Steininger's avatar
Theo Steininger committed
218
219

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
223
224
225
226
        Returns
        -------
        delta : List
            List of the new scalar coefficients (deltas).
227

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
228
        """
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
        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]

Martin Reinecke's avatar
fix  
Martin Reinecke committed
245
        for j in xrange(m):
246
247
248
249
250
251
            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

Theo Steininger's avatar
Theo Steininger committed
252
    def add_new_point(self, x, gradient):
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
253
254
        """Updates the s list and y list.

255
256
257
        Calculates the new position and gradient differences and adds them to
        the respective list.

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
258
        """
Martin Reinecke's avatar
fixes  
Martin Reinecke committed
259
260
        self.s[self.k] = x - self.last_x
        self.y[self.k] = gradient - self.last_gradient
Theo Steininger's avatar
Theo Steininger committed
261

262
263
        self.last_x = x.copy()
        self.last_gradient = gradient.copy()
Theo Steininger's avatar
Theo Steininger committed
264

Martin Reinecke's avatar
fixes  
Martin Reinecke committed
265
266
        self.k += 1

Theo Steininger's avatar
Theo Steininger committed
267
268

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

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
271
272
273
274
    Parameters
    ----------
    history_length : integer
        Maximum number of stored past updates.
275

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
276
277
    Attributes
    ----------
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
278
279
280
    _pos : integer
        Number of items that have been added so far.
    _storage : list of length history_length
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
281
        List where input values are stored.
282

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
283
    """
Theo Steininger's avatar
Theo Steininger committed
284
    def __init__(self, history_length):
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
285
        self._storage = [None]*history_length
Theo Steininger's avatar
Theo Steininger committed
286
287

    def __getitem__(self, index):
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
288
        """Returns the element with index [index].
289

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
290
291
292
293
        Parameters
        ----------
        index : integer
            Index of the selected element.
294

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
295
296
297
298
        Returns
        -------
            selected element
        """
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
299
        return self._storage[index%len(self._storage)]
300

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
301
302
    def __setitem__(self, index, value):
        self._storage[index%len(self._storage)] = value