Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Weighted TSNE/3.4/GaussianProcess/GaussianProcessModel.cs @ 15568

Last change on this file since 15568 was 15165, checked in by gkronber, 7 years ago

#2782: fixed calculation of log pseudo-likelihood by adding the noise term to the covariance function

File size: 16.5 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2016 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
22using System;
23using System.Collections.Generic;
24using System.Linq;
25using HeuristicLab.Common;
26using HeuristicLab.Core;
27using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
28using HeuristicLab.Problems.DataAnalysis;
29
30namespace HeuristicLab.Algorithms.DataAnalysis {
31  /// <summary>
32  /// Represents a Gaussian process model.
33  /// </summary>
34  [StorableClass]
35  [Item("GaussianProcessModel", "Represents a Gaussian process posterior.")]
36  public sealed class GaussianProcessModel : RegressionModel, IGaussianProcessModel {
37    public override IEnumerable<string> VariablesUsedForPrediction {
38      get { return allowedInputVariables; }
39    }
40
41    [Storable]
42    private double negativeLogLikelihood;
43    public double NegativeLogLikelihood {
44      get { return negativeLogLikelihood; }
45    }
46
47    [Storable]
48    private double loocvNegLogPseudoLikelihood;
49    public double LooCvNegativeLogPseudoLikelihood {
50      get { return loocvNegLogPseudoLikelihood; }
51    }
52
53    [Storable]
54    private double[] hyperparameterGradients;
55    public double[] HyperparameterGradients {
56      get {
57        var copy = new double[hyperparameterGradients.Length];
58        Array.Copy(hyperparameterGradients, copy, copy.Length);
59        return copy;
60      }
61    }
62
63    [Storable]
64    private ICovarianceFunction covarianceFunction;
65    public ICovarianceFunction CovarianceFunction {
66      get { return covarianceFunction; }
67    }
68    [Storable]
69    private IMeanFunction meanFunction;
70    public IMeanFunction MeanFunction {
71      get { return meanFunction; }
72    }
73
74    [Storable]
75    private string[] allowedInputVariables;
76    public string[] AllowedInputVariables {
77      get { return allowedInputVariables; }
78    }
79
80    [Storable]
81    private double[] alpha;
82    [Storable]
83    private double sqrSigmaNoise;
84    public double SigmaNoise {
85      get { return Math.Sqrt(sqrSigmaNoise); }
86    }
87
88    [Storable]
89    private double[] meanParameter;
90    [Storable]
91    private double[] covarianceParameter;
92
93    private double[,] l; // used to be storable in previous versions (is calculated lazily now)
94    private double[,] x; // scaled training dataset, used to be storable in previous versions (is calculated lazily now)
95
96    // BackwardsCompatibility3.4
97    #region Backwards compatible code, remove with 3.5
98    [Storable(Name = "l")] // restore if available but don't store anymore
99    private double[,] l_storable {
100      set { this.l = value; }
101      get {
102        if (trainingDataset == null) return l; // this model has been created with an old version
103        else return null; // if the training dataset is available l should not be serialized
104      }
105    }
106    [Storable(Name = "x")] // restore if available but don't store anymore
107    private double[,] x_storable {
108      set { this.x = value; }
109      get {
110        if (trainingDataset == null) return x; // this model has been created with an old version
111        else return null; // if the training dataset is available x should not be serialized
112      }
113    }
114    #endregion
115
116
117    [Storable]
118    private IDataset trainingDataset; // it is better to store the original training dataset completely because this is more efficient in persistence
119    [Storable]
120    private int[] trainingRows;
121
122    [Storable]
123    private Scaling inputScaling;
124
125
126    [StorableConstructor]
127    private GaussianProcessModel(bool deserializing) : base(deserializing) { }
128    private GaussianProcessModel(GaussianProcessModel original, Cloner cloner)
129      : base(original, cloner) {
130      this.meanFunction = cloner.Clone(original.meanFunction);
131      this.covarianceFunction = cloner.Clone(original.covarianceFunction);
132      if (original.inputScaling != null)
133        this.inputScaling = cloner.Clone(original.inputScaling);
134      this.trainingDataset = cloner.Clone(original.trainingDataset);
135      this.negativeLogLikelihood = original.negativeLogLikelihood;
136      this.loocvNegLogPseudoLikelihood = original.loocvNegLogPseudoLikelihood;
137      this.sqrSigmaNoise = original.sqrSigmaNoise;
138      if (original.meanParameter != null) {
139        this.meanParameter = (double[])original.meanParameter.Clone();
140      }
141      if (original.covarianceParameter != null) {
142        this.covarianceParameter = (double[])original.covarianceParameter.Clone();
143      }
144
145      // shallow copies of arrays because they cannot be modified
146      this.trainingRows = original.trainingRows;
147      this.allowedInputVariables = original.allowedInputVariables;
148      this.alpha = original.alpha;
149      this.l = original.l;
150      this.x = original.x;
151    }
152    public GaussianProcessModel(IDataset ds, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
153      IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction,
154      bool scaleInputs = true)
155      : base(targetVariable) {
156      this.name = ItemName;
157      this.description = ItemDescription;
158      this.meanFunction = (IMeanFunction)meanFunction.Clone();
159      this.covarianceFunction = (ICovarianceFunction)covarianceFunction.Clone();
160      this.allowedInputVariables = allowedInputVariables.ToArray();
161
162
163      int nVariables = this.allowedInputVariables.Length;
164      meanParameter = hyp
165        .Take(this.meanFunction.GetNumberOfParameters(nVariables))
166        .ToArray();
167
168      covarianceParameter = hyp.Skip(this.meanFunction.GetNumberOfParameters(nVariables))
169                                             .Take(this.covarianceFunction.GetNumberOfParameters(nVariables))
170                                             .ToArray();
171      sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
172      try {
173        CalculateModel(ds, rows, scaleInputs);
174      } catch (alglib.alglibexception ae) {
175        // wrap exception so that calling code doesn't have to know about alglib implementation
176        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
177      }
178    }
179
180    private void CalculateModel(IDataset ds, IEnumerable<int> rows, bool scaleInputs = true) {
181      this.trainingDataset = (IDataset)ds.Clone();
182      this.trainingRows = rows.ToArray();
183      this.inputScaling = scaleInputs ? new Scaling(ds, allowedInputVariables, rows) : null;
184
185      x = GetData(ds, this.allowedInputVariables, this.trainingRows, this.inputScaling);
186
187      IEnumerable<double> y;
188      y = ds.GetDoubleValues(TargetVariable, rows);
189
190      int n = x.GetLength(0);
191
192      var columns = Enumerable.Range(0, x.GetLength(1)).ToArray();
193      // calculate cholesky decomposed (lower triangular) covariance matrix
194      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
195      this.l = CalculateL(x, cov, sqrSigmaNoise);
196
197      // calculate mean
198      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
199      double[] m = Enumerable.Range(0, x.GetLength(0))
200        .Select(r => mean.Mean(x, r))
201        .ToArray();
202
203      // calculate sum of diagonal elements for likelihood
204      double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
205
206      // solve for alpha
207      double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
208
209      int info;
210      alglib.densesolverreport denseSolveRep;
211
212      alglib.spdmatrixcholeskysolve(l, n, false, ym, out info, out denseSolveRep, out alpha);
213      for (int i = 0; i < alpha.Length; i++)
214        alpha[i] = alpha[i] / sqrSigmaNoise;
215      negativeLogLikelihood = 0.5 * Util.ScalarProd(ym, alpha) + diagSum + (n / 2.0) * Math.Log(2.0 * Math.PI * sqrSigmaNoise);
216
217      // derivatives
218      int nAllowedVariables = x.GetLength(1);
219
220      alglib.matinvreport matInvRep;
221      double[,] lCopy = new double[l.GetLength(0), l.GetLength(1)];
222      Array.Copy(l, lCopy, lCopy.Length);
223
224      alglib.spdmatrixcholeskyinverse(ref lCopy, n, false, out info, out matInvRep);
225      if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");
226
227      // LOOCV log pseudo-likelihood (or log predictive probability) (GPML page 116 and 117)
228      var sumLoo = 0.0;
229      var ki = new double[n];
230      for (int i = 0; i < n; i++) {
231        for (int j = 0; j < n; j++) ki[j] = cov.Covariance(x, i, j);
232        ki[i] += sqrSigmaNoise;
233
234        var yi = Util.ScalarProd(ki, alpha);
235        var yi_loo = yi - alpha[i] / (lCopy[i, i] / sqrSigmaNoise);
236        var s2_loo = 1.0 / (lCopy[i, i] / sqrSigmaNoise);
237        var err = ym[i] - yi_loo;
238        var nll_loo = 0.5 * Math.Log(2 * Math.PI * s2_loo) + 0.5 * err * err / s2_loo;
239        sumLoo += nll_loo;
240      }
241      loocvNegLogPseudoLikelihood = sumLoo;
242
243      for (int i = 0; i < n; i++) {
244        for (int j = 0; j <= i; j++)
245          lCopy[i, j] = lCopy[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
246      }
247
248      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => lCopy[i, i]).Sum();
249
250      double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
251      for (int k = 0; k < meanGradients.Length; k++) {
252        var meanGrad = new double[alpha.Length];
253        for (int g = 0; g < meanGrad.Length; g++)
254          meanGrad[g] = mean.Gradient(x, g, k);
255        meanGradients[k] = -Util.ScalarProd(meanGrad, alpha);
256      }
257
258      double[] covGradients = new double[covarianceFunction.GetNumberOfParameters(nAllowedVariables)];
259      if (covGradients.Length > 0) {
260        for (int i = 0; i < n; i++) {
261          for (int j = 0; j < i; j++) {
262            var g = cov.CovarianceGradient(x, i, j);
263            for (int k = 0; k < covGradients.Length; k++) {
264              covGradients[k] += lCopy[i, j] * g[k];
265            }
266          }
267
268          var gDiag = cov.CovarianceGradient(x, i, i);
269          for (int k = 0; k < covGradients.Length; k++) {
270            // diag
271            covGradients[k] += 0.5 * lCopy[i, i] * gDiag[k];
272          }
273        }
274      }
275
276      hyperparameterGradients =
277        meanGradients
278        .Concat(covGradients)
279        .Concat(new double[] { noiseGradient }).ToArray();
280
281    }
282
283    private static double[,] GetData(IDataset ds, IEnumerable<string> allowedInputs, IEnumerable<int> rows, Scaling scaling) {
284      if (scaling != null) {
285        // BackwardsCompatibility3.3
286        #region Backwards compatible code, remove with 3.4
287        // TODO: completely remove Scaling class
288        List<string> variablesList = allowedInputs.ToList();
289        List<int> rowsList = rows.ToList();
290
291        double[,] matrix = new double[rowsList.Count, variablesList.Count];
292
293        int col = 0;
294        foreach (string column in variablesList) {
295          var values = scaling.GetScaledValues(ds, column, rowsList);
296          int row = 0;
297          foreach (var value in values) {
298            matrix[row, col] = value;
299            row++;
300          }
301          col++;
302        }
303        return matrix;
304        #endregion
305      } else {
306        return ds.ToArray(allowedInputs, rows);
307      }
308    }
309
310    private static double[,] CalculateL(double[,] x, ParameterizedCovarianceFunction cov, double sqrSigmaNoise) {
311      int n = x.GetLength(0);
312      var l = new double[n, n];
313
314      // calculate covariances
315      for (int i = 0; i < n; i++) {
316        for (int j = i; j < n; j++) {
317          l[j, i] = cov.Covariance(x, i, j) / sqrSigmaNoise;
318          if (j == i) l[j, i] += 1.0;
319        }
320      }
321
322      // cholesky decomposition
323      var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
324      if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
325      return l;
326    }
327
328
329    public override IDeepCloneable Clone(Cloner cloner) {
330      return new GaussianProcessModel(this, cloner);
331    }
332
333    // is called by the solution creator to set all parameter values of the covariance and mean function
334    // to the optimized values (necessary to make the values visible in the GUI)
335    public void FixParameters() {
336      covarianceFunction.SetParameter(covarianceParameter);
337      meanFunction.SetParameter(meanParameter);
338      covarianceParameter = new double[0];
339      meanParameter = new double[0];
340    }
341
342    #region IRegressionModel Members
343    public override IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
344      return GetEstimatedValuesHelper(dataset, rows);
345    }
346    public override IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
347      return new GaussianProcessRegressionSolution(this, new RegressionProblemData(problemData));
348    }
349    #endregion
350
351
352    private IEnumerable<double> GetEstimatedValuesHelper(IDataset dataset, IEnumerable<int> rows) {
353      try {
354        if (x == null) {
355          x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
356        }
357        int n = x.GetLength(0);
358
359        double[,] newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
360        int newN = newX.GetLength(0);
361
362        var Ks = new double[newN][];
363        var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
364        var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
365        var ms = Enumerable.Range(0, newX.GetLength(0))
366        .Select(r => mean.Mean(newX, r))
367        .ToArray();
368        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
369        for (int i = 0; i < newN; i++) {
370          Ks[i] = new double[n];
371          for (int j = 0; j < n; j++) {
372            Ks[i][j] = cov.CrossCovariance(x, newX, j, i);
373          }
374        }
375
376        return Enumerable.Range(0, newN)
377          .Select(i => ms[i] + Util.ScalarProd(Ks[i], alpha));
378      } catch (alglib.alglibexception ae) {
379        // wrap exception so that calling code doesn't have to know about alglib implementation
380        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
381      }
382    }
383
384    public IEnumerable<double> GetEstimatedVariances(IDataset dataset, IEnumerable<int> rows) {
385      try {
386        if (x == null) {
387          x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
388        }
389        int n = x.GetLength(0);
390
391        var newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
392        int newN = newX.GetLength(0);
393
394        var kss = new double[newN];
395        double[,] sWKs = new double[n, newN];
396        var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
397        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
398
399        if (l == null) {
400          l = CalculateL(x, cov, sqrSigmaNoise);
401        }
402
403        // for stddev
404        for (int i = 0; i < newN; i++)
405          kss[i] = cov.Covariance(newX, i, i);
406
407        for (int i = 0; i < newN; i++) {
408          for (int j = 0; j < n; j++) {
409            sWKs[j, i] = cov.CrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);
410          }
411        }
412
413        // for stddev
414        alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
415
416        for (int i = 0; i < newN; i++) {
417          var col = Util.GetCol(sWKs, i).ToArray();
418          var sumV = Util.ScalarProd(col, col);
419          kss[i] += sqrSigmaNoise; // kss is V(f), add noise variance of predictive distibution to get V(y)
420          kss[i] -= sumV;
421          if (kss[i] < 0) kss[i] = 0;
422        }
423        return kss;
424      } catch (alglib.alglibexception ae) {
425        // wrap exception so that calling code doesn't have to know about alglib implementation
426        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
427      }
428    }
429
430  }
431}
Note: See TracBrowser for help on using the repository browser.