Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs @ 8484

Last change on this file since 8484 was 8484, checked in by gkronber, 12 years ago

#1902 changed interface for covariance functions to improve readability, fixed several bugs in the covariance functions and in the line chart for Gaussian process models.

File size: 10.1 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2012 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 : NamedItem, IGaussianProcessModel {
37    [Storable]
38    private double negativeLogLikelihood;
39    public double NegativeLogLikelihood {
40      get { return negativeLogLikelihood; }
41    }
42
43    [Storable]
44    private double[] hyperparameterGradients;
45    public double[] HyperparameterGradients {
46      get {
47        var copy = new double[hyperparameterGradients.Length];
48        Array.Copy(hyperparameterGradients, copy, copy.Length);
49        return copy;
50      }
51    }
52
53    [Storable]
54    private ICovarianceFunction covarianceFunction;
55    public ICovarianceFunction CovarianceFunction {
56      get { return covarianceFunction; }
57    }
58    [Storable]
59    private IMeanFunction meanFunction;
60    public IMeanFunction MeanFunction {
61      get { return meanFunction; }
62    }
63    [Storable]
64    private string targetVariable;
65    public string TargetVariable {
66      get { return targetVariable; }
67    }
68    [Storable]
69    private string[] allowedInputVariables;
70    public string[] AllowedInputVariables {
71      get { return allowedInputVariables; }
72    }
73
74    [Storable]
75    private double[] alpha;
76    [Storable]
77    private double sqrSigmaNoise;
78
79    [Storable]
80    private double[,] l;
81
82    [Storable]
83    private double[,] x;
84    [Storable]
85    private Scaling inputScaling;
86
87
88    [StorableConstructor]
89    private GaussianProcessModel(bool deserializing) : base(deserializing) { }
90    private GaussianProcessModel(GaussianProcessModel original, Cloner cloner)
91      : base(original, cloner) {
92      this.meanFunction = cloner.Clone(original.meanFunction);
93      this.covarianceFunction = cloner.Clone(original.covarianceFunction);
94      this.inputScaling = cloner.Clone(original.inputScaling);
95      this.negativeLogLikelihood = original.negativeLogLikelihood;
96      this.targetVariable = original.targetVariable;
97      this.sqrSigmaNoise = original.sqrSigmaNoise;
98
99      // shallow copies of arrays because they cannot be modified
100      this.allowedInputVariables = original.allowedInputVariables;
101      this.alpha = original.alpha;
102      this.l = original.l;
103      this.x = original.x;
104    }
105    public GaussianProcessModel(Dataset ds, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
106      IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction)
107      : base() {
108      this.name = ItemName;
109      this.description = ItemDescription;
110      this.meanFunction = (IMeanFunction)meanFunction.Clone();
111      this.covarianceFunction = (ICovarianceFunction)covarianceFunction.Clone();
112      this.targetVariable = targetVariable;
113      this.allowedInputVariables = allowedInputVariables.ToArray();
114
115
116      int nVariables = this.allowedInputVariables.Length;
117      this.meanFunction.SetParameter(hyp
118        .Take(this.meanFunction.GetNumberOfParameters(nVariables))
119        .ToArray());
120      this.covarianceFunction.SetParameter(hyp.Skip(this.meanFunction.GetNumberOfParameters(nVariables))
121        .Take(this.covarianceFunction.GetNumberOfParameters(nVariables))
122        .ToArray());
123      sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
124
125      CalculateModel(ds, rows);
126    }
127
128    private void CalculateModel(Dataset ds, IEnumerable<int> rows) {
129      inputScaling = new Scaling(ds, allowedInputVariables, rows);
130      x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
131      var y = ds.GetDoubleValues(targetVariable, rows);
132
133      int n = x.GetLength(0);
134      l = new double[n, n];
135
136      meanFunction.SetData(x);
137
138      // calculate means and covariances
139      double[] m = meanFunction.GetMean(x);
140      for (int i = 0; i < n; i++) {
141        for (int j = i; j < n; j++) {
142          l[j, i] = covarianceFunction.GetCovariance(x, i, j) / sqrSigmaNoise;
143          if (j == i) l[j, i] += 1.0;
144        }
145      }
146
147      // cholesky decomposition
148      int info;
149      alglib.densesolverreport denseSolveRep;
150
151      var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
152      if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
153
154      // calculate sum of diagonal elements for likelihood
155      double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
156
157      // solve for alpha
158      double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
159
160      alglib.spdmatrixcholeskysolve(l, n, false, ym, out info, out denseSolveRep, out alpha);
161      for (int i = 0; i < alpha.Length; i++)
162        alpha[i] = alpha[i] / sqrSigmaNoise;
163      negativeLogLikelihood = 0.5 * Util.ScalarProd(ym, alpha) + diagSum + (n / 2.0) * Math.Log(2.0 * Math.PI * sqrSigmaNoise);
164
165      // derivatives
166      int nAllowedVariables = x.GetLength(1);
167
168      alglib.matinvreport matInvRep;
169      double[,] lCopy = new double[l.GetLength(0), l.GetLength(1)];
170      Array.Copy(l, lCopy, lCopy.Length);
171
172      alglib.spdmatrixcholeskyinverse(ref lCopy, n, false, out info, out matInvRep);
173      if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");
174      for (int i = 0; i < n; i++) {
175        for (int j = 0; j <= i; j++)
176          lCopy[i, j] = lCopy[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
177      }
178
179      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => lCopy[i, i]).Sum();
180
181      double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
182      for (int i = 0; i < meanGradients.Length; i++) {
183        var meanGrad = meanFunction.GetGradients(i, x);
184        meanGradients[i] = -Util.ScalarProd(meanGrad, alpha);
185      }
186
187      double[] covGradients = new double[covarianceFunction.GetNumberOfParameters(nAllowedVariables)];
188      if (covGradients.Length > 0) {
189        for (int i = 0; i < n; i++) {
190          for (int j = 0; j < i; j++) {
191            var g = covarianceFunction.GetGradient(x, i, j).ToArray();
192            for (int k = 0; k < covGradients.Length; k++) {
193              covGradients[k] += lCopy[i, j] * g[k];
194            }
195          }
196
197          var gDiag = covarianceFunction.GetGradient(x, i, i).ToArray();
198          for (int k = 0; k < covGradients.Length; k++) {
199            // diag
200            covGradients[k] += 0.5 * lCopy[i, i] * gDiag[k];
201          }
202        }
203      }
204
205      hyperparameterGradients =
206        meanGradients
207        .Concat(covGradients)
208        .Concat(new double[] { noiseGradient }).ToArray();
209
210    }
211
212
213    public override IDeepCloneable Clone(Cloner cloner) {
214      return new GaussianProcessModel(this, cloner);
215    }
216
217    #region IRegressionModel Members
218    public IEnumerable<double> GetEstimatedValues(Dataset dataset, IEnumerable<int> rows) {
219      return GetEstimatedValuesHelper(dataset, rows);
220    }
221    public GaussianProcessRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
222      return new GaussianProcessRegressionSolution(this, problemData);
223    }
224    IRegressionSolution IRegressionModel.CreateRegressionSolution(IRegressionProblemData problemData) {
225      return CreateRegressionSolution(problemData);
226    }
227    #endregion
228
229    private IEnumerable<double> GetEstimatedValuesHelper(Dataset dataset, IEnumerable<int> rows) {
230      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
231      int newN = newX.GetLength(0);
232      int n = x.GetLength(0);
233      var Ks = new double[newN, n];
234      meanFunction.SetData(newX);
235      var ms = meanFunction.GetMean(newX);
236      for (int i = 0; i < newN; i++) {
237        for (int j = 0; j < n; j++) {
238          Ks[i, j] = covarianceFunction.GetCrossCovariance(x, newX, j, i);
239        }
240      }
241
242      return Enumerable.Range(0, newN)
243        .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha));
244    }
245
246    public IEnumerable<double> GetEstimatedVariance(Dataset dataset, IEnumerable<int> rows) {
247      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
248      int newN = newX.GetLength(0);
249      int n = x.GetLength(0);
250
251      var kss = new double[newN];
252      double[,] sWKs = new double[n, newN];
253
254      // for stddev
255      for (int i = 0; i < newN; i++)
256        kss[i] = covarianceFunction.GetCovariance(newX, i, i);
257
258      for (int i = 0; i < newN; i++) {
259        for (int j = 0; j < n; j++) {
260          sWKs[j, i] = covarianceFunction.GetCrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);
261        }
262      }
263
264      // for stddev
265      alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
266
267      for (int i = 0; i < newN; i++) {
268        var sumV = Util.ScalarProd(Util.GetCol(sWKs, i), Util.GetCol(sWKs, i));
269        kss[i] -= sumV;
270        if (kss[i] < 0) kss[i] = 0;
271      }
272      return kss;
273    }
274  }
275}
Note: See TracBrowser for help on using the repository browser.