Changeset 15442

Timestamp:

10/31/17 20:01:02 (7 years ago)

Author:

gkronber

Message:

#2789: added cross-validation and automatic determination of regularization parameter (naive CV is slow!).

Research produced ACM Algorithm 642.f (Fortran) for univariate cubic spline smoothing with generalized cross-validation in O(n) described in HUTCHINSON, M. F., AND DE HOOG, F.R. Smoothing noisy data with spline functions. Numer. Math. 4 7 (1985), 99-106. (TODO: translate to C#)

Location:

branches/MathNetNumerics-Exploration-2789

Files:

• HeuristicLab.Algorithms.DataAnalysis.Experimental/642.f (added)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs (modified) (1 diff)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs (modified) (13 diffs)
• Test/Test.csproj (modified) (1 diff)
• Test/UnitTest1.cs (modified) (2 diffs)

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs

@@ -218 +218 @@
           product = product.Zip(problemData.Dataset.GetDoubleValues(inputVars[i], problemData.TrainingIndices), (pi, vi) => pi * vi).ToArray();
         }
-        // Umständlich!
-        return Splines.CalculatePenalizedRegressionSpline(
+        double optTolerance, looRMSE;
+        return Splines.CalculateSmoothingSplineReinsch(
           product,
-          (double[])target.Clone(), lambda,
-          problemData.TargetVariable, inputVars
+          (double[])target.Clone(), inputVars,
+          problemData.TargetVariable, 
+          out optTolerance, out looRMSE
           );
       } else return new ConstantModel(target.Average(), problemData.TargetVariable);

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs

@@ -21 +21 @@
 
 using System;
-using System.Collections.Concurrent;
 using System.Collections.Generic;
 using System.Linq;
 using System.Threading;
-using System.Threading.Tasks;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
 using HeuristicLab.Data;
-using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
 using HeuristicLab.Optimization;
 using HeuristicLab.Parameters;
 using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
 using HeuristicLab.Problems.DataAnalysis;
-using HeuristicLab.Problems.DataAnalysis.Symbolic;
-using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
 using HeuristicLab.Random;
 using MathNet.Numerics.Interpolation;
@@ -71 +66 @@
           "Smoothing Spline (Mine)",
           "Smoothing Spline (Reinsch)",
-          "Smoothing Spline (Reinsch with automatic tolerance determination using CV)",
+          "Smoothing Spline (Reinsch with automatic tolerance determination using LOOCV)",
           "B-Spline Smoothing",
           "Penalized Regression Spline (alglib)"
@@ -150 +145 @@
               break;
             }
-          case "Smoothing Spline (Reinsch with automatic tolerance determination using CV)": {
-              var model = CalculateSmoothingSplineReinschWithAutomaticTolerance(x, y, inputVars, Problem.ProblemData.TargetVariable);
-              Results.Add(new Result("Solution", new RegressionSolution(model, (IRegressionProblemData)Problem.ProblemData.Clone())));
+          case "Smoothing Spline (Reinsch with automatic tolerance determination using LOOCV)": {
+              double optTol, looRMSE;
+              var model = CalculateSmoothingSplineReinsch(x, y, inputVars, Problem.ProblemData.TargetVariable, out optTol, out looRMSE);
+              Results.Add(new Result("Solution", model.CreateRegressionSolution((IRegressionProblemData)Problem.ProblemData.Clone())));
+              Results.Add(new Result("Optimal tolerance", new DoubleValue(optTol)));
+              Results.Add(new Result("RMSE (LOO-CV)", new DoubleValue(looRMSE)));
               break;
             }
@@ -159 +157 @@
               break;
             }
-          case "Penalized Regression Spline (alglib)":
-          {
+          case "Penalized Regression Spline (alglib)": {
               var lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value; // controls extent of smoothing
               var model = CalculatePenalizedRegressionSpline(x, y, lambda, Problem.ProblemData.TargetVariable, inputVars);
@@ -168 +165 @@
 
               break;
-          }
+            }
           default: throw new NotSupportedException();
         }
