python中五对角矩阵的Cholesky优化

Question

所以我正在处理这个五对角矩阵A，大小n：

（这里还有五角矩阵的一般信息： https://en.wikipedia.org/wiki/Pentadiagonal_matrix) 我正在使用 Cholesky 分解 来获得矩阵 L 的矩阵 A ，其中 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 的 band 矩阵并在其上应用 cholesky ？

Answer 1

肯定可以做更多的事情来优化带状矩阵的 Cholesky（在您的情况下是五对角矩阵）。

特别是，我会向您指出 Python 基础架构中的现有解决方案：scipy.linalg.cholesky_banded。

使用它可以让您利用已经实现的优化，同时查看该子例程背后的实际代码，您可以找到背后的想法。

稀疏矩阵分解和保持稀疏性的研究非常活跃，并且有很多解决方案可用，因此根据您的需求（编码/研究新算法与分解特定矩阵），您有不同的行动方案。