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))

• examples/direct_solve.py

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[0], 3*v[1]])
...
>>> 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[0], 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()    
  

  • examples/pyamg_with_lobpcg.py

• example by Nils Wagner:

  • examples/lobpcg_sakurai.py

• 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