krypy.linsys - Linear Algebraic Systems Solver

The linsys module provides functions for the solution of linear algebraic systems.

class krypy.linsys.LinearSystem(A, b, M=None, Minv=None, Ml=None, Mr=None, ip_B=None, normal=None, self_adjoint=False, positive_definite=False, exact_solution=None)

Bases: object

Representation of a (preconditioned) linear system.

Represents a linear system

\[Ax=b\]

or a preconditioned linear system

\[M M_l A M_r y = M M_l b \quad\text{with}\quad x=M_r y.\]
Parameters:
  • A – a linear operator on \(\mathbb{C}^N\) (has to be compatible with get_linearoperator()).
  • b – the right hand side in \(\mathbb{C}^N\), i.e., b.shape == (N, 1).
  • M – (optional) a self-adjoint and positive definite preconditioner, linear operator on \(\mathbb{C}^N\) with respect to the inner product defined by ip_B. This preconditioner changes the inner product to \(\langle x,y\rangle_M = \langle Mx,y\rangle\) where \(\langle \cdot,\cdot\rangle\) is the inner product defined by the parameter ip_B. Defaults to the identity.
  • Minv – (optional) the inverse of the preconditioner provided by M. This operator is needed, e.g., for orthonormalizing vectors for the computation of Ritz vectors in deflated methods.
  • Ml – (optional) left preconditioner, linear operator on \(\mathbb{C}^N\). Defaults to the identity.
  • Mr – (optional) right preconditioner, linear operator on \(\mathbb{C}^N\). Defaults to the identity.
  • ip_B – (optional) defines the inner product, see inner().
  • normal – (bool, optional) Is \(M_l A M_r\) normal in the inner product defined by ip_B? Defaults to False.
  • self_adjoint – (bool, optional) Is \(M_l A M_r\) self-adjoint in the inner product defined by ip_B? self_adjoint=True also sets normal=True. Defaults to False.
  • positive_definite – (bool, optional) Is \(M_l A M_r\) positive (semi-)definite with respect to the inner product defined by ip_B? Defaults to False.
  • exact_solution – (optional) If an exact solution \(x\) is known, it can be provided as a numpy.array with exact_solution.shape == (N,1). Then error norms can be computed (for debugging or research purposes). Defaults to None.
MMlb_norm = None

Norm of the right hand side.

\[\|M M_l b\|_{M^{-1}}\]
N = None

Dimension \(N\) of the space \(\mathbb{C}^N\) where the linear system is defined.

get_ip_Minv_B()

Returns the inner product that is implicitly used with the positive definite preconditioner M.

get_residual(z, compute_norm=False)

Compute residual.

For a given \(z\in\mathbb{C}^N\), the residual

\[r = M M_l ( b - A z )\]

is computed. If compute_norm == True, then also the absolute residual norm

\[\| M M_l (b-Az)\|_{M^{-1}}\]

is computed.

Parameters:
  • z – approximate solution with z.shape == (N, 1).
  • compute_norm – (bool, optional) pass True if also the norm of the residual should be computed.
class krypy.linsys.Cg(linear_system, **kwargs)

Bases: krypy.linsys._KrylovSolver

Preconditioned CG method.

The preconditioned conjugate gradient method can be used to solve a system of linear algebraic equations where the linear operator is self-adjoint and positive definite. Let the following linear algebraic system be given:

\[M M_l A M_r y = M M_l b,\]

where \(x=M_r y\) and \(M_l A M_r\) is self-adjoint and positive definite with respect to the inner product \(\langle \cdot,\cdot \rangle\) defined by ip_B. The preconditioned CG method then computes (in exact arithmetics!) iterates \(x_k \in x_0 + M_r K_k\) with \(K_k:= K_k(M M_l A M_r, r_0)\) such that

\[\|x - x_k\|_A = \min_{z \in x_0 + M_r K_k} \|x - z\|_A.\]

The Lanczos alorithm is used with the operator \(M M_l A M_r\) and the inner product defined by \(\langle x,y \rangle_{M^{-1}} = \langle M^{-1}x,y \rangle\). The initial vector for Lanczos is \(r_0 = M M_l (b - Ax_0)\) - note that \(M_r\) is not used for the initial vector.

Memory consumption is:

  • if store_arnoldi==False: 3 vectors or 6 vectors if \(M\) is used.
  • if store_arnoldi==True: about maxiter+1 vectors for the Lanczos basis. If \(M\) is used the memory consumption is 2*(maxiter+1).

Caution: CG’s convergence may be delayed significantly due to round-off errors, cf. chapter 5.9 in [LieS13].

All parameters of _KrylovSolver are valid in this solver. Note the restrictions on M, Ml, A, Mr and ip_B above.

static operations(nsteps)

Returns the number of operations needed for nsteps of CG

class krypy.linsys.Minres(linear_system, ortho='lanczos', **kwargs)

Bases: krypy.linsys._KrylovSolver

Preconditioned MINRES method.

The preconditioned minimal residual method can be used to solve a system of linear algebraic equations where the linear operator is self-adjoint. Let the following linear algebraic system be given:

\[M M_l A M_r y = M M_l b,\]

where \(x=M_r y\) and \(M_l A M_r\) is self-adjoint with respect to the inner product \(\langle \cdot,\cdot \rangle\) defined by inner_product. The preconditioned MINRES method then computes (in exact arithmetics!) iterates \(x_k \in x_0 + M_r K_k\) with \(K_k:= K_k(M M_l A M_r, r_0)\) such that

