Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
07/30/15 18:09:15 (9 years ago)
Author:
gkronber
Message:

#2449: improved persistence of GaussianProcessModel (Cholesky decomposed covariance matrix is not stored and recalculated lazily)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

    r12817 r12819  
    8585    private double[] covarianceParameter;
    8686
    87     [Storable]
    88     private double[,] l;
    89 
    90     [Storable]
    91     private double[,] x;
     87    private double[,] l; // used to be storable in previous versions (is calculated lazily now)
     88    private double[,] x; // scaled training dataset, used to be storable in previous versions (is calculated lazily now)
     89
     90    // BackwardsCompatibility3.4
     91    #region Backwards compatible code, remove with 3.5
     92    [Storable(Name = "l")] // restore if available but don't store anymore
     93    private double[,] l_storable {
     94      set { this.l = value; }
     95      get {
     96        if (trainingDataset == null) return l; // this model has been created with an old version
     97        else return null; // if the training dataset is available l should not be serialized
     98      }
     99    }
     100    [Storable(Name = "x")] // restore if available but don't store anymore
     101    private double[,] x_storable {
     102      set { this.x = value; }
     103      get {
     104        if (trainingDataset == null) return x; // this model has been created with an old version
     105        else return null; // if the training dataset is available x should not be serialized
     106      }
     107    }
     108    #endregion
     109
     110
     111    [Storable]
     112    private IDataset trainingDataset; // it is better to store the original training dataset completely because this is more efficient in persistence
     113    [Storable]
     114    private int[] trainingRows;
     115
    92116    [Storable]
    93117    private Scaling inputScaling;
     
    101125      this.covarianceFunction = cloner.Clone(original.covarianceFunction);
    102126      this.inputScaling = cloner.Clone(original.inputScaling);
     127      this.trainingDataset = cloner.Clone(original.trainingDataset);
    103128      this.negativeLogLikelihood = original.negativeLogLikelihood;
    104129      this.targetVariable = original.targetVariable;
     
    112137
    113138      // shallow copies of arrays because they cannot be modified
     139      this.trainingRows = original.trainingRows;
    114140      this.allowedInputVariables = original.allowedInputVariables;
    115141      this.alpha = original.alpha;
     
    137163                                             .ToArray();
    138164      sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
    139 
    140165      CalculateModel(ds, rows);
    141166    }
    142167
    143168    private void CalculateModel(IDataset ds, IEnumerable<int> rows) {
    144       inputScaling = new Scaling(ds, allowedInputVariables, rows);
    145       x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
     169      this.trainingDataset = (IDataset)ds.Clone();
     170      this.trainingRows = rows.ToArray();
     171      this.inputScaling = new Scaling(trainingDataset, allowedInputVariables, rows);
     172      this.x = CalculateX(trainingDataset, allowedInputVariables, rows, inputScaling);
    146173      var y = ds.GetDoubleValues(targetVariable, rows);
    147174
    148175      int n = x.GetLength(0);
    149       l = new double[n, n];
    150 
    151       // calculate means and covariances
     176
     177      // calculate cholesky decomposed (lower triangular) covariance matrix
     178      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
     179      this.l = CalculateL(x, cov, sqrSigmaNoise);
     180
     181      // calculate mean
    152182      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, x.GetLength(1)));
    153183      double[] m = Enumerable.Range(0, x.GetLength(0))
     
    155185        .ToArray();
    156186
    157       var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
    158       for (int i = 0; i < n; i++) {
    159         for (int j = i; j < n; j++) {
    160           l[j, i] = cov.Covariance(x, i, j) / sqrSigmaNoise;
    161           if (j == i) l[j, i] += 1.0;
    162         }
    163       }
    164 
    165 
    166       // cholesky decomposition
     187
     188
     189      // calculate sum of diagonal elements for likelihood
     190      double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
     191
     192      // solve for alpha
     193      double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
     194
    167195      int info;
    168196      alglib.densesolverreport denseSolveRep;
    169 
    170       var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
    171       if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
    172 
    173       // calculate sum of diagonal elements for likelihood
    174       double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
    175 
    176       // solve for alpha
    177       double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
    178197
    179198      alglib.spdmatrixcholeskysolve(l, n, false, ym, out info, out denseSolveRep, out alpha);
     
    230249    }
    231250
     251    private static double[,] CalculateX(IDataset ds, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows, Scaling inputScaling) {
     252      return AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
     253    }
     254
     255    private static double[,] CalculateL(double[,] x, ParameterizedCovarianceFunction cov, double sqrSigmaNoise) {
     256      int n = x.GetLength(0);
     257      var l = new double[n, n];
     258
     259      // calculate covariances
     260      for (int i = 0; i < n; i++) {
     261        for (int j = i; j < n; j++) {
     262          l[j, i] = cov.Covariance(x, i, j) / sqrSigmaNoise;
     263          if (j == i) l[j, i] += 1.0;
     264        }
     265      }
     266
     267      // cholesky decomposition
     268      var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
     269      if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
     270      return l;
     271    }
     272
    232273
    233274    public override IDeepCloneable Clone(Cloner cloner) {
     
    258299
    259300    private IEnumerable<double> GetEstimatedValuesHelper(IDataset dataset, IEnumerable<int> rows) {
     301      if (x == null) {
     302        this.x = CalculateX(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
     303      }
     304      int n = x.GetLength(0);
     305
    260306      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
    261307      int newN = newX.GetLength(0);
    262       int n = x.GetLength(0);
     308
    263309      var Ks = new double[newN, n];
    264310      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, newX.GetLength(1)));
     
    278324
    279325    public IEnumerable<double> GetEstimatedVariance(IDataset dataset, IEnumerable<int> rows) {
     326      if (x == null) {
     327        this.x = CalculateX(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
     328      }
     329      int n = x.GetLength(0);
     330
    280331      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
    281332      int newN = newX.GetLength(0);
    282       int n = x.GetLength(0);
    283333
    284334      var kss = new double[newN];
    285335      double[,] sWKs = new double[n, newN];
    286336      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
     337
     338      if (l == null) {
     339        l = CalculateL(x, cov, sqrSigmaNoise);
     340      }
    287341
    288342      // for stddev
Note: See TracChangeset for help on using the changeset viewer.