python的简化--Cholesky分解

正在学习Cholesky分解
需要return一个上三角矩阵
因为是第一次接触python,感觉自己的代码写的很复杂,想问问大家该怎么变得更简洁
希望得到类似于这个帖子的第一个代码
https://blog.csdn.net/ReDreamme/article/details/108745255?spm=1001.2014.3001.5501
以下是我自己的代码

def Cholesky(A):
    len=A.shape[0]
    R=np.zeros((len,len))
    for i in range(0,len):
        for j in range(i, len):
            temp=A[i,j]
            k=i-1
            while(k>=0):
                temp-=R[k,j]*R[k,i]
                k=k-1
            if(i==j):
                R[i,j]=np.sqrt(temp)
            else:
                R[i,j]=temp/R[i,i]
    return R

这样?使用np.dot函数代替循环:代码二中使用了np.dot函数来计算矩阵内积,这样可以使代码简洁,运算效率也更高。

利用对称性:Cholesky分解生成的矩阵是上三角矩阵,因此只需要存储一半的数据。代码二中G是一个对称矩阵,但只使用了一半的空间。

使用矩阵列索引:代码二中使用了矩阵列索引,G[j,:i]的意思是G矩阵的第j行的前i列。

优化循环:代码二中的循环只需要从i+1开始,因为G的对称性已经保证了前面的值已经计算过。


import numpy as np

def Cholesky(matrix):
    n = matrix.shape[0]
    L = np.zeros((n,n))
    for i in range(n):
        L[i,i] = (matrix[i,i] - np.dot(L[i,:i],L[i,:i].T))**0.5
        for j in range(i+1,n):
            L[j,i] = (matrix[j,i] - np.dot(L[j,:i],L[i,:i].T))/L[i,i]
    return L

下面是一种更简洁的代码:

def Cholesky(A):
    len=A.shape[0]
    R=np.zeros((len,len))
    for i in range(len):
        for j in range(i, len):
            temp=A[i,j]
            for k in range(i):
                temp-=R[k,j]*R[k,i]
            if(i==j):
                R[i,j]=np.sqrt(temp)
            else:
                R[i,j]=temp/R[i,i]
    return R