Restarted and block eigensolvers

Let $A\in {\mathbb{C}}^{n\times n}$, let $v\mathrm{\ne }0$, and choose integers $1\le r<m\le n$, a residual tolerance $\tau \ge 0$, a restart budget $N_{\max }\ge 1$, and a restart mode thick or implicit. All operations below use the exact arithmetic model and projected-solver input specified next. The number $k$ of locked pairs is the number of linearly independent locked vectors; it starts at zero. If a restart is needed, a target rule may be supplied for thick mode and is required for an implicit-mode restart after one or more pairs are locked. When invoked, it must return $r-k$ current Ritz vectors independent modulo the span of all locked vectors, or report unavailable; no rule, an unavailable selection, or a wrong-size or dependent selection gives status rank-loss. For an implicit-mode restart with no locked pairs, a shift rule supplies exactly $p=m-r$ complex shifts for the current cycle or reports unavailable; an unavailable list or a list of any other length gives status shift-count-unavailable. The shifts may be chosen from unwanted Ritz values, but the procedure allows any supplied list and asserts no convergence for an arbitrary schedule.

Projected-solver input and exact operations. Supply a projected eigensolver $S$. For every requested exact square matrix $H\in {\mathbb{C}}^{d\times d}$ with $1\le d\le m$, it must return in finite work either unavailable or success with an ordered list of exactly $d$ pairs $({\theta }_j,y_j)$, $0\le j<d$, with ${\theta }_j\in \mathbb{C}$, $0\mathrm{\ne }{y}_j\in {\mathbb{C}}^d$, $H{y}_j={\theta }_j{y}_j$, and $\det (z{I}_d-H)={\prod }_{j=0}^{d-1}(z-{\theta }_j)\mathrm{.}$ Thus eigenvalues are represented with their algebraic multiplicities; projected eigenvectors need not be independent when $H$ is defective. The procedure does not assert existence of an always-successful exact eigensolver in any fixed arithmetic or encoding. A provider may report unavailable; absent provider data give the same result without invoking a solver. Every original successful exact projected-eigensolve instance is retained by supplying its provider.

Outside $S$, use exact complex field operations, conjugation, nonnegative real square roots, and exact equality/order comparisons on the real quantities compared. Complex operations use their real and imaginary component formulas. The finite-return target and shift rules have the output contracts already stated. Exact rank tests can be performed by Gaussian elimination, scanning columns and choosing the first remaining row with nonzero pivot; normalization divides only by a checked positive Euclidean norm. This is an exact-operation model, not a bit-cost or machine floating-point claim.

At every projected solve, invoke $S$ on the current $H_d$. An absent provider or unavailable return gives status projected-eigensolve-unavailable. A provider status other than success or unavailable, wrong list/vector shape, zero vector, nonzero exact residual $H_d{y}_j-{\theta }_j{y}_j$, or failed characteristic-polynomial identity gives status projected-eigensolve-invalid. Check that identity by exact coefficient arithmetic: expand $\det (z{I}_d-H_d)$ by the finite permutation formula of For $n\ge 1$, the determinant over a commutative ring by the Leibniz formula, and $\mathrm{\mid }\det A\mathrm{\mid }$ for a real matrix, in lexicographic permutation order, multiplying each term's $d$ affine-or-constant entries in increasing column order and using the inversion parity for its sign; form the right side by successive multiplication of its $d$ linear factors in returned list order. Polynomial addition and multiplication mean finite coefficient-array sums and convolutions, and compare all coefficients through degree $d$ exactly. This finite check extracts no roots. All these failures preserve the locked outputs and the most recently completed valid unlocked Ritz list (initially empty), with the failing cycle and projected dimension. Only after the entire successful list passes these checks normalize the vectors, lift them and replace the current Ritz list. The supplied finite-return promise applies to $S$ itself; no timeout or generic eigensolver implementation is silently provided.

