Simultaneous Hessenberg Process

• Simultaneous Hessenberg Process is a Krylov method that constructs a single non-orthogonal basis to address multiple shifted systems through shift-invariance.
• It leverages a restarted shifted formulation and an extended block approach to handle both single and multiple right-hand sides efficiently.
• The method minimizes computational cost by reusing basis generation across shifts, while pivoting ensures numerical robustness.

The simultaneous Hessenberg process denotes a family of Hessenberg-type Krylov procedures in which a single non-orthogonal basis is constructed once and then reused across a coupled set of problems. In the restarted shifted formulation, the central task is the family of nonsymmetric systems

$(A+\sigma_i I)x_i=b,\qquad i=1,\ldots,\nu,$

with a common right-hand side, where shift-invariance makes one Krylov basis sufficient for all shifts when the initial guesses are chosen consistently. In the extended block formulation, a single unit lower trapezoidal block basis simultaneously spans positive and negative powers of $A$ applied to a block $V$, so that the same projected matrix can be reused for both $f(A)V$ and shifted systems with multiple right-hand sides. These two formulations give the term “simultaneous” two precise meanings: simultaneous treatment of many shifts, and simultaneous treatment of many right-hand sides together with both $A$ and $A^{-1}$ actions (Gu et al., 2015, Bentbib et al., 2024).

1. Scope and historical placement

The modern simultaneous Hessenberg viewpoint lies at the intersection of generalized Hessenberg reduction, restarted shifted Krylov methods, and block or extended projection methods. The foundational source identified for the Hessenberg procedure is M. Heyouni’s 1996 thesis. The restarted shifted framework is explicitly connected with V. Simoncini’s restarted full orthogonalization method for shifted systems, while the CMRH line of work of H. Sadok and the later analysis by H. Sadok and D.B. Szyld provide further background. The stability discussion in the literature also invokes J.H. Wilkinson and Howell–Stephens.

In this literature, simultaneous treatment is not tied to orthonormal bases. On the contrary, the distinctive feature of the Hessenberg approach is that it replaces Arnoldi orthonormalization by a pivoted non-orthogonal construction. In the single-vector restarted setting, this yields a Krylov basis $V_m=[v_1,\ldots,v_m]$ and an upper Hessenberg matrix $H_m$. In the block-extended setting, it yields a unit lower trapezoidal block basis $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$ and a block upper Hessenberg projected matrix $\mathbb{T}_{2m}$ (Gu et al., 2015, Bentbib et al., 2024).

Formulation	Shared basis object	Small projected problem
Restarted shifted Hessenberg	$A$0 for all $A$1	$A$2
Extended block Hessenberg	$A$3 for all shifts and block RHS	$A$4 or $A$5

A common misunderstanding is that simultaneous treatment requires building separate Krylov spaces for each shifted problem. The shifted Hessenberg framework uses the opposite fact: when the right-hand side is the same across shifts and initial guesses are chosen consistently, the Krylov subspaces generated by $A$6 and by $A$7 are identical.

2. Core Hessenberg relations and shift-invariance

The basic Hessenberg procedure builds a non-orthogonal Krylov basis with pivoting and produces the decomposition

$A$8

This is the analogue of the Arnoldi relation, but it is associated with an oblique projection rather than an orthogonal one. The $A$9-th basis vector has its first $V$0 entries zero by construction, and pivoting is used to maintain numerical robustness (Gu et al., 2015).

For shifted problems, the same decomposition immediately yields

$V$1

Because

$V$2

for any $V$3, a single basis $V$4 generated from the seed system is sufficient for all shifts. For each $V$5, the projected Galerkin system is

$V$6

with $V$7, and the approximation is updated by

$V$8

The residual formula is the essential structural fact: $V$9 Hence every shifted residual is a scalar multiple of the same next basis vector $f(A)V$0. This makes the simultaneous process viable under restart, because the new cycle can be seeded from a common residual direction rather than from $f(A)V$1 unrelated residual vectors (Gu et al., 2015).

3. Restarted shifted Hessenberg method

The restarted shifted Hessenberg method is designed for practical settings in which the Krylov dimension must be capped by a restart length $f(A)V$2. A cycle begins from $f(A)V$3, typically $f(A)V$4, with $f(A)V$5. The Hessenberg process with pivoting is applied once to build $f(A)V$6 and $f(A)V$7. Then, for each unconverged shift, one solves the reduced system $f(A)V$8, updates $f(A)V$9, computes the scalar

