krypy.utils
- Krylov Subspace Utilities¶
The utils module provides helper functions for common tasks in the process of solving linear algebraic systems.
Collection of standard functions.
This method provides functions like inner products, norms, …
-
exception
krypy.utils.
ArgumentError
¶ Bases:
Exception
Raised when an argument is invalid.
Analogue to
ValueError
which is not used here in order to be able to distinguish between built-in errors andkrypy
errors.
-
exception
krypy.utils.
AssumptionError
¶ Bases:
Exception
Raised when an assumption is not satisfied.
Differs from
ArgumentError
in that all passed arguments are valid but computations reveal that assumptions are not satisfied and the result cannot be computed.
-
exception
krypy.utils.
ConvergenceError
(msg, solver)¶ Bases:
Exception
Raised when a method did not converge.
The
ConvergenceError
holds a message describing the error and the attributesolver
through which the last approximation and other relevant information can be retrieved.
-
exception
krypy.utils.
LinearOperatorError
¶ Bases:
Exception
Raised when a
LinearOperator
cannot be applied.
-
exception
krypy.utils.
InnerProductError
¶ Bases:
Exception
Raised when the inner product is indefinite.
-
exception
krypy.utils.
RuntimeError
¶ Bases:
Exception
Raised for errors that do not fit in any other exception.
-
class
krypy.utils.
Arnoldi
(A, v, maxiter=None, ortho='mgs', M=None, Mv=None, Mv_norm=None, ip_B=None)¶ Bases:
object
Arnoldi algorithm.
Computes V and H such that \(AV_n=V_{n+1}\underline{H}_n\). If the Krylov subspace becomes A-invariant then V and H are truncated such that \(AV_n = V_n H_n\).
Parameters: - A – a linear operator that can be used with scipy’s
aslinearoperator with
shape==(N,N)
. - v – the initial vector with
shape==(N,1)
. - maxiter – (optional) maximal number of iterations. Default: N.
- ortho –
(optional) orthogonalization algorithm: may be one of
'mgs'
: modified Gram-Schmidt (default).'dmgs'
: double Modified Gram-Schmidt.'lanczos'
: Lanczos short recurrence.'house'
: Householder.
- M – (optional) a self-adjoint and positive definite
preconditioner. If
M
is provided, then also a second basis \(P_n\) is constructed such that \(V_n=MP_n\). This is of importance in preconditioned methods.M
has to beNone
ifortho=='house'
(seeB
). - ip_B –
(optional) defines the inner product to use. See
inner()
.ip_B
has to beNone
ifortho=='house'
. It’s unclear to me (andrenarchy), how a variant of the Householder QR algorithm can be used with a non-Euclidean inner product. Compare http://math.stackexchange.com/questions/433644/is-householder-orthogonalization-qr-practicable-for-non-euclidean-inner-products
-
advance
()¶ Carry out one iteration of Arnoldi.
-
get
()¶
-
get_last
()¶
- A – a linear operator that can be used with scipy’s
aslinearoperator with
-
class
krypy.utils.
BoundCG
(evals, exclude_zeros=False)¶ Bases:
object
CG residual norm bound.
Computes the \(\kappa\)-bound for the CG error \(A\)-norm when the eigenvalues of the operator are given, see [LieS13].
Parameters: - evals – an array of eigenvalues \(\lambda_1,\ldots,\lambda_N\in\mathbb{R}\). The eigenvalues will be sorted internally such that \(0=\lambda_1=\ldots=\lambda_{t-1}<\lambda_t\leq\ldots\lambda_N\) for \(t\in\mathbb{N}\).
- steps – (optional) the number of steps \(k\) to compute the bound
for. If steps is
None
(default), then \(k=N\) is used.
Returns: array \([\eta_0,\ldots,\eta_k]\) with
\[\eta_n = 2 \left( \frac{\sqrt{\kappa_{\text{eff}}} - 1} {\sqrt{\kappa_{\text{eff}}} + 1} \right)^n \quad\text{for}\quad n\in\{0,\ldots,k\}\]where \(\kappa_{\text{eff}}=\frac{\lambda_N}{\lambda_t}\).
Initialize with array/list of eigenvalues or Intervals object.
-
eval_step
(step)¶ Evaluate bound for given step.
-
get_step
(tol)¶ Return step at which bound falls below tolerance.
-
class
krypy.utils.
BoundMinres
(evals)¶ Bases:
object
MINRES residual norm bound.
Computes a bound for the MINRES residual norm when the eigenvalues of the operator are given, see [Gre97].
Parameters: - evals – an array of eigenvalues \(\lambda_1,\ldots,\lambda_N\in\mathbb{R}\). The eigenvalues will be sorted internally such that \(\lambda_1\leq\ldots\lambda_s<0=\lambda_{s+1}=\ldots=\lambda_{s+t-1}<\lambda_t\leq\ldots\lambda_N\) for \(s,t\in\mathbb{N}\) and \(s<t\).
- steps – (optional) the number of steps \(k\) to compute the bound
for. If steps is
None
(default), then \(k=N\) is used.
Returns: array \([\eta_0,\ldots,\eta_k]\) with
\[\eta_n = 2 \left( \frac{ \sqrt{|\lambda_1\lambda_N|} - \sqrt{|\lambda_s\lambda_t|}} { \sqrt{|\lambda_1\lambda_N|} + \sqrt{|\lambda_s\lambda_t|}} \right)^{\left[\frac{n}{2}\right]} \quad\text{for}\quad n\in\{0,\ldots,k\}\]if \(s>0\). If \(s=0\), i.e., if the eigenvalues are non-negative, then the result of
bound_cg()
is returned.Initialize with array/list of eigenvalues or Intervals object.
-
static
__new__
(cls, evals)¶ Use BoundCG if all eigenvalues are non-negative.
-
eval_step
(step)¶ Evaluate bound for given step.
-
get_step
(tol)¶ Return step at which bound falls below tolerance.
-
exception
krypy.utils.
ConvergenceError
(msg, solver) Bases:
Exception
Raised when a method did not converge.
The
ConvergenceError
holds a message describing the error and the attributesolver
through which the last approximation and other relevant information can be retrieved.
-
class
krypy.utils.
Givens
(x)¶ Bases:
object
Compute Givens rotation for provided vector x.
Computes Givens rotation \(G=\begin{bmatrix}c&s\\-\overline{s}&c\end{bmatrix}\) such that \(Gx=\begin{bmatrix}r\\0\end{bmatrix}\).
-
apply
(x)¶ Apply Givens rotation to vector x.
-
-
class
krypy.utils.
House
(x)¶ Bases:
object
Compute Householder transformation for given vector.
Initialize Householder transformation \(H\) such that \(Hx = \alpha \|x\|_2 e_1\) with \(|\alpha|=1\)
The algorithm is a combination of Algorithm 5.1.1 on page 236 and the treatment of the complex case in Section 5.1.13 on page 243 in Golub, Van Loan. Matrix computations. Fourth Edition. 2013.
-
apply
(x)¶ Apply Householder transformation to vector x.
Applies the Householder transformation efficiently to the given vector.
-
matrix
()¶ Build matrix representation of Householder transformation.
Builds the matrix representation \(H = I - \beta vv^*\).
Use with care! This routine may be helpful for testing purposes but should not be used in production codes for high dimensions since the resulting matrix is dense.
-
-
class
krypy.utils.
IdentityLinearOperator
(shape)¶ Bases:
krypy.utils.LinearOperator
-
class
krypy.utils.
LinearOperator
(shape, dtype, dot=None, dot_adj=None)¶ Bases:
object
Linear operator.
Is partly based on the LinearOperator from scipy (BSD License).
-
__add__
(X)¶
-
__mul__
(X)¶
-
__neg__
()¶
-
__pow__
(X)¶
-
__repr__
()¶ Return repr(self).
-
__rmul__
(X)¶
-
__sub__
(X)¶
-
adj
¶
-
dot
(X)¶
-
dot_adj
(X)¶
-
-
class
krypy.utils.
MatrixLinearOperator
(A)¶ Bases:
krypy.utils.LinearOperator
-
__repr__
()¶ Return repr(self).
-
-
class
krypy.utils.
NormalizedRootsPolynomial
(roots)¶ Bases:
object
A polynomial with specified roots and p(0)=1.
Represents the polynomial
\[p(\lambda) = \prod_{i=1}^n \left(1-\frac{\lambda}{\theta_i}\right).\]Parameters: roots – array with roots \(\theta_1,\dots,\theta_n\) of the polynomial and roots.shape==(n,)
.-
__call__
(points)¶ Evaluate polyonmial at given points.
Parameters: points – a point \(x\) or array of points \(x_1,\dots,x_m\) with points.shape==(m,)
.Returns: \(p(x)\) or array of shape (m,)
with \(p(x_1),\dots,p(x_m)\).
-
minmax_candidates
()¶ Get points where derivative is zero.
Useful for computing the extrema of the polynomial over an interval if the polynomial has real roots. In this case, the maximum is attained for one of the interval endpoints or a point from the result of this function that is contained in the interval.
-
-
class
krypy.utils.
Projection
(X, Y=None, ip_B=None, orthogonalize=True, iterations=2)¶ Bases:
object
Generic projection.
This class can represent any projection (orthogonal and oblique) on a N-dimensional Hilbert space. A projection is a linear operator \(P\) with \(P^2=P\). A projection is uniquely defined by its range \(\mathcal{V}:=\operatorname{range}(P)\) and its kernel \(\mathcal{W}:=\operatorname{ker}(P)\); this projection is called \(P_{\mathcal{V},\mathcal{W}}\).
Let X and Y be two full rank arrays with
shape==(N,k)
and let \(\mathcal{X}\oplus\mathcal{Y}^\perp=\mathbb{C}^N\) where \(\mathcal{X}:=\operatorname{colspan}(X)\) and \(\mathcal{Y}:=\operatorname{colspan}(Y)\). Then this class constructs the projection \(P_{\mathcal{X},\mathcal{Y}^\perp}\). The requirement \(\mathcal{X}\oplus\mathcal{Y}^\perp=\mathbb{C}^N\) is equivalent to\langle X,Y\rangle
being nonsingular.Parameters: - X – array with
shape==(N,k)
and \(\operatorname{rank}(X)=k\). - Y – (optional)
None
or array withshape==(N,k)
and \(\operatorname{rank}(X)=k\). If Y isNone
then Y is set to X which means that the resulting projection is orthogonal. - ip_B – (optional) inner product, see
inner()
.None
, anumpy.array
or aLinearOperator
is preferred due to the applicability of the proposed algorithms in [Ste11], see below. - orthogonalize – (optional) True orthogonalizes the bases provided in X and Y with respect to the inner product defined by ip_B. Defaults to True as the orthogonalization is suggested by the round-off error analysis in [Ste11].
- iterations – (optional) number of applications of the projection. It was suggested in [Ste11] to use 2 iterations (default) in order to achieve high accuracy (“twice is enough” as in the orthogonal case).
This projection class makes use of the round-off error analysis of oblique projections in the work of Stewart [Ste11] and implements the algorithms that are considered as the most stable ones (e.g., the XQRY representation in [Ste11]).
-
apply
(a, return_Ya=False)¶ Apply the projection to an array.
The computation is carried out without explicitly forming the matrix corresponding to the projection (which would be an array with
shape==(N,N)
).See also
_apply()
.
-
apply_adj
(a)¶
-
apply_complement
(a, return_Ya=False)¶ Apply the complementary projection to an array.
Parameters: z – array with shape==(N,m)
.Returns: \(P_{\mathcal{Y}^\perp,\mathcal{X}}z = z - P_{\mathcal{X},\mathcal{Y}^\perp} z\).
-
apply_complement_adj
(a)¶
-
matrix
()¶ Builds matrix representation of projection.
Builds the matrix representation \(P = X \langle Y,X\rangle^{-1} \langle Y, I_N\rangle\).
Use with care! This routine may be helpful for testing purposes but should not be used in production codes for high dimensions since the resulting matrix is dense.
-
operator
()¶ Get a
LinearOperator
corresponding to apply().Returns: a LinearOperator that calls apply().
-
operator_complement
()¶ Get a
LinearOperator
corresponding to apply_complement().Returns: a LinearOperator that calls apply_complement().
- X – array with
-
class
krypy.utils.
Timer
¶ Bases:
list
Measure execution time of multiple code blocks with
with
.Example:
t = Timer() with t: print('time me!') print('don\'t time me!') with t: print('time me, too!') print(t)
Result:
time me! don't time me! time me, too! [6.389617919921875e-05, 6.008148193359375e-05]
-
__enter__
()¶
-
__exit__
(a, b, c)¶
-
-
krypy.utils.
angles
(F, G, ip_B=None, compute_vectors=False)¶ Principal angles between two subspaces.
This algorithm is based on algorithm 6.2 in Knyazev, Argentati. Principal angles between subspaces in an A-based scalar product: algorithms and perturbation estimates. 2002. This algorithm can also handle small angles (in contrast to the naive cosine-based svd algorithm).
Parameters: - F – array with
shape==(N,k)
. - G – array with
shape==(N,l)
. - ip_B – (optional) angles are computed with respect to this
inner product. See
inner()
. - compute_vectors – (optional) if set to
False
then only the angles are returned (default). If set toTrue
then also the principal vectors are returned.
Returns: theta
ifcompute_vectors==False
theta, U, V
ifcompute_vectors==True
where
theta
is the array withshape==(max(k,l),)
containing the principal angles \(0\leq\theta_1\leq\ldots\leq\theta_{\max\{k,l\}}\leq \frac{\pi}{2}\).U
are the principal vectors from F with \(\langle U,U \rangle=I_k\).V
are the principal vectors from G with \(\langle V,V \rangle=I_l\).
The principal angles and vectors fulfill the relation \(\langle U,V \rangle = \begin{bmatrix} \cos(\Theta) & 0_{m,l-m} \\ 0_{k-m,m} & 0_{k-m,l-m} \end{bmatrix}\) where \(m=\min\{k,l\}\) and \(\cos(\Theta)=\operatorname{diag}(\cos(\theta_1),\ldots,\cos(\theta_m))\). Furthermore, \(\theta_{m+1}=\ldots=\theta_{\max\{k,l\}}=\frac{\pi}{2}\).
- F – array with
-
krypy.utils.
arnoldi
(*args, **kwargs)¶
-
krypy.utils.
arnoldi_res
(A, V, H, ip_B=None)¶ Measure Arnoldi residual.
Parameters: - A – a linear operator that can be used with scipy’s aslinearoperator
with
shape==(N,N)
. - V – Arnoldi basis matrix with
shape==(N,n)
. - H – Hessenberg matrix: either \(\underline{H}_{n-1}\) with
shape==(n,n-1)
or \(H_n\) withshape==(n,n)
(if the Arnoldi basis spans an A-invariant subspace). - ip_B – (optional) the inner product to use, see
inner()
.
Returns: either \(\|AV_{n-1} - V_n \underline{H}_{n-1}\|\) or \(\|A V_n - V_n H_n\|\) (in the invariant case).
- A – a linear operator that can be used with scipy’s aslinearoperator
with
-
krypy.utils.
arnoldi_projected
(H, P, k, ortho='mgs')¶ Compute (perturbed) Arnoldi relation for projected operator.
Assume that you have computed an Arnoldi relation
\[A V_n = V_{n+1} \underline{H}_n\]where \(V_{n+1}\in\mathbb{C}^{N,n+1}\) has orthogonal columns (with respect to an inner product \(\langle\cdot,\cdot\rangle\)) and \(\underline{H}_n\in\mathbb{C}^{n+1,n}\) is an extended upper Hessenberg matrix.
For \(k<n\) you choose full rank matrices \(X\in\mathbb{C}^{n-1,k}\) and \(Y\in\mathbb{C}^{n,k}\) and define \(\tilde{X}:=A V_{n_1}X = V_n \underline{H}_{n-1} X\) and \(\tilde{Y}:=V_n Y\) such that \(\langle \tilde{Y}, \tilde{X} \rangle = Y^*\underline{H}_{n-1} X\) is invertible. Then the projections \(P\) and \(\tilde{P}\) characterized by
- \(\tilde{P}x = x - \tilde{X} \langle \tilde{Y},\tilde{X} \rangle^{-1} \langle\tilde{Y},x\rangle\)
- \(P = I - \underline{H}_{n-1}X (Y^*\underline{H}_{n-1}X)^{-1}Y^*\)
are well defined and \(\tilde{P}V_{n+1} = [V_n P, v_{n+1}]\) holds.
This method computes for \(i<n-k\) the Arnoldi relation
\[(\tilde{P}A + E_i) W_i = W_{i+1} \underline{G}_i\]where \(W_{i+1}=V_n U_{i+1}\) has orthogonal columns with respect to \(\langle\cdot,\cdot\rangle\), \(\underline{G}_i\) is an extended upper Hessenberg matrix and \(E_i x = v_{n+1} F_i \langle W_i,x\rangle\) with \(F_i=[f_1,\ldots,f_i]\in\mathbb{C}^{1,i}\).
The perturbed Arnoldi relation can also be generated with the operator \(P_{V_n} \tilde{P} A\):
\[P_{V_n} \tilde{P} A W_i = W_{i+1} \underline{G}_i.\]In a sense the perturbed Arnoldi relation is the best prediction for the behavior of the Krylov subspace \(K_i(\tilde{P}A,\tilde{P}v_1)\) that can be generated only with the data from \(K_{n+1}(A,v_1)\) and without carrying out further matrix-vector multiplications with A.
Parameters: - H – the extended upper Hessenberg matrix
\(\underline{H}_n\) with
shape==(n+1,n)
. - P – the projection
\(P:\mathbb{C}^n\longrightarrow\mathbb{C}^n\) (has to be
compatible with
get_linearoperator()
). - k – the dimension of the null space of P.
Returns: U, G, F where
- U is the coefficient matrix \(U_{i+1}\) with
shape==(n,i+1)
, - G is the extended upper Hessenberg matrix \(\underline{G}_i\)
with
shape==(i+1,i)
, - F is the error matrix \(F_i\) with
shape==(1,i)
.
-
krypy.utils.
bound_perturbed_gmres
(pseudo, p, epsilon, deltas)¶ Compute GMRES perturbation bound based on pseudospectrum
Computes the GMRES bound from [SifEM13].
-
krypy.utils.
gap
(lamda, sigma, mode='individual')¶ Compute spectral gap.
Useful for eigenvalue/eigenvector bounds. Computes the gap \(\delta\geq 0\) between two sets of real numbers
lamda
andsigma
. The gap can be computed in several ways and may not exist, see themode
parameter.Parameters: - lamda – a non-empty set
\(\Lambda=\{\lambda_1,\ldots,\lambda_n\}\) given as a single real
number or a list or
numpy.array
with real numbers. - sigma – a non-empty set \(\Sigma=\{\sigma_1,\ldots,\sigma_m\}\).
See
lamda
. - mode –
(optional). Defines how the gap should be computed. May be one of
'individual'
(default): \(\delta=\min_{\substack{i\in\{1,\ldots,n\}\\j\in\{1,\ldots,m\}}} |\lambda_i - \sigma_j|\). With this mode, the gap is always be defined.'interval'
: determine the maximal \(\delta\) such that \(\Sigma\subset\mathbb{R}\setminus[\min_{\lambda\in\Lambda}\lambda-\delta,\max_{\lambda\in\Lambda}\lambda+\delta]\). If the gap does not exists,None
is returned.
Returns: \(\delta\) or
None
.- lamda – a non-empty set
\(\Lambda=\{\lambda_1,\ldots,\lambda_n\}\) given as a single real
number or a list or
-
krypy.utils.
get_linearoperator
(shape, A, timer=None)¶ Enhances aslinearoperator if A is None.
-
krypy.utils.
hegedus
(A, b, x0, M=None, Ml=None, ip_B=None)¶ Rescale initial guess appropriately (Hegedüs trick).
The Hegedüs trick rescales the initial guess to \(\gamma_{\min} x_0\) such that
\[\|r_0\|_{M^{-1}} = \| M M_l (b - A \gamma_{\min} x_0) \|_{M^{-1}} = \min_{\gamma\in\mathbb{C}} \| M M_l (b - A \gamma x_0) \|_{M^{-1}} \leq \| M M_l b \|_{M^{-1}}.\]This is achieved by \(\gamma_{\min} = \frac{\langle z, M M_l b \rangle_{M^{-1}}}{\|z\|_{M^{-1}}^2}\) for \(z=M M_l A x_0\) because then \(r_0=P_{z^\perp}b\). (Note that the right hand side of formula (5.8.16) in [LieS13] has to be complex conjugated.)
The parameters are the parameters you want to pass to
gmres()
,minres()
orcg()
.Returns: the adapted initial guess with the above property.
-
krypy.utils.
inner
(X, Y, ip_B=None)¶ Euclidean and non-Euclidean inner product.
numpy.vdot only works for vectors and numpy.dot does not use the conjugate transpose.
Parameters: - X – numpy array with
shape==(N,m)
- Y – numpy array with
shape==(N,n)
- ip_B –
(optional) May be one of the following
None
: Euclidean inner product.- a self-adjoint and positive definite operator \(B\) (as
numpy.array
orLinearOperator
). Then \(X^*B Y\) is returned. - a callable which takes 2 arguments X and Y and returns \(\langle X,Y\rangle\).
Caution: a callable should only be used if necessary. The choice potentially has an impact on the round-off behavior, e.g. of projections.
Returns: numpy array \(\langle X,Y\rangle\) with shape==(m,n)
.- X – numpy array with
-
krypy.utils.
ip_euclid
(X, Y)¶ Euclidean inner product.
numpy.vdot only works for vectors and numpy.dot does not use the conjugate transpose.
Parameters: - X – numpy array with
shape==(N,m)
- Y – numpy array with
shape==(N,n)
Returns: numpy array \(X^* Y\) with
shape==(m,n)
.- X – numpy array with
-
krypy.utils.
norm
(x, y=None, ip_B=None)¶ Compute norm (Euclidean and non-Euclidean).
Parameters: - x – a 2-dimensional
numpy.array
. - y – a 2-dimensional
numpy.array
. - ip_B – see
inner()
.
Compute \(\sqrt{\langle x,y\rangle}\) where the inner product is defined via
ip_B
.- x – a 2-dimensional
-
krypy.utils.
norm_MMlr
(M, Ml, A, Mr, b, x0, yk, inner_product=<function ip_euclid>)¶
-
krypy.utils.
norm_squared
(x, Mx=None, inner_product=<function ip_euclid>)¶ Compute the norm^2 w.r.t. to a given scalar product.
-
krypy.utils.
orthonormality
(V, ip_B=None)¶ Measure orthonormality of given basis.
Parameters: - V – a matrix \(V=[v_1,\ldots,v_n]\) with
shape==(N,n)
. - ip_B – (optional) the inner product to use, see
inner()
.
Returns: \(\| I_n - \langle V,V \rangle \|_2\).
- V – a matrix \(V=[v_1,\ldots,v_n]\) with
-
krypy.utils.
qr
(X, ip_B=None, reorthos=1)¶ QR factorization with customizable inner product.
Parameters: - X – array with
shape==(N,k)
- ip_B – (optional) inner product, see
inner()
. - reorthos – (optional) numer of reorthogonalizations. Defaults to 1 (i.e. 2 runs of modified Gram-Schmidt) which should be enough in most cases (TODO: add reference).
Returns: Q, R where \(X=QR\) with \(\langle Q,Q \rangle=I_k\) and R upper triangular.
- X – array with
-
krypy.utils.
ritz
(H, V=None, hermitian=False, type='ritz')¶ Compute several kinds of Ritz pairs from an Arnoldi/Lanczos relation.
This function computes Ritz, harmonic Ritz or improved harmonic Ritz values and vectors with respect to the Krylov subspace \(K_n(A,v)\) from the extended Hessenberg matrix \(\underline{H}_n\) generated with n iterations the Arnoldi algorithm applied to A and v.
Parameters: - H – Hessenberg matrix from Arnoldi/Lanczos algorithm.
- V –
(optional) Arnoldi/Lanczos vectors, \(V\in\mathbb{C}^{N,n+1}\). If provided, the Ritz vectors are also returned. The Arnoldi vectors have to form an orthonormal basis with respect to an inner product.
Caution: if you are using the Lanzcos or Gram-Schmidt Arnoldi algorithm without reorthogonalization, then the orthonormality of the basis is usually lost. For accurate results it is advisable to use the Householder Arnoldi (
ortho='house'
) or modified Gram-Schmidt with reorthogonalization (ortho='dmgs'
). - hermitian – (optional) if set to
True
the matrix \(H_n\) must be Hermitian. A Hermitian matrix \(H_n\) allows for faster and often more accurate computation of Ritz pairs. - type –
(optional) type of Ritz pairs, may be one of
'ritz'
,'harmonic'
or'harmonic_like'
. Two choices of Ritz pairs fit in the following description:Given two n-dimensional subspaces \(X,Y\subseteq \mathbb{C}^N\), find a basis \(z_1,\ldots,z_n\) of \(X\) and \(\theta_1,\ldots,\theta_n\in\mathbb{C}\) such that \(A z_i - \theta_i z_i \perp Y\) for all \(i\in\{1,\ldots,n\}\).
In this setting the choices are
'ritz'
: regular Ritz pairs, i.e. \(X=Y=K_n(A,v)\).'harmonic'
: harmonic Ritz pairs, i.e.- \(X=K_n(A,v)\) and \(Y=AK_n(A,v)\).
'harmonic_improved'
: the returned vectorsU
(andV
, if requested) are the same as withtype='harmonic'
. Thetheta
array contains the improved Ritz values \(\theta_i = u_i^* H_n u_i\), cf. section 2 in Morgan, Zeng. Harmonic Projection Methods for Large Non-symmetric Eigenvalue Problems. 1998. It can be shown that the residual norm of improved Ritz pairs is always less than or equal to the residual norm of the harmonic Ritz pairs. However, the improved Ritz pairs do not fit into the framework above since the orthogonality condition is lost.
Returns: - If V is not
None
thentheta, U, resnorm, Z
is returned. - If V is
None
thentheta, U, resnorm
is returned.
Where
theta
are the Ritz values \([\theta_1,\ldots,\theta_n]\).U
are the coefficients of the Ritz vectors in the Arnoldi basis, i.e. \(z_i=Vu_i\) where \(u_i\) is the i-th column of U.resnorm
is a residual norm vector.Z
are the actual Ritz vectors, i.e.Z=dot(V,U)
.
-
krypy.utils.
shape_vec
(x)¶ Take a (n,) ndarray and return it as (n,1) ndarray.
-
krypy.utils.
shape_vecs
(*args)¶ Reshape all ndarrays with
shape==(n,)
toshape==(n,1)
.Recognizes ndarrays and ignores all others.