Changeset 13287 for stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess

Timestamp:

11/19/15 11:59:51 (9 years ago)

Author:

gkronber

Message:

#1967: partial merge of r13160 from trunk to stable

Location:

stable/HeuristicLab.Algorithms.DataAnalysis

Files:

• . (modified) (1 prop)
• 3.4/GaussianProcess/GaussianProcessCovarianceOptimizationProblem.cs (modified) (3 diffs)
• 3.4/GaussianProcess/GaussianProcessModel.cs (modified) (2 diffs)

stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessCovarianceOptimizationProblem.cs

@@ -26 +26 @@
 using HeuristicLab.Data;
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
+using HeuristicLab.Optimization;
 using HeuristicLab.Parameters;
 using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
@@ -233 +234 @@
     }
 
-    public void ObjectiveFunction(double[] x, ref double func, double[] grad, object obj) {
+    public override void Analyze(ISymbolicExpressionTree[] trees, double[] qualities, ResultCollection results, IRandom random) {
+      if (!results.ContainsKey("Best Solution Quality")) {
+        results.Add(new Result("Best Solution Quality", typeof(DoubleValue)));
+      }
+      if (!results.ContainsKey("Best Tree")) {
+        results.Add(new Result("Best Tree", typeof(ISymbolicExpressionTree)));
+      }
+      if (!results.ContainsKey("Best Solution")) {
+        results.Add(new Result("Best Solution", typeof(GaussianProcessRegressionSolution)));
+      }
+
+      var bestQuality = qualities.Max();
+
+      if (results["Best Solution Quality"].Value == null || bestQuality > ((DoubleValue)results["Best Solution Quality"].Value).Value) {
+        var bestIdx = Array.IndexOf(qualities, bestQuality);
+        var bestClone = (ISymbolicExpressionTree)trees[bestIdx].Clone();
+        results["Best Tree"].Value = bestClone;
+        results["Best Solution Quality"].Value = new DoubleValue(bestQuality);
+        results["Best Solution"].Value = CreateSolution(bestClone, random);
+      }
+    }
+
+    private IItem CreateSolution(ISymbolicExpressionTree tree, IRandom random) {
+      // again tune the hyper-parameters.
+      // this is suboptimal because 1) more effort and 2) we cannot be sure to find the same local optimum
+      var meanFunction = new MeanConst();
+      var problemData = ProblemData;
+      var ds = problemData.Dataset;
+      var targetVariable = problemData.TargetVariable;
+      var allowedInputVariables = problemData.AllowedInputVariables.ToArray();
+      var nVars = allowedInputVariables.Length;
+      var trainingRows = problemData.TrainingIndices.ToArray();
+      var bestObjValue = new double[1] { double.MinValue };
+
+      // use the same covariance function for each restart
+      var covarianceFunction = TreeToCovarianceFunction(tree);
+      // data that is necessary for the objective function
+      var data = Tuple.Create(ds, targetVariable, allowedInputVariables, trainingRows, (IMeanFunction)meanFunction, covarianceFunction, bestObjValue);
+
+      // allocate hyperparameters
+      var hyperParameters = new double[meanFunction.GetNumberOfParameters(nVars) + covarianceFunction.GetNumberOfParameters(nVars) + 1]; // mean + cov + noise
+
+      // initialize hyperparameters
+      hyperParameters[0] = ds.GetDoubleValues(targetVariable).Average(); // mean const
+
+      for (int i = 0; i < covarianceFunction.GetNumberOfParameters(nVars); i++) {
+        hyperParameters[1 + i] = random.NextDouble() * 2.0 - 1.0;
+      }
+      hyperParameters[hyperParameters.Length - 1] = 1.0; // s² = exp(2), TODO: other inits better?
+
+      // use alglib.bfgs for hyper-parameter optimization ...
+      double epsg = 0;
+      double epsf = 0.00001;
+      double epsx = 0;
+      double stpmax = 1;
+      int maxits = ConstantOptIterations;
+      alglib.mincgstate state;
+      alglib.mincgreport rep;
+
+      alglib.mincgcreate(hyperParameters, out state);
+      alglib.mincgsetcond(state, epsg, epsf, epsx, maxits);
+      alglib.mincgsetstpmax(state, stpmax);
+      alglib.mincgoptimize(state, ObjectiveFunction, null, data);
+
+      alglib.mincgresults(state, out hyperParameters, out rep);
+
+      if (rep.terminationtype >= 0) {
+
+        var model = new GaussianProcessModel(ds, targetVariable, allowedInputVariables, trainingRows, hyperParameters, meanFunction, covarianceFunction);
+        return model.CreateRegressionSolution(ProblemData);
+      } else return null;
+    }
+
+    private void ObjectiveFunction(double[] x, ref double func, double[] grad, object obj) {
       // we want to optimize the model likelihood by changing the hyperparameters and also return the gradient for each hyperparameter
       var data = (Tuple<IDataset, string, string[], int[], IMeanFunction, ICovarianceFunction, double[]>)obj;
@@ -252 +326 @@
         var gradients = model.HyperparameterGradients;
         Array.Copy(gradients, grad, gradients.Length);
-      } catch (Exception) {
+      } catch (ArgumentException) {
         // building the GaussianProcessModel might fail, in this case we return the worst possible objective value
         func = 1.0E+300;

stable/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

@@ -165 +165 @@
                                              .ToArray();
       sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
-      CalculateModel(ds, rows, scaleInputs);
+      try {
+        CalculateModel(ds, rows, scaleInputs);
+      } catch (alglib.alglibexception ae) {
+        // wrap exception so that calling code doesn't have to know about alglib implementation
+        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
+      }
     }
 
@@ -306 +311 @@
 
     private IEnumerable<double> GetEstimatedValuesHelper(IDataset dataset, IEnumerable<int> rows) {
-      if (x == null) {
-        x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
-      }
-      int n = x.GetLength(0);
-
-      double[,] newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
-      int newN = newX.GetLength(0);
-
-      var Ks = new double[newN, n];
-      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, newX.GetLength(1)));
-      var ms = Enumerable.Range(0, newX.GetLength(0))
-      .Select(r => mean.Mean(newX, r))
-      .ToArray();
-      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, newX.GetLength(1)));
-      for (int i = 0; i < newN; i++) {
-        for (int j = 0; j < n; j++) {
-          Ks[i, j] = cov.CrossCovariance(x, newX, j, i);
-        }
-      }
-
-      return Enumerable.Range(0, newN)
-        .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha));
+      try {
+        if (x == null) {
+          x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
+        }
+        int n = x.GetLength(0);
+
+        double[,] newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
+        int newN = newX.GetLength(0);
+
+        var Ks = new double[newN, n];
+        var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, newX.GetLength(1)));
+        var ms = Enumerable.Range(0, newX.GetLength(0))
+        .Select(r => mean.Mean(newX, r))
+        .ToArray();
+        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, newX.GetLength(1)));
+        for (int i = 0; i < newN; i++) {
+          for (int j = 0; j < n; j++) {
+            Ks[i, j] = cov.CrossCovariance(x, newX, j, i);
+          }
+        }
+
+        return Enumerable.Range(0, newN)
+          .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha));
+      } catch (alglib.alglibexception ae) {
+        // wrap exception so that calling code doesn't have to know about alglib implementation
+        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
+      }
     }
 
     public IEnumerable<double> GetEstimatedVariance(IDataset dataset, IEnumerable<int> rows) {
-      if (x == null) {
-        x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
-      }
-      int n = x.GetLength(0);
-
-      var newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
-      int newN = newX.GetLength(0);
-
-      var kss = new double[newN];
-      double[,] sWKs = new double[n, newN];
-      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
-
-      if (l == null) {
-        l = CalculateL(x, cov, sqrSigmaNoise);
-      }
-
-      // for stddev 
-      for (int i = 0; i < newN; i++)
-        kss[i] = cov.Covariance(newX, i, i);
-
-      for (int i = 0; i < newN; i++) {
-        for (int j = 0; j < n; j++) {
-          sWKs[j, i] = cov.CrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);
-        }
-      }
-
-      // for stddev 
-      alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
-
-      for (int i = 0; i < newN; i++) {
-        var sumV = Util.ScalarProd(Util.GetCol(sWKs, i), Util.GetCol(sWKs, i));
-        kss[i] += sqrSigmaNoise; // kss is V(f), add noise variance of predictive distibution to get V(y)
-        kss[i] -= sumV;
-        if (kss[i] < 0) kss[i] = 0;
-      }
-      return kss;
+      try {
+        if (x == null) {
+          x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
+        }
+        int n = x.GetLength(0);
+
+        var newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
+        int newN = newX.GetLength(0);
+
+        var kss = new double[newN];
+        double[,] sWKs = new double[n, newN];
+        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
+
+        if (l == null) {
+          l = CalculateL(x, cov, sqrSigmaNoise);
+        }
+
+        // for stddev 
+        for (int i = 0; i < newN; i++)
+          kss[i] = cov.Covariance(newX, i, i);
+
+        for (int i = 0; i < newN; i++) {
+          for (int j = 0; j < n; j++) {
+            sWKs[j, i] = cov.CrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);
+          }
+        }
+
+        // for stddev 
+        alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
+
+        for (int i = 0; i < newN; i++) {
+          var sumV = Util.ScalarProd(Util.GetCol(sWKs, i), Util.GetCol(sWKs, i));
+          kss[i] += sqrSigmaNoise; // kss is V(f), add noise variance of predictive distibution to get V(y)
+          kss[i] -= sumV;
+          if (kss[i] < 0) kss[i] = 0;
+        }
+        return kss;
+      } catch (alglib.alglibexception ae) {
+        // wrap exception so that calling code doesn't have to know about alglib implementation
+        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
+      }
     }
   }