python - Python 中的 Cholesky 分解，递归

Question

我正在尝试制作一个递归程序来计算 Cholesky factorization ，但输出不是下三角形式。我怎样才能改变它以正确计算它？

def cholesky(A):

    n = np.shape(A)[0]
    A[0,0] = math.sqrt(abs(A[0,0]))

    if n == 1:
        return A
    else:
        A[1:,0] /= A[0,0]
        A[1:,1:] -= np.dot(A[1:,0], (A[1:,0]).T)
        cholesky(A[1:,1:])

Answer 1

我确信有一种更漂亮的方法(特别是numpy.linalg.cholesky)，但这看起来在合理的测试数据上确实可行。我遵循维基百科文章中的符号，并使用他们的示例作为测试数据。

import numpy as np

def cholesky_reduce(A):
    pivot = A[0, 0]
    b = np.mat(A[1:, 0])
    B = A[1:, 1:]
    return B - (b.T * b) / pivot


def L(A):
    n = A.shape[0]
    if n == 1:
        return np.sqrt(A)
    b = np.mat(A[1:, 0])
    pivot = np.sqrt(A[0, 0])
    return np.bmat([
        [np.mat(pivot), np.zeros((1, n - 1))],
        [b.T / pivot, L(cholesky_reduce(A))]
    ])


def __main():
    A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]])
    print(L(A))


if __name__ == '__main__':
    __main()