Changeset 8463 for trunk/sources/HeuristicLab.Algorithms.DataAnalysis

Timestamp:

08/09/12 16:32:44 (13 years ago)

Author:

gkronber

Message:

#1902 improved GPR implementation

Location:

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4

Files:

• GaussianProcess/CovariancePeriodic.cs (modified) (1 diff)
• GaussianProcess/CovarianceProd.cs (modified) (5 diffs)
• GaussianProcess/CovarianceSEiso.cs (modified) (2 diffs)
• GaussianProcess/CovarianceSum.cs (modified) (4 diffs)
• GaussianProcess/GaussianProcessHyperparameterInitializer.cs (modified) (1 diff)
• GaussianProcess/GaussianProcessModel.cs (modified) (7 diffs)
• GaussianProcess/GaussianProcessRegressionSolutionCreator.cs (modified) (1 diff)
• GaussianProcess/MeanProd.cs (modified) (1 diff)
• GaussianProcess/Util.cs (modified) (2 diffs)
• Linear/Scaling.cs (modified) (1 diff)

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovariancePeriodic.cs ¶

@@ -102 +102 @@
     }
 
-
-    public double[] GetDiagonalCovariances() {
-      if (x != xt) throw new InvalidOperationException();
-      int rows = x.GetLength(0);
-      var cov = new double[rows];
-      for (int i = 0; i < rows; i++) {
-        double k = Math.Sqrt(Util.SqrDist(Util.GetRow(x, i), Util.GetRow(xt, i)));
-        k = Math.PI * k / p;
-        k = Math.Sin(k) / l;
-        k = k * k;
-        cov[i] = sf2 * Math.Exp(-2.0 * k);
-      }
-      return cov;
-    }
-
     public double GetGradient(int i, int j, int k) {
       double v = Math.PI * sd[i, j] / p;

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceProd.cs ¶

@@ -20 +20 @@
 #endregion
 
+using System;
+using System.Collections.Generic;
 using System.Linq;
 using HeuristicLab.Common;
@@ -48 +50 @@
       this.factors = cloner.Clone(original.factors);
       this.numberOfVariables = original.numberOfVariables;
+      AttachEventHandlers();
     }
 
@@ -53 +56 @@
       : base() {
       this.factors = new ItemList<ICovarianceFunction>();
+      AttachEventHandlers();
+    }
+
+    private void AttachEventHandlers() {
+      this.factors.CollectionReset += (sender, args) => ClearCache();
+      this.factors.ItemsAdded += (sender, args) => ClearCache();
+      this.factors.ItemsRemoved += (sender, args) => ClearCache();
+      this.factors.ItemsReplaced += (sender, args) => ClearCache();
+      this.factors.ItemsMoved += (sender, args) => ClearCache();
     }
 
@@ -86 +98 @@
     }
 
