Changeset 14888 for branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/KernelRidgeRegression/KernelRidgeRegressionModel.cs

Timestamp:

04/25/17 09:42:18 (7 years ago)

Author:

gkronber

Message:

#2699: re-added calculation of leave one out cv estimate

File:

• branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/KernelRidgeRegression/KernelRidgeRegressionModel.cs (modified) (6 diffs)

branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/KernelRidgeRegression/KernelRidgeRegressionModel.cs

@@ -21 +21 @@
 
 using System;
-using System.Collections.Generic;    
+using System.Collections.Generic;
 using System.Linq;
 using HeuristicLab.Common;
@@ -42 +42 @@
     }
 
+
+    [Storable]
+    public double LooCvRMSE { get; private set; }
+
     [Storable]
     private readonly double[] alpha;
@@ -58 +62 @@
 
     [Storable]
-    private readonly double yOffset; // implementation works for zero-mean target variables
+    private readonly double yOffset; // implementation works for zero-mean, unit-variance target variables
 
     [Storable]
@@ -73 +77 @@
       scaling = original.scaling;
       lambda = original.lambda;
+      LooCvRMSE = original.LooCvRMSE;
 
       yOffset = original.yOffset;
@@ -104 +109 @@
         }
         int info;
+        int n = trainX.GetLength(0);
         alglib.densesolverreport denseSolveRep;
-        var gram = BuildGramMatrix(trainX, lambda); 
-        int n = trainX.GetLength(0);
+        var gram = BuildGramMatrix(trainX, lambda);
+        var l = new double[n, n]; Array.Copy(gram, l, l.Length);
 
         // cholesky decomposition
-        var res = alglib.trfac.spdmatrixcholesky(ref gram, n, false);
-        if(res == false) throw new ArgumentException("Could not decompose matrix. Is it quadratic symmetric positive definite?");
-
-        alglib.spdmatrixcholeskysolve(gram, n, false, y, out info, out denseSolveRep, out alpha);
+        var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
+        if (res == false) throw new ArgumentException("Could not decompose matrix. Is it quadratic symmetric positive definite?");
+
+        alglib.spdmatrixcholeskysolve(l, n, false, y, out info, out denseSolveRep, out alpha);
         if (info != 1) throw new ArgumentException("Could not create model.");
+
+
+        {
+          // for LOO-CV we need to build the inverse of the gram matrix
+          alglib.matinvreport rep;
+          var invG = l;   // rename 
+          alglib.spdmatrixcholeskyinverse(ref invG, n, false, out info, out rep);          
+          if (info != 1) throw new ArgumentException("Could not invert Gram matrix.");
+          var ssqLooError = 0.0;
+          for (int i = 0; i < n; i++) {
+            var pred_i = Util.ScalarProd(Util.GetRow(gram, i).ToArray(), alpha);
+            var looPred_i = pred_i - alpha[i] / invG[i, i];
+            var error = (y[i] - looPred_i) / yScale;
+            ssqLooError += error * error;
+          }
+          LooCvRMSE = Math.Sqrt(ssqLooError / n);
+        }
       } catch (alglib.alglibexception ae) {
         // wrap exception so that calling code doesn't have to know about alglib implementation
@@ -151 +174 @@
       for (var i = 0; i < n; i++) {
         for (var j = i; j < n; j++) {
-          gram[j, i] = cov.Covariance(data, i, j); // symmetric matrix --> fill only lower triangle
+          gram[i, j] = gram[j, i] = cov.Covariance(data, i, j); // symmetric matrix 
         }
         gram[i, i] += lambda;