不完全的Cholesky因式分解很慢

Question

背景：我正在为我的数值线性代数课程做一个项目。对于这个项目，我决定尝试使用半精度算术进行不完全的cholesky因式分解，并将结果用作迭代方法的前提。我首先尝试实现此Matlab 2019b（具有半精度数据类型），但它不支持半精度sparse矩阵，因此我不得不使用完整矩阵。但是半精度的算术在Matlab中要慢得多，我发现花了大约20分钟才能分解出500 x 500的矩阵（而我想提高到1000 x 1000）。但是，以单精度/双精度，500 x 500矩阵花费的时间不到一秒钟。

我以为如果我可以真正利用矩阵的稀疏性，我将有更好的运气缩放到更高的矩阵。我记得numpy / scipy有一个float 16数据类型，所以我决定尝试在python中实现它。所以我写了这个

from scipy.io import loadmat
def icholesky(a):
    n = a.shape[0]
    for k in tqdm(list(range(n))): 
        a[k,k] = np.sqrt(a[k,k])
        #for i in range(k+1,n):
        #    if (a[i,k] !=0):
        #        a[i,k] = a[i,k]/a[k,k]
        i,_= a[:,k].nonzero()
        if len(i) > 0:
            a[i,k] = a[i,k]/a[k,k]
        for j in range(k+1,n):
            #for i in range(j,n):
            #    if (a[i,j]!=0):
            #        a[i,j] = a[i,j]-a[i,k]*a[j,k]  
            i,_ = a[j:,j].nonzero()
            if len(i) > 0: 
                a[i,j]  = a[i,j] - a[i,k]*a[j,k]     
    return a

bus = loadmat('494_bus.mat') #From University of Florida's Sparse Matrix Collection
A = bus['Problem'][0,0][1]
H = A.copy()
icholesky(H)

其中'a'是具有CSC格式的稀疏矩阵。 （注释掉的代码只是完全写出的算法，没有尝试利用稀疏性）。我发现运行此过程大约需要6分钟，这比使用半精度浮点数的MATLAB代码要快得多，但仍比使用单精度/双精度浮点数（不到一秒）的matlab代码要慢得多。 ，即使MATLAB使用的是完整矩阵。

总是有可能我只是在某个地方的代码中犯了一个错误，而我实际上并没有获得正确的运行时间，因此我将再次进行研究。但是我想知道是否还有更习惯于scipy / numpy的人看到我选择实现上述代码的方式有什么毛病。

我还有另一种理论，为什么python代码可能这么慢。我正在学校的高性能计算机上运行此程序，可能是将matlab设置为自动利用并行性的情况，而python没有。这看起来像是合理的假设吗？如果是这样，您是否对如何正确并行化算法有任何建议？

Answer 1

我想出了一种减轻问题的方法。我不需要j遍历从k + 1到n的所有数字。我只需要检查a [j，k]为非零的j。这是我计算i_时刚刚计算的结果。

def icholesky(a):
    n = a.shape[0]
    for k in range(n): 
        a[k,k] = np.sqrt(a[k,k])
        i_,_ = a[k+1:,k].nonzero() 
        if len(i_) > 0:
            i_= i_ + (k+1)
            a[i_,k] = a[i_,k]/a[k,k]
        for j in i_: 
            i2_,_ = a[j:n,j].nonzero()
            if len(i2_) > 0:
                i2_ = i2_ + j
                a[i2_,j]  = a[i2_,j] - a[i2_,k]*a[j,k]   


    return a