+    private Dictionary<int, Tuple<int, int>> cachedParameterMap;
     public double GetGradient(int i, int j, int k) {
-      // map from parameter index to factor
-      var vi = factors.Select((f, idx) => Enumerable.Repeat(idx, f.GetNumberOfParameters(numberOfVariables))).SelectMany(x => x).ToArray();
+      if (cachedParameterMap == null) {
+        CalculateParameterMap();
+      }
+      int ti = cachedParameterMap[k].Item1;
+      k = cachedParameterMap[k].Item2;
       double res = 1.0;
-      int jj = Enumerable.Range(0, k).Count(e => vi[e] == vi[k]);
       for (int ii = 0; ii < factors.Count; ii++) {
         var f = factors[ii];
-        if (ii == vi[k]) {
-          res *= f.GetGradient(i, j, jj);
+        if (ii == ti) {
+          res *= f.GetGradient(i, j, k);
         } else {
           res *= f.GetCovariance(i, j);
@@ -101 +116 @@
       return res;
     }
+
+    private void ClearCache() {
+      cachedParameterMap = null;
+    }
+
+    private void CalculateParameterMap() {
+      cachedParameterMap = new Dictionary<int, Tuple<int, int>>();
+      int k = 0;
+      for (int ti = 0; ti < factors.Count; ti++) {
+        for (int ti_k = 0; ti_k < factors[ti].GetNumberOfParameters(numberOfVariables); ti_k++) {
+          cachedParameterMap[k++] = Tuple.Create(ti, ti_k);
+        }
+      }
+    }
   }
 }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceSEiso.cs ¶

@@ -114 +114 @@
       int cols = xt.GetLength(0);
       sd = new double[rows, cols];
+      double lInv = 1.0 / l;
       if (symmetric) {
         for (int i = 0; i < rows; i++) {
           for (int j = i; j < rows; j++) {
-            sd[i, j] = Util.SqrDist(Util.GetRow(x, i).Select(e => e / l), Util.GetRow(xt, j).Select(e => e / l));
+            sd[i, j] = Util.SqrDist(Util.GetRow(x, i).Select(e => e * lInv), Util.GetRow(xt, j).Select(e => e * lInv));
             sd[j, i] = sd[i, j];
           }
@@ -124 +125 @@
         for (int i = 0; i < rows; i++) {
           for (int j = 0; j < cols; j++) {
-            sd[i, j] = Util.SqrDist(Util.GetRow(x, i).Select(e => e / l), Util.GetRow(xt, j).Select(e => e / l));
+            sd[i, j] = Util.SqrDist(Util.GetRow(x, i).Select(e => e * lInv), Util.GetRow(xt, j).Select(e => e * lInv));
           }
         }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceSum.cs ¶

@@ -20 +20 @@
 #endregion
 
+using System;
+using System.Collections.Generic;
 using System.Linq;
 using HeuristicLab.Common;
@@ -48 +50 @@
       this.terms = cloner.Clone(original.terms);
       this.numberOfVariables = original.numberOfVariables;
+      AttachEventHandlers();
     }
 
@@ -53 +56 @@
       : base() {
       this.terms = new ItemList<ICovarianceFunction>();
+      AttachEventHandlers();
+    }
+
+    private void AttachEventHandlers() {
+      this.terms.CollectionReset += (sender, args) => ClearCache();
+      this.terms.ItemsAdded += (sender, args) => ClearCache();
+      this.terms.ItemsRemoved += (sender, args) => ClearCache();
+      this.terms.ItemsReplaced += (sender, args) => ClearCache();
+      this.terms.ItemsMoved += (sender, args) => ClearCache();
     }
 
@@ -86 +98 @@
     }
 
+    private Dictionary<int, Tuple<int, int>> cachedParameterMap;
     public double GetGradient(int i, int j, int k) {
-      int ii = 0;
-      while (k > terms[ii].GetNumberOfParameters(numberOfVariables)) {
-        k -= terms[ii].GetNumberOfParameters(numberOfVariables);
+      if (cachedParameterMap == null) {
+        CalculateParameterMap();
       }
-      return terms[ii].GetGradient(i, j, k);
+      int ti = cachedParameterMap[k].Item1;
+      k = cachedParameterMap[k].Item2;
+      return terms[ti].GetGradient(i, j, k);
+    }
+    private void ClearCache() {
+      cachedParameterMap = null;
+    }
+
+    private void CalculateParameterMap() {
+      cachedParameterMap = new Dictionary<int, Tuple<int, int>>();
+      int k = 0;
+      for (int ti = 0; ti < terms.Count; ti++) {
+        for (int ti_k = 0; ti_k < terms[ti].GetNumberOfParameters(numberOfVariables); ti_k++) {
+          cachedParameterMap[k++] = Tuple.Create(ti, ti_k);
+        }
+      }
     }
   }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessHyperparameterInitializer.cs ¶

@@ -91 +91 @@
       var rand = RandomParameter.ActualValue;
       for (int i = 0; i < r.Length; i++)
-        r[i] = rand.NextDouble() * 4 - 2;
+        r[i] = rand.NextDouble() * 2 - 1;
 
       HyperparameterParameter.ActualValue = r;

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

@@ -73 +73 @@
     private double[,] x;
     [Storable]
-    private Scaling scaling;
+    private Scaling inputScaling;
+    [Storable]
+    private Scaling targetScaling;
 
 
@@ -82 +84 @@
       this.meanFunction = cloner.Clone(original.meanFunction);
       this.covarianceFunction = cloner.Clone(original.covarianceFunction);
-      this.scaling = cloner.Clone(original.scaling);
+      this.inputScaling = cloner.Clone(original.inputScaling);
+      this.targetScaling = cloner.Clone(original.targetScaling);
       this.negativeLogLikelihood = original.negativeLogLikelihood;
       this.targetVariable = original.targetVariable;
@@ -118 +121 @@
 
     private void CalculateModel(Dataset ds, IEnumerable<int> rows) {
-      scaling = new Scaling(ds, allowedInputVariables, rows);
-      x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, scaling);
-
-      var y = ds.GetDoubleValues(targetVariable, rows).ToArray();
+      inputScaling = new Scaling(ds, allowedInputVariables, rows);
+      x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
+
+
+      targetScaling = new Scaling(ds, new string[] { targetVariable }, rows);
+      var y = targetScaling.GetScaledValues(ds, targetVariable, rows);
 
       int n = x.GetLength(0);
@@ -162 +167 @@
       int n = x.GetLength(0);
       int nAllowedVariables = x.GetLength(1);
-      double[,] q = new double[n, n];
-      double[,] eye = new double[n, n];
-      for (int i = 0; i < n; i++) eye[i, i] = 1.0;
 
       int info;
-      alglib.densesolverreport denseSolveRep;
-
-      alglib.spdmatrixcholeskysolvem(l, n, false, eye, n, out info, out denseSolveRep, out q);
-      // double[,] a2 = outerProd(alpha, alpha);
+      alglib.matinvreport matInvRep;
+
+      alglib.spdmatrixcholeskyinverse(ref l, n, false, out info, out matInvRep);
+      if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");
       for (int i = 0; i < n; i++) {
-        for (int j = 0; j < n; j++)
-          q[i, j] = q[i, j] / sqrSigmaNoise - alpha[i] * alpha[j]; // a2[i, j];
-      }
-
-      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => q[i, i]).Sum();
+        for (int j = 0; j <= i; j++)
+          l[i, j] = l[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
+      }
+
+      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => l[i, i]).Sum();
 
       double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
