source: branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/CMAOperators/CMAUpdater.cs @ 9298

Last change on this file since 9298 was 9298, checked in by abeham, 7 years ago

#1961: fixed bug introduced in last commit regarding initialization of cCov

File size: 18.1 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2013 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using HeuristicLab.Common;
23using HeuristicLab.Core;
24using HeuristicLab.Data;
25using HeuristicLab.Encodings.RealVectorEncoding;
26using HeuristicLab.Operators;
27using HeuristicLab.Optimization;
28using HeuristicLab.Parameters;
29using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
30using System;
31using System.Linq;
32
33namespace HeuristicLab.Algorithms.CMAEvolutionStrategy {
34  [Item("CMAUpdater", "Updates the covariance matrix and strategy parameters of CMA-ES.")]
35  [StorableClass]
36  public class CMAUpdater : SingleSuccessorOperator, ICMAUpdater, IIterationBasedOperator {
37
38    public Type CMAType {
39      get { return typeof(CMAParameters); }
40    }
41
42    #region Parameter Properties
43    public ILookupParameter<CMAParameters> StrategyParametersParameter {
44      get { return (ILookupParameter<CMAParameters>)Parameters["StrategyParameters"]; }
45    }
46
47    public ILookupParameter<RealVector> MeanParameter {
48      get { return (ILookupParameter<RealVector>)Parameters["Mean"]; }
49    }
50
51    public ILookupParameter<RealVector> OldMeanParameter {
52      get { return (ILookupParameter<RealVector>)Parameters["OldMean"]; }
53    }
54
55    public IScopeTreeLookupParameter<RealVector> OffspringParameter {
56      get { return (IScopeTreeLookupParameter<RealVector>)Parameters["Offspring"]; }
57    }
58
59    public IScopeTreeLookupParameter<DoubleValue> QualityParameter {
60      get { return (IScopeTreeLookupParameter<DoubleValue>)Parameters["Quality"]; }
61    }
62
63    public ILookupParameter<IntValue> IterationsParameter {
64      get { return (ILookupParameter<IntValue>)Parameters["Iterations"]; }
65    }
66
67    public IValueLookupParameter<IntValue> MaximumIterationsParameter {
68      get { return (IValueLookupParameter<IntValue>)Parameters["MaximumIterations"]; }
69    }
70
71    public IValueLookupParameter<IntValue> MaximumEvaluatedSolutionsParameter {
72      get { return (IValueLookupParameter<IntValue>)Parameters["MaximumEvaluatedSolutions"]; }
73    }
74
75    public ILookupParameter<BoolValue> DegenerateStateParameter {
76      get { return (ILookupParameter<BoolValue>)Parameters["DegenerateState"]; }
77    }
78    #endregion
79
80    [StorableConstructor]
81    protected CMAUpdater(bool deserializing) : base(deserializing) { }
82    protected CMAUpdater(CMAUpdater original, Cloner cloner) : base(original, cloner) { }
83    public CMAUpdater()
84      : base() {
85      Parameters.Add(new LookupParameter<CMAParameters>("StrategyParameters", "The strategy parameters of CMA-ES."));
86      Parameters.Add(new LookupParameter<RealVector>("Mean", "The new mean."));
87      Parameters.Add(new LookupParameter<RealVector>("OldMean", "The old mean."));
88      Parameters.Add(new ScopeTreeLookupParameter<RealVector>("Offspring", "The created offspring solutions."));
89      Parameters.Add(new ScopeTreeLookupParameter<DoubleValue>("Quality", "The quality of the offspring."));
90      Parameters.Add(new LookupParameter<IntValue>("Iterations", "The number of iterations passed."));
91      Parameters.Add(new ValueLookupParameter<IntValue>("MaximumIterations", "The maximum number of iterations."));
92      Parameters.Add(new ValueLookupParameter<IntValue>("MaximumEvaluatedSolutions", "The maximum number of evaluated solutions."));
93      Parameters.Add(new LookupParameter<BoolValue>("DegenerateState", "Whether the algorithm state has degenerated and should be terminated."));
94      MeanParameter.ActualName = "XMean";
95      OldMeanParameter.ActualName = "XOld";
96      OffspringParameter.ActualName = "RealVector";
97    }
98
99    public override IDeepCloneable Clone(Cloner cloner) {
100      return new CMAUpdater(this, cloner);
101    }
102
103    public override IOperation Apply() {
104      var iterations = IterationsParameter.ActualValue.Value;
105
106      var xold = OldMeanParameter.ActualValue;
107      var xmean = MeanParameter.ActualValue;
108      var offspring = OffspringParameter.ActualValue;
109      var quality = QualityParameter.ActualValue;
110      var lambda = offspring.Length;
111
112      var N = xmean.Length;
113      var sp = StrategyParametersParameter.ActualValue;
114
115      #region Initialize default values for strategy parameter adjustment
116      if (sp.ChiN == 0) sp.ChiN = Math.Sqrt(N) * (1.0 - 1.0 / (4.0 * N) + 1.0 / (21.0 * N * N));
117      if (sp.MuEff == 0) sp.MuEff = sp.Weights.Sum() * sp.Weights.Sum() / sp.Weights.Sum(x => x * x);
118      if (sp.CS == 0) sp.CS = (sp.MuEff + 2) / (N + sp.MuEff + 3);
119      if (sp.Damps == 0) {
120        var maxIterations = MaximumIterationsParameter.ActualValue.Value;
121        var maxEvals = MaximumEvaluatedSolutionsParameter.ActualValue.Value;
122        sp.Damps = 2 * Math.Max(0, Math.Sqrt((sp.MuEff - 1) / (N + 1)) - 1)
123                                * Math.Max(0.3, 1 - N / (1e-6 + Math.Min(maxIterations, maxEvals / lambda))) + sp.CS + 1;
124      }
125      if (sp.CC == 0) sp.CC = 4.0 / (N + 4);
126      if (sp.MuCov == 0) sp.MuCov = sp.MuEff;
127      if (sp.CCov == 0) sp.CCov = 2.0 / ((N + 1.41) * (N + 1.41) * sp.MuCov)
128                             + (1 - (1.0 / sp.MuCov)) * Math.Min(1, (2 * sp.MuEff - 1) / (sp.MuEff + (N + 2) * (N + 2)));
129      if (sp.CCovSep == 0) sp.CCovSep = Math.Min(1, sp.CCov * (N + 1.5) / 3);
130      #endregion
131
132      sp.QualityHistory.Enqueue(quality[0].Value);
133      while (sp.QualityHistory.Count > sp.QualityHistorySize && sp.QualityHistorySize >= 0)
134        sp.QualityHistory.Dequeue();
135
136      for (int i = 0; i < N; i++) {
137        sp.BDz[i] = Math.Sqrt(sp.MuEff) * (xmean[i] - xold[i]) / sp.Sigma;
138      }
139
140      if (sp.InitialIterations >= iterations) {
141        for (int i = 0; i < N; i++) {
142          sp.PS[i] = (1 - sp.CS) * sp.PS[i]
143                     + Math.Sqrt(sp.CS * (2 - sp.CS)) * sp.BDz[i] / sp.D[i];
144        }
145      } else {
146        var artmp = new double[N];
147        for (int i = 0; i < N; i++) {
148          var sum = 0.0;
149          for (int j = 0; j < N; j++) {
150            sum += sp.B[j, i] * sp.BDz[j];
151          }
152          artmp[i] = sum / sp.D[i];
153        }
154        for (int i = 0; i < N; i++) {
155          var sum = 0.0;
156          for (int j = 0; j < N; j++) {
157            sum += sp.B[i, j] * artmp[j];
158          }
159          sp.PS[i] = (1 - sp.CS) * sp.PS[i] + Math.Sqrt(sp.CS * (2 - sp.CS)) * sum;
160        }
161      }
162      var normPS = Math.Sqrt(sp.PS.Select(x => x * x).Sum());
163      var hsig = normPS / Math.Sqrt(1 - Math.Pow(1 - sp.CS, 2 * iterations)) / sp.ChiN < 1.4 + 2.0 / (N + 1) ? 1.0 : 0.0;
164      for (int i = 0; i < sp.PC.Length; i++) {
165        sp.PC[i] = (1 - sp.CC) * sp.PC[i]
166                   + hsig * Math.Sqrt(sp.CC * (2 - sp.CC)) * sp.BDz[i];
167      }
168
169      if (sp.CCov > 0) {
170        if (sp.InitialIterations >= iterations) {
171          for (int i = 0; i < N; i++) {
172            sp.C[i, i] = (1 - sp.CCovSep) * sp.C[i, i]
173                         + sp.CCov * (1 / sp.MuCov)
174                         * (sp.PC[i] * sp.PC[i] + (1 - hsig) * sp.CC * (2 - sp.CC) * sp.C[i, i]);
175            for (int k = 0; k < sp.Mu; k++) {
176              sp.C[i, i] += sp.CCov * (1 - 1 / sp.MuCov) * sp.Weights[k] * (offspring[k][i] - xold[i]) *
177                            (offspring[k][i] - xold[i]) / (sp.Sigma * sp.Sigma);
178            }
179          }
180        } else {
181          for (int i = 0; i < N; i++) {
182            for (int j = 0; j < N; j++) {
183              sp.C[i, j] = (1 - sp.CCov) * sp.C[i, j]
184                           + sp.CCov * (1 / sp.MuCov)
185                           * (sp.PC[i] * sp.PC[j] + (1 - hsig) * sp.CC * (2 - sp.CC) * sp.C[i, j]);
186              for (int k = 0; k < sp.Mu; k++) {
187                sp.C[i, j] += sp.CCov * (1 - 1 / sp.MuCov) * sp.Weights[k] * (offspring[k][i] - xold[i]) *
188                              (offspring[k][j] - xold[j]) / (sp.Sigma * sp.Sigma);
189              }
190            }
191          }
192        }
193      }
194      sp.Sigma *= Math.Exp((sp.CS / sp.Damps) * (normPS / sp.ChiN - 1));
195
196      double minSqrtdiagC = int.MaxValue, maxSqrtdiagC = int.MinValue;
197      for (int i = 0; i < N; i++) {
198        if (Math.Sqrt(sp.C[i, i]) < minSqrtdiagC) minSqrtdiagC = Math.Sqrt(sp.C[i, i]);
199        if (Math.Sqrt(sp.C[i, i]) > maxSqrtdiagC) maxSqrtdiagC = Math.Sqrt(sp.C[i, i]);
200      }
201
202      // ensure maximal and minimal standard deviations
203      if (sp.SigmaBounds != null && sp.SigmaBounds.GetLength(0) > 0) {
204        for (int i = 0; i < N; i++) {
205          var d = sp.SigmaBounds[Math.Min(i, sp.SigmaBounds.GetLength(0) - 1), 0];
206          if (d > sp.Sigma * minSqrtdiagC) sp.Sigma = d / minSqrtdiagC;
207        }
208        for (int i = 0; i < N; i++) {
209          var d = sp.SigmaBounds[Math.Min(i, sp.SigmaBounds.GetLength(0) - 1), 1];
210          if (d > sp.Sigma * maxSqrtdiagC) sp.Sigma = d / maxSqrtdiagC;
211        }
212      }
213      // end ensure ...
214
215      // testAndCorrectNumerics
216      double fac = 1;
217      if (sp.D.Max() < 1e-6)
218        fac = 1.0 / sp.D.Max();
219      else if (sp.D.Min() > 1e4)
220        fac = 1.0 / sp.D.Min();
221
222      if (fac != 1.0) {
223        sp.Sigma /= fac;
224        for (int i = 0; i < N; i++) {
225          sp.PC[i] *= fac;
226          sp.D[i] *= fac;
227          for (int j = 0; j < N; j++)
228            sp.C[i, j] *= fac * fac;
229        }
230      }
231      // end testAndCorrectNumerics
232
233
234      if (sp.InitialIterations >= iterations) {
235        for (int i = 0; i < N; i++)
236          sp.D[i] = Math.Sqrt(sp.C[i, i]);
237        DegenerateStateParameter.ActualValue = new BoolValue(false);
238      } else {
239
240        // set B <- C
241        for (int i = 0; i < N; i++) {
242          for (int j = 0; j < N; j++) {
243            sp.B[i, j] = sp.C[i, j];
244          }
245        }
246        /*double[] eVal;
247        double[,] eVec;
248        var success = alglib.smatrixevd(sp.B, N, 1, false, out eVal, out eVec);
249        sp.B = eVec;
250        sp.D = eVal;*/
251        var success = Eigendecomposition(N, sp.B, sp.D);
252
253        DegenerateStateParameter.ActualValue = new BoolValue(!success);
254
255        // assign D to eigenvalue square roots
256        for (int i = 0; i < N; i++) {
257          if (sp.D[i] < 0) { // numerical problem?
258            DegenerateStateParameter.ActualValue.Value = true;
259            sp.D[i] = 0;
260          } else sp.D[i] = Math.Sqrt(sp.D[i]);
261        }
262
263        if (sp.D.Min() == 0.0) sp.AxisRatio = double.PositiveInfinity;
264        else sp.AxisRatio = sp.D.Max() / sp.D.Min();
265      }
266      return base.Apply();
267    }
268
269    private bool Eigendecomposition(int N, double[,] B, double[] diagD) {
270      bool result = true;
271      // eigendecomposition
272      var offdiag = new double[N];
273      try {
274        tred2(N, B, diagD, offdiag);
275        tql2(N, diagD, offdiag, B);
276      } catch { result = false; }
277
278      return result;
279    } // eigendecomposition
280
281
282    // Symmetric Householder reduction to tridiagonal form, taken from JAMA package.
283    private void tred2(int n, double[,] V, double[] d, double[] e) {
284
285      //  This is derived from the Algol procedures tred2 by
286      //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
287      //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
288      //  Fortran subroutine in EISPACK.
289
290      for (int j = 0; j < n; j++) {
291        d[j] = V[n - 1, j];
292      }
293
294      // Householder reduction to tridiagonal form.
295
296      for (int i = n - 1; i > 0; i--) {
297
298        // Scale to avoid under/overflow.
299
300        double scale = 0.0;
301        double h = 0.0;
302        for (int k = 0; k < i; k++) {
303          scale = scale + Math.Abs(d[k]);
304        }
305        if (scale == 0.0) {
306          e[i] = d[i - 1];
307          for (int j = 0; j < i; j++) {
308            d[j] = V[i - 1, j];
309            V[i, j] = 0.0;
310            V[j, i] = 0.0;
311          }
312        } else {
313
314          // Generate Householder vector.
315
316          for (int k = 0; k < i; k++) {
317            d[k] /= scale;
318            h += d[k] * d[k];
319          }
320          double f = d[i - 1];
321          double g = Math.Sqrt(h);
322          if (f > 0) {
323            g = -g;
324          }
325          e[i] = scale * g;
326          h = h - f * g;
327          d[i - 1] = f - g;
328          for (int j = 0; j < i; j++) {
329            e[j] = 0.0;
330          }
331
332          // Apply similarity transformation to remaining columns.
333
334          for (int j = 0; j < i; j++) {
335            f = d[j];
336            V[j, i] = f;
337            g = e[j] + V[j, j] * f;
338            for (int k = j + 1; k <= i - 1; k++) {
339              g += V[k, j] * d[k];
340              e[k] += V[k, j] * f;
341            }
342            e[j] = g;
343          }
344          f = 0.0;
345          for (int j = 0; j < i; j++) {
346            e[j] /= h;
347            f += e[j] * d[j];
348          }
349          double hh = f / (h + h);
350          for (int j = 0; j < i; j++) {
351            e[j] -= hh * d[j];
352          }
353          for (int j = 0; j < i; j++) {
354            f = d[j];
355            g = e[j];
356            for (int k = j; k <= i - 1; k++) {
357              V[k, j] -= (f * e[k] + g * d[k]);
358            }
359            d[j] = V[i - 1, j];
360            V[i, j] = 0.0;
361          }
362        }
363        d[i] = h;
364      }
365
366      // Accumulate transformations.
367
368      for (int i = 0; i < n - 1; i++) {
369        V[n - 1, i] = V[i, i];
370        V[i, i] = 1.0;
371        double h = d[i + 1];
372        if (h != 0.0) {
373          for (int k = 0; k <= i; k++) {
374            d[k] = V[k, i + 1] / h;
375          }
376          for (int j = 0; j <= i; j++) {
377            double g = 0.0;
378            for (int k = 0; k <= i; k++) {
379              g += V[k, i + 1] * V[k, j];
380            }
381            for (int k = 0; k <= i; k++) {
382              V[k, j] -= g * d[k];
383            }
384          }
385        }
386        for (int k = 0; k <= i; k++) {
387          V[k, i + 1] = 0.0;
388        }
389      }
390      for (int j = 0; j < n; j++) {
391        d[j] = V[n - 1, j];
392        V[n - 1, j] = 0.0;
393      }
394      V[n - 1, n - 1] = 1.0;
395      e[0] = 0.0;
396    }
397
398    // Symmetric tridiagonal QL algorithm, taken from JAMA package.
399    private void tql2(int n, double[] d, double[] e, double[,] V) {
400
401      //  This is derived from the Algol procedures tql2, by
402      //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
403      //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
404      //  Fortran subroutine in EISPACK.
405
406      for (int i = 1; i < n; i++) {
407        e[i - 1] = e[i];
408      }
409      e[n - 1] = 0.0;
410
411      double f = 0.0;
412      double tst1 = 0.0;
413      double eps = Math.Pow(2.0, -52.0);
414      for (int l = 0; l < n; l++) {
415
416        // Find small subdiagonal element
417
418        tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l]));
419        int m = l;
420        while (m < n) {
421          if (Math.Abs(e[m]) <= eps * tst1) {
422            break;
423          }
424          m++;
425        }
426
427        // If m == l, d[l] is an eigenvalue,
428        // otherwise, iterate.
429
430        if (m > l) {
431          int iter = 0;
432          do {
433            iter = iter + 1;  // (Could check iteration count here.)
434
435            // Compute implicit shift
436
437            double g = d[l];
438            double p = (d[l + 1] - g) / (2.0 * e[l]);
439            double r = hypot(p, 1.0);
440            if (p < 0) {
441              r = -r;
442            }
443            d[l] = e[l] / (p + r);
444            d[l + 1] = e[l] * (p + r);
445            double dl1 = d[l + 1];
446            double h = g - d[l];
447            for (int i = l + 2; i < n; i++) {
448              d[i] -= h;
449            }
450            f = f + h;
451
452            // Implicit QL transformation.
453
454            p = d[m];
455            double c = 1.0;
456            double c2 = c;
457            double c3 = c;
458            double el1 = e[l + 1];
459            double s = 0.0;
460            double s2 = 0.0;
461            for (int i = m - 1; i >= l; i--) {
462              c3 = c2;
463              c2 = c;
464              s2 = s;
465              g = c * e[i];
466              h = c * p;
467              r = hypot(p, e[i]);
468              e[i + 1] = s * r;
469              s = e[i] / r;
470              c = p / r;
471              p = c * d[i] - s * g;
472              d[i + 1] = h + s * (c * g + s * d[i]);
473
474              // Accumulate transformation.
475
476              for (int k = 0; k < n; k++) {
477                h = V[k, i + 1];
478                V[k, i + 1] = s * V[k, i] + c * h;
479                V[k, i] = c * V[k, i] - s * h;
480              }
481            }
482            p = -s * s2 * c3 * el1 * e[l] / dl1;
483            e[l] = s * p;
484            d[l] = c * p;
485
486            // Check for convergence.
487
488          } while (Math.Abs(e[l]) > eps * tst1);
489        }
490        d[l] = d[l] + f;
491        e[l] = 0.0;
492      }
493
494      // Sort eigenvalues and corresponding vectors.
495
496      for (int i = 0; i < n - 1; i++) {
497        int k = i;
498        double p = d[i];
499        for (int j = i + 1; j < n; j++) {
500          if (d[j] < p) { // NH find smallest k>i
501            k = j;
502            p = d[j];
503          }
504        }
505        if (k != i) {
506          d[k] = d[i]; // swap k and i
507          d[i] = p;
508          for (int j = 0; j < n; j++) {
509            p = V[j, i];
510            V[j, i] = V[j, k];
511            V[j, k] = p;
512          }
513        }
514      }
515    }
516
517    /** sqrt(a^2 + b^2) without under/overflow. **/
518    private double hypot(double a, double b) {
519      double r = 0;
520      if (Math.Abs(a) > Math.Abs(b)) {
521        r = b / a;
522        r = Math.Abs(a) * Math.Sqrt(1 + r * r);
523      } else if (b != 0) {
524        r = a / b;
525        r = Math.Abs(b) * Math.Sqrt(1 + r * r);
526      }
527      return r;
528    }
529  }
530}
Note: See TracBrowser for help on using the repository browser.