python中五对角矩阵的Cholesky优化

Question

所以我在研究五对角矩阵A，大小n：

（这里还有五对角矩阵的一般信息： https://en.wikipedia.org/wiki/Pentadiagonal_matrix） 我使用Cholesky分解来得到矩阵A的L，其中L*L.T=A，（L.T是L的转置）根据算法。因此numpy的标准算法是：

def mycholesky(A):                                                                                                 
    """Performs a Cholesky decomposition of A, which must                                                          
    be a symmetric and positive definite matrix. The function                                                      
    returns the lower variant triangular matrix, L."""                                                             
    n = len(A)                                                                                                     

    # Create zero matrix for L                                                                                     
    #L = [[0.0] * n for i in range(n)]                                                                             

    n = len(A)                                                                                                     
    L = np.zeros((n, n), dtype=float)                                                                              





# Perform the Cholesky decomposition                                                                           

for i in range(n):                                                                                             
    for k in range(i+1):                                                                                       
        tmp_sum = sum(L[i][j] * L[k][j] for j in range(k))                                                 

        if (i == k):  # Diagonal elements                                                                      
            # LaTeX: l_{kk} = \sqrt{ a_{kk} - \sum^{k-1}_{j=1} l^2_{kj}}                                       
            L[i][k] = sqrt(A[i][i] - tmp_sum)                                                                  
        else:                                                                                                  
            # LaTeX: l_{ik} = \frac{1}{l_{kk}} \left( a_{ik} - \sum^{k-1}_{j=1} l_{ij} l_{kj} \right)          
            L[i][k] = (1.0 / L[k][k] * (A[i][k] - tmp_sum))                                                    
return L         

您可以在这里看到页面，其中也有数学公式。我对算法进行了一些修改，以供Python3和numpy使用： https://www.quantstart.com/articles/Cholesky-Decomposition-in-Python-and-NumPy

我想优化算法，因为我正在处理的A矩阵是一个稀疏的矩阵，我想测试非常大的n（即n=10000）。经典的cholesky没有优化，因为有很多零不需要访问。到目前为止，我尝试的是改变代码行的范围

tmp_sum = sum(L[i][j] * L[k][j] for j in range(k))

收件人：

tmp_sum = sum(L[i][j] * L[k][j] for j in range(k-2,k))

为了避免每次也计算零的sum。能否进一步优化？因为仍然可以访问零，并且不需要进行计算。 或者另一个解决方案是，取原始A的带矩阵，并在其上应用cholesky？你知道吗

Answer 1

对于带状矩阵的Cholesky（在您的例子中是五对角矩阵）的优化，肯定还有更多的工作可以做。你知道吗

特别是，我要向您指出Python基础结构中的一个现有解决方案：^{}。你知道吗

使用它可以让您利用已经实现的优化，同时查看此子例程背后的实际代码，您可以找到其背后的思想。你知道吗

稀疏矩阵分解和保持稀疏性的研究非常活跃，有很多解决方案，因此根据您的需要（编码/研究新算法与获得特定矩阵分解），您有不同的行动方案。你知道吗