Initialize the restart counter $c=0$ before the first cycle. The initial cycle is an $m$-step fully reorthogonalized Arnoldi run from $v_1=v\mathrm{/}\mathrm{\parallel }v{\mathrm{\parallel }}_2$: each step computes a matrix-vector product, performs two complete orthogonalization passes against the current basis, and accumulates the coefficients in the projected matrix. If a residual vanishes at step $j$, use the $j$-column basis and mark its span invariant. If the residual after step $m$ is nonzero, retain its normalized direction and coefficient as the Arnoldi residual data, but use only the first $m$ columns in the projected eigensolve. Every later cycle uses the restarted basis and projected matrix constructed in step 3 or 4. The output is the locked Ritz pairs, the most recently computed unlocked Ritz pairs, and one of converged, budget-exhausted, invariant-subspace-breakdown, rank-loss, shift-count-unavailable, or filtered-breakdown, projected-eigensolve-unavailable, or projected-eigensolve-invalid.

1. For the initial cycle, start with $v_1=v\mathrm{/}\mathrm{\parallel }v{\mathrm{\parallel }}_2$, no locked pairs, and run the $m$-step Arnoldi procedure described above. Let $d$ be the projected basis dimension, so $d=m$ unless exact breakdown occurs earlier; on a nonbreakdown cycle also retain the residual direction and coefficient from step $m$. For every cycle, invoke and validate the supplied projected eigensolver as specified above, returning its failure status before processing any pair if it fails. For each returned projected eigenvector $y$, normalize it so $\mathrm{\parallel }y{\mathrm{\parallel }}_2=1$ and lift it to $u=V_{d}y$. Process the Ritz pairs in the order returned by the projected eigensolve. For each pair $(\theta ,u)$, compute the original residual $r_u=Au-\theta u$ directly (on a standard Arnoldi cycle without breakdown it may be evaluated by An Arnoldi Ritz pair has residual norm controlled by the last Hessenberg subdiagonal entry). If $\mathrm{\parallel }{r}_u{\mathrm{\parallel }}_2\le \tau$ and adding $u$ increases the exact rank of the locked-vector matrix, append $(\theta ,u)$ to the locked outputs and increment $k$. A residual-accepted vector in the span of the locked vectors is not counted or appended again. Locked vectors are retained as a subspace at every later restart. If $k\ge r$, return converged.

2. If the current projected subspace is marked invariant and fewer than $r$ pairs are locked, return invariant-subspace-breakdown. If the restart counter satisfies $c=N_{\max }$, return budget-exhausted.

3. For a thick restart, or for an implicit-mode restart when $k>0$, invoke the target rule. Form a target block from every locked vector and the selected $r-k$ wanted Ritz vectors. For either block to be rank-revealing, use the following exact column-pivoted orthogonalization: start with an empty orthonormal block $Q$ and all input columns remaining. For each remaining column $y_j$, compute $c_j=y_j-Q(Q^{\ast }{y}_j)$. Choose a column with largest $\mathrm{\parallel }{c}_j{\mathrm{\parallel }}_2$, breaking ties by original column order. If the largest norm is zero, stop; otherwise append $c_j\mathrm{/}\mathrm{\parallel }{c}_j{\mathrm{\parallel }}_2$ to $Q$ and remove that pivot column. The columns appended to $Q$ form an ordered orthonormal basis of the input range. Apply this procedure first to the target block; if it fails to retain all locked and selected directions, return rank-loss. With its resulting orthonormal basis $U$, compute the residual block $R_U=(I-U{U}^{\ast })AU$. Apply the same orthogonalization to $R_U$ to obtain an ordered orthonormal basis $F$ of its range. If $U$ has $s$ columns, retain only the first $\min (\dim \mathrm{range}(R_U),m-s)$ columns of $F$, so $[U,F]$ has at most $m$ columns. Set the ordered basis to $[U,F]$ and a queue to its columns in order. While the basis has fewer than $m$ columns and the queue is nonempty, remove the first queued column $q$, compute $w=Aq$, and orthogonalize $w$ against every current basis column in two complete passes. If the residual is nonzero, normalize it, append it to the basis, and append the new column to the queue; if it is zero, continue with the remaining queue. Set $V_d$ to the resulting basis and compute the projected matrix $H_d=V_d^{\ast }A{V}_d$ exactly. If the queue empties before the basis has dimension $m$, or if $d=n$, mark this subspace invariant; first process its Ritz pairs in step 5, then return invariant-subspace-breakdown at step 2 if fewer than $r$ pairs have been locked. This projection is required because a thick restart can begin from several retained and residual directions rather than from a single Arnoldi starting vector.

