scipy slow sparse matrix solver

Question

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

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?

In case of a dense matrix

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

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.

Convergence also depends on the initial guess you choose. If you happen to pick the correct answer as the starting point you'll converge in one iteration; make a poor selection and you may require many more iterations. It's not just about the matrix and its conditioning. — duffymo
A (1000, 1000) matrix is trivially small. Try something bigger and sparser. — ali_m
Adding x0=x to the bicg call in my second example does'nt fix the problem. When reducing the size to 10, I get a valid result with bicg without starting guess. Adding x0=x in gives me a wrong solution. What am I doing wrong? — zeus300

Warren Weckesser Warren Weckesser · Accepted Answer · 2016-07-19T02:47:56

You create A using

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

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.

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

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)

Output:

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 slow sparse matrix solver

1 Answers