Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Util.cs @ 10662

Last change on this file since 10662 was 10422, checked in by gkronber, 11 years ago

#1967 added Cholesky decomposition for Toeplitz matrices to allow sampling from one-dimensional Gaussian processes

File size: 7.6 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.Algorithms.DataAnalysis;
26using HeuristicLab.Core;
27using HeuristicLab.Random;
28
29namespace HeuristicLab.Problems.Instances.DataAnalysis {
30  public class Util {
31
32
33    public static List<double> SampleGaussianProcess(IRandom random, ParameterizedCovarianceFunction covFunction, List<List<double>> data) {
34      int n = data[0].Count;
35
36      var normalRand = new NormalDistributedRandom(random, 0, 1);
37      var alpha = (from i in Enumerable.Range(0, n)
38                   select normalRand.NextDouble()).ToArray();
39      return SampleGaussianProcess(random, covFunction, data, alpha);
40    }
41
42    public static List<double> SampleGaussianProcess(IRandom random, ParameterizedCovarianceFunction covFunction, List<List<double>> data, double[] alpha) {
43      if (alpha.Length != data[0].Count) throw new ArgumentException();
44
45      double[,] x = new double[data[0].Count, data.Count];
46      for (int i = 0; i < x.GetLength(0); i++)
47        for (int j = 0; j < x.GetLength(1); j++)
48          x[i, j] = data[j][i];
49      double[,] K = new double[x.GetLength(0), x.GetLength(0)];
50      for (int i = 0; i < K.GetLength(0); i++)
51        for (int j = i; j < K.GetLength(1); j++)
52          K[i, j] = covFunction.Covariance(x, i, j);
53
54      // if (!alglib.spdmatrixcholesky(ref K, K.GetLength(0), true)) throw new ArgumentException();
55      K = toeplitz_cholesky_lower(K.GetLength(0), K);
56      List<double> target = new List<double>(K.GetLength(0));
57      for (int i = 0; i < K.GetLength(0); i++) {
58        double s = 0.0;
59        for (int j = K.GetLength(0) - 1; j >= 0; j--) {
60
61          s += K[j, i] * alpha[j];
62        }
63
64        target.Add(s);
65      }
66
67      return target;
68    }
69
70    //****************************************************************************80
71
72    public static double[,] toeplitz_cholesky_lower(int n, double[,] a)
73
74    //****************************************************************************80
75      //
76      //  Purpose:
77      //
78      //    TOEPLITZ_CHOLESKY_LOWER: lower Cholesky factor of a Toeplitz matrix.
79      //
80      //  Discussion:
81      //
82      //    The Toeplitz matrix must be positive semi-definite.
83      //
84      //    After factorization, A = L * L'.
85      //
86      //  Licensing:
87      //
88      //    This code is distributed under the GNU LGPL license.
89      //
90      //  Modified:
91      //
92      //    13 November 2012
93      //    29 January  2014: adapted to C# by Gabriel Kronberger
94      //  Author:
95      //
96      //    John Burkardt
97      //
98      //  Reference:
99      //
100      //    Michael Stewart,
101      //    Cholesky factorization of semi-definite Toeplitz matrices.
102      //
103      //  Parameters:
104      //
105      //    Input, int N, the order of the matrix.
106      //
107      //    Input, double A[N,N], the Toeplitz matrix.
108      //
109      //    Output, double TOEPLITZ_CHOLESKY_LOWER[N,N], the lower Cholesky factor.
110      //
111    {
112      double div;
113      double[] g;
114      double g1j;
115      double g2j;
116      int i;
117      int j;
118      double[,] l;
119      double rho;
120
121      l = new double[n, n];
122
123      for (j = 0; j < n; j++) {
124        for (i = 0; i < n; i++) {
125          l[i, j] = 0.0;
126        }
127      }
128      g = new double[2 * n];
129
130      for (j = 0; j < n; j++) {
131        g[0 + j * 2] = a[0, j];
132      }
133      g[1 + 0 * 2] = 0.0;
134      for (j = 1; j < n; j++) {
135        g[1 + j * 2] = a[j, 0];
136      }
137
138      for (i = 0; i < n; i++) {
139        l[i, 0] = g[0 + i * 2];
140      }
141      for (j = n - 1; 1 <= j; j--) {
142        g[0 + j * 2] = g[0 + (j - 1) * 2];
143      }
144      g[0 + 0 * 2] = 0.0;
145
146      for (i = 1; i < n; i++) {
147        rho = -g[1 + i * 2] / g[0 + i * 2];
148        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));
149
150        for (j = i; j < n; j++) {
151          g1j = g[0 + j * 2];
152          g2j = g[1 + j * 2];
153          g[0 + j * 2] = (g1j + rho * g2j) / div;
154          g[1 + j * 2] = (rho * g1j + g2j) / div;
155        }
156
157        for (j = i; j < n; j++) {
158          l[j, i] = g[0 + j * 2];
159        }
160        for (j = n - 1; i < j; j--) {
161          g[0 + j * 2] = g[0 + (j - 1) * 2];
162        }
163        g[0 + i * 2] = 0.0;
164      }
165
166
167      return l;
168    }
169    //****************************************************************************80
170
171    public static double[,] toeplitz_cholesky_upper(int n, double[,] a)
172
173    //****************************************************************************80
174      //
175      //  Purpose:
176      //
177      //    TOEPLITZ_CHOLESKY_UPPER: upper Cholesky factor of a Toeplitz matrix.
178      //
179      //  Discussion:
180      //
181      //    The Toeplitz matrix must be positive semi-definite.
182      //
183      //    After factorization, A = R' * R.
184      //
185      //  Licensing:
186      //
187      //    This code is distributed under the GNU LGPL license.
188      //
189      //  Modified:
190      //
191      //    14 November 2012
192      //    29 January  2014: adapted to C# by Gabriel Kronberger
193      //
194      //  Author:
195      //
196      //    John Burkardt
197      //
198      //  Reference:
199      //
200      //    Michael Stewart,
201      //    Cholesky factorization of semi-definite Toeplitz matrices.
202      //
203      //  Parameters:
204      //
205      //    Input, int N, the order of the matrix.
206      //
207      //    Input, double A[N,N], the Toeplitz matrix.
208      //
209      //    Output, double TOEPLITZ_CHOLESKY_UPPER[N,N], the upper Cholesky factor.
210      //
211    {
212      double div;
213      double[] g;
214      double g1j;
215      double g2j;
216      int i;
217      int j;
218      double[,] r;
219      double rho;
220
221      r = new double[n, n];
222
223      for (j = 0; j < n; j++) {
224        for (i = 0; i < n; i++) {
225          r[i, j] = 0.0;
226        }
227      }
228      g = new double[2 * n];
229
230      for (j = 0; j < n; j++) {
231        g[0 + j * 2] = a[0, j];
232      }
233
234      g[1 + 0 * 2] = 0.0;
235      for (j = 1; j < n; j++) {
236        g[1 + j * 2] = a[j, 0];
237      }
238      for (j = 0; j < n; j++) {
239        r[0, j] = g[0 + j * 2];
240      }
241      for (j = n - 1; 1 <= j; j--) {
242        g[0 + j * 2] = g[0 + (j - 1) * 2];
243      }
244      g[0 + 0 * 2] = 0.0;
245
246      for (i = 1; i < n; i++) {
247        rho = -g[1 + i * 2] / g[0 + i * 2];
248        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));
249        for (j = i; j < n; j++) {
250          g1j = g[0 + j * 2];
251          g2j = g[1 + j * 2];
252          g[0 + j * 2] = (g1j + rho * g2j) / div;
253          g[1 + j * 2] = (rho * g1j + g2j) / div;
254        }
255        for (j = i; j < n; j++) {
256          r[i, j] = g[0 + j * 2];
257        }
258        for (j = n - 1; i < j; j--) {
259          g[0 + j * 2] = g[0 + (j - 1) * 2];
260        }
261        g[0 + i * 2] = 0.0;
262      }
263
264      return r;
265    }
266  }
267}
Note: See TracBrowser for help on using the repository browser.