@@ -187 +189 @@
       if (covGradients.Length > 0) {
         for (int i = 0; i < n; i++) {
-          for (int j = 0; j < n; j++) {
-            for (int k = 0; k < covGradients.Length; k++) {
-              covGradients[k] += q[i, j] * covarianceFunction.GetGradient(i, j, k);
+          for (int k = 0; k < covGradients.Length; k++) {
+            for (int j = 0; j < i; j++) {
+              covGradients[k] += l[i, j] * covarianceFunction.GetGradient(i, j, k);
             }
+            covGradients[k] += 0.5 * l[i, i] * covarianceFunction.GetGradient(i, i, k);
           }
         }
-        covGradients = covGradients.Select(g => g / 2.0).ToArray();
       }
 
@@ -219 +221 @@
 
     private IEnumerable<double> GetEstimatedValuesHelper(Dataset dataset, IEnumerable<int> rows) {
-      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, scaling);
+      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
       int newN = newX.GetLength(0);
       int n = x.GetLength(0);
@@ -251 +253 @@
       // alglib.rmatrixsolvem(l, n, sWKs, newN, true, out info, out denseSolveRep, out v);
 
-
-      for (int i = 0; i < newN; i++) {
-        // predMean[i] = ms[i] + prod(GetRow(Ks, i), alpha);
-        yield return ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha);
-        // var sumV2 = prod(GetCol(v, i), GetCol(v, i));
-        // predVar[i] = kss[i] - sumV2;
-      }
+      double targetScaleMin, targetScaleMax;
+      targetScaling.GetScalingParameters(targetVariable, out targetScaleMin, out targetScaleMax);
+      return Enumerable.Range(0, newN)
+        .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha))
+        .Select(m => m * (targetScaleMax - targetScaleMin) + targetScaleMin);
+      //for (int i = 0; i < newN; i++) {
+      //  // predMean[i] = ms[i] + prod(GetRow(Ks, i), alpha);
+      //  // var sumV2 = prod(GetCol(v, i), GetCol(v, i));
+      //  // predVar[i] = kss[i] - sumV2;
+      //}
 
     }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessRegressionSolutionCreator.cs ¶

@@ -75 +75 @@
 
     public override IOperation Apply() {
-      var m = ModelParameter.ActualValue;
-      var data = ProblemDataParameter.ActualValue;
+      var m = (IGaussianProcessModel)ModelParameter.ActualValue.Clone();
+      var data = (IRegressionProblemData)ProblemDataParameter.ActualValue.Clone();
       var s = new GaussianProcessRegressionSolution(m, data);
 

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/MeanProd.cs ¶

@@ -82 +82 @@
     public double[] GetGradients(int k, double[,] x) {
       double[] res = Enumerable.Repeat(1.0, x.GetLength(0)).ToArray();
-      foreach (var f in factors) {
-        var numParam = f.GetNumberOfParameters(numberOfVariables);
-        if (k >= 0 && k < numParam) {
+      // find index of factor for the given k
+      int j = 0;
+      while (k >= factors[j].GetNumberOfParameters(numberOfVariables)) {
+        k -= factors[j].GetNumberOfParameters(numberOfVariables);
+        j++;
+      }
+      for (int i = 0; i < factors.Count; i++) {
+        var f = factors[i];
+        if (i == j) {
           // multiply gradient
           var g = f.GetGradients(k, x);
-          for (int i = 0; i < res.Length; i++) res[i] *= g[i];
-          k -= numParam;
+          for (int ii = 0; ii < res.Length; ii++) res[ii] *= g[ii];
         } else {
           // multiply mean
           var m = f.GetMean(x);
-          for (int i = 0; i < res.Length; i++) res[i] *= m[i];
-          k -= numParam;
+          for (int ii = 0; ii < res.Length; ii++) res[ii] *= m[ii];
         }
       }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/Util.cs ¶

@@ -20 +20 @@
 #endregion
 
-using System;
 using System.Collections.Generic;
 using System.Linq;
@@ -32 +31 @@
     public static double SqrDist(double x, double y) {
       double d = x - y;
-      return Math.Max(d * d, 0.0);
+      return d * d;
     }
 
     public static double SqrDist(IEnumerable<double> x, IEnumerable<double> y) {
-      return x.Zip(y, SqrDist).Sum();
+      return x.Zip(y, (a, b) => (a - b) * (a - b)).Sum();
     }
 

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/Linear/Scaling.cs ¶

@@ -59 +59 @@
       return ds.GetDoubleValues(variable, rows).Select(x => (x - min) / (max - min));
     }
+
+    public void GetScalingParameters(string variable, out double min, out double max) {
+      min = scalingParameters[variable].Item1;
+      max = scalingParameters[variable].Item2;
+    }
   }
 }