Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Async/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.4/CMAOperators/CMAUpdater.cs @ 14544

Last change on this file since 14544 was 12987, checked in by mkommend, 9 years ago

#2487: Extended check for degenerate states in CMAUpdater.

File size: 10.8 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 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, ISingleObjectiveOperator {
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        double[] d;
241        double[,] b;
242        var success = alglib.smatrixevd(sp.C, N, 1, true, out d, out b);
243        sp.D = d;
244        sp.B = b;
245
246        DegenerateStateParameter.ActualValue = new BoolValue(!success);
247
248        // assign D to eigenvalue square roots
249        for (int i = 0; i < N; i++) {
250          if (sp.D[i] <= 0) { // numerical problem?
251            DegenerateStateParameter.ActualValue.Value = true;
252            sp.D[i] = 0;
253          } else sp.D[i] = Math.Sqrt(sp.D[i]);
254        }
255
256        if (sp.D.Min() == 0.0) sp.AxisRatio = double.PositiveInfinity;
257        else sp.AxisRatio = sp.D.Max() / sp.D.Min();
258      }
259      return base.Apply();
260    }
261  }
262}
Note: See TracBrowser for help on using the repository browser.