---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.18.0-dev
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

# Linear System Solvers

- sparse matrix/eigenvalue problem solvers live in {mod}`scipy.sparse.linalg`

- the submodules:
  - {mod}`dsolve`: direct factorization methods for solving linear systems
  - {mod}`isolve`: iterative methods for solving linear systems
  - {mod}`eigen`: sparse eigenvalue problem solvers

All solvers are accessible from:

```{code-cell}
import scipy as sp
sp.sparse.linalg.__all__
```

## Sparse Direct Solvers

- default solver: SuperLU
  - 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 {mod}`scikits.umfpack`
  - check-out the new {mod}`scikits.suitesparse` by Nathaniel Smith

### Examples

Import the whole module, and see its docstring:

```{code-cell}
help(sp.sparse.linalg.spsolve)
```

Both superlu and umfpack can be used (if the latter is installed) as follows.

Prepare a linear system:

```{code-cell}
import numpy as np
mtx = sp.sparse.spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]], [0, 1], 5, 5, "csc")
mtx.toarray()
```

```{code-cell}
rhs = np.array([1, 2, 3, 4, 5], dtype=np.float32)
```

Solve as single precision real:

```{code-cell}
mtx1 = mtx.astype(np.float32)
x = sp.sparse.linalg.spsolve(mtx1, rhs, use_umfpack=False)
print(x)
```

```{code-cell}
print("Error: %s" % (mtx1 * x - rhs))
```

Solve as double precision real:

```{code-cell}
mtx2 = mtx.astype(np.float64)
x = sp.sparse.linalg.spsolve(mtx2, rhs, use_umfpack=True)
print(x)
```

```{code-cell}
print("Error: %s" % (mtx2 * x - rhs))
```

Solve as single precision complex:

```{code-cell}
mtx1 = mtx.astype(np.complex64)
x = sp.sparse.linalg.spsolve(mtx1, rhs, use_umfpack=False)
print(x)
```

```{code-cell}
print("Error: %s" % (mtx1 * x - rhs))
```

Solve as double precision complex:

```{code-cell}
mtx2 = mtx.astype(np.complex128)
x = sp.sparse.linalg.spsolve(mtx2, rhs, use_umfpack=True)
print(x)
```

```{code-cell}
print("Error: %s" % (mtx2 * x - rhs))
```

{download}`examples/direct_solve.py`

+++

## Iterative Solvers

- the {mod}`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` : The N-by-N matrix of the linear system.
  - `b`: Right hand side of the linear system. Has shape (N,) or (N,1).

- optional:
  - `x0`: Starting guess for the solution.
  - `tol` : Relative tolerance to achieve before terminating.
  - `maxiter` : Maximum number of iterations. Iteration will stop after maxiter
    steps even if the specified tolerance has not been achieved.
  - `M` : 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` : User-supplied function to call after each iteration. It is
    called as `callback(xk)`, where `xk` is the current solution vector.

+++

### LinearOperator Class

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

Here is an example:

```{code-cell}
import numpy as np
import scipy as sp

def mv(v):
    return np.array([2 * v[0], 3 * v[1]])
```

```{code-cell}
A = sp.sparse.linalg.LinearOperator((2, 2), matvec=mv)
A
```

```{code-cell}
A.matvec(np.ones(2))
```

```{code-cell}
A * np.ones(2)
```

### A Few Notes on Preconditioning

- problem specific
- often hard to develop
- if not sure, try ILU
  - available in {mod}`scipy.sparse.linalg` as {func}`spilu()`

## Eigenvalue Problem Solvers

### The {mod}`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](https://github.com/pyamg/pyamg)
  - example by Nathan Bell:

    {download}`examples/pyamg_with_lobpcg.py`

Another example by Nils Wagner:

{download}`examples/lobpcg_sakurai.py`

Output:

```bash
$ 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
```

![](figures/lobpcg_eigenvalues.png)