0

背景：我正在为我的数值线性代数课程做一个项目。对于这个项目，我决定尝试使用半精度算术进行不完全 Cholesky 分解，并将结果用作迭代方法的预处理器。我首先尝试实现这个 Matlab 2019b（它具有半精度数据类型），但它不支持半精度稀疏矩阵，所以我不得不使用全矩阵。但是半精度的算术在 Matlab 中要慢得多，我发现像 500 x 500 矩阵这样的因子需要 20 分钟（我想达到 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 格式的 scipy 稀疏矩阵。（注释掉的代码只是完全写出的算法，而不是试图利用稀疏性）。我发现这大约需要 6 分钟才能运行，这比我使用半精度浮点数时的 MATLAB 代码快得多，但当我使用单/双精度浮点数时（耗时不到一秒），它仍然比 matlab 代码慢很多，即使 MATLAB 使用的是完整矩阵。

总是有可能我只是在某个地方的代码中犯了一个错误，而我实际上并没有获得正确的运行时间，所以我会再看一遍。但我想知道是否有人更习惯于 scipy / numpy 看到我选择实现上述代码的方式有什么不妥之处。

我确实有另一种理论来解释为什么 python 代码可能这么慢。我在我学校的高性能计算机上运行它，可能是 matlab 设置为自动利用并行性，但 python 不是。这看起来像是一个合理的假设吗？如果是这样，您对我如何正确并行化我的算法有什么建议吗？

4

1 回答 1

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
于 2020-04-30T03:55:12.923 回答