@@ -176 +173 @@
 
 
-
-    public static IRegressionModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, string targetVar, string[] inputVars)
-    {
+    public static IRegressionModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, string targetVar, string[] inputVars) {
       int info;
       alglib.spline1dinterpolant interpolant;
@@ -189 +184 @@
     }
 
+    public static IRegressionEnsembleModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, string targetVar, string[] inputVars,
+      out double avgTrainRMSE,
+      out double loocvRMSE) {
+      int info;
+      alglib.spline1dinterpolant interpolant;
+      alglib.spline1dfitreport rep;
+      int n = x.Length;
+      n = (int)Math.Max(50, 3 * Math.Sqrt(n));
+
+      double[] x_loo = new double[x.Length - 1];
+      double[] y_loo = new double[y.Length - 1];
+
+      var models = new List<IRegressionModel>();
+      var sumTrainRmse = 0.0;
+      var sse = 0.0;
+
+      for (int i = 0; i < x.Length; i++) {
+        Array.Copy(x, 0, x_loo, 0, i);
+        Array.Copy(x, i + 1, x_loo, i, x.Length - (i + 1));
+        Array.Copy(y, 0, y_loo, 0, i);
+        Array.Copy(y, i + 1, y_loo, i, y.Length - (i + 1));
+        alglib.spline1dfitpenalized(x_loo, y_loo, n, lambda, out info, out interpolant, out rep);
+        var y_pred = alglib.spline1dcalc(interpolant, x[i]);
+        double res = y[i] - y_pred;
+        sse += res * res;
+        sumTrainRmse += rep.rmserror;
+
+        var model = new Spline1dModel(interpolant, targetVar, inputVars);
+        models.Add(model);
+      }
+
+      loocvRMSE = Math.Sqrt(sse / x.Length);
+      avgTrainRMSE = sumTrainRmse / x.Length;
+
+      return new RegressionEnsembleModel(models);
+    }
+
     public void BSplineSmoothing(double[] x, double[] y, string[] inputVars) {
       Array.Sort(x, y);
-
       CubicSpline2(x, y, inputVars);
-
-      //  // see Elements of Statistical Learning Appendix: Computations for Splines
-      //  int M = 4; // order of spline (cubic)
-      //  int K = x.Length;
-      //  int n = x.Length;
-      // 
-      //  double[] t = new double[K + 2 * M];
-      //  for (int i = 0; i < M; i++) t[i] = x[0];
-      //  for (int i = K; i < K + 2 * M; i++) t[i] = x[K - 1];
-      // 
-      //  double[,] B = new double[n, n + 4]; // basis function matrix
-      //  double[,] W = new double[n + 4, n + 4]; // penalty matrix (W[j,k] = 
-      //  for (int i = 0; i < M; i++) B[i] = new double[K + 2 * M];
-      //  for (int j = 0; j < K; j++) {
-      //    for (int i = 0; i < K + 2M - 1; i++) {
-      //      if (t[i] <= x[j] && x[j] < t[i + 1])
-      //        B[1][i] = 1.0;
-      //    }
-      //  }
-
-
     }
 
@@ -312 +322 @@
     }
 
