scipy 慢速稀疏矩阵求解器 [英] scipy slow sparse matrix solver

问题描述

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

I can't seem to get a benefit out of scipy's CG and sparse matrix algorithms.

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

When I try to solve this banded matrix equation

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 次迭代才能达到结果.因此，由于有关条带的先验知识，我预计它至少会一样快，但实际上要快得多.是不是和条件数这么差有关?

it takes about 20 times as long as the LU solver. From what I understood, CG isn't that much more costly than LU, even if it has to use all N iterations to reach the result. So I expected it to be at least as fast, but actually a lot faster, due to the prior knowledge about the bandedness. Does it have to do with the condition number to be so bad?

在密集矩阵的情况下

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的条件数似乎没有那么大，但是我对这种事情没有太多经验.

I don't get convergence at all. In this case, the condition number of A doesn't seem that large, but then I don't have a lot of experience with this kind of thing.

推荐答案

You create A using

You create A using

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

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

That call to .toarray() converts the sparse matrix to a regular numpy array. So you are passing a regular array to a sparse solver, which means the sparse solver can't take advantage of any sparsity structure. If you pass the original sparse matrix to the solver, it is much faster.

为了解决带状系统，一个快速的替代方法是 scipy.linalg.solve_banded.(还有 scipy.linalg.solveh_banded 用于 Hermitian 系统.)

For solving a banded system, a fast alternative is scipy.linalg.solve_banded. (There is also scipy.linalg.solveh_banded for Hermitian systems.)

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

Here's your example, but with the sparse matrix passed to the sparse solver. Also included is a solution computed using scipy.linalg.solve_banded, which turns out to be much faster than the other two methods.

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

这篇关于scipy 慢速稀疏矩阵求解器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！