4. For an implicit restart with no locked pairs, set $p=m-r$ and request a shift list ${\mu }_1,\dots ,{\mu }_p$ from the shift rule. If it reports unavailable or returns a list of the wrong length, return shift-count-unavailable. Starting from the no-breakdown $m$-step Arnoldi relation $A{V}_m=V_m{H}_m+{\beta }_m{v}_{m+1}{e}_m^{\ast }$, apply exact shifted QR steps to $H_m$ using the transformations of Real Givens rotations and complex Givens transformations. Set $H_0=H_m$ and $q(z)={\prod }_{j=1}^p(z-{\mu }_j)$. For $j=1,\dots ,p$, initialize $R=H_{j-1}-{\mu }_{j}I$ and $Q_j=I_m$. For $\mathrm{\ell }=1,\dots ,m-1$, generate a left factor from $(a,b)=(R_{\mathrm{\ell },\mathrm{\ell }},R_{\mathrm{\ell }+1,\mathrm{\ell }})$ and apply it to rows $\mathrm{\ell },\mathrm{\ell }+1$ by setting $R\leftarrow LR$, and update $Q_j\leftarrow Q_j{L}^{\ast }$. Use the identity if $b=0$; if $b\mathrm{\ne }0$, put $s_0=\sqrt{\mathrm{\mid }a{\mathrm{\mid }}^2+\mathrm{\mid }b{\mathrm{\mid }}^2}$ and choose the exact coefficients $(c,s)=(a\mathrm{/}{s}_0,b\mathrm{/}{s}_0)$ over $\mathbb{R}$, while over $\mathbb{C}$ choose $(c,s)=(0,\bar{b}\mathrm{/}\mathrm{\mid }b\mathrm{\mid })$ if $a=0$ and $(c,s)=(\mathrm{\mid }a\mathrm{\mid }\mathrm{/}{s}_0,(a\mathrm{/}\mathrm{\mid }a\mathrm{\mid })\bar{b}\mathrm{/}{s}_0)$ otherwise. Use $L=\left(\begin{array}{cc}c& s\\ -\bar{s}& c\end{array}\right)$. At the end set $R_j=R$ and $H_j=R_j{Q}_j+{\mu }_{j}I=Q_j^{\ast }{H}_{j-1}{Q}_j$. Let $H^{+}=H_p$ and put $Q=Q_1\cdots Q_p$. Test the filtered Krylov dimension exactly by forming $w=q(A)v_1$ and the columns $w,Aw,\dots ,A^{r-1}w$. If their rank is below $r$, return filtered-breakdown. Otherwise retain $V_r^{+}=V_{m}Q(:,1:r)$ and set $H_r^{+}=H^{+}(1:r,1:r)$. Compute $g=A(V_r^{+}{)}_{:r}-V_r^{+}{H}_r^{+}{e}_r$, the residual of the last retained column; the first $r-1$ retained columns have zero residual by the upper-Hessenberg relation and shift count $p=m-r$. If $g=0$, use $V_r^{+}$ as the current invariant basis with dimension $r$. If $g\mathrm{\ne }0$, set $v_{r+1}^{+}=g\mathrm{/}\mathrm{\parallel }g{\mathrm{\parallel }}_2$ and ${\beta }_r^{+}=\mathrm{\parallel }g{\mathrm{\parallel }}_2$, initialize the continued Arnoldi relation with $(V_r^{+},H_r^{+},v_{r+1}^{+},{\beta }_r^{+})$, and continue Arnoldi with two complete orthogonalization passes until dimension $m$ or exact breakdown. On a nonbreakdown continuation, also compute the residual direction and coefficient after its $m$th Arnoldi step for the next implicit restart. Mark the current subspace invariant when this continuation breaks down.

