Free cookie consent management tool by TermsFeed Policy Generator

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

Last change on this file since 13200 was 13200, checked in by gkronber, 7 years ago

#1967: also added the Gaussian process solution as a result

File size: 18.4 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 = bestHyperParameters;
262          this.meanFunc = meanFunc;
263          this.covFunc = covFunc;
264        }
265      }
266    }
267
268    public override void Analyze(ISymbolicExpressionTree[] trees, double[] qualities, ResultCollection results, IRandom random) {
269      if (!results.ContainsKey("Best Solution Quality")) {
270        results.Add(new Result("Best Solution Quality", typeof(DoubleValue)));
271      }
272      if (!results.ContainsKey("Best Tree")) {
273        results.Add(new Result("Best Tree", typeof(ISymbolicExpressionTree)));
274      }
275      if (!results.ContainsKey("Best Solution")) {
276        results.Add(new Result("Best Solution", typeof(GaussianProcessRegressionSolution)));
277      }
278
279      var bestQuality = qualities.Max();
280
281      if (results["Best Solution Quality"].Value == null || bestQuality > ((DoubleValue)results["Best Solution Quality"].Value).Value) {
282        var bestIdx = Array.IndexOf(qualities, bestQuality);
283        var bestClone = (ISymbolicExpressionTree)trees[bestIdx].Clone();
284        results["Best Tree"].Value = bestClone;
285        results["Best Solution Quality"].Value = new DoubleValue(bestQuality);
286        results["Best Solution"].Value = CreateSolution();
287      }
288    }
289
290    private IItem CreateSolution() {
291      var problemData = ProblemData;
292      var ds = problemData.Dataset;
293      var targetVariable = problemData.TargetVariable;
294      var allowedInputVariables = problemData.AllowedInputVariables.ToArray();
295      var trainingRows = problemData.TrainingIndices.ToArray();
296
297      lock (problemStateLocker) {
298        var model = new GaussianProcessModel(ds, targetVariable, allowedInputVariables, trainingRows, bestHyperParameters, (IMeanFunction)meanFunc.Clone(), (ICovarianceFunction)covFunc.Clone());
299        model.FixParameters();
300        return model.CreateRegressionSolution((IRegressionProblemData)ProblemData.Clone());
301      }
302    }
303
304    private void ObjectiveFunction(double[] x, ref double func, double[] grad, object obj) {
305      // we want to optimize the model likelihood by changing the hyperparameters and also return the gradient for each hyperparameter
306      var data = (Tuple<IDataset, string, string[], int[], IMeanFunction, ICovarianceFunction, double[]>)obj;
307      var ds = data.Item1;
308      var targetVariable = data.Item2;
309      var allowedInputVariables = data.Item3;
310      var trainingRows = data.Item4;
311      var meanFunction = data.Item5;
312      var covarianceFunction = data.Item6;
313      var bestObjValue = data.Item7;
314      var hyperParameters = x; // the decision variable vector
315
316      try {
317        var model = new GaussianProcessModel(ds, targetVariable, allowedInputVariables, trainingRows, hyperParameters, meanFunction, covarianceFunction);
318
319        func = model.NegativeLogLikelihood; // mincgoptimize, so we return negative likelihood
320        bestObjValue[0] = Math.Max(bestObjValue[0], -func); // problem itself is a maximization problem
321        var gradients = model.HyperparameterGradients;
322        Array.Copy(gradients, grad, gradients.Length);
323      } catch (ArgumentException) {
324        // building the GaussianProcessModel might fail, in this case we return the worst possible objective value
325        func = 1.0E+300;
326        Array.Clear(grad, 0, grad.Length);
327      }
328    }
329
330    private ICovarianceFunction TreeToCovarianceFunction(ISymbolicExpressionTree tree) {
331      return TreeToCovarianceFunction(tree.Root.GetSubtree(0).GetSubtree(0)); // skip programroot and startsymbol
332    }
333
334    private ICovarianceFunction TreeToCovarianceFunction(ISymbolicExpressionTreeNode node) {
335      switch (node.Symbol.Name) {
336        case "Sum": {
337            var sum = new CovarianceSum();
338            sum.Terms.Add(TreeToCovarianceFunction(node.GetSubtree(0)));
339            sum.Terms.Add(TreeToCovarianceFunction(node.GetSubtree(1)));
340            return sum;
341          }
342        case "Product": {
343            var prod = new CovarianceProduct();
344            prod.Factors.Add(TreeToCovarianceFunction(node.GetSubtree(0)));
345            prod.Factors.Add(TreeToCovarianceFunction(node.GetSubtree(1)));
346            return prod;
347          }
348        // covFunction is cloned by the model so we can reuse instances of terminal covariance functions
349        case "Linear": return linear;
350        case "LinearArd": return linearArd;
351        case "MaternIso1": return maternIso1;
352        case "MaternIso3": return maternIso3;
353        case "MaternIso5": return maternIso5;
354        case "NeuralNetwork": return neuralNetwork;
355        case "Periodic": return periodic;
356        case "PiecewisePolynomial0": return piecewisePoly0;
357        case "PiecewisePolynomial1": return piecewisePoly1;
358        case "PiecewisePolynomial2": return piecewisePoly2;
359        case "PiecewisePolynomial3": return piecewisePoly3;
360        case "Polynomial2": return poly2;
361        case "Polynomial3": return poly3;
362        case "RationalQuadraticArd": return ratQuadraticArd;
363        case "RationalQuadraticIso": return ratQuadraticIso;
364        case "SpectralMixture1": return spectralMixture1;
365        case "SpectralMixture3": return spectralMixture3;
366        case "SpectralMixture5": return spectralMixture5;
367        case "SquaredExponentialArd": return sqrExpArd;
368        case "SquaredExponentialIso": return sqrExpIso;
369        default: throw new InvalidProgramException(string.Format("Found invalid symbol {0}", node.Symbol.Name));
370      }
371    }
372
373
374    // persistence
375    [StorableConstructor]
376    private GaussianProcessCovarianceOptimizationProblem(bool deserializing) : base(deserializing) { }
377    [StorableHook(HookType.AfterDeserialization)]
378    private void AfterDeserialization() {
379    }
380
381    // cloning
382    private GaussianProcessCovarianceOptimizationProblem(GaussianProcessCovarianceOptimizationProblem original, Cloner cloner)
383      : base(original, cloner) {
384      bestQ = original.bestQ;
385      meanFunc = cloner.Clone(original.meanFunc);
386      covFunc = cloner.Clone(original.covFunc);
387      if (bestHyperParameters != null) {
388        bestHyperParameters = new double[original.bestHyperParameters.Length];
389        Array.Copy(original.bestHyperParameters, bestHyperParameters, bestHyperParameters.Length);
390      }
391    }
392    public override IDeepCloneable Clone(Cloner cloner) {
393      return new GaussianProcessCovarianceOptimizationProblem(this, cloner);
394    }
395
396    public void Load(IRegressionProblemData data) {
397      this.ProblemData = data;
398      OnProblemDataChanged();
399    }
400
401    public IRegressionProblemData Export() {
402      return ProblemData;
403    }
404
405    #region events
406    public event EventHandler ProblemDataChanged;
407
408
409    private void OnProblemDataChanged() {
410      var handler = ProblemDataChanged;
411      if (handler != null)
412        handler(this, EventArgs.Empty);
413    }
414    #endregion
415
416  }
417}
Note: See TracBrowser for help on using the repository browser.