Fixed bug in XNES covariance matrix update.

pints/_optimisers/_xnes.py

@@ -14,6 +14,9 @@
 import warnings
 
 
+count = 0
+
+
 class XNES(pints.PopulationBasedOptimiser):
     """
     Finds the best parameters using the xNES method described in [1]_, [2]_.
@@ -47,7 +50,7 @@ def __init__(self, x0, sigma0=None, boundaries=None):
         self._bounded_ids = None  # Indices of those xs
 
         # Normalisation / distribution
-        self._mu = np.array(self._x0)   # Mean
+        self._mu = pints.vector(x0)     # Mean
         self._A = None                  # Covariance
 
         # Best solution seen
@@ -73,6 +76,36 @@ def ask(self):
         self._xs = np.array([self._mu + np.dot(self._A, self._zs[i])
                              for i in range(self._population_size)])
 
+        global count
+        if count > 0:
+            print('> ' * 40)
+        print('Local center', self._mu)
+        print('Local A', self._A)
+
+
+        if count == 0:
+            self._xs = np.array([
+                [ 1.82917497, -0.56082756],
+                [-0.03493942,  0.41688281],
+                [-0.4400481 ,  0.60255597],
+                [1.04353976, 1.41187322],
+                [1.43655514, 1.11339303],
+                [1.14385767, 0.58171564]])
+        else:
+            self._xs = np.array([
+                [-0.7315636 ,  0.36498419],
+                [-0.44829245, -0.59588781],
+                [ 0.79268259, -0.64598044],
+                [0.52818196, 2.17088402],
+                [-0.95416713,  0.62276164],
+                [-0.09167635,  0.81058063]])
+        count += 1
+
+        zs = self._xs - self._mu
+        Ai = np.linalg.inv(self._A)
+        self._zs = np.array([np.dot(Ai, z) for z in zs])
+
+
         # Boundaries? Then only pass user xs that are within bounds
         if self._boundaries is not None:
             self._bounded_ids = np.nonzero(
@@ -84,6 +117,9 @@ def ask(self):
         else:
             self._bounded_xs = self._xs
 
+        print('Local population', self._xs)
+        print('Local normalised', self._zs)
+
         # Set as read-only and return
         self._bounded_xs.setflags(write=False)
         return self._bounded_xs
@@ -106,13 +142,13 @@ def _initialise(self):
         d = self._n_parameters
         n = self._population_size
 
-        # Learning rates
+        # Learning rates, see Table 1
         # TODO Allow changing before run() with method call
         self._eta_mu = 1
         # TODO Allow changing before run() with method call
         self._eta_A = 0.6 * (3 + np.log(d)) * d ** -1.5
 
-        # Pre-calculated utilities
+        # Pre-calculated utilities, see Table 1
         self._us = np.maximum(0, np.log(n / 2 + 1) - np.log(1 + np.arange(n)))
         self._us /= np.sum(self._us)
         self._us -= 1 / n
@@ -147,6 +183,8 @@ def tell(self, fx):
             raise Exception('ask() not called before tell()')
         self._ready_for_tell = False
 
+        print('Local fx', fx)
+
         # Boundaries? Then reconstruct full fx vector
         if self._boundaries is not None and len(fx) < self._population_size:
             bounded_fx = fx
@@ -157,15 +195,46 @@ def tell(self, fx):
         order = np.argsort(fx)
         self._zs = self._zs[order]
 
+        '''
+        I = eye(self.numParameters)
+        self._produceSamples()
+        utilities = self.shapingFunction(self._currentEvaluations)
+        utilities /= sum(utilities)  # make the utilities sum to 1
+        if self.uniformBaseline:
+            utilities -= 1./self.batchSize
+        samples = array(list(map(self._base2sample, self._population)))
+
+        dCenter = dot(samples.T, utilities)
+        covGradient = dot(array([outer(s,s) - I for s in samples]).T, utilities)
+        covTrace = trace(covGradient)
+        covGradient -= covTrace/self.numParameters * I
+        dA = 0.5 * (self.scaleLearningRate * covTrace/self.numParameters * I
+                    +self.covLearningRate * covGradient)
+
+        self._lastLogDetA = self._logDetA
+        self._lastInvA = self._invA
+
+        self._center += self.centerLearningRate * dot(self._A, dCenter)
+        self._A = dot(self._A, expm(dA))
+        self._invA = dot(expm(-dA), self._invA)
+        self._logDetA += 0.5 * self.scaleLearningRate * covTrace
+        if self.storeAllDistributions:
+            self._allDistributions.append((self._center.copy(), self._A.copy()))
+        '''
+
         # Update center
         Gd = np.dot(self._us, self._zs)
         self._mu += self._eta_mu * np.dot(self._A, Gd)
+        print('Local dCenter', Gd)
 
         # Update root of covariance matrix
-        Gm = np.dot(
-            np.array([np.outer(z, z).T - self._I for z in self._zs]).T,
-            self._us)
-        self._A *= scipy.linalg.expm(np.dot(0.5 * self._eta_A, Gm))
+        # Note that this is equation 11 (for the eta-sigma=eta-B case), not the
+        # more general equations 9&10 version given in Algorithm 1
+        Gm = 0.5 * np.sum([
+            u * (np.outer(z, z.T) - self._I)
+            for u, z in zip(self._us, self._zs)], axis=0)
+        self._A = np.dot(self._A, scipy.linalg.expm(self._eta_A * Gm))
+        print('Local new A', self._A)
 
         # Update f_guessed on the assumption that the lowest value in our
         # sample approximates f(mu)