source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/CovarianceFunctions/CovarianceSpectralMixture.cs @ 13784

Last change on this file since 13784 was 13784, checked in by pfleck, 5 years ago

#2591 Made the creation of a GaussianProcessModel faster by avoiding additional iterators during calculation of the hyperparameter gradients.
The gradients of the hyperparameters are now calculated in one sweep and returned as IList, instead of returning an iterator (with yield return).
This avoids a large amount of Move-calls of the iterator, especially for covariance functions with a lot of hyperparameters.
Besides, the signature of the CovarianceGradientFunctionDelegate is changed, to return an IList instead of an IEnumerable to avoid unnececary ToList or ToArray calls.

File size: 10.3 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 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.Data;
28using HeuristicLab.Parameters;
29using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
30
31namespace HeuristicLab.Algorithms.DataAnalysis {
32  [StorableClass]
33  [Item(Name = "CovarianceSpectralMixture",
34    Description = "The spectral mixture kernel described in Wilson A. G. and Adams R.P., Gaussian Process Kernels for Pattern Discovery and Exptrapolation, ICML 2013.")]
35  public sealed class CovarianceSpectralMixture : ParameterizedNamedItem, ICovarianceFunction {
36    public const string QParameterName = "Number of components (Q)";
37    public const string WeightParameterName = "Weight";
38    public const string FrequencyParameterName = "Component frequency (mu)";
39    public const string LengthScaleParameterName = "Length scale (nu)";
40    public IValueParameter<IntValue> QParameter {
41      get { return (IValueParameter<IntValue>)Parameters[QParameterName]; }
42    }
43
44    public IValueParameter<DoubleArray> WeightParameter {
45      get { return (IValueParameter<DoubleArray>)Parameters[WeightParameterName]; }
46    }
47    public IValueParameter<DoubleArray> FrequencyParameter {
48      get { return (IValueParameter<DoubleArray>)Parameters[FrequencyParameterName]; }
49    }
50
51    public IValueParameter<DoubleArray> LengthScaleParameter {
52      get { return (IValueParameter<DoubleArray>)Parameters[LengthScaleParameterName]; }
53    }
54
55    private bool HasFixedWeightParameter {
56      get { return WeightParameter.Value != null; }
57    }
58    private bool HasFixedFrequencyParameter {
59      get { return FrequencyParameter.Value != null; }
60    }
61    private bool HasFixedLengthScaleParameter {
62      get { return LengthScaleParameter.Value != null; }
63    }
64
65    [StorableConstructor]
66    private CovarianceSpectralMixture(bool deserializing)
67      : base(deserializing) {
68    }
69
70    private CovarianceSpectralMixture(CovarianceSpectralMixture original, Cloner cloner)
71      : base(original, cloner) {
72    }
73
74    public CovarianceSpectralMixture()
75      : base() {
76      Name = ItemName;
77      Description = ItemDescription;
78      Parameters.Add(new ValueParameter<IntValue>(QParameterName, "The number of Gaussians (Q) to use for the spectral mixture.", new IntValue(10)));
79      Parameters.Add(new OptionalValueParameter<DoubleArray>(WeightParameterName, "The weight of the component w (peak height of the Gaussian in spectrum)."));
80      Parameters.Add(new OptionalValueParameter<DoubleArray>(FrequencyParameterName, "The inverse component period parameter mu_q (location of the Gaussian in spectrum)."));
81      Parameters.Add(new OptionalValueParameter<DoubleArray>(LengthScaleParameterName, "The length scale parameter (nu_q) (variance of the Gaussian in the spectrum)."));
82    }
83
84    public override IDeepCloneable Clone(Cloner cloner) {
85      return new CovarianceSpectralMixture(this, cloner);
86    }
87
88    public int GetNumberOfParameters(int numberOfVariables) {
89      var q = QParameter.Value.Value;
90      return
91        (HasFixedWeightParameter ? 0 : q) +
92        (HasFixedFrequencyParameter ? 0 : q * numberOfVariables) +
93        (HasFixedLengthScaleParameter ? 0 : q * numberOfVariables);
94    }
95
96    public void SetParameter(double[] p) {
97      double[] weight, frequency, lengthScale;
98      GetParameterValues(p, out weight, out frequency, out lengthScale);
99      WeightParameter.Value = new DoubleArray(weight);
100      FrequencyParameter.Value = new DoubleArray(frequency);
101      LengthScaleParameter.Value = new DoubleArray(lengthScale);
102    }
103
104
105    private void GetParameterValues(double[] p, out double[] weight, out double[] frequency, out double[] lengthScale) {
106      // gather parameter values
107      int c = 0;
108      int q = QParameter.Value.Value;
109      // guess number of elements for frequency and length (=q * numberOfVariables)
110      int n = WeightParameter.Value == null ? ((p.Length - q) / 2) : (p.Length / 2);
111      if (HasFixedWeightParameter) {
112        weight = WeightParameter.Value.ToArray();
113      } else {
114        weight = p.Skip(c).Select(Math.Exp).Take(q).ToArray();
115        c += q;
116      }
117      if (HasFixedFrequencyParameter) {
118        frequency = FrequencyParameter.Value.ToArray();
119      } else {
120        frequency = p.Skip(c).Select(Math.Exp).Take(n).ToArray();
121        c += n;
122      }
123      if (HasFixedLengthScaleParameter) {
124        lengthScale = LengthScaleParameter.Value.ToArray();
125      } else {
126        lengthScale = p.Skip(c).Select(Math.Exp).Take(n).ToArray();
127        c += n;
128      }
129      if (p.Length != c) throw new ArgumentException("The length of the parameter vector does not match the number of free parameters for CovarianceSpectralMixture", "p");
130    }
131
132    public ParameterizedCovarianceFunction GetParameterizedCovarianceFunction(double[] p, int[] columnIndices) {
133      double[] weight, frequency, lengthScale;
134      GetParameterValues(p, out weight, out frequency, out lengthScale);
135      var fixedWeight = HasFixedWeightParameter;
136      var fixedFrequency = HasFixedFrequencyParameter;
137      var fixedLengthScale = HasFixedLengthScaleParameter;
138      // create functions
139      var cov = new ParameterizedCovarianceFunction();
140      cov.Covariance = (x, i, j) => {
141        return GetCovariance(x, x, i, j, QParameter.Value.Value, weight, frequency,
142                             lengthScale, columnIndices);
143      };
144      cov.CrossCovariance = (x, xt, i, j) => {
145        return GetCovariance(x, xt, i, j, QParameter.Value.Value, weight, frequency,
146                             lengthScale, columnIndices);
147      };
148      cov.CovarianceGradient = (x, i, j) => GetGradient(x, i, j, QParameter.Value.Value, weight, frequency,
149                             lengthScale, columnIndices, fixedWeight, fixedFrequency, fixedLengthScale);
150      return cov;
151    }
152
153    private static double GetCovariance(double[,] x, double[,] xt, int i, int j, int maxQ, double[] weight, double[] frequency, double[] lengthScale, int[] columnIndices) {
154      // tau = x - x' (only for selected variables)
155      double[] tau =
156        Util.GetRow(x, i, columnIndices).Zip(Util.GetRow(xt, j, columnIndices), (xi, xj) => xi - xj).ToArray();
157      int numberOfVariables = lengthScale.Length / maxQ;
158      double k = 0;
159      // for each component
160      for (int q = 0; q < maxQ; q++) {
161        double kc = weight[q]; // weighted kernel component
162
163        int idx = 0; // helper index for tau
164        // for each selected variable
165        for (int c = 0; c < columnIndices.Length; c++) {
166          var col = columnIndices[c];
167          kc *= f1(tau[idx], lengthScale[q * numberOfVariables + col]) * f2(tau[idx], frequency[q * numberOfVariables + col]);
168          idx++;
169        }
170        k += kc;
171      }
172      return k;
173    }
174
175    public static double f1(double tau, double lengthScale) {
176      return Math.Exp(-2 * Math.PI * Math.PI * tau * tau * lengthScale);
177    }
178    public static double f2(double tau, double frequency) {
179      return Math.Cos(2 * Math.PI * tau * frequency);
180    }
181
182    // order of returned gradients must match the order in GetParameterValues!
183    private static IList<double> GetGradient(double[,] x, int i, int j, int maxQ, double[] weight, double[] frequency, double[] lengthScale, int[] columnIndices,
184      bool fixedWeight, bool fixedFrequency, bool fixedLengthScale) {
185      double[] tau = Util.GetRow(x, i, columnIndices).Zip(Util.GetRow(x, j, columnIndices), (xi, xj) => xi - xj).ToArray();
186      int numberOfVariables = lengthScale.Length / maxQ;
187
188      var g = new List<double>((!fixedWeight ? maxQ : 0) + (!fixedFrequency ? maxQ * columnIndices.Length : 0) + (!fixedLengthScale ? maxQ * columnIndices.Length : 0));
189      if (!fixedWeight) {
190        // weight
191        // for each component
192        for (int q = 0; q < maxQ; q++) {
193          double k = weight[q];
194          int idx = 0; // helper index for tau
195          // for each selected variable
196          for (int c = 0; c < columnIndices.Length; c++) {
197            var col = columnIndices[c];
198            k *= f1(tau[idx], lengthScale[q * numberOfVariables + col]) * f2(tau[idx], frequency[q * numberOfVariables + col]);
199            idx++;
200          }
201          g.Add(k);
202        }
203      }
204
205      if (!fixedFrequency) {
206        // frequency
207        // for each component
208        for (int q = 0; q < maxQ; q++) {
209          int idx = 0; // helper index for tau
210          // for each selected variable
211          foreach (var c in columnIndices) {
212            double k = f1(tau[idx], lengthScale[q * numberOfVariables + c]) *
213                       -2 * Math.PI * tau[idx] * frequency[q * numberOfVariables + c] *
214                       Math.Sin(2 * Math.PI * tau[idx] * frequency[q * numberOfVariables + c]);
215            idx++;
216            g.Add(weight[q] * k);
217          }
218        }
219      }
220
221      if (!fixedLengthScale) {
222        // length scale
223        // for each component
224        for (int q = 0; q < maxQ; q++) {
225          int idx = 0; // helper index for tau
226          // for each selected variable
227          foreach (var c in columnIndices) {
228            double k = -2 * Math.PI * Math.PI * tau[idx] * tau[idx] * lengthScale[q * numberOfVariables + c] *
229                       f1(tau[idx], lengthScale[q * numberOfVariables + c]) *
230                       f2(tau[idx], frequency[q * numberOfVariables + c]);
231            idx++;
232            g.Add(weight[q] * k);
233          }
234        }
235      }
236
237      return g;
238    }
239  }
240}
Note: See TracBrowser for help on using the repository browser.