5. Increment the restart counter by setting $c\leftarrow c+1$, invoke and validate $S$ on the new projected matrix as in step 1, and return the stated projected-solve failure status if it fails. On success, recompute the lifted Ritz pairs and their original residuals, apply the acceptance rule in step 1, and return to step 2 unless a return status applies. No convergence is asserted for an arbitrary shift schedule. The residual acceptance test certifies only the computed residual size; it does not by itself give a forward eigenvalue bound for a nonnormal matrix.

Let $A\in {\mathbb{C}}^{n\times n}$, and suppose an exact $m$-step Arnoldi factorization without breakdown is $A{V}_m=V_m{H}_m+{\beta }_m{v}_{m+1}{e}_m^{\ast },V_m^{\ast }{V}_m=I_m,$ where $H_m$ is upper Hessenberg and ${\beta }_m>0$. Let $1\le r<m$, $p=m-r$, and let $H_{j-1}-{\mu }_{j}I=Q_j{R}_j$ be exact shifted QR factorizations using the ordered top-down adjacent-row Givens eliminations on pairs $(1,2),(2,3),\dots ,(m-1,m)$, with $R_j$ upper triangular. Put $Q=Q_1\cdots Q_p,H^{+}=Q^{\ast }{H}_{m}Q,q(z)={\prod }_{j=1}^p(z-{\mu }_j)\mathrm{.}$ Assume $q(A)v_1\mathrm{\ne }0$, the exact filtered Krylov space $K_r(A,q(A)v_1)$ has dimension $r$, and $\rho \mathrm{\ne }0$ in the factorization $q(H_m)=QR$ with $R$ upper triangular. Then $\mathrm{range}(V_{m}Q(:,1:r))=K_r(A,q(A)v_1)\mathrm{.}$ Writing $H_r^{+}=H^{+}(1:r,1:r)$, $h^{+}=H_{r+1,r}^{+}$, and $g=h^{+}{V}_{m}Q{e}_{r+1}+{\beta }_m(e_m^{\ast }Q{e}_r)v_{m+1},$ the truncated restart has the exact residual identity $A{V}_{m}Q(:,1:r)=V_{m}Q(:,1:r)H_r^{+}+g{e}_r^{\ast }\mathrm{.}$ If $g\mathrm{\ne }0$, then $v_{r+1}^{+}=g\mathrm{/}\mathrm{\parallel }g{\mathrm{\parallel }}_2$ is orthogonal to the retained columns and the new Arnoldi residual coefficient is ${\beta }_r^{+}=\mathrm{\parallel }g{\mathrm{\parallel }}_2$; if $g=0$, the retained subspace is exactly invariant.

Let $A=A^{\ast }\in {\mathbb{C}}^{n\times n}$, let $B=B^{\ast }>0$, and let $T=T^{\ast }>0$ be a fixed preconditioning operator. Here $B>0$ and $T>0$ mean $x^{\ast }Bx>0$ and $x^{\ast }Tx>0$ for every nonzero $x$. The $B$-inner product and norm are those of The B-inner product, B-norm, and dual residual norm. Choose a block width $1\le q\le n$, a residual tolerance $\tau \ge 0$, a finite iteration budget $N_{\max }\ge 0$, and an input matrix $X_0\in {\mathbb{C}}^{n\times q}$. The output is a block of Ritz pairs and one of converged, budget-exhausted, rank-loss, or invalid-initial-block.

