Changeset 10422 for branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression

Timestamp:

01/29/14 14:34:44 (11 years ago)

Author:

gkronber

Message:

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

Location:

branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression

Files:

• Instances.DataAnalysis.GaussianProcessRegression-3.3.csproj (modified) (1 diff)
• Plugin.cs.frame (modified) (1 diff)
• Util.cs (modified) (2 diffs)

branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Instances.DataAnalysis.GaussianProcessRegression-3.3.csproj

@@ -98 +98 @@
   </PropertyGroup>
   <ItemGroup>
-    <Reference Include="ALGLIB-3.6.0, Version=3.6.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
-      <SpecificVersion>False</SpecificVersion>
-      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.6.0.dll</HintPath>
-      <Private>False</Private>
+    <Reference Include="ALGLIB-3.7.0, Version=3.7.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.7.0.dll</HintPath>
     </Reference>
     <Reference Include="HeuristicLab.Algorithms.DataAnalysis-3.4, Version=3.4.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">

branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Plugin.cs.frame

@@ -25 +25 @@
   [Plugin("HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression", "3.3.7.$WCREV$")]
   [PluginFile("HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression-3.3.dll", PluginFileType.Assembly)]
-  [PluginDependency("HeuristicLab.ALGLIB", "3.6.0")]
+  [PluginDependency("HeuristicLab.ALGLIB", "3.7.0")]
   [PluginDependency("HeuristicLab.Algorithms.DataAnalysis", "3.4")]
   [PluginDependency("HeuristicLab.Collections", "3.3")]

branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.Instances.DataAnalysis.GaussianProcessRegression/Util.cs

@@ -52 +52 @@
           K[i, j] = covFunction.Covariance(x, i, j);
 
-      if (!alglib.spdmatrixcholesky(ref K, K.GetLength(0), true)) throw new ArgumentException();
-
+      // if (!alglib.spdmatrixcholesky(ref K, K.GetLength(0), true)) throw new ArgumentException();
+      K = toeplitz_cholesky_lower(K.GetLength(0), K);
       List<double> target = new List<double>(K.GetLength(0));
       for (int i = 0; i < K.GetLength(0); i++) {
@@ -67 +67 @@
       return target;
     }
