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

Timestamp:

04/26/17 11:06:51 (7 years ago)

Author:

bwerth

Message:

#2699 reworked kenel functions (beta is always a scaling factor now), added LU-Decomposition as a fall-back if Cholesky-decomposition fails

File:

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

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

@@ -32 +32 @@
   [Item("KernelRidgeRegressionModel", "A kernel ridge regression model")]
   public sealed class KernelRidgeRegressionModel : RegressionModel {
-    public override IEnumerable<string> VariablesUsedForPrediction {
+    public override IEnumerable<string> VariablesUsedForPrediction
+    {
       get { return allowedInputVariables; }
     }
@@ -38 +39 @@
     [Storable]
     private readonly string[] allowedInputVariables;
-    public string[] AllowedInputVariables {
+    public string[] AllowedInputVariables
+    {
       get { return allowedInputVariables; }
     }
@@ -114 +116 @@
         var l = new double[n, n]; Array.Copy(gram, l, l.Length);
 
+        double[,] invG;
         // cholesky decomposition
         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.");
-
-
-        {
+        if (res == false) { //throw new ArgumentException("Could not decompose matrix. Is it quadratic symmetric positive definite?");
+          int[] pivots;
+          var lua = new double[n, n];
+          Array.Copy(gram, lua, lua.Length);
+          alglib.rmatrixlu(ref lua, n, n, out pivots);
+          alglib.rmatrixlusolve(lua, pivots, n, y, out info, out denseSolveRep, out alpha);
+          if (info != 1) throw new ArgumentException("Could not create model.");
+          alglib.matinvreport rep;
+          invG = lua;  // rename
+          alglib.rmatrixluinverse(ref invG, pivots, n, out info, out rep);
+          if (info != 1) throw new ArgumentException("Could not invert Gram matrix.");
+        } else {
+          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);          
+          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) {
+        }
+
+        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
         throw new ArgumentException("There was a problem in the calculation of the kernel ridge regression model", ae);