1. If $X_0$ does not have full column rank or does not satisfy $X_0^{\ast }B{X}_0=I_q$, return invalid-initial-block. Diagonalize the Hermitian compression $X_0^{\ast }A{X}_0=Z_0{\mathrm{\Theta }}_0{Z}_0^{\ast }$ with $Z_0$ unitary and the $q$ Ritz values in increasing order, then replace $X_0$ by $X_0{Z}_0$. Thus the columns used in the first residual block are the actual B-orthonormal Ritz vectors paired with the diagonal entries of ${\mathrm{\Theta }}_0$. Set the previous search block $D_0$ empty and initialize the iteration counter $k=0$.

2. At iteration $k$, form the original residual block $R_k=A{X}_k-B{X}_k{\mathrm{\Theta }}_k$. If every column has Euclidean norm at most $\tau$, return converged. If $k=N_{\max }$, return budget-exhausted with $({\mathrm{\Theta }}_k,X_k)$.

3. Compute the preconditioned residual block $W_k=T{R}_k$ and set $Y_k=[X_k,W_k,D_k]$, omitting the empty block at $k=0$. Compute an exact rank-revealing $B$-QR with column pivoting as follows. Starting with no columns in $Q_k$, for each remaining column $y_j$ of $Y_k$ compute its residual $c_j=y_j-Q_k(Q_k^{\ast }B{y}_j)$ and choose a column with largest $B$-norm $\mathrm{\parallel }{c}_j{\mathrm{\parallel }}_B$ (break ties by original column order). If the largest residual norm is zero, stop; otherwise append the normalized residual to $Q_k$ and remove that pivot column. Let $s_k$ be the number of columns appended and let ${\mathrm{\Pi }}_k$ record the pivot order followed by the unused columns. If $s_k<q$, return rank-loss with the current Ritz pairs; no projected solve is attempted on that branch. Otherwise set $S_k=Q_k^{\ast }B{Y}_k{\mathrm{\Pi }}_k(:,1:s_k)\mathrm{.}$ The pivot construction gives $Y_k{\mathrm{\Pi }}_k(:,1:s_k)=Q_k{S}_k,Q_k^{\ast }B{Q}_k=I_{s_k},$ and $S_k$ is invertible because its selected columns are independent. Dependent trial columns are discarded and their count is recorded as a rank-compression diagnostic.

4. The projected matrix $Q_k^{\ast }A{Q}_k$ is Hermitian, hence normal, so its orthonormal eigenbasis and real eigenvalues are supplied by Complex spectral theorem: a normal endomorphism of a finite-dimensional complex inner product space has an orthonormal eigenbasis, and conversely. Solve $Q_k^{\ast }A{Q}_k{Z}_k=Z_k{\mathrm{\Theta }}_{k+1}$ and choose its $q$ smallest eigenvalues. Set $C_k={\mathrm{\Pi }}_k(:,1:s_k)S_k^{-1}{Z}_k$; then $Y_k{C}_k=Q_k{Z}_k$. Partition $C_k=\left[\begin{array}{c}{C}_X\\ C_W\\ C_D\end{array}\right]$, omitting the $C_D$ block when $D_k$ is empty. The projected metric is positive definite because for nonzero $z$, $(Q_{k}z{)}^{\ast }B(Q_{k}z)=\mathrm{\parallel }z{\mathrm{\parallel }}_2^2>0$. Set $X_{k+1}=X_k{C}_X+W_k{C}_W+D_k{C}_D,D_{k+1}=W_k{C}_W+D_k{C}_D,$ again omitting the empty-block term. Let ${\mathrm{\Theta }}_{k+1}$ be the selected diagonal Ritz-value matrix and increment $k$.

5. Repeat from step 2. No convergence, rate, or accuracy is asserted when the budget or rank-loss branch is returned. The exact SPD conditions on $B$ and $T$ define the permitted metric and preconditioner; no particular incomplete-factorization procedure is assumed.

