Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessCovarianceOptimizationProblem.cs @ 13201

Last change on this file since 13201 was 13201, checked in by gkronber, 8 years ago

#1967: create a copy of the hyper-parameter array

File size: 18.6 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 System;
23using System.Linq;
24using HeuristicLab.Common;
25using HeuristicLab.Core;
26using HeuristicLab.Data;
27using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
28using HeuristicLab.Optimization;
29using HeuristicLab.Parameters;
30using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
31using HeuristicLab.Problems.DataAnalysis;
32using HeuristicLab.Problems.Instances;
33
34
35namespace HeuristicLab.Algorithms.DataAnalysis {
36  [Item("Gaussian Process Covariance Optimization Problem", "")]
37  [Creatable(CreatableAttribute.Categories.GeneticProgrammingProblems, Priority = 300)]
38  [StorableClass]
39  public sealed class GaussianProcessCovarianceOptimizationProblem : SymbolicExpressionTreeProblem, IRegressionProblem, IProblemInstanceConsumer<IRegressionProblemData>, IProblemInstanceExporter<IRegressionProblemData> {
40    #region static variables and ctor
41    private static readonly CovarianceMaternIso maternIso1;
42    private static readonly CovarianceMaternIso maternIso3;
43    private static readonly CovarianceMaternIso maternIso5;
44    private static readonly CovariancePiecewisePolynomial piecewisePoly0;
45    private static readonly CovariancePiecewisePolynomial piecewisePoly1;
46    private static readonly CovariancePiecewisePolynomial piecewisePoly2;
47    private static readonly CovariancePiecewisePolynomial piecewisePoly3;
48    private static readonly CovariancePolynomial poly2;
49    private static readonly CovariancePolynomial poly3;
50    private static readonly CovarianceSpectralMixture spectralMixture1;
51    private static readonly CovarianceSpectralMixture spectralMixture3;
52    private static readonly CovarianceSpectralMixture spectralMixture5;
53    private static readonly CovarianceLinear linear;
54    private static readonly CovarianceLinearArd linearArd;
55    private static readonly CovarianceNeuralNetwork neuralNetwork;
56    private static readonly CovariancePeriodic periodic;
57    private static readonly CovarianceRationalQuadraticIso ratQuadraticIso;
58    private static readonly CovarianceRationalQuadraticArd ratQuadraticArd;
59    private static readonly CovarianceSquaredExponentialArd sqrExpArd;
60    private static readonly CovarianceSquaredExponentialIso sqrExpIso;
61
62    static GaussianProcessCovarianceOptimizationProblem() {
63      // cumbersome initialization because of ConstrainedValueParameters
64      maternIso1 = new CovarianceMaternIso(); SetConstrainedValueParameter(maternIso1.DParameter, 1);
65      maternIso3 = new CovarianceMaternIso(); SetConstrainedValueParameter(maternIso3.DParameter, 3);
66      maternIso5 = new CovarianceMaternIso(); SetConstrainedValueParameter(maternIso5.DParameter, 5);
67
68      piecewisePoly0 = new CovariancePiecewisePolynomial(); SetConstrainedValueParameter(piecewisePoly0.VParameter, 0);
69      piecewisePoly1 = new CovariancePiecewisePolynomial(); SetConstrainedValueParameter(piecewisePoly1.VParameter, 1);
70      piecewisePoly2 = new CovariancePiecewisePolynomial(); SetConstrainedValueParameter(piecewisePoly2.VParameter, 2);
71      piecewisePoly3 = new CovariancePiecewisePolynomial(); SetConstrainedValueParameter(piecewisePoly3.VParameter, 3);
72
73      poly2 = new CovariancePolynomial(); poly2.DegreeParameter.Value.Value = 2;
74      poly3 = new CovariancePolynomial(); poly3.DegreeParameter.Value.Value = 3;
75
76      spectralMixture1 = new CovarianceSpectralMixture(); spectralMixture1.QParameter.Value.Value = 1;
77      spectralMixture3 = new CovarianceSpectralMixture(); spectralMixture3.QParameter.Value.Value = 3;
78      spectralMixture5 = new CovarianceSpectralMixture(); spectralMixture5.QParameter.Value.Value = 5;
79
80      linear = new CovarianceLinear();
81      linearArd = new CovarianceLinearArd();
82      neuralNetwork = new CovarianceNeuralNetwork();
83      periodic = new CovariancePeriodic();
84      ratQuadraticArd = new CovarianceRationalQuadraticArd();
85      ratQuadraticIso = new CovarianceRationalQuadraticIso();
86      sqrExpArd = new CovarianceSquaredExponentialArd();
87      sqrExpIso = new CovarianceSquaredExponentialIso();
88    }
89
90    private static void SetConstrainedValueParameter(IConstrainedValueParameter<IntValue> param, int val) {
91      param.Value = param.ValidValues.Single(v => v.Value == val);
92    }
93
94    #endregion
95
96    #region parameter names
97
98    private const string ProblemDataParameterName = "ProblemData";
99    private const string ConstantOptIterationsParameterName = "Constant optimization steps";
100    private const string RestartsParameterName = "Restarts";
101
102    #endregion
103
104    #region Parameter Properties
105    IParameter IDataAnalysisProblem.ProblemDataParameter { get { return ProblemDataParameter; } }
106
107    public IValueParameter<IRegressionProblemData> ProblemDataParameter {
108      get { return (IValueParameter<IRegressionProblemData>)Parameters[ProblemDataParameterName]; }
109    }
110    public IFixedValueParameter<IntValue> ConstantOptIterationsParameter {
111      get { return (IFixedValueParameter<IntValue>)Parameters[ConstantOptIterationsParameterName]; }
112    }
113    public IFixedValueParameter<IntValue> RestartsParameter {
114      get { return (IFixedValueParameter<IntValue>)Parameters[RestartsParameterName]; }
115    }
116    #endregion
117
118    #region Properties
119
120    public IRegressionProblemData ProblemData {
121      get { return ProblemDataParameter.Value; }
122      set { ProblemDataParameter.Value = value; }
123    }
124    IDataAnalysisProblemData IDataAnalysisProblem.ProblemData { get { return ProblemData; } }
125
126    public int ConstantOptIterations {
127      get { return ConstantOptIterationsParameter.Value.Value; }
128      set { ConstantOptIterationsParameter.Value.Value = value; }
129    }
130
131    public int Restarts {
132      get { return RestartsParameter.Value.Value; }
133      set { RestartsParameter.Value.Value = value; }
134    }
135    #endregion
136
137    public override bool Maximization {
138      get { return true; } // return log likelihood (instead of negative log likelihood as in GPR
139    }
140
141    // problem stores a few variables for information exchange from Evaluate() to Analyze()
142    private object problemStateLocker = new object();
143    [Storable]
144    private double bestQ;
145    [Storable]
146    private double[] bestHyperParameters;
147    [Storable]
148    private IMeanFunction meanFunc;
149    [Storable]
150    private ICovarianceFunction covFunc;
151
152    public GaussianProcessCovarianceOptimizationProblem()
153      : base() {
154      Parameters.Add(new ValueParameter<IRegressionProblemData>(ProblemDataParameterName, "The data for the regression problem", new RegressionProblemData()));
155      Parameters.Add(new FixedValueParameter<IntValue>(ConstantOptIterationsParameterName, "Number of optimization steps for hyperparameter values", new IntValue(50)));
156      Parameters.Add(new FixedValueParameter<IntValue>(RestartsParameterName, "The number of random restarts for constant optimization.", new IntValue(10)));
157      Parameters["Restarts"].Hidden = true;
158      var g = new SimpleSymbolicExpressionGrammar();
159      g.AddSymbols(new string[] { "Sum", "Product" }, 2, 2);
160      g.AddTerminalSymbols(new string[]
161      {
162        "Linear",
163        "LinearArd",
164        "MaternIso1",
165        "MaternIso3",
166        "MaternIso5",
167        "NeuralNetwork",
168        "Periodic",
169        "PiecewisePolynomial0",
170        "PiecewisePolynomial1",
171        "PiecewisePolynomial2",
172        "PiecewisePolynomial3",
173        "Polynomial2",
174        "Polynomial3",
175        "RationalQuadraticArd",
176        "RationalQuadraticIso",
177        "SpectralMixture1",
178        "SpectralMixture3",
179        "SpectralMixture5",
180        "SquaredExponentialArd",
181        "SquaredExponentialIso"
182      });
183      base.Encoding = new SymbolicExpressionTreeEncoding(g, 10, 5);
184    }
185
186    protected override void OnReset() {
187      base.OnReset();
188      meanFunc = null;
189      covFunc = null;
190      bestQ = double.NegativeInfinity;
191      bestHyperParameters = null;
192    }
193
194    public override double Evaluate(ISymbolicExpressionTree tree, IRandom random) {
195
196      var meanFunction = new MeanConst();
197      var problemData = ProblemData;
198      var ds = problemData.Dataset;
199      var targetVariable = problemData.TargetVariable;
200      var allowedInputVariables = problemData.AllowedInputVariables.ToArray();
201      var nVars = allowedInputVariables.Length;
202      var trainingRows = problemData.TrainingIndices.ToArray();
203
204      // use the same covariance function for each restart
205      var covarianceFunction = TreeToCovarianceFunction(tree);
206
207      // allocate hyperparameters
208      var hyperParameters = new double[meanFunction.GetNumberOfParameters(nVars) + covarianceFunction.GetNumberOfParameters(nVars) + 1]; // mean + cov + noise
209      double[] bestHyperParameters = new double[hyperParameters.Length];
210      var bestObjValue = new double[1] { double.MinValue };
211
212      // data that is necessary for the objective function
213      var data = Tuple.Create(ds, targetVariable, allowedInputVariables, trainingRows, (IMeanFunction)meanFunction, covarianceFunction, bestObjValue);
214
215      for (int t = 0; t < Restarts; t++) {
216        var prevBest = bestObjValue[0];
217        var prevBestHyperParameters = new double[hyperParameters.Length];
218        Array.Copy(bestHyperParameters, prevBestHyperParameters, bestHyperParameters.Length);
219
220        // initialize hyperparameters
221        hyperParameters[0] = ds.GetDoubleValues(targetVariable).Average(); // mean const
222
223        for (int i = 0; i < covarianceFunction.GetNumberOfParameters(nVars); i++) {
224          hyperParameters[1 + i] = random.NextDouble() * 2.0 - 1.0;
225        }
226        hyperParameters[hyperParameters.Length - 1] = 1.0; // s² = exp(2), TODO: other inits better?
227
228        // use alglib.bfgs for hyper-parameter optimization ...
229        double epsg = 0;
230        double epsf = 0.00001;
231        double epsx = 0;
232        double stpmax = 1;
233        int maxits = ConstantOptIterations;
234        alglib.mincgstate state;
235        alglib.mincgreport rep;
236
237        alglib.mincgcreate(hyperParameters, out state);
238        alglib.mincgsetcond(state, epsg, epsf, epsx, maxits);
239        alglib.mincgsetstpmax(state, stpmax);
240        alglib.mincgoptimize(state, ObjectiveFunction, null, data);
241
242        alglib.mincgresults(state, out bestHyperParameters, out rep);
243
244        if (rep.terminationtype < 0) {
245          // error -> restore previous best quality
246          bestObjValue[0] = prevBest;
247          Array.Copy(prevBestHyperParameters, bestHyperParameters, prevBestHyperParameters.Length);
248        }
249      }
250
251      UpdateBestSoFar(bestObjValue[0], bestHyperParameters, meanFunction, covarianceFunction);
252
253      return bestObjValue[0];
254    }
255
256    // updates the overall best quality and overall best model for Analyze()
257    private void UpdateBestSoFar(double bestQ, double[] bestHyperParameters, IMeanFunction meanFunc, ICovarianceFunction covFunc) {
258      lock (problemStateLocker) {
259        if (bestQ > this.bestQ) {
260          this.bestQ = bestQ;
261          this.bestHyperParameters = new double[bestHyperParameters.Length];
262          Array.Copy(bestHyperParameters, this.bestHyperParameters, this.bestHyperParameters.Length);
263          this.meanFunc = meanFunc;
264          this.covFunc = covFunc;
265        }
266      }
267    }
268
269    public override void Analyze(ISymbolicExpressionTree[] trees, double[] qualities, ResultCollection results, IRandom random) {
270      if (!results.ContainsKey("Best Solution Quality")) {
271        results.Add(new Result("Best Solution Quality", typeof(DoubleValue)));
272      }
273      if (!results.ContainsKey("Best Tree")) {
274        results.Add(new Result("Best Tree", typeof(ISymbolicExpressionTree)));
275      }
276      if (!results.ContainsKey("Best Solution")) {
277        results.Add(new Result("Best Solution", typeof(GaussianProcessRegressionSolution)));
278      }
279
280      var bestQuality = qualities.Max();
281
282      if (results["Best Solution Quality"].Value == null || bestQuality > ((DoubleValue)results["Best Solution Quality"].Value).Value) {
283        var bestIdx = Array.IndexOf(qualities, bestQuality);
284        var bestClone = (ISymbolicExpressionTree)trees[bestIdx].Clone();
285        results["Best Tree"].Value = bestClone;
286        results["Best Solution Quality"].Value = new DoubleValue(bestQuality);
287        results["Best Solution"].Value = CreateSolution();
288      }
289    }
290
291    private IItem CreateSolution() {
292      var problemData = ProblemData;
293      var ds = problemData.Dataset;
294      var targetVariable = problemData.TargetVariable;
295      var allowedInputVariables = problemData.AllowedInputVariables.ToArray();
296      var trainingRows = problemData.TrainingIndices.ToArray();
297
298      lock (problemStateLocker) {
299        var model = new GaussianProcessModel(ds, targetVariable, allowedInputVariables, trainingRows, bestHyperParameters, (IMeanFunction)meanFunc.Clone(), (ICovarianceFunction)covFunc.Clone());
300        model.FixParameters();
301        return model.CreateRegressionSolution((IRegressionProblemData)ProblemData.Clone());
302      }
303    }
304
305    private void ObjectiveFunction(double[] x, ref double func, double[] grad, object obj) {
306      // we want to optimize the model likelihood by changing the hyperparameters and also return the gradient for each hyperparameter
307      var data = (Tuple<IDataset, string, string[], int[], IMeanFunction, ICovarianceFunction, double[]>)obj;
308      var ds = data.Item1;
309      var targetVariable = data.Item2;
310      var allowedInputVariables = data.Item3;
311      var trainingRows = data.Item4;
312      var meanFunction = data.Item5;
313      var covarianceFunction = data.Item6;
314      var bestObjValue = data.Item7;
315      var hyperParameters = x; // the decision variable vector
316
317      try {
318        var model = new GaussianProcessModel(ds, targetVariable, allowedInputVariables, trainingRows, hyperParameters, meanFunction, covarianceFunction);
319
320        func = model.NegativeLogLikelihood; // mincgoptimize, so we return negative likelihood
321        bestObjValue[0] = Math.Max(bestObjValue[0], -func); // problem itself is a maximization problem
322        var gradients = model.HyperparameterGradients;
323        Array.Copy(gradients, grad, gradients.Length);
324      } catch (ArgumentException) {
325        // building the GaussianProcessModel might fail, in this case we return the worst possible objective value
326        func = 1.0E+300;
327        Array.Clear(grad, 0, grad.Length);
328      }
329    }
330
331    private ICovarianceFunction TreeToCovarianceFunction(ISymbolicExpressionTree tree) {
332      return TreeToCovarianceFunction(tree.Root.GetSubtree(0).GetSubtree(0)); // skip programroot and startsymbol
333    }
334
335    private ICovarianceFunction TreeToCovarianceFunction(ISymbolicExpressionTreeNode node) {
336      switch (node.Symbol.Name) {
337        case "Sum": {
338            var sum = new CovarianceSum();
339            sum.Terms.Add(TreeToCovarianceFunction(node.GetSubtree(0)));
340            sum.Terms.Add(TreeToCovarianceFunction(node.GetSubtree(1)));
341            return sum;
342          }
343        case "Product": {
344            var prod = new CovarianceProduct();
345            prod.Factors.Add(TreeToCovarianceFunction(node.GetSubtree(0)));
346            prod.Factors.Add(TreeToCovarianceFunction(node.GetSubtree(1)));
347            return prod;
348          }
349        // covFunction is cloned by the model so we can reuse instances of terminal covariance functions
350        case "Linear": return linear;
351        case "LinearArd": return linearArd;
352        case "MaternIso1": return maternIso1;
353        case "MaternIso3": return maternIso3;
354        case "MaternIso5": return maternIso5;
355        case "NeuralNetwork": return neuralNetwork;
356        case "Periodic": return periodic;
357        case "PiecewisePolynomial0": return piecewisePoly0;
358        case "PiecewisePolynomial1": return piecewisePoly1;
359        case "PiecewisePolynomial2": return piecewisePoly2;
360        case "PiecewisePolynomial3": return piecewisePoly3;
361        case "Polynomial2": return poly2;
362        case "Polynomial3": return poly3;
363        case "RationalQuadraticArd": return ratQuadraticArd;
364        case "RationalQuadraticIso": return ratQuadraticIso;
365        case "SpectralMixture1": return spectralMixture1;
366        case "SpectralMixture3": return spectralMixture3;
367        case "SpectralMixture5": return spectralMixture5;
368        case "SquaredExponentialArd": return sqrExpArd;
369        case "SquaredExponentialIso": return sqrExpIso;
370        default: throw new InvalidProgramException(string.Format("Found invalid symbol {0}", node.Symbol.Name));
371      }
372    }
373
374
375    // persistence
376    [StorableConstructor]
377    private GaussianProcessCovarianceOptimizationProblem(bool deserializing) : base(deserializing) { }
378    [StorableHook(HookType.AfterDeserialization)]
379    private void AfterDeserialization() {
380    }
381
382    // cloning
383    private GaussianProcessCovarianceOptimizationProblem(GaussianProcessCovarianceOptimizationProblem original, Cloner cloner)
384      : base(original, cloner) {
385      bestQ = original.bestQ;
386      meanFunc = cloner.Clone(original.meanFunc);
387      covFunc = cloner.Clone(original.covFunc);
388      if (bestHyperParameters != null) {
389        bestHyperParameters = new double[original.bestHyperParameters.Length];
390        Array.Copy(original.bestHyperParameters, bestHyperParameters, bestHyperParameters.Length);
391      }
392    }
393    public override IDeepCloneable Clone(Cloner cloner) {
394      return new GaussianProcessCovarianceOptimizationProblem(this, cloner);
395    }
396
397    public void Load(IRegressionProblemData data) {
398      this.ProblemData = data;
399      OnProblemDataChanged();
400    }
401
402    public IRegressionProblemData Export() {
403      return ProblemData;
404    }
405
406    #region events
407    public event EventHandler ProblemDataChanged;
408
409
410    private void OnProblemDataChanged() {
411      var handler = ProblemDataChanged;
412      if (handler != null)
413        handler(this, EventArgs.Empty);
414    }
415    #endregion
416
417  }
418}
Note: See TracBrowser for help on using the repository browser.