$A$0

and sets

$A$1

Converged systems are removed, and the next cycle is started from the common direction $A$2, appropriately scaled (Gu et al., 2015).

The method’s main theorem states that after each cycle all residuals are collinear with $A$3 and therefore with each other. In the unshifted seed system, the Galerkin condition with oblique test space gives

$A$4

For shifted systems the same argument yields the shifted residual formula above, so the restart preserves a single residual direction across all systems. The proof is constructive: inserting the Hessenberg relation into the residual expression eliminates the $A$5 contribution and leaves only the next-basis direction.

A polynomial interpretation clarifies restart propagation. Since $A$6, one has

$A$7

for a monic polynomial $A$8 of degree $A$9. Under restart,

$A^{-1}$0

Because the basis for the $A^{-1}$1-th cycle is generated from the seed residual $A^{-1}$2, both the seed and shifted residuals at the end of the cycle remain in the same one-dimensional direction. This is the algebraic reason that simultaneous treatment survives restarting (Gu et al., 2015).

This residual-collinearity mechanism sharply distinguishes restarted shifted Hessenberg and restarted shifted FOM from restarted GMRES. GMRES loses collinearity after restart unless additional constraints are imposed, and those constraints compromise minimal residual properties for non-seed shifts and increase cost. FOM preserves collinearity, but its Arnoldi orthonormalization is more expensive per cycle.

4. Block and extended simultaneous Hessenberg constructions

The extended block Hessenberg process generalizes the simultaneous idea to a block $A^{-1}$3, with $A^{-1}$4, and to an extended Krylov space that contains both positive and negative powers of $A^{-1}$5. The space used is

$A^{-1}$6

This is the sum of two block Krylov subspaces, one for $A^{-1}$7 and one for $A^{-1}$8. The basis is not orthonormal; instead, it is unit lower trapezoidal and is obtained by PLU factorizations with pivoting (Bentbib et al., 2024).

Initialization is already block-structured: $A^{-1}$9 The second block is formed from $V_m=[v_1,\ldots,v_m]$0 after removing coupling with the first block, followed by another PLU factorization. Thereafter, for each step $V_m=[v_1,\ldots,v_m]$1, the process generates a pair of blocks: $V_m=[v_1,\ldots,v_m]$2 from the positive direction $V_m=[v_1,\ldots,v_m]$3, and $V_m=[v_1,\ldots,v_m]$4 from the negative direction $V_m=[v_1,\ldots,v_m]$5. Galerkin-like orthogonality is imposed not by Gram–Schmidt but against selector matrices $V_m=[v_1,\ldots,v_m]$6 that correspond to pivot rows.

The resulting block decomposition can be written through the left-inverse

$V_m=[v_1,\ldots,v_m]$7

and the projected matrix

$V_m=[v_1,\ldots,v_m]$8

The simultaneous Hessenberg relation is then

$V_m=[v_1,\ldots,v_m]$9

A corresponding relation also holds for $H_m$0.

This construction supports two major tasks. For matrix functions, the approximation is

$H_m$1

The approximation is exact for Laurent polynomials of positive degree $H_m$2 and negative degree $H_m$3, and for $H_m$4 with $H_m$5 dissipative the paper gives an a priori bound involving the residual block $H_m$6. For shifted systems with multiple right-hand sides, one seeks

$H_m$7

with reduced problem

$H_m$8

and residual

$H_m$9

Thus one block basis and one projected matrix can be reused across many shifts $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$0 and all $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$1 right-hand sides simultaneously (Bentbib et al., 2024).

5. Arithmetic cost, storage, and numerical behavior

The principal advantage of the Hessenberg process over Arnoldi is lower orthogonalization cost. In the restarted shifted single-vector setting, the per-cycle basis-construction cost is

$\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$2

operations, compared with

$\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$3

for Arnoldi/FOM/GMRES. Storage is $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$4 for $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$5 plus $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$6 for $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$7. For $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$8 shifts, the projected stage requires $\mathbb{V}_{2m}=[V_1,\ldots,V_{2m}]$9 independent $\mathbb{T}_{2m}$0 solves; if factored independently, this is $\mathbb{T}_{2m}$1, which is cheap when $\mathbb{T}_{2m}$2 is moderate. Residual norms are available without forming full residual vectors: $\mathbb{T}_{2m}$3 and analogously for shifted systems (Gu et al., 2015).

In the extended block case, the paper enumerates block costs explicitly: $\mathbb{T}_{2m}$4 costs $\mathbb{T}_{2m}$5; applying $\mathbb{T}_{2m}$6 to a block via LU/backslash costs $\mathbb{T}_{2m}$7; LU factorization of an $\mathbb{T}_{2m}$8 matrix costs $\mathbb{T}_{2m}$9; computing each $A$00 via LU on $A$01 costs $A$02; and multiplying $A$03 costs $A$04. The resulting total flop estimate for EBHA is given explicitly in the paper, while storage is $A$05 for the $A$06 basis blocks together with the block upper Hessenberg matrix (Bentbib et al., 2024).

The price of these savings is non-orthogonality. The Hessenberg process with pivoting is not provably backward stable in finite precision, but the reported backward errors are usually small in practice, and pivoting is used specifically to mitigate near-breakdowns. If $A$07 before $A$08, the process terminates early and the current subspace may already deliver convergence. In the block-extended case, breakdown is described as “very rare” for small $A$09. Numerical safeguards mentioned in the literature include partial or mild reorthogonalization, re-normalization of $A$10, monitoring of the $A$11, and increased pivoting vigilance (Gu et al., 2015, Bentbib et al., 2024).

Preconditioning is delicate. In shifted systems, preserving shift-invariance requires a preconditioner $A$12 such that $A$13 across shifts. Polynomial preconditioners or nested Krylov strategies can be effective, but flexible preconditioning may destroy collinearity unless it is specially designed. In the extended block setting, repeated $A$14 applications are intrinsic; the paper uses LU/backslash and notes that factorizing $A$15 once and reusing sparse triangular solves is essential to efficiency.

6. Applications and empirical profile

The restarted shifted Hessenberg method was tested on UF sparse matrices with real negative shifts, on QCD Dirac operator discretizations with complex shifts, on a $A$16D heat equation matrix-exponential action computed through the rational approximation

$A$17

and on fractional convection–diffusion–reaction problems using the Carathéodory–Fejér rational approximation for $A$18. In these experiments, $A$19Hessen$A$20 consistently required similar or fewer matrix-vector products than $A$21FOM$A$22 and notably less CPU time, while weighted $A$23FOM frequently failed to converge. In the heat-equation experiment, the setup used $A$24 complex shifts and $A$25. In the fractional problem, $A$26Hessen$A$27 outperformed the competing methods, and several short-recurrence methods failed on some instances. Typical settings were $A$28, maximum MVP counts capped at values such as $A$29, and restart lengths $A$30–$A$31 (Gu et al., 2015).

The extended block Hessenberg paper reports matrix-function experiments with $A$32 and $A$33 for

$A$34

For $A$35 and $A$36, MF-EBH required less CPU time than MF-EBA while producing relative errors in the $A$37–$A$38 range, depending on the function and $A$39. In a stiff tridiagonal example with $A$40, MF-EBH reached the target relative error $A$41 in fewer seconds than MF-EBA, although MF-EBA converged in slightly fewer iterations. For shifted systems with $A$42 shifts in $A$43 and $A$44, restarted-EBH required fewer seconds than restarted-EBA on PDE-like discretizations and the SuiteSparse matrices add32, epb1, and memplus, with comparable numbers of restarts and residuals near $A$45. Both iterative methods significantly outperformed direct Gaussian elimination over all shifts (Bentbib et al., 2024).

Taken together, these results identify the simultaneous Hessenberg process as a computational strategy rather than a single narrow algorithm. Its defining principle is the reuse of one Hessenberg-type basis across a structured family of problems. In the shifted single-vector case, the structure is residual collinearity across shifts and across restarts. In the extended block case, the structure is simultaneous representation of multiple right-hand sides together with positive and negative powers of $A$46. A plausible implication is that the method is most attractive when basis construction and orthogonalization dominate wall-clock time, when many shifts or right-hand sides are present, and when the reduced problems remain small enough that their independent solution is negligible relative to sparse matvecs or sparse triangular solves.