Consider any execution of Finite-budget thick and implicitly restarted Arnoldi with residual acceptance and locking or Finite-budget block LOBPCG for a Hermitian-definite pencil with explicit rank-loss status on inputs satisfying that procedure's stated hypotheses. Assume every invoked target rule returns in finite work either a valid requested block or unavailable, every invoked shift rule returns in finite work either a shift list of the required length or unavailable, and every projected eigensolve, finite-dimensional arithmetic, orthogonalization, and rank computation returns its exact value in finite work. Then the execution returns one of the statuses listed by its procedure after finitely many operations. In particular, the Arnoldi procedure performs at most $N_{\max }$ restarts and at most $N_{\max }+1$ Arnoldi cycles, and the LOBPCG procedure performs at most $N_{\max }$ updates.

Let $A=A^{\ast }\in {\mathbb{C}}^{n\times n}$ and $B=B^{\ast }>0$. Write the generalized eigenvalues of $(A,B)$ in increasing order ${\lambda }_1\le \cdots \le {\lambda }_n$. For a subspace $S$ with a full-column-rank $B$-orthonormal basis matrix $Q_S$, let ${\theta }_1(S)\le \cdots \le {\theta }_{\dim S}(S)$ be the eigenvalues of $Q_S^{\ast }A{Q}_S$. These values are independent of the chosen $B$-orthonormal basis. Then, for $1\le j\le \dim S$, ${\lambda }_j\le {\theta }_j(S)\mathrm{.}$ If $S\subseteq T$ and $j\le \dim S$, then ${\theta }_j(T)\le {\theta }_j(S)\mathrm{.}$

Let $A=A^{\ast }\in {\mathbb{C}}^{n\times n}$, $B=B^{\ast }>0$, and ${\lambda }_1\le \cdots \le {\lambda }_n$ be the generalized eigenvalues of $(A,B)$. Fix $1\le q<n$ with the strict gap ${\lambda }_q<{\lambda }_{q+1}$ and let $U=\mathrm{span}\{u_1,\dots ,u_q\}$ be the corresponding sum of the first $q$ eigenspaces, with $(u_i)$ $B$-orthonormal. For a $q$-dimensional trial subspace $S$ with $B$-orthonormal basis $Q$, define $P_{U}x={\sum }_{i=1}^q{u}_i(u_i^{\ast }Bx),\sin \vartheta ={\sup }_{\begin{array}{c}x\in S\\ \mathrm{\parallel }x{\mathrm{\parallel }}_B=1\end{array}}\mathrm{\parallel }(I-P_U)x{\mathrm{\parallel }}_B\mathrm{.}$ Let ${\theta }_1(S)\le \cdots \le {\theta }_q(S)$ be the generalized Ritz values from Generalized Ritz values obey min-max bounds and are monotone under trial-space inclusion. Then ${\lambda }_j\le {\theta }_j(S)(1\le j\le q),{\theta }_q(S)\le {\lambda }_q+({\lambda }_n-{\lambda }_q){\sin }^2\vartheta \mathrm{.}$ The strict gap identifies $U$ uniquely; the estimate is a subspace enclosure, not an iteration convergence rate.

Let $A=A^{\ast }$ and $B=B^{\ast }>0$. Suppose one successful update of Finite-budget block LOBPCG for a Hermitian-definite pencil with explicit rank-loss status starts from a $q$-dimensional $B$-orthonormal block $X_k$ with selected Ritz values ${\theta }_1(X_k)\le \cdots \le {\theta }_q(X_k)$. The update forms a trial space $S_k=\mathrm{span}(X_k,W_k,D_k)$, retains the old block, and selects the $q$ smallest generalized Ritz values on $S_k$. Then ${\theta }_j(X_{k+1})\le {\theta }_j(X_k)(1\le j\le q)\mathrm{.}$ This assertion applies only to an update that reaches the Rayleigh-Ritz step; the algorithm's rank-loss and budget-exhausted exits make no update claim.

