Cholesky 分解在 Java 中生成 NaN

Question

我不确定这是 maths.se 还是 SO 问题，但我会选择 SO，因为我认为它与我的 Java 有关。

我正在关注一本关于高斯过程 (R&W) 的教科书，并在 Java 中实现了一些示例。几个示例的一个常见步骤是生成协方差矩阵的 Cholesky 分解。在我的尝试中，我可以获得最大尺寸有限（33x33）的矩阵的成功结果。然而，对于任何更大的 NaN 出现在对角线（在 32,32 处），因此矩阵中的所有后续值同样是 NaN。

代码如下所示，在cholesky方法中注明了NaN的来源。本质上，协方差元素 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第201页：

在实践中，可能需要添加一个小的倍数 单位矩阵 $\epsilon I$ 到数值的协方差矩阵 原因。这是因为矩阵 K 的特征值可以衰减 非常迅速 [...] 并且没有这种稳定，Cholesky 分解失败。对生成的样本的影响是添加 额外的独立方差噪声$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 编写了自己的 Cholesky 分解例程版本。它不是计算 L，而是计算 U，但我很想用导致 NaN 错误的矩阵对其进行测试。您能在此处发布矩阵吗，或与我联系（个人资料中的信息。）