scipy慢稀疏矩阵求解器

Question

我似乎无法从 scipy 的 CG 和稀疏矩阵算法中获益。

当我尝试求解这个带状矩阵方程时

import time
import scipy.sparse.linalg as ssla
import scipy.sparse as ss
import numpy as np

size = 3000
x = np.ones(size)
A = ss.diags([1, -2, 1], [-1, 0, 1], shape=[size, size]).toarray()
print 'cond. nr.:{}'.format(np.linalg.cond(A))
b = np.dot(A, x)

start_time = time.clock()
sol = np.linalg.solve(A, b)
elapsed_time = time.clock() - start_time
error = np.sum(np.abs(sol - x))
print 'LU time, error:', elapsed_time, error

start_time = time.clock()
sol, info = ssla.bicg(A, b, tol=1e-12)
elapsed_time = time.clock() - start_time
error = np.sum(np.abs(sol - x))
print 'CG time, ret code, error:', elapsed_time, info, error

它所花费的时间大约是 LU 求解器的 20 倍。据我了解，CG 并不比 LU 贵多少，即使它必须使用所有 N 次迭代才能达到结果。所以我预计它至少会一样快，但实际上要快得多，因为先验知识是关于带状的。是不是跟条件数这么差有关系？

在密集矩阵的情况下

import time
import scipy.sparse.linalg as ssla
import numpy as np
import matplotlib.pyplot as plt

size = 1000
x = np.ones(size)
np.random.seed(1)
A = np.random.random_integers(-size, size,
                              size=(size, size))
print 'cond. nr.:{}'.format(np.linalg.cond(A))
b = np.dot(A, x)

start_time = time.clock()
sol = np.linalg.solve(A, b)
elapsed_time = time.clock() - start_time
error = np.sum(np.abs(sol - x))
print 'LU time, error:', elapsed_time, error

start_time = time.clock()
sol, info = ssla.bicg(A, b, tol=1e-12)
elapsed_time = time.clock() - start_time
error = np.sum(np.abs(sol - x))
print 'CG time, ret code, error:', elapsed_time, info, error

我根本没有收敛。在这种情况下，A的条件数似乎并没有那么大，但是我对这种事情没有太多经验。

Answer 1

您使用

创建 A

A = ss.diags([1, -2, 1], [-1, 0, 1], shape=[size, size]).toarray()

对.toarray() 的调用将稀疏矩阵转换为常规的numpy 数组。因此，您将常规数组传递给稀疏求解器，这意味着稀疏求解器无法利用任何稀疏结构。 如果将原始稀疏矩阵传递给求解器，它会快得多。

为了解决带状系统，一个快速的替代方案是scipy.linalg.solve_banded。 （对于 Hermitian 系统，还有 scipy.linalg.solveh_banded。）

这是您的示例，但将稀疏矩阵传递给稀疏求解器。还包括一个使用scipy.linalg.solve_banded 计算的解决方案，结果证明它比其他两种方法快很多。

import time
import scipy.sparse.linalg as ssla
import scipy.sparse as ss
from scipy.linalg import solve_banded
import numpy as np

size = 3000
x = np.ones(size)
S = ss.diags([1, -2, 1], [-1, 0, 1], shape=[size, size])
A = S.toarray()
print 'cond. nr.:{}'.format(np.linalg.cond(A))
b = np.dot(A, x)

start_time = time.clock()
sol = np.linalg.solve(A, b)
elapsed_time = time.clock() - start_time
error = np.sum(np.abs(sol - x))
print 'LU time, error          :  %9.6f    %g' % (elapsed_time, error)

start_time = time.clock()
sol, info = ssla.bicg(S, b, tol=1e-12)
elapsed_time = time.clock() - start_time
error = np.sum(np.abs(sol - x))
print 'CG time, ret code, error:  %9.6f %2d %g' % (elapsed_time, info, error)

B = np.empty((3, size))
B[0, :] = 1
B[1, :] = -2
B[2, :] = 1
start_time = time.clock()
sol = solve_banded((1, 1), B, b)
elapsed_time = time.clock() - start_time
error = np.sum(np.abs(sol - x))
print 'solve_banded time, error:  %9.6f    %g' % (elapsed_time, error)

输出：

cond. nr.:3649994.05818
LU time, error          :   0.858295    4.05262e-09
CG time, ret code, error:   0.552952  0 6.7263e-11
solve_banded time, error:   0.000750    4.05262e-09