Cholesky分解在python，遞歸

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()