Changeset 12819 for trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

Timestamp:

07/30/15 18:09:15 (9 years ago)

Author:

gkronber

Message:

#2449: improved persistence of GaussianProcessModel (Cholesky decomposed covariance matrix is not stored and recalculated lazily)

File:

• trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs (modified) (8 diffs)

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

@@ -85 +85 @@
     private double[] covarianceParameter;
 
-    [Storable]
-    private double[,] l;
-
-    [Storable]
-    private double[,] x;
+    private double[,] l; // used to be storable in previous versions (is calculated lazily now)
+    private double[,] x; // scaled training dataset, used to be storable in previous versions (is calculated lazily now)
+
+    // BackwardsCompatibility3.4
+    #region Backwards compatible code, remove with 3.5
+    [Storable(Name = "l")] // restore if available but don't store anymore
+    private double[,] l_storable {
+      set { this.l = value; }
+      get {
+        if (trainingDataset == null) return l; // this model has been created with an old version 
+        else return null; // if the training dataset is available l should not be serialized
+      }
+    }
+    [Storable(Name = "x")] // restore if available but don't store anymore
+    private double[,] x_storable {
+      set { this.x = value; }
+      get {
+        if (trainingDataset == null) return x; // this model has been created with an old version 
+        else return null; // if the training dataset is available x should not be serialized
+      }
+    }
+    #endregion
+
+
+    [Storable]
+    private IDataset trainingDataset; // it is better to store the original training dataset completely because this is more efficient in persistence
+    [Storable]
+    private int[] trainingRows;
+
     [Storable]
     private Scaling inputScaling;
@@ -101 +125 @@
       this.covarianceFunction = cloner.Clone(original.covarianceFunction);
       this.inputScaling = cloner.Clone(original.inputScaling);
+      this.trainingDataset = cloner.Clone(original.trainingDataset);
       this.negativeLogLikelihood = original.negativeLogLikelihood;
       this.targetVariable = original.targetVariable;
@@ -112 +137 @@
 
       // shallow copies of arrays because they cannot be modified
+      this.trainingRows = original.trainingRows;
       this.allowedInputVariables = original.allowedInputVariables;
       this.alpha = original.alpha;
@@ -137 +163 @@
                                              .ToArray();
       sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
-
       CalculateModel(ds, rows);
     }
 
     private void CalculateModel(IDataset ds, IEnumerable<int> rows) {
-      inputScaling = new Scaling(ds, allowedInputVariables, rows);
-      x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
+      this.trainingDataset = (IDataset)ds.Clone();
+      this.trainingRows = rows.ToArray();
+      this.inputScaling = new Scaling(trainingDataset, allowedInputVariables, rows);
+      this.x = CalculateX(trainingDataset, allowedInputVariables, rows, inputScaling);
       var y = ds.GetDoubleValues(targetVariable, rows);
 
       int n = x.GetLength(0);
-      l = new double[n, n];
-
-      // calculate means and covariances
+
+      // calculate cholesky decomposed (lower triangular) covariance matrix
+      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
+      this.l = CalculateL(x, cov, sqrSigmaNoise);
+
+      // calculate mean
       var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, x.GetLength(1)));
       double[] m = Enumerable.Range(0, x.GetLength(0))
@@ -155 +185 @@
         .ToArray();
 
-      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
-      for (int i = 0; i < n; i++) {
-        for (int j = i; j < n; j++) {
-          l[j, i] = cov.Covariance(x, i, j) / sqrSigmaNoise;
-          if (j == i) l[j, i] += 1.0;
-        }
-      }
-
-
-      // cholesky decomposition
+
+
+      // calculate sum of diagonal elements for likelihood
+      double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
+
+      // solve for alpha
+      double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
+
       int info;
       alglib.densesolverreport denseSolveRep;
-
-      var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
-      if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
-
-      // calculate sum of diagonal elements for likelihood
-      double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
-
-      // solve for alpha
-      double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
 
       alglib.spdmatrixcholeskysolve(l, n, false, ym, out info, out denseSolveRep, out alpha);
@@ -230 +249 @@
     }
 
+    private static double[,] CalculateX(IDataset ds, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows, Scaling inputScaling) {
+      return AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
+    }
+
+    private static double[,] CalculateL(double[,] x, ParameterizedCovarianceFunction cov, double sqrSigmaNoise) {
+      int n = x.GetLength(0);
+      var l = new double[n, n];
+
+      // calculate covariances
+      for (int i = 0; i < n; i++) {
+        for (int j = i; j < n; j++) {
+          l[j, i] = cov.Covariance(x, i, j) / sqrSigmaNoise;
+          if (j == i) l[j, i] += 1.0;
+        }
+      }
+
+      // cholesky decomposition
+      var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
+      if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
+      return l;
+    }
+
 
     public override IDeepCloneable Clone(Cloner cloner) {
@@ -258 +299 @@
 
     private IEnumerable<double> GetEstimatedValuesHelper(IDataset dataset, IEnumerable<int> rows) {
+      if (x == null) {
+        this.x = CalculateX(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
+      }
+      int n = x.GetLength(0);
+
       var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
       int newN = newX.GetLength(0);
-      int n = x.GetLength(0);
+
       var Ks = new double[newN, n];
       var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, newX.GetLength(1)));
@@ -278 +324 @@
 
     public IEnumerable<double> GetEstimatedVariance(IDataset dataset, IEnumerable<int> rows) {
+      if (x == null) {
+        this.x = CalculateX(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
+      }
+      int n = x.GetLength(0);
+
       var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
       int newN = newX.GetLength(0);
-      int n = x.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