为什么 cholesky 分解对于大尺寸的帕斯卡矩阵会失败？答案

【问题标题】：Why does cholesky factorization fail for pascal matrices of big size?为什么 cholesky 分解对于大尺寸的帕斯卡矩阵会失败？
【发布时间】：2022-01-06 05:54:01
【问题描述】：

我想对大小为 50 的帕斯卡矩阵执行 Cholesky 分解。值变得太大，所以 scipy.linalg.pascal 返回 object 类型的矩阵。

A_scipy = scipy.linalg.pascal(50)
A_scipy.dtype
>dtype('O')

如果手工构建：

def construct_pascal_triangle(n):
    L = np.zeros((n, n), dtype=np.float64)
    L[:, 0] = 1
    for i in range(1,n):
        for j in range(1,i+1):
            L[i][j] = L[i-1][j] + L[i-1][j-1]
    return L
L = construct_pascal_triangle(n)
A = L @ L.T

那么它不同于A_scipy。我假设np.float64 也无法处理它，所以当我在函数construct_pascal_triangle 中将dtype 转换为object 时，A 和A_scipy 重合。 np.linalg.cholesky 无法处理 object 类型矩阵。所以我写了自己的函数

def cholesky(A):
    n = A.shape[0]
    M = A.copy()
    L = np.zeros_like(A)
    for i in range(n):
        L[i, i] = M[i, i] ** 0.5
        L[i, i + 1:] = M[i, i + 1:] / L[i, i]
        for j in range(i + 1, n):
            M[j, j:] = M[j, j:] - L[i, j] * L[i, j:]

但它也失败了，因为M[i, i] 在某些时候变成了负数。我想知道为什么会这样。帕斯卡矩阵对于任何大小都是正定的，因此 cholesky 分解总是存在的。是不是类型已经有问题了，而且数字太大了，甚至对象都无法处理它们？或者这是别的什么？

【问题讨论】：

标签： python numpy floating-point

【解决方案1】：

你认为M[i,i] ** 0.5 应用于长整数时会做什么？

如果你发现它将结果转换为最接近的 float64，那么你会得到一个不精确的解释......

float64 使用不精确的算术（算术运算的结果必须四舍五入到最接近的 float64）。

将数字与精确算术（至少对于 + - *）和不精确算术混合会导致算术不精确。所以这个** 步骤失去了额外的精度。

即使M[i,i] ** 0.5本身是一个精确运算（应该是n=50），像L[i, j] * L[i, j:]这样的下面的运算也可能会超过float64的53bits精度，并且不精确。

当平方根是精确整数时，您必须找到（或编写）保留整数结果的 sqrt 版本。（例如，Squeak Smalltalk sqrt 就是这样做的）。

【讨论】：