+
+    //****************************************************************************80
+
+    public static double[,] toeplitz_cholesky_lower(int n, double[,] a)
+
+    //****************************************************************************80
+      //
+      //  Purpose:
+      //
+      //    TOEPLITZ_CHOLESKY_LOWER: lower Cholesky factor of a Toeplitz matrix.
+      //
+      //  Discussion:
+      //
+      //    The Toeplitz matrix must be positive semi-definite.
+      //
+      //    After factorization, A = L * L'.
+      //
+      //  Licensing:
+      //
+      //    This code is distributed under the GNU LGPL license.
+      //
+      //  Modified:
+      //
+      //    13 November 2012
+      //    29 January  2014: adapted to C# by Gabriel Kronberger
+      //  Author:
+      //
+      //    John Burkardt
+      //
+      //  Reference:
+      //
+      //    Michael Stewart,
+      //    Cholesky factorization of semi-definite Toeplitz matrices.
+      //
+      //  Parameters:
+      //
+      //    Input, int N, the order of the matrix.
+      //
+      //    Input, double A[N,N], the Toeplitz matrix.
+      //
+      //    Output, double TOEPLITZ_CHOLESKY_LOWER[N,N], the lower Cholesky factor.
+      //
+    {
+      double div;
+      double[] g;
+      double g1j;
+      double g2j;
+      int i;
+      int j;
+      double[,] l;
+      double rho;
+
+      l = new double[n, n];
+
+      for (j = 0; j < n; j++) {
+        for (i = 0; i < n; i++) {
+          l[i, j] = 0.0;
+        }
+      }
+      g = new double[2 * n];
+
+      for (j = 0; j < n; j++) {
+        g[0 + j * 2] = a[0, j];
+      }
+      g[1 + 0 * 2] = 0.0;
+      for (j = 1; j < n; j++) {
+        g[1 + j * 2] = a[j, 0];
+      }
+
+      for (i = 0; i < n; i++) {
+        l[i, 0] = g[0 + i * 2];
+      }
+      for (j = n - 1; 1 <= j; j--) {
+        g[0 + j * 2] = g[0 + (j - 1) * 2];
+      }
+      g[0 + 0 * 2] = 0.0;
+
+      for (i = 1; i < n; i++) {
+        rho = -g[1 + i * 2] / g[0 + i * 2];
+        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));
+
+        for (j = i; j < n; j++) {
+          g1j = g[0 + j * 2];
+          g2j = g[1 + j * 2];
+          g[0 + j * 2] = (g1j + rho * g2j) / div;
+          g[1 + j * 2] = (rho * g1j + g2j) / div;
+        }
+
+        for (j = i; j < n; j++) {
+          l[j, i] = g[0 + j * 2];
+        }
+        for (j = n - 1; i < j; j--) {
+          g[0 + j * 2] = g[0 + (j - 1) * 2];
+        }
+        g[0 + i * 2] = 0.0;
+      }
+
+
+      return l;
+    }
+    //****************************************************************************80
+
+    public static double[,] toeplitz_cholesky_upper(int n, double[,] a)
+
+    //****************************************************************************80
+      //
+      //  Purpose:
+      //
+      //    TOEPLITZ_CHOLESKY_UPPER: upper Cholesky factor of a Toeplitz matrix.
+      //
+      //  Discussion:
+      //
+      //    The Toeplitz matrix must be positive semi-definite.
+      //
+      //    After factorization, A = R' * R.
+      //
+      //  Licensing:
+      //
+      //    This code is distributed under the GNU LGPL license.
+      //
+      //  Modified:
+      //
+      //    14 November 2012
+      //    29 January  2014: adapted to C# by Gabriel Kronberger
+      //
+      //  Author:
+      //
+      //    John Burkardt
+      //
+      //  Reference:
+      //
+      //    Michael Stewart,
+      //    Cholesky factorization of semi-definite Toeplitz matrices.
+      //
+      //  Parameters:
+      //
+      //    Input, int N, the order of the matrix.
+      //
+      //    Input, double A[N,N], the Toeplitz matrix.
+      //
+      //    Output, double TOEPLITZ_CHOLESKY_UPPER[N,N], the upper Cholesky factor.
+      //
+    {
+      double div;
+      double[] g;
+      double g1j;
+      double g2j;
+      int i;
+      int j;
+      double[,] r;
+      double rho;
+
+      r = new double[n, n];
+
+      for (j = 0; j < n; j++) {
+        for (i = 0; i < n; i++) {
+          r[i, j] = 0.0;
+        }
+      }
+      g = new double[2 * n];
+
+      for (j = 0; j < n; j++) {
+        g[0 + j * 2] = a[0, j];
+      }
+
+      g[1 + 0 * 2] = 0.0;
+      for (j = 1; j < n; j++) {
+        g[1 + j * 2] = a[j, 0];
+      }
+      for (j = 0; j < n; j++) {
+        r[0, j] = g[0 + j * 2];
+      }
+      for (j = n - 1; 1 <= j; j--) {
+        g[0 + j * 2] = g[0 + (j - 1) * 2];
+      }
+      g[0 + 0 * 2] = 0.0;
+
+      for (i = 1; i < n; i++) {
+        rho = -g[1 + i * 2] / g[0 + i * 2];
+        div = Math.Sqrt((1.0 - rho) * (1.0 + rho));
+        for (j = i; j < n; j++) {
+          g1j = g[0 + j * 2];
+          g2j = g[1 + j * 2];
+          g[0 + j * 2] = (g1j + rho * g2j) / div;
+          g[1 + j * 2] = (rho * g1j + g2j) / div;
+        }
+        for (j = i; j < n; j++) {
+          r[i, j] = g[0 + j * 2];
+        }
+        for (j = n - 1; i < j; j--) {
+          g[0 + j * 2] = g[0 + (j - 1) * 2];
+        }
+        g[0 + i * 2] = 0.0;
+      }
+
+      return r;
+    }
   }
 }