Let $A\in {\mathbb{C}}^{n\times n}$, $n\ge 1$, choose a shift $\sigma \in \mathbb{C}$, and put $M=A-\sigma I$. This procedure returns candidate values and their residual data; it asserts no convergence or forward-error bound. Its branch statuses are available, arnoldi-breakdown, linear-solve-unavailable, zero-transformed-Ritz-value, infinite-harmonic-value, projected-pencil-unresolved, qz-budget-exhausted, and qz-certificate-inconclusive.

Shift-invert branch. This branch requires $M$ to be invertible, equivalently that $\sigma$ is not an eigenvalue of $A$ by For a finite-dimensional space, $\lambda$ is an eigenvalue of $T$ if and only if $T-\lambda I$ is not invertible. Its inputs also include a nonzero starting vector $v$, an integer Arnoldi dimension $1\le m\le n$, a residual tolerance $\tau \ge 0$, and a solve routine. Normalize $v_1=v\mathrm{/}\mathrm{\parallel }v{\mathrm{\parallel }}_2$. For $j=1,\dots ,m$, request either the exact solution $y_j=M^{-1}{v}_j$ or a computed vector $y_j$ together with its actual equation residual $s_j=M{y}_j-v_j$. If a requested solve is unavailable, return linear-solve-unavailable with the residuals already recorded. Set $w=y_j$. For $i=1,\dots ,j$, compute $h_{ij}=⟨w,v_i⟩$ and replace $w\leftarrow w-h_{ij}{v}_i$. Repeat this pass once, adding each second-pass coefficient to the corresponding $h_{ij}$. Set $h_{j+1,j}=\mathrm{\parallel }w{\mathrm{\parallel }}_2$. If this norm is zero, set the current dimension to $k=j$ and stop with status arnoldi-breakdown; otherwise set $v_{j+1}=w\mathrm{/}{h}_{j+1,j}$. After $m$ columns, set $k=m$ and stop with status available. When every solve is exact, this is Arnoldi applied to $M^{-1}$. When any solve is inexact, the projected matrix is only the one produced from the computed vectors; the procedure makes no exact-Arnoldi or convergence claim for that branch.

Compute the Ritz pairs $(\theta ,z)$ of the current $k\times k$ matrix $H$ as in Ritz values and Ritz vectors extracted from the Arnoldi Hessenberg reduction, normalize $z$ to unit norm, and lift $u=V_{k}z$. For each $\theta \mathrm{\ne }0$, report the transformed candidate $\lambda =\sigma +1\mathrm{/}\theta$, the lifted vector $u$, and the directly computed original residual $r_A=Au-\lambda u$. Mark this candidate accepted exactly when $\mathrm{\parallel }{r}_A{\mathrm{\parallel }}_2\le \tau$. If $\theta =0$, assign no finite inverse map and mark that candidate zero-transformed-Ritz-value. To record the solve-error transfer for a candidate with $\theta \mathrm{\ne }0$, request one additional solve $My=u$; if it is available, store $s_0=My-u$ and $r_T=y-\theta u$ alongside the directly computed residual $r_A$. If this additional solve is unavailable, return linear-solve-unavailable with the candidate and accumulated ledger.

Harmonic Ritz branch. This branch takes a full-column-rank trial basis $U\in {\mathbb{C}}^{n\times k}$, $k\ge 1$; it does not require $M$ to be invertible. Form $W=MU$ and $G=W^{\ast }W=U^{\ast }{M}^{\ast }MU,C=W^{\ast }U=U^{\ast }{M}^{\ast }U\mathrm{.}$ Compute an exact column-pivoted QR of $W$ as in Full, reduced, and column-pivoted computational QR factorisations and let $r_W$ be its exact rank. If $r_W<k$, choose $z\mathrm{\ne }0$ with $Wz=0$. Then $z^{\ast }G=(Wz{)}^{\ast }W=0$ and $z^{\ast }C=(Wz{)}^{\ast }U=0$. Thus the nonzero row $z^{\ast }$ annihilates $G-\eta C$ for every $\eta$; if this matrix were invertible, right multiplication by its inverse would force $z^{\ast }=0$, a contradiction. Hence $G-\eta C$ is singular for every $\eta$; return projected-pencil-unresolved. If $W$ has full column rank, then $z^{\ast }Gz=\mathrm{\parallel }Wz{\mathrm{\parallel }}_2^2>0$ for every $z\mathrm{\ne }0$. Thus $G$ is positive definite and the pencil $G-\eta C$ is regular in the convention of Regular matrix pencils and homogeneous generalized eigenvalues, since $\det (G-0C)=\det G\mathrm{\ne }0$.

