# 2.5.3. Linear System Solvers¶

• sparse matrix/eigenvalue problem solvers live in `scipy.sparse.linalg`

• the submodules:
• `dsolve`: direct factorization methods for solving linear systems
• `isolve`: iterative methods for solving linear systems
• `eigen`: sparse eigenvalue problem solvers
• all solvers are accessible from:

```>>> import scipy.sparse.linalg as spla
>>> spla.__all__
['LinearOperator', 'Tester', 'arpack', 'aslinearoperator', 'bicg',
'bicgstab', 'cg', 'cgs', 'csc_matrix', 'csr_matrix', 'dsolve',
'eigen', 'eigen_symmetric', 'factorized', 'gmres', 'interface',
'isolve', 'iterative', 'lgmres', 'linsolve', 'lobpcg', 'lsqr',
'minres', 'np', 'qmr', 'speigs', 'spilu', 'splu', 'spsolve', 'svd',
'test', 'umfpack', 'use_solver', 'utils', 'warnings']
```

## 2.5.3.1. Sparse Direct Solvers¶

• default solver: SuperLU 4.0
• included in SciPy
• real and complex systems
• both single and double precision
• optional: umfpack
• real and complex systems
• double precision only
• recommended for performance
• wrappers now live in `scikits.umfpack`
• check-out the new `scikits.suitesparse` by Nathaniel Smith

### Examples¶

• import the whole module, and see its docstring:

```>>> from scipy.sparse.linalg import dsolve
>>> help(dsolve)
```
• both superlu and umfpack can be used (if the latter is installed) as follows:

• prepare a linear system:

```>>> import numpy as np
>>> from scipy import sparse
>>> mtx = sparse.spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]], [0, 1], 5, 5)
>>> mtx.todense()
matrix([[ 1,  5,  0,  0,  0],
[ 0,  2,  8,  0,  0],
[ 0,  0,  3,  9,  0],
[ 0,  0,  0,  4, 10],
[ 0,  0,  0,  0,  5]])
>>> rhs = np.array([1, 2, 3, 4, 5], dtype=np.float32)
```
• solve as single precision real:

```>>> mtx1 = mtx.astype(np.float32)
>>> x = dsolve.spsolve(mtx1, rhs, use_umfpack=False)
>>> print(x)
[106.   -21.     5.5   -1.5    1. ]
>>> print("Error: %s" % (mtx1 * x - rhs))
Error:  [0.  0.  0.  0.  0.]
```
• solve as double precision real:

```>>> mtx2 = mtx.astype(np.float64)
>>> x = dsolve.spsolve(mtx2, rhs, use_umfpack=True)
>>> print(x)
[106.   -21.     5.5   -1.5    1. ]
>>> print("Error: %s" % (mtx2 * x - rhs))
Error:  [0.  0.  0.  0.  0.]
```
• solve as single precision complex:

```>>> mtx1 = mtx.astype(np.complex64)
>>> x = dsolve.spsolve(mtx1, rhs, use_umfpack=False)
>>> print(x)
[106. +0.j  -21. +0.j    5.5+0.j   -1.5+0.j    1. +0.j]
>>> print("Error: %s" % (mtx1 * x - rhs))
Error:  [0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
```
• solve as double precision complex:

```>>> mtx2 = mtx.astype(np.complex128)
>>> x = dsolve.spsolve(mtx2, rhs, use_umfpack=True)
>>> print(x)
[106. +0.j  -21. +0.j    5.5+0.j   -1.5+0.j    1. +0.j]
>>> print("Error: %s" % (mtx2 * x - rhs))
Error:  [0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
```
```"""
Solve a linear system
=======================

Construct a 1000x1000 lil_matrix and add some values to it, convert it
to CSR format and solve A x = b for x:and solve a linear system with a
direct solver.
"""
import numpy as np
import scipy.sparse as sps
from matplotlib import pyplot as plt
from scipy.sparse.linalg.dsolve import linsolve

rand = np.random.rand

mtx = sps.lil_matrix((1000, 1000), dtype=np.float64)
mtx[0, :100] = rand(100)
mtx[1, 100:200] = mtx[0, :100]
mtx.setdiag(rand(1000))

plt.clf()
plt.spy(mtx, marker='.', markersize=2)
plt.show()

mtx = mtx.tocsr()
rhs = rand(1000)

x = linsolve.spsolve(mtx, rhs)

print('rezidual: %r' % np.linalg.norm(mtx * x - rhs))
```

