1

所以我正在处理这个五角矩阵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（即 for 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？

4

1 回答 1

3

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

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

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

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

于 2019-12-15T21:34:05.883 回答