\[\|M M_l(b - A x_k)\|_{M^{-1}} = \min_{z \in x_0 + M_r K_k} \|M M_l (b - A z)\|_{M^{-1}}.\]

The Lanczos alorithm is used with the operator \(M M_l A M_r\) and the inner product defined by \(\langle x,y \rangle_{M^{-1}} = \langle M^{-1}x,y \rangle\). The initial vector for Lanczos is \(r_0 = M M_l (b - Ax_0)\) - note that \(M_r\) is not used for the initial vector.

Memory consumption is:

  • if store_arnoldi==False: 3 vectors or 6 vectors if \(M\) is used.
  • if store_arnoldi==True: about maxiter+1 vectors for the Lanczos basis. If \(M\) is used the memory consumption is 2*(maxiter+1).

Caution: MINRES’ convergence may be delayed significantly or even stagnate due to round-off errors, cf. chapter 5.9 in [LieS13].

In addition to the attributes described in _KrylovSolver, the following attributes are available in an instance of this solver:

  • lanczos: the Lanczos relation (an instance of Arnoldi).

All parameters of _KrylovSolver are valid in this solver. Note the restrictions on M, Ml, A, Mr and ip_B above.

static operations(nsteps)

Returns the number of operations needed for nsteps of MINRES

class krypy.linsys.Gmres(linear_system, ortho='mgs', **kwargs)

Bases: krypy.linsys._KrylovSolver

Preconditioned GMRES method.

The preconditioned generalized minimal residual method can be used to solve a system of linear algebraic equations. Let the following linear algebraic system be given:

\[M M_l A M_r y = M M_l b,\]

where \(x=M_r y\). The preconditioned GMRES method then computes (in exact arithmetics!) iterates \(x_k \in x_0 + M_r K_k\) with \(K_k:= K_k(M M_l A M_r, r_0)\) such that

\[\|M M_l(b - A x_k)\|_{M^{-1}} = \min_{z \in x_0 + M_r K_k} \|M M_l (b - A z)\|_{M^{-1}}.\]

The Arnoldi alorithm is used with the operator \(M M_l A M_r\) and the inner product defined by \(\langle x,y \rangle_{M^{-1}} = \langle M^{-1}x,y \rangle\). The initial vector for Arnoldi is \(r_0 = M M_l (b - Ax_0)\) - note that \(M_r\) is not used for the initial vector.

Memory consumption is about maxiter+1 vectors for the Arnoldi basis. If \(M\) is used the memory consumption is 2*(maxiter+1).

If the operator \(M_l A M_r\) is self-adjoint then consider using the MINRES method Minres.

All parameters of _KrylovSolver are valid in this solver.

static operations(nsteps)

Returns the number of operations needed for nsteps of GMRES

class krypy.linsys._KrylovSolver(linear_system, x0=None, tol=1e-05, maxiter=None, explicit_residual=False, store_arnoldi=False, dtype=None)

Prototype of a Krylov subspace method for linear systems.

Init standard attributes and perform checks.

All Krylov subspace solvers in this module are applied to a LinearSystem. The specific methods may impose further restrictions on the operators

Parameters:
  • linear_system – a LinearSystem.
  • x0 – (optional) the initial guess to use. Defaults to zero vector. Unless you have a good reason to use a nonzero initial guess you should use the zero vector, cf. chapter 5.8.3 in Liesen, Strakos. Krylov subspace methods. 2013. See also hegedus().
  • tol –

    (optional) the tolerance for the stopping criterion with respect to the relative residual norm:

    \[\frac{ \| M M_l (b-A (x_0+M_r y_k))\|_{M^{-1}} } { \|M M_l b\|_{M^{-1}}} \leq \text{tol}\]
  • maxiter – (optional) maximum number of iterations. Defaults to N.
  • explicit_residual – (optional) if set to False (default), the updated residual norm from the used method is used in each iteration. If set to True, the residual is computed explicitly in each iteration and thus requires an additional application of M, Ml, A and Mr in each iteration.
  • store_arnoldi – (optional) if set to True then the computed Arnoldi basis and the Hessenberg matrix are set as attributes V and H on the returned object. If M is not None, then also P is set where V=M*P. Defaults to False. If the method is based on the Lanczos method (e.g., Cg or Minres), then H is real, symmetric and tridiagonal.
  • dtype – (optional) an optional dtype that is used to determine the dtype for the Arnoldi/Lanczos basis and matrix.

Upon convergence, the instance contains the following attributes:

  • xk: the approximate solution \(x_k\).
  • resnorms: relative residual norms of all iterations, see parameter tol.
  • errnorms: the error norms of all iterations if exact_solution was provided.
  • V, H and P if store_arnoldi==True, see store_arnoldi

If the solver does not converge, a ConvergenceError is thrown which can be used to examine the misconvergence.

errnorms = None

Error norms.

iter = None

Iteration number.

static operations(nsteps)

Returns the number of operations needed for nsteps of the solver.

Parameters:nsteps – number of steps.
Returns:a dictionary with the same keys as the timings parameter. Each value is the number of operations of the corresponding type for nsteps iterations of the method.
resnorms = None

Relative residual norms as described for parameter tol.

xk = None

Approximate solution.