## 2.5.3.2. Iterative Solvers¶

• the `isolve` module contains the following solvers:
• `bicg` (BIConjugate Gradient)
• `bicgstab` (BIConjugate Gradient STABilized)
• `cg` (Conjugate Gradient) - symmetric positive definite matrices only
• `cgs` (Conjugate Gradient Squared)
• `gmres` (Generalized Minimal RESidual)
• `minres` (MINimum RESidual)
• `qmr` (Quasi-Minimal Residual)

### Common Parameters¶

• mandatory:

A : {sparse matrix, dense matrix, LinearOperator}

The N-by-N matrix of the linear system.

b : {array, matrix}

Right hand side of the linear system. Has shape (N,) or (N,1).

• optional:

x0 : {array, matrix}

Starting guess for the solution.

tol : float

Relative tolerance to achieve before terminating.

maxiter : integer

Maximum number of iterations. Iteration will stop after maxiter steps even if the specified tolerance has not been achieved.

M : {sparse matrix, dense matrix, LinearOperator}

Preconditioner for A. The preconditioner should approximate the inverse of A. Effective preconditioning dramatically improves the rate of convergence, which implies that fewer iterations are needed to reach a given error tolerance.

callback : function

User-supplied function to call after each iteration. It is called as callback(xk), where xk is the current solution vector.

### LinearOperator Class¶

```from scipy.sparse.linalg.interface import LinearOperator
```
• common interface for performing matrix vector products
• useful abstraction that enables using dense and sparse matrices within the solvers, as well as matrix-free solutions
• has shape and matvec() (+ some optional parameters)
• example:
```>>> import numpy as np
>>> from scipy.sparse.linalg import LinearOperator
>>> def mv(v):
...     return np.array([2*v, 3*v])
...
>>> A = LinearOperator((2, 2), matvec=mv)
>>> A
<2x2 _CustomLinearOperator with dtype=float64>
>>> A.matvec(np.ones(2))
array([2.,  3.])
>>> A * np.ones(2)
array([2.,  3.])
```

### A Few Notes on Preconditioning¶

• problem specific
• often hard to develop
• if not sure, try ILU
• available in `dsolve` as `spilu()`

## 2.5.3.3. Eigenvalue Problem Solvers¶

### The `eigen` module¶

• `arpack` * a collection of Fortran77 subroutines designed to solve large scale eigenvalue problems

• `lobpcg` (Locally Optimal Block Preconditioned Conjugate Gradient Method) * works very well in combination with PyAMG * example by Nathan Bell:

```"""
Compute eigenvectors and eigenvalues using a preconditioned eigensolver
========================================================================

In this example Smoothed Aggregation (SA) is used to precondition
the LOBPCG eigensolver on a two-dimensional Poisson problem with
Dirichlet boundary conditions.
"""

import scipy
from scipy.sparse.linalg import lobpcg
import matplotlib.pyplot as plt

from pyamg import smoothed_aggregation_solver
from pyamg.gallery import poisson

N = 100
K = 9
A = poisson((N,N), format='csr')

# create the AMG hierarchy
ml = smoothed_aggregation_solver(A)

# initial approximation to the K eigenvectors
X = scipy.rand(A.shape, K)

# preconditioner based on ml
M = ml.aspreconditioner()

# compute eigenvalues and eigenvectors with LOBPCG
W,V = lobpcg(A, X, M=M, tol=1e-8, largest=False)

#plot the eigenvectors
plt.figure(figsize=(9,9))

for i in range(K):
plt.subplot(3, 3, i+1)
plt.title('Eigenvector %d' % i)
plt.pcolor(V[:,i].reshape(N,N))
plt.axis('equal')
plt.axis('off')
plt.show()
```
• example by Nils Wagner:

• output:

```\$ python examples/lobpcg_sakurai.py
Results by LOBPCG for n=2500

[ 0.06250083  0.06250028  0.06250007]

Exact eigenvalues

[ 0.06250005  0.0625002   0.06250044]

Elapsed time 7.01
``` 