正在学习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