JavaでNaNを生成するコレスキー分解

Question

これはmaths.seかSOのどちらの質問か分かりませんが、私のJavaに関連していると思います

私はGaussian Processes（R&W）のテキストブックに従い、Javaでいくつかの例を実装しています。いくつかの例の1つの共通ステップは、共分散行列のコレスキー分解を生成することです。私の試みでは、限られたサイズ（33x33）までのマトリックスの成功を得ることができます。しかし、より大きなものについては、NaNが対角線（32,32）に現れます。したがって、その後の行列の値も同様にNaNです。

コードは以下に示され、NaNのソースはcholeskyメソッドで示されます。本質的に共分散要素a[32][32]は1.0ですが、sumの値はこれを少し上回っています（1.0000001423291431）ので、平方根は虚数です。だから、私の質問は以下のとおりです。これは、線形代数から期待される結果、または、例えば、私の実装の アーチファクト

1. ですか？
2. どのようにこの問題を回避するのが最も効果的ですか？

私は使用するライブラリの推奨を探しているわけではありません。これは私の理解のためのものです。

長さについて謝罪が、私は完全にMWEを提供しようとしました：

import static org.junit.Assert.assertFalse; 

import org.junit.Test; 

public class CholeskyTest { 

    @Test 
    public void testCovCholesky() { 
     final int n = 34; // Test passes for n<34 
     final double[] xData = getSpread(-5, 5, n); 
     double[][] cov = covarianceSE(xData); 
     double[][] lower = cholesky(cov); 
     for(int i=0; i<n; ++i) { 
      for(int j=0; j<n; ++j) { 
       assertFalse("NaN at " + i + "," + j, Double.isNaN(lower[i][j])); 
      } 
     } 
    } 

    /** 
    * Generate n evenly space values from min to max inclusive 
    */ 
    private static double[] getSpread(final double min, final double max, final int n) { 
     final double[] values = new double[n]; 
     final double delta = (max - min)/(n - 1); 
     for(int i=0; i<n; ++i) { 
      values[i] = min + i*delta; 
     } 
     return values; 
    } 

    /** 
    * Calculate the covariance matrix for the given observations using 
    * the squared exponential (SE) covariance function. 
    */ 
    private static double[][] covarianceSE (double[] v) { 
     final int m = v.length; 
     double[][] k = new double[m][]; 
     for(int i=0; i<m; ++i) { 
      double vi = v[i]; 
      double row[] = new double[m]; 
      for(int j=0; j<m; ++j) { 
       double dist = vi - v[j]; 
       row[j] = Math.exp(-0.5*dist*dist); 
      } 
      k[i] = row; 
     } 
     return k; 
    } 

    /** 
    * Calculate lower triangular matrix L such that LL^T = A 
    * Using Cholesky decomposition from 
    * https://rosettacode.org/wiki/Cholesky_decomposition#Java 
    */ 
    private static double[][] cholesky(double[][] a) { 
     final int m = a.length; 
     double[][] l = new double[m][m]; 
     for(int i = 0; i< m;i++){ 
      for(int k = 0; k < (i+1); k++){ 
       double sum = 0; 
       for(int j = 0; j < k; j++){ 
        sum += l[i][j] * l[k][j]; 
       } 
       l[i][k] = (i == k) ? Math.sqrt(a[i][i] - sum) : // Source of NaN at 32,32 
        (1.0/l[k][k] * (a[i][k] - sum)); 
      } 
     } 
     return l; 
    } 
} 

Answer 1

うーん、私は、私は、次のされたのと同じ教科書から、自分の質問への答えを見つけたと思います。 R&Wからp.201：

実際には、私が数値 の理由から共分散行列に$ 単位行列の$ \イプシロンの小さな倍数を追加する必要があるかもしれません。これは、行列Kの固有値が を非常に急速に減衰させることができ、この安定化がなければ、コレスキー 分解が失敗するからです。生成されたサンプルへの影響は、 追加の分散ノイズ$ epsilon $のノイズを追加することです。

したがって、次の変更は十分であると思わ：

private static double[][] cholesky(double[][] a) { 
    final int m = a.length; 
    double epsilon = 0.000001; // Small extra noise value 
    double[][] l = new double[m][m]; 
    for(int i = 0; i< m;i++){ 
     for(int k = 0; k < (i+1); k++){ 
      double sum = 0; 
      for(int j = 0; j < k; j++){ 
       sum += l[i][j] * l[k][j]; 
      } 
      l[i][k] = (i == k) ? Math.sqrt(a[i][i]+epsilon - sum) : // Add noise to diagonal values 
       (1.0/l[k][k] * (a[i][k] - sum)); 
     } 
    } 
    return l; 
} 

Answer 2

私はC++とJavaScriptでコレスキー分解ルーチンの私自身のバージョンを書き終えました。 Lを計算する代わりに、Uを計算しますが、NaNエラーの原因となる行列を使ってテストするのは興味があります。