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 parameterip_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 toFalse
. - self_adjoint – (bool, optional) Is \(M_l A M_r\) self-adjoint
in the inner product defined by
ip_B
?self_adjoint=True
also setsnormal=True
. Defaults toFalse
. - positive_definite – (bool, optional) Is \(M_l A M_r\)
positive (semi-)definite with respect to the inner product defined by
ip_B
? Defaults toFalse
. - exact_solution – (optional) If an exact solution \(x\) is
known, it can be provided as a
numpy.array
withexact_solution.shape == (N,1)
. Then error norms can be computed (for debugging or research purposes). Defaults toNone
.
-
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.
- z – approximate solution with
- A – a linear operator on \(\mathbb{C}^N\) (has to be
compatible with
-
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 onM
,Ml
,A
,Mr
andip_B
above.-
static
operations
(nsteps)¶ Returns the number of operations needed for nsteps of CG
- if
-
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 ofArnoldi
).
All parameters of
_KrylovSolver
are valid in this solver. Note the restrictions onM
,Ml
,A
,Mr
andip_B
above.-
static
operations
(nsteps)¶ Returns the number of operations needed for nsteps of MINRES
- if
-
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
-
static
-
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 operatorsParameters: - 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 toTrue
, the residual is computed explicitly in each iteration and thus requires an additional application ofM
,Ml
,A
andMr
in each iteration. - store_arnoldi – (optional)
if set to
True
then the computed Arnoldi basis and the Hessenberg matrix are set as attributesV
andH
on the returned object. IfM
is notNone
, then alsoP
is set whereV=M*P
. Defaults toFalse
. If the method is based on the Lanczos method (e.g.,Cg
orMinres
), thenH
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 parametertol
.errnorms
: the error norms of all iterations ifexact_solution
was provided.V
,H
andP
ifstore_arnoldi==True
, seestore_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.
- linear_system – a