1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21 


22  using System;


23  using System.Collections.Generic;


24  using System.Linq;


25  using HeuristicLab.Common;


26  using HeuristicLab.Core;


27  using HEAL.Attic;


28  using HeuristicLab.Problems.DataAnalysis;


29  using HeuristicLab.Data;


30 


31  namespace HeuristicLab.Algorithms.DataAnalysis {


32  [StorableType("4148D88C60814D84B718C949CA5AA766")]


33  [Item("KernelRidgeRegressionModel", "A kernel ridge regression model")]


34  public sealed class KernelRidgeRegressionModel : RegressionModel {


35  public override IEnumerable<string> VariablesUsedForPrediction {


36  get { return allowedInputVariables; }


37  }


38 


39  [Storable]


40  private readonly string[] allowedInputVariables;


41  public string[] AllowedInputVariables {


42  get { return allowedInputVariables.ToArray(); }


43  }


44 


45 


46  [Storable]


47  public double LooCvRMSE { get; private set; }


48 


49  [Storable]


50  private readonly double[] alpha;


51 


52  [Storable]


53  private readonly double[,] trainX; // it is better to store the original training dataset completely because this is more efficient in persistence


54 


55  [Storable]


56  private readonly ITransformation<double>[] scaling;


57 


58  [Storable]


59  private readonly ICovarianceFunction kernel;


60 


61  [Storable]


62  private readonly double lambda;


63 


64  [Storable]


65  private readonly double yOffset; // implementation works for zeromean, unitvariance target variables


66 


67  [Storable]


68  private readonly double yScale;


69 


70  [StorableConstructor]


71  private KernelRidgeRegressionModel(StorableConstructorFlag _) : base(_) { }


72  private KernelRidgeRegressionModel(KernelRidgeRegressionModel original, Cloner cloner)


73  : base(original, cloner) {


74  // shallow copies of arrays because they cannot be modified


75  allowedInputVariables = original.allowedInputVariables;


76  alpha = original.alpha;


77  trainX = original.trainX;


78  scaling = original.scaling;


79  lambda = original.lambda;


80  LooCvRMSE = original.LooCvRMSE;


81 


82  yOffset = original.yOffset;


83  yScale = original.yScale;


84  kernel = original.kernel;


85  }


86  public override IDeepCloneable Clone(Cloner cloner) {


87  return new KernelRidgeRegressionModel(this, cloner);


88  }


89 


90  public static KernelRidgeRegressionModel Create(IDataset dataset, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,


91  bool scaleInputs, ICovarianceFunction kernel, double lambda = 0.1) {


92  var trainingRows = rows.ToArray();


93  var model = new KernelRidgeRegressionModel(dataset, targetVariable, allowedInputVariables, trainingRows, scaleInputs, kernel, lambda);


94 


95  try {


96  int info;


97  int n = model.trainX.GetLength(0);


98  alglib.densesolverreport denseSolveRep;


99  var gram = BuildGramMatrix(model.trainX, lambda, kernel);


100  var l = new double[n, n];


101  Array.Copy(gram, l, l.Length);


102 


103  double[] alpha = new double[n];


104  double[,] invG;


105  var y = dataset.GetDoubleValues(targetVariable, trainingRows).ToArray();


106  for (int i = 0; i < y.Length; i++) {


107  y[i] = model.yOffset;


108  y[i] *= model.yScale;


109  }


110  // cholesky decomposition


111  var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);


112  if (res == false) { //try lua decomposition if cholesky faild


113  int[] pivots;


114  var lua = new double[n, n];


115  Array.Copy(gram, lua, lua.Length);


116  alglib.rmatrixlu(ref lua, n, n, out pivots);


117  alglib.rmatrixlusolve(lua, pivots, n, y, out info, out denseSolveRep, out alpha);


118  if (info != 1) throw new ArgumentException("Could not create model.");


119  alglib.matinvreport rep;


120  invG = lua; // rename


121  alglib.rmatrixluinverse(ref invG, pivots, n, out info, out rep);


122  } else {


123  alglib.spdmatrixcholeskysolve(l, n, false, y, out info, out denseSolveRep, out alpha);


124  if (info != 1) throw new ArgumentException("Could not create model.");


125  // for LOOCV we need to build the inverse of the gram matrix


126  alglib.matinvreport rep;


127  invG = l; // rename


128  alglib.spdmatrixcholeskyinverse(ref invG, n, false, out info, out rep);


129  }


130  if (info != 1) throw new ArgumentException("Could not invert Gram matrix.");


131 


132  var ssqLooError = 0.0;


133  for (int i = 0; i < n; i++) {


134  var pred_i = Util.ScalarProd(Util.GetRow(gram, i).ToArray(), alpha);


135  var looPred_i = pred_i  alpha[i] / invG[i, i];


136  var error = (y[i]  looPred_i) / model.yScale;


137  ssqLooError += error * error;


138  }


139 


140  Array.Copy(alpha, model.alpha, n);


141  model.LooCvRMSE = Math.Sqrt(ssqLooError / n);


142  } catch (alglib.alglibexception ae) {


143  // wrap exception so that calling code doesn't have to know about alglib implementation


144  throw new ArgumentException("There was a problem in the calculation of the kernel ridge regression model", ae);


145  }


146  return model;


147  }


148 


149  private KernelRidgeRegressionModel(IDataset dataset, string targetVariable, IEnumerable<string> allowedInputVariables, int[] rows,


150  bool scaleInputs, ICovarianceFunction kernel, double lambda = 0.1) : base(targetVariable) {


151  this.allowedInputVariables = allowedInputVariables.ToArray();


152  if (kernel.GetNumberOfParameters(this.allowedInputVariables.Length) > 0) throw new ArgumentException("All parameters in the kernel function must be specified.");


153  name = ItemName;


154  description = ItemDescription;


155 


156  this.kernel = (ICovarianceFunction)kernel.Clone();


157  this.lambda = lambda;


158  if (scaleInputs) {


159  var trans = new ShiftToRangeTransformation(allowedInputVariables);


160  trans.Range.Start = 0;


161  trans.Range.End = 1;


162  scaling = Transformation.CreateTransformations(trans, dataset, rows, this.allowedInputVariables);


163  }


164 


165  trainX = ExtractData(dataset, rows, this.allowedInputVariables, scaling);


166  var y = dataset.GetDoubleValues(targetVariable, rows).ToArray();


167  yOffset = y.Average();


168  yScale = 1.0 / y.StandardDeviation();


169  alpha = new double[trainX.GetLength(0)];


170  }


171 


172 


173  #region IRegressionModel Members


174  public override IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {


175  var newX = ExtractData(dataset, rows, allowedInputVariables, scaling);


176  var dim = newX.GetLength(1);


177  var cov = kernel.GetParameterizedCovarianceFunction(new double[0], Enumerable.Range(0, dim).ToArray());


178 


179  var pred = new double[newX.GetLength(0)];


180  for (int i = 0; i < pred.Length; i++) {


181  double sum = 0.0;


182  for (int j = 0; j < alpha.Length; j++) {


183  sum += alpha[j] * cov.CrossCovariance(trainX, newX, j, i);


184  }


185  pred[i] = sum / yScale + yOffset;


186  }


187  return pred;


188  }


189  public override IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {


190  return new RegressionSolution(this, new RegressionProblemData(problemData));


191  }


192  #endregion


193 


194  #region helpers


195  private static double[,] BuildGramMatrix(double[,] data, double lambda, ICovarianceFunction kernel) {


196  var n = data.GetLength(0);


197  var dim = data.GetLength(1);


198  var cov = kernel.GetParameterizedCovarianceFunction(new double[0], Enumerable.Range(0, dim).ToArray());


199  var gram = new double[n, n];


200  // G = (K + λ I)


201  for (var i = 0; i < n; i++) {


202  for (var j = i; j < n; j++) {


203  gram[i, j] = gram[j, i] = cov.Covariance(data, i, j); // symmetric matrix


204  }


205  gram[i, i] += lambda;


206  }


207  return gram;


208  }


209 


210  private static double[,] ExtractData(IDataset dataset, IEnumerable<int> rows, IReadOnlyCollection<string> allowedInputVariables, ITransformation<double>[] scaling = null) {


211  double[][] variables;


212  if (scaling != null) {


213  variables =


214  allowedInputVariables.Select((var, i) => scaling[i].Apply(dataset.GetDoubleValues(var, rows)).ToArray())


215  .ToArray();


216  } else {


217  variables =


218  allowedInputVariables.Select(var => dataset.GetDoubleValues(var, rows).ToArray()).ToArray();


219  }


220  int n = variables.First().Length;


221  var res = new double[n, variables.Length];


222  for (int r = 0; r < n; r++)


223  for (int c = 0; c < variables.Length; c++) {


224  res[r, c] = variables[c][r];


225  }


226  return res;


227  }


228  #endregion


229  }


230  }

