Changeset 14872 for branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/RadialBasisFunctions/RadialBasisFunctionModel.cs

Timestamp:

04/14/17 17:53:30 (7 years ago)

Author:

gkronber

Message:

#2699: made a number of changes mainly to RBF regression model

File:

• branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/RadialBasisFunctions/RadialBasisFunctionModel.cs (modified) (5 diffs)

branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/RadialBasisFunctions/RadialBasisFunctionModel.cs

@@ -1 @@
-#region License Information
+#region License Information
 /* HeuristicLab
  * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
@@ -22 +22 @@
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
-using HeuristicLab.Data;
 using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
 using HeuristicLab.Problems.DataAnalysis;
@@ -31 +31 @@
 namespace HeuristicLab.Algorithms.DataAnalysis {
   /// <summary>
-  /// Represents a Radial Basis Function regression model.
+  /// Represents an RBF regression model.
   /// </summary>
   [StorableClass]
-  [Item("RBFModel", "Represents a Gaussian process posterior.")]
+  [Item("RBFModel", "An RBF regression model")]
   public sealed class RadialBasisFunctionModel : RegressionModel, IConfidenceRegressionModel {
-    public override IEnumerable<string> VariablesUsedForPrediction
-    {
+    public override IEnumerable<string> VariablesUsedForPrediction {
       get { return allowedInputVariables; }
     }
 
     [Storable]
-    private string[] allowedInputVariables;
-    public string[] AllowedInputVariables
-    {
+    private readonly string[] allowedInputVariables;
+    public string[] AllowedInputVariables {
       get { return allowedInputVariables; }
     }
 
     [Storable]
-    private double[] alpha;
-    [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 IKernelFunction<double[]> kernel;
-    [Storable]
-    private DoubleMatrix gramInvert;
-
+    private readonly double[] alpha;
+
+    [Storable]
+    private readonly double[,] trainX; // it is better to store the original training dataset completely because this is more efficient in persistence
+
+    [Storable]
+    private readonly ITransformation<double>[] scaling;
+
+    [Storable]
+    private readonly ICovarianceFunction kernel;
+
+    private double[,] gramInvert; // not storable as it can be large (recreate after deserialization as required)
+
+    [Storable]
+    private readonly double meanOffset; // implementation works for zero-mean target variables
 
     [StorableConstructor]
@@ -65 +68 @@
       : base(original, cloner) {
       // shallow copies of arrays because they cannot be modified
-      trainingRows = original.trainingRows;
       allowedInputVariables = original.allowedInputVariables;
       alpha = original.alpha;
-      trainingDataset = original.trainingDataset;
-      kernel = original.kernel;
-    }
-    public RadialBasisFunctionModel(IDataset dataset, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows, IKernelFunction<double[]> kernel)
+      trainX = original.trainX;
+      gramInvert = original.gramInvert;
+      scaling = original.scaling;
+
+      meanOffset = original.meanOffset;
+      if (original.kernel != null)
+        kernel = cloner.Clone(original.kernel);
+    }
+    public RadialBasisFunctionModel(IDataset dataset, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
+      bool scaleInputs, ICovarianceFunction kernel)
       : base(targetVariable) {
+      if (kernel.GetNumberOfParameters(allowedInputVariables.Count()) > 0) throw new ArgumentException("All parameters in the kernel function must be specified.");
       name = ItemName;
       description = ItemDescription;
       this.allowedInputVariables = allowedInputVariables.ToArray();
-      trainingRows = rows.ToArray();
-      trainingDataset = dataset;
-      this.kernel = (IKernelFunction<double[]>)kernel.Clone();
+      var trainingRows = rows.ToArray();
+      this.kernel = (ICovarianceFunction)kernel.Clone();
       try {
-        var data = ExtractData(dataset, trainingRows);
-        var qualities = dataset.GetDoubleValues(targetVariable, trainingRows).ToArray();
+        if (scaleInputs)
+          scaling = CreateScaling(dataset, trainingRows);
+        trainX = ExtractData(dataset, trainingRows, scaling);
+        var y = dataset.GetDoubleValues(targetVariable, trainingRows).ToArray();
+        meanOffset = y.Average();
+        for (int i = 0; i < y.Length; i++) y[i] -= meanOffset;
         int info;
+        // TODO: improve efficiency by decomposing matrix once instead of solving the system and then inverting the matrix
         alglib.densesolverlsreport denseSolveRep;
-        var gr = BuildGramMatrix(data);
-        alglib.rmatrixsolvels(gr, data.Length + 1, data.Length + 1, qualities.Concat(new[] { 0.0 }).ToArray(), 0.0, out info, out denseSolveRep, out alpha);
-        if (info != 1) throw new ArgumentException("Could not create Model.");
-        gramInvert = new DoubleMatrix(gr).Invert();
-      }
-      catch (alglib.alglibexception ae) {
+        gramInvert = BuildGramMatrix(trainX);
+        int n = trainX.GetLength(0);
+        alglib.rmatrixsolvels(gramInvert, n, n, y, 0.0, out info, out denseSolveRep, out alpha);
+        if (info != 1) throw new ArgumentException("Could not create model.");
+
+        alglib.matinvreport report;
+        alglib.rmatrixinverse(ref gramInvert, out info, out report);
+        if (info != 1) throw new ArgumentException("Could not invert matrix. Is it quadratic symmetric positive definite?");
+
+      } 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 RBF process model", ae);
-      }
-    }
-    private double[][] ExtractData(IDataset dataset, IEnumerable<int> rows) {
-      return rows.Select(r => allowedInputVariables.Select(v => dataset.GetDoubleValue(v, r)).ToArray()).ToArray();
+        throw new ArgumentException("There was a problem in the calculation of the RBF model", ae);
+      }
+    }
+
+    private ITransformation<double>[] CreateScaling(IDataset dataset, int[] rows) {
+      var trans = new ITransformation<double>[allowedInputVariables.Length];
+      int i = 0;
+      foreach (var variable in allowedInputVariables) {
+        var lin = new LinearTransformation(allowedInputVariables);
+        var max = dataset.GetDoubleValues(variable, rows).Max();
+        var min = dataset.GetDoubleValues(variable, rows).Min();
+        lin.Multiplier = 1.0 / (max - min);
+        lin.Addend = -min / (max - min);
+        trans[i] = lin;
+        i++;
+      }
+      return trans;
+    }
+
+    private double[,] ExtractData(IDataset dataset, IEnumerable<int> rows, ITransformation<double>[] scaling = null) {
+      double[][] variables;
+      if (scaling != null) {
+        variables =
+          allowedInputVariables.Select((var, i) => scaling[i].Apply(dataset.GetDoubleValues(var, rows)).ToArray())
+            .ToArray();
+      } else {
+        variables =
+        allowedInputVariables.Select(var => dataset.GetDoubleValues(var, rows).ToArray()).ToArray();
+      }
+      int n = variables.First().Length;
+      var res = new double[n, variables.Length];
+      for (int r = 0; r < n; r++)
+        for (int c = 0; c < variables.Length; c++) {
+          res[r, c] = variables[c][r];
+        }
+      return res;
     }
 
@@ -104 +152 @@
     #region IRegressionModel Members
     public override IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
-      var solutions = ExtractData(dataset, rows);
-      var data = ExtractData(trainingDataset, trainingRows);
-      return solutions.Select(solution => alpha.Zip(data, (a, d) => a * kernel.Get(solution, d)).Sum() + 1 * alpha[alpha.Length - 1]).ToArray();
+      var newX = ExtractData(dataset, rows, scaling);
+      var dim = newX.GetLength(1);
+      var cov = kernel.GetParameterizedCovarianceFunction(new double[0], Enumerable.Range(0, dim).ToArray());
+
+      var pred = new double[newX.GetLength(0)];
+      for (int i = 0; i < pred.Length; i++) {
+        double sum = meanOffset;
+        for (int j = 0; j < alpha.Length; j++) {
+          sum += alpha[j] * cov.CrossCovariance(trainX, newX, j, i);
+        }
+        pred[i] = sum;
+      }
+      return pred;
     }
     public override IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
-      return new RadialBasisFunctionRegressionSolution(this, new RegressionProblemData(problemData));
+      return new ConfidenceRegressionSolution(this, new RegressionProblemData(problemData));
     }
     #endregion
 
     public IEnumerable<double> GetEstimatedVariances(IDataset dataset, IEnumerable<int> rows) {
-      var data = ExtractData(trainingDataset, trainingRows);
-      return ExtractData(dataset, rows).Select(x => GetVariance(x, data));
+      if (gramInvert == null) CreateAndInvertGramMatrix();
+      int n = gramInvert.GetLength(0);
+      var newData = ExtractData(dataset, rows, scaling);
+      var dim = newData.GetLength(1);
+      var cov = kernel.GetParameterizedCovarianceFunction(new double[0], Enumerable.Range(0, dim).ToArray());
+
+      // TODO perf (matrix matrix multiplication)
+      for (int i = 0; i < newData.GetLength(0); i++) {
+        double[] p = new double[n];
+
+        for (int j = 0; j < trainX.GetLength(0); j++) {
+          p[j] = cov.CrossCovariance(trainX, newData, j, i);
+        }
+
+        var Ap = new double[n];
+        alglib.ablas.rmatrixmv(n, n, gramInvert, 0, 0, 0, p, 0, ref Ap, 0);
+        var res = 0.0;
+        // dot product
+        for (int j = 0; j < p.Length; j++) res += p[j] * Ap[j];
+        yield return res > 0 ? res : 0;
+      }
     }
     public double LeaveOneOutCrossValidationRootMeanSquaredError() {
-      return Math.Sqrt(alpha.Select((t, i) => t / gramInvert[i, i]).Sum(d => d * d) / gramInvert.Rows);
-    }
-
+      if (gramInvert == null) CreateAndInvertGramMatrix();
+      var n = gramInvert.GetLength(0);
+      var s = 1.0 / n;
+
+      var sum = 0.0;
+      for (int i = 0; i < alpha.Length; i++) {
+        var x = alpha[i] / gramInvert[i, i];
+        sum += x * x;
+      }
+      sum *= s;
+      return Math.Sqrt(sum);
+    }
+
+    private void CreateAndInvertGramMatrix() {
+      try {
+        gramInvert = BuildGramMatrix(trainX);
+        int info = 0;
+        alglib.matinvreport report;
+        alglib.rmatrixinverse(ref gramInvert, out info, out report);
+        if (info != 1)
+          throw new ArgumentException("Could not invert matrix. Is it quadratic symmetric positive definite?");
+      } catch (alglib.alglibexception) {
+        // wrap exception so that calling code doesn't have to know about alglib implementation
+        throw new ArgumentException("Could not invert matrix. Is it quadratic symmetric positive definite?");
+      }
+    }
     #region helpers
-    private double[,] BuildGramMatrix(double[][] data) {
-      var size = data.Length + 1;
-      var gram = new double[size, size];
-      for (var i = 0; i < size; i++)
-        for (var j = i; j < size; j++) {
-          if (j == size - 1 && i == size - 1) gram[i, j] = 0;
-          else if (j == size - 1 || i == size - 1) gram[j, i] = gram[i, j] = 1;
-          else gram[j, i] = gram[i, j] = kernel.Get(data[i], data[j]); //symmteric Matrix --> half of the work
+    private double[,] BuildGramMatrix(double[,] data) {
+      var n = data.GetLength(0);
+      var dim = data.GetLength(1);
+      var cov = kernel.GetParameterizedCovarianceFunction(new double[0], Enumerable.Range(0, dim).ToArray());
+      var gram = new double[n, n];
+      for (var i = 0; i < n; i++)
+        for (var j = i; j < n; j++) {
+          gram[j, i] = gram[i, j] = cov.Covariance(data, i, j); // symmetric matrix --> half of the work
         }
       return gram;
     }
-    private double GetVariance(double[] solution, IEnumerable<double[]> data) {
-      var phiT = data.Select(x => kernel.Get(x, solution)).Concat(new[] { 1.0 }).ToColumnVector();
-      var res = phiT.Transpose().Mul(gramInvert.Mul(phiT))[0, 0];
-      return res > 0 ? res : 0;
-    }
+
     #endregion
   }