-    public void CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars) {
-      double s = xOrig.Length + 1;
+    public void CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars) {
+      double s = x.Length + 1;
       double stdTol = ((IValueParameter<DoubleValue>)Parameters["StdDev (noise)"]).Value.Value; // controls extent of smoothing
 
-      var model = CalculateSmoothingSplineReinsch(xOrig, yOrig, inputVars, stdTol, s, Problem.ProblemData.TargetVariable);
-
-      Results.Add(new Result("Solution", new RegressionSolution(model, (IRegressionProblemData)Problem.ProblemData.Clone())));
-
-
-    }
-
-
-
-    public static IRegressionModel CalculateSmoothingSplineReinschWithAutomaticTolerance(double[] xOrig, double[] yOrig, string[] inputVars, string targetVar) {
-      int n = xOrig.Length;
-      var rows = Enumerable.Range(0, n);
-      var rand = new FastRandom(1234);
-      var trainRows = rows.Shuffle(rand).Take((int)(n * 0.66)).ToArray();
-      var testRows = rows.Except(trainRows).ToArray();
-
-      var trainX = trainRows.Select(i => xOrig[i]).ToArray();
-      var trainY = trainRows.Select(i => yOrig[i]).ToArray();
-      var testX = testRows.Select(i => xOrig[i]).ToArray();
-      var testY = testRows.Select(i => yOrig[i]).ToArray();
-
-      var shuffledDs = new Dataset(inputVars.Concat(new string[] { targetVar }), new[] { trainX.Concat(testX).ToArray(), trainY.Concat(testY).ToArray() });
-      var shuffeledProblemData = new RegressionProblemData(shuffledDs, inputVars, targetVar);
-      shuffeledProblemData.TrainingPartition.Start = 0;
-      shuffeledProblemData.TrainingPartition.End = trainX.Length;
-      shuffeledProblemData.TestPartition.Start = trainX.Length;
-      shuffeledProblemData.TestPartition.End = shuffledDs.Rows + 1;
-
-      double curTol = trainY.StandardDeviation();
-      double prevTestRMSE = double.PositiveInfinity;
-      double testRMSE = double.PositiveInfinity;
-      IRegressionModel prevModel = null;
-      IRegressionModel model = null;
-      do {                                     
-        prevModel = model;
-        prevTestRMSE = testRMSE;
-        model = CalculateSmoothingSplineReinsch(trainX, trainY, inputVars, curTol, s: trainX.Length + 1, targetVar: targetVar);
-
-        var sol = model.CreateRegressionSolution(shuffeledProblemData);
-        var trainRMSE = sol.TrainingRootMeanSquaredError;
-        testRMSE = sol.TestRootMeanSquaredError;
-        curTol = Math.Max(curTol * 0.5, 1e-12 * trainY.StandardDeviation());
-      } while (testRMSE < prevTestRMSE);
-      return prevModel;
-    }
-
-
-    public static IRegressionModel CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars, double stdTol, double s, string targetVar) {
+      var model = CalculateSmoothingSplineReinsch(x, y, inputVars, stdTol, s, Problem.ProblemData.TargetVariable);
+
+      Results.Add(new Result("Solution", model.CreateRegressionSolution((IRegressionProblemData)Problem.ProblemData.Clone())));
+    }
+
+
+    public static IRegressionEnsembleModel CalculateSmoothingSplineReinsch(double[] x, double[] y, double tol, string targetVar, string[] inputVars,
+      out double avgTrainRMSE,
+      out double loocvRMSE) {
+      int info;
+
+      double[] x_loo = new double[x.Length - 1];
+      double[] y_loo = new double[y.Length - 1];
+
+      var models = new List<IRegressionModel>();
+      var sumTrainRmse = 0.0;
+      var sse = 0.0;
+
+      for (int i = 0; i < x.Length; i++) {
+        Array.Copy(x, i + 1, x_loo, i, x.Length - (i + 1));
+        Array.Copy(y, i + 1, y_loo, i, y.Length - (i + 1));
+
+        var model = CalculateSmoothingSplineReinsch(x_loo, y_loo, inputVars, tol, x_loo.Length + 1, targetVar);
+
+        var y_pred = model.GetEstimatedValue(x[i]);
+        double res = y[i] - y_pred;
+        sse += res * res;
+
+        models.Add(model);
+
+        // add current row for next iteration
+        if (i < x_loo.Length) {
+          x_loo[i] = x[i];
+          y_loo[i] = y[i];
+        }
+      }
+
+      loocvRMSE = Math.Sqrt(sse / x.Length);
+      avgTrainRMSE = sumTrainRmse / x.Length;
+
+      return new RegressionEnsembleModel(models);
+    }
+
+
+    // automatically determines tolerance using loo cv
+    public static IRegressionModel CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars, string targetVar,
+      out double optTolerance, out double looRMSE) {
+      var maxTolerance = y.StandardDeviation();
+      var minTolerance = maxTolerance * 1E-6;
+      optTolerance = maxTolerance;
+      var tol = maxTolerance;
+      looRMSE = double.PositiveInfinity;
+      IRegressionModel bestModel = new ConstantModel(y.Average(), targetVar);
+      while (tol > minTolerance) {
+        double loo_rmse;
+        double avgTrainRmse;
+        var model = Splines.CalculateSmoothingSplineReinsch(
+          x, y, tol, targetVar, inputVars, out avgTrainRmse, out loo_rmse);
+
+        if (loo_rmse < looRMSE) {
+          looRMSE = loo_rmse;
+          optTolerance = tol;
+          bestModel = model;
+        }
+        tol *= 0.85;
+      }
+      return bestModel;
+    }
+
+    public static ReinschSmoothingSplineModel CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars, double stdTol, double s, string targetVar) {
       var minX = xOrig.Min();
       var maxX = xOrig.Max();
@@ -605 +637 @@
 
   // UNFINISHED
-  internal class ReinschSmoothingSplineModel : NamedItem, IRegressionModel {
+  public class ReinschSmoothingSplineModel : NamedItem, IRegressionModel {
     private MyArray<double> a;
     private MyArray<double> b;
@@ -631 +663 @@
       this.offset = orig.offset;
     }
-    public ReinschSmoothingSplineModel(
+    internal ReinschSmoothingSplineModel(
       MyArray<double> a,
       MyArray<double> b,
@@ -650 +682 @@
       // extrapolate for xx > x[n2]
       b[b.Length] = b[b.Length - 1];
-      d[b.Length] = d[d.Length - 1];
+      d[1] = 0;
+      // d[b.Length] = -d[d.Length - 1];
     }
 
@@ -661 +694 @@
     }
 
+    public double GetEstimatedValue(double xx) {
+      int n = x.Length;
+      if (xx <= x[1]) {
+        double h = xx - x[1];
+        return a[1] + h * (b[1] + h * (c[1] + h * d[1]));
+      } else if (xx >= x[n]) {
+        double h = xx - x[n];
+        return a[n] + h * (b[n] + h * (c[n] + h * d[n]));
+      } else {
+        // binary search
+        int lower = 1;
+        int upper = n;
+        while (true) {
+          if (upper < lower) throw new InvalidProgramException();
+          int i = lower + (upper - lower) / 2;
+          if (x[i] <= xx && xx < x[i + 1]) {
+            double h = xx - x[i];
+            return a[i] + h * (b[i] + h * (c[i] + h * d[i]));
+          } else if (xx < x[i]) {
+            upper = i - 1;
+          } else {
+            lower = i + 1;
+          }
+        }
+      }
+    }
+
     public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
       int n = x.Length;
-      foreach (var xx in dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows).Select(xi => (xi - offset) * scale)) {
-        if (xx <= x[1]) {
-          double h = xx - x[1];
-          yield return a[1] + h * (b[1] + h * (c[1] + h * d[1]));
-        } else if (xx >= x[n]) {
-          double h = xx - x[n];
-          yield return a[n] + h * (b[n] + h * (c[n] + h * d[n]));
-        } else {
-          // binary search
-          int lower = 1;
-          int upper = n;
-          while (true) {
-            if (upper < lower) throw new InvalidProgramException();
-            int i = lower + (upper - lower) / 2;
-            if (x[i] <= xx && xx < x[i + 1]) {
-              double h = xx - x[i];
-              yield return a[i] + h * (b[i] + h * (c[i] + h * d[i]));
-              break;
-            } else if (xx < x[i]) {
-              upper = i - 1;
-            } else {
-              lower = i + 1;
-            }
-          }
-        }
+      var product = dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows).ToArray();
+      foreach (var factor in VariablesUsedForPrediction.Skip(1)) {
+        product = product.Zip(dataset.GetDoubleValues(factor, rows), (pi, fi) => pi * fi).ToArray();
+      }
+      foreach (var xx in product.Select(xi => (xi - offset) * scale)) {
+        yield return GetEstimatedValue(xx);
       }
     }
@@ -764 +805 @@
     public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
       var product = dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows).ToArray();
-      foreach(var factor in VariablesUsedForPrediction.Skip(1)) {
+      foreach (var factor in VariablesUsedForPrediction.Skip(1)) {
         product = product.Zip(dataset.GetDoubleValues(factor, rows), (pi, fi) => pi * fi).ToArray();
       }

branches/MathNetNumerics-Exploration-2789/Test/Test.csproj

@@ -36 +36 @@
   </PropertyGroup>
   <ItemGroup>
+    <Reference Include="ALGLIB-3.7.0, Version=3.7.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.7.0.dll</HintPath>
+    </Reference>
+    <Reference Include="HeuristicLab.Algorithms.DataAnalysis-3.4, Version=3.4.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec" />
     <Reference Include="HeuristicLab.Algorithms.GeneticAlgorithm-3.3">
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Algorithms.GeneticAlgorithm-3.3.dll</HintPath>

branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs

@@ -1 +1 @@
 using System;
+using System.Diagnostics;
+using System.Linq;
 using HeuristicLab.Algorithms.DataAnalysis.Experimental;
 using HeuristicLab.Algorithms.GeneticAlgorithm;
+using HeuristicLab.Common;
+using HeuristicLab.Problems.DataAnalysis;
 using HeuristicLab.SequentialEngine;
 using Microsoft.VisualStudio.TestTools.UnitTesting;
@@ -22 +26 @@
     }
   }
+
+  [TestClass]
+  public class PenalizedRegressionSplineTests {
+    [TestMethod]
+    public void TestPenalizedRegressionSplinesCrossValidation() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToArray();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToArray();
+
+      double bestRho = -15;
+      double best_loo_rmse = double.PositiveInfinity;
+      var clock = new Stopwatch();
+      clock.Start();
+      int iters = 0;
+      for (int rho = -15; rho <= 15; rho++) {
+        double loo_rmse;
+        double avgTrainRmse;
+        Splines.CalculatePenalizedRegressionSpline(xs, ys_noise, rho, "y", new string[] { "x" }, out avgTrainRmse, out loo_rmse);
+        iters++;
+        Console.WriteLine("{0} {1} {2}", rho, avgTrainRmse, loo_rmse);
+        if (loo_rmse < best_loo_rmse) {
+          best_loo_rmse = loo_rmse;
+          bestRho = rho;
+        }
+      }
+      clock.Stop();
+      Console.WriteLine("Best rho {0}, RMSE (LOO): {1}, ms/run {2}", bestRho, best_loo_rmse, clock.ElapsedMilliseconds / (double)iters);
+    }
+
+    [TestMethod]
+    public void TestReinschSplineCrossValidation() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToList();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToList();
+
+      var tol = ys_noise.StandardDeviation();
+
+      var ds = new Dataset(new string[] { "x", "y" }, new[] { xs, ys_noise });
+      var rows = Enumerable.Range(0, xs.Count);
+
+      var bestTol = double.PositiveInfinity;
+      var bestLooRmse = double.PositiveInfinity;
+      var clock = new Stopwatch();
+      clock.Start();
+      var iters = 0;
+      while (tol > 0.0001) {
+        double loo_rmse;
+        double avgTrainRmse;
+        var model = Splines.CalculateSmoothingSplineReinsch(
+          xs.ToArray(), ys_noise.ToArray(), tol, "y", new string[] { "x" }, out avgTrainRmse, out loo_rmse);
+        var y_pred = model.GetEstimatedValues(ds, rows);
+
+        Console.WriteLine("{0} {1} {2}", tol, avgTrainRmse, loo_rmse);
+        if (loo_rmse < bestLooRmse) {
+          bestLooRmse = loo_rmse;
+          bestTol = tol;
+        }
+        tol *= 0.8;
+        iters++;
+      }
+      clock.Stop();
+      Console.WriteLine("Best tolerance {0}, RMSE (LOO): {1}, ms/run {2}", bestTol, bestLooRmse, clock.ElapsedMilliseconds / (double)iters);
+    }
+
+    [TestMethod]
+    public void TestReinschSplineAutomaticTolerance() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToList();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToList();
+
+      double optTol;
+      double looRMSE;
+      Splines.CalculateSmoothingSplineReinsch(xs.ToArray(), ys_noise.ToArray(), new string[] { "x" }, "y", out optTol, out looRMSE);
+
+      Console.WriteLine("Best tolerance {0}, RMSE (LOO): {1}", optTol, looRMSE);
+
+    }
+  }
 }