Cholesky 分解在 Java 中生成 NaN

Question

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

有关

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

代码如下，NaN的来源在cholesky方法中注明。本质上协方差元素a[32][32]是1.0，但是sum的值稍微超过这个（1.0000001423291431），所以平方根是虚数。所以我的问题是：

这是线性代数的预期结果吗，或者，例如，我实施的人工制品？
实践中如何最好地避免这个问题？

请注意，我不是在寻找要使用的库的建议。这只是我自己的理解。

抱歉篇幅太长，但我已尝试提供完整的 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 页：

In practice it may be necessary to add a small multiple of the identity matrix $\epsilon I$ to the covariance matrix for numerical reasons. This is because the eigenvalues of the matrix K can decay very rapidly [...] and without this stabilization the Cholesky decomposition fails. The effect on the generated samples is to add additional independent noise of variance $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 错误的矩阵对其进行测试。你能post这里的矩阵，或者联系我（个人资料中的信息。）

Cholesky 分解在 Java 中生成 NaN

Cholesky decomposition generating NaNs in Java

java

linear-algebra

matrix-factorization