Supply a finite QZ sweep budget $N\ge 0$ and tolerance ${\tau }_{QZ}\ge 0$, and run Finite-budget QZ iteration with homogeneous status and residual output on $(G,C)$. Compute $G,C$ exactly in the mathematical procedure; a finite implementation must include their formation errors in the perturbation ledger or return qz-certificate-inconclusive. On qz-computed, require its exact-factor and local-update certificates; if either is unavailable, return qz-certificate-inconclusive. With these certificates, Every finite-budget QZ run returns a status and records unitary transforms and explicit deflations then records equivalence to a perturbed projected pencil. For each diagonal index $i$, set ${\alpha }_i=S_{ii}$ and ${\beta }_i=T_{ii}$. If $({\alpha }_i,{\beta }_i)=(0,0)$, return projected-pencil-unresolved; if ${\beta }_i\mathrm{\ne }0$, report the computed harmonic Ritz value ${\lambda }_i=\sigma +{\alpha }_i\mathrm{/}{\beta }_i$; if ${\beta }_i=0$, mark that candidate infinite-harmonic-value. These are values of the computed perturbed pencil, with the QZ residual and perturbation ledger retained; no exact-root, eigenvector, or forward-error claim is made. After processing all diagonal candidates, return available with the candidate list, its per-candidate marks, and the residual and perturbation ledger. Map a QZ budget-exhausted status to qz-budget-exhausted; map a QZ certificate-inconclusive or range-failure to qz-certificate-inconclusive; map singular-pencil-suspected to projected-pencil-unresolved.

Consider a valid shift-invert or harmonic-extraction execution of Shift-invert and harmonic Ritz extraction with solve and pencil status. Assume each invoked solve, projected eigensolve, finite-dimensional arithmetic, QR, and orthogonalization operation returns its value in finite work, except that a solve may return unavailable and a harmonic-branch QZ call may return one of its listed statuses. On the harmonic branch, also assume the finite-budget QZ run satisfies the operationwise arithmetic hypotheses of Every finite-budget QZ run returns a status and records unitary transforms and explicit deflations. Then the execution returns one of the statuses listed by the algorithm after finitely many operations.

Let $A\in {\mathbb{C}}^{n\times n}$, $\sigma \in \mathbb{C}$, and $M=A-\sigma I$ be invertible. Then a nonzero vector $u$ is an eigenvector of $A$ with eigenvalue $\lambda$ if and only if it is an eigenvector of $M^{-1}$ with eigenvalue $\theta =(\lambda -\sigma {)}^{-1}$. In particular, $\lambda =\sigma +1\mathrm{/}\theta$ and $\theta \mathrm{\ne }0$.

For any nonzero candidate $u$, any $\theta \mathrm{\ne }0$, and any computed vector $y$ with solve residual $s_0=My-u$, define the transformed residual $r_T=y-\theta u$ and $\lambda =\sigma +1\mathrm{/}\theta$. Then $(A-\lambda I)u={\theta }^{-1}(s_0-M{r}_T),$ and consequently $\mathrm{\parallel }(A-\lambda I)u{\mathrm{\parallel }}_2\le \frac{\mathrm{\parallel }{s}_0{\mathrm{\parallel }}_2+\mathrm{\parallel }M{r}_T{\mathrm{\parallel }}_2}{\mathrm{\mid }\theta \mathrm{\mid }}\mathrm{.}$