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 builtin 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 Ainvariant 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 GramSchmidt (default).'dmgs'
: double Modified GramSchmidt.'lanczos'
: Lanczos short recurrence.'house'
: Householder.
 M – (optional) a selfadjoint 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 nonEuclidean inner product. Compare http://math.stackexchange.com/questions/433644/ishouseholderorthogonalizationqrpracticablefornoneuclideaninnerproducts

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_{t1}<\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+t1}<\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 nonnegative, then the result of
bound_cg()
is returned.Initialize with array/list of eigenvalues or Intervals object.

static
__new__
(evals)¶ Use BoundCG if all eigenvalues are nonnegative.

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__
()¶

__rmul__
(X)¶

__sub__
(X)¶

adj
¶

dot
(X)¶

dot_adj
(X)¶


class
krypy.utils.
MatrixLinearOperator
(A)¶ Bases:
krypy.utils.LinearOperator

__repr__
()¶


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 Ndimensional 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 roundoff 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 roundoff 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.389617919921875e05, 6.008148193359375e05]

__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 Abased scalar product: algorithms and perturbation estimates. 2002. This algorithm can also handle small angles (in contrast to the naive cosinebased 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,lm} \\ 0_{km,m} & 0_{km,lm} \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}_{n1}\) with
shape==(n,n1)
or \(H_n\) withshape==(n,n)
(if the Arnoldi basis spans an Ainvariant subspace).  ip_B – (optional) the inner product to use, see
inner()
.
Returns: either \(\AV_{n1}  V_n \underline{H}_{n1}\\) 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}^{n1,k}\) and \(Y\in\mathbb{C}^{n,k}\) and define \(\tilde{X}:=A V_{n_1}X = V_n \underline{H}_{n1} X\) and \(\tilde{Y}:=V_n Y\) such that \(\langle \tilde{Y}, \tilde{X} \rangle = Y^*\underline{H}_{n1} 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}_{n1}X (Y^*\underline{H}_{n1}X)^{1}Y^*\)
are well defined and \(\tilde{P}V_{n+1} = [V_n P, v_{n+1}]\) holds.
This method computes for \(i<nk\) 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 matrixvector 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 nonempty set
\(\Lambda=\{\lambda_1,\ldots,\lambda_n\}\) given as a single real
number or a list or
numpy.array
with real numbers.  sigma – a nonempty 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 nonempty 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 nonEuclidean 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 selfadjoint 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 roundoff 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 nonEuclidean).
Parameters:  x – a 2dimensional
numpy.array
.  y – a 2dimensional
numpy.array
.  ip_B – see
inner()
.
Compute \(\sqrt{\langle x,y\rangle}\) where the inner product is defined via
ip_B
. x – a 2dimensional

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 GramSchmidt) 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 GramSchmidt 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 GramSchmidt 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 ndimensional 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 Nonsymmetric 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 ith 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.