Bi-Lanczos Algorithm Overview

• Bi-Lanczos algorithm is a Krylov subspace method for non-Hermitian operators that builds paired biorthogonal bases and a tridiagonal projection for efficient computation.
• It employs coupled three-term recurrences to generate the right and left Krylov bases, enabling matrix function approximations and eigenvalue extraction.
• Extensions include tensor generalizations, look-ahead variants, and reduced-order modeling to address breakdowns and improve numerical stability.

The Bi-Lanczos algorithm is a foundational Krylov subspace method for non-Hermitian (nonsymmetric) linear operators, generalizing the classical Hermitian Lanczos algorithm via the construction of biorthogonal sequences. It is central to iterative solvers, matrix function approximation, quantum chaos diagnostics, tensor computations, and algorithms for nonlinear eigenvalue problems in both mathematical and physical sciences. The bi-Lanczos process formulates a pair of coupled three-term recurrences that generate two mutually biorthogonal Krylov bases—respecting the intrinsic structure of non-Hermitian operators—and yields a tridiagonal (or block-tridiagonal) projection that underpins both theoretical analysis and efficient computation.

1. Bi-Lanczos Recursion and Biorthogonal Krylov Bases

The bi-Lanczos algorithm constructs two sequences of vectors, $\{|q_n\rangle\}$ and $\{|p_n\rangle\}$, which form biorthonormal bases for the right and left Krylov subspaces associated with a non-Hermitian operator $H$ (or, more generally, matrix $A$): $\langle p_n | q_m \rangle = \delta_{nm} \ .$ These bases satisfy recursive relations involving complex coefficients $\{a_n, b_n, c_n\}$: $$\begin{aligned} |r_{n+1}\rangle &= (H - a_n) |q_n\rangle - b_n |q_{n-1}\rangle,\quad |q_n\rangle = c_n^{-1} |r_n\rangle \ |l_{n+1}\rangle &= (H^\dagger - a_n^*) |p_n\rangle - c_n^* |p_{n-1}\rangle, \quad |p_n\rangle = (b_n^*)^{-1} |l_n\rangle \ . \end{aligned}$$ The tridiagonal structure emerges as the representation of $H$ in the biorthogonal Krylov basis: $$T = \begin{pmatrix} a_0 & b_1 & 0 & \cdots & 0 \ c_1 & a_1 & b_2 & \ddots & \vdots \ 0 & c_2 & a_2 & \ddots & 0 \ \vdots & \ddots & \ddots & \ddots & b_{n-1} \ 0 & \cdots & 0 & c_{n-1} & a_{n-1} \end{pmatrix} .$$ This structure underlies projection methods, matrix functions, and spectral diagnostics.

2. Algorithmic Properties, Stability, and Generalizations

The standard bi-Lanczos process, as described above, forms the basis for iterative solvers such as BiCG, QMR, BiLQ, and extensions to nonlinear eigenproblems and time-dependent operator equations.

Algorithmic workflow:

1. Initialization: Choose starting vectors $|q_0\rangle$ and $|p_0\rangle$ such that $\langle p_0|q_0 \rangle = 1$.
2. Recursion: Iteratively compute $|q_{n+1}\rangle$, $|p_{n+1}\rangle$, and coefficients $a_n$, $b_n$, $c_n$ via biorthogonalization.
3. Breakdown: The process can suffer breakdowns (e.g., vanishing biorthogonality), but look-ahead or block extensions circumvent most practical obstacles (Pozza et al., 2019).
4. Extension to tensors: The formalism generalizes to higher-order tensors, with biorthogonal bases built through tensor contractions, extending applicability to operator-valued ODEs and bilinear forms on time-ordered exponentials (Cipolla et al., 2022).

Novel algorithmic families can be constructed by leveraging higher-degree recurrence relations involving orthogonal polynomials with degree differences beyond one. New recurrences of "A-type" and "B-type" allow for combinatorial flexibility and can be leveraged to respond to breakdowns or tailor the method to specific matrix structures, though with trade-offs between storage and computational cost (Farooq et al., 2014).

3. Connections to Gauss Quadrature and Realization Theory

The bi-Lanczos process serves as the finite-dimensional backbone of general Gauss quadrature and minimal partial realization for arbitrary (possibly indefinite or complex) linear functionals (Pozza et al., 2019). For a linear functional $\mathcal{L}(p) = w^* p(A) v$ with moments $m_j = \mathcal{L}(\lambda^j)$, the block-tridiagonal matrix $T_n$ produced by bi-Lanczos (with possible look-ahead) encodes the quadrature: $$\mathcal{G}_n(f) = \mu m_{\nu(1)-1} e_1^T f(T_n) e_{\nu(1)}$$ and satisfies the "matching moment property": $$\mu m_{\nu(1)-1} e_1^T T_n^k e_{\nu(1)} = m_k, \quad k=0, \dots$$ In cases of incurable breakdown, the spectrum of $T_n$ coincides with a subset of $A$'s spectrum ("mismatch theorem"). This linkage also forms the minimal partial realization in systems theory, enabling reduced-order models grounded in moment matching.

4. Krylov Complexity, Quantum Chaos, and Non-Hermitian Systems

For diagnostics of quantum chaos in open quantum systems, the bi-Lanczos algorithm supplies a principled platform for computing Krylov complexity (KC) when the generator is non-Hermitian due to dissipative effects or environmental coupling (Baggioli et al., 19 Aug 2025, Bhattacharya et al., 2023).

Key aspects:

• Operator growth in the biorthogonal Krylov basis quantifies spreading under nonunitary evolution.
• Krylov complexity for time-evolved states, using bi-Lanczos amplitudes $\Phi_n^{p}(t), \Phi_n^{q}(t)$:

$$C(t) = \sum_n n \left| \tilde{\Phi}_n^{p*}(t) \tilde{\Phi}_n^q(t) \right|$$

• Chaos diagnostics: KC displays a pronounced early-time peak in chaotic regimes but not integrable ones, matching universal spectral statistics (Ginibre vs. Poisson), with the universality further encoded in coefficient relations (e.g., $1/\sqrt{2}|a_n| \approx |b_n| = c_n$).
• Superiority over SVD-based approaches: Only the bi-Lanczos method, which maintains the full biorthogonal structure, provides a robust diagnostic in non-Hermitian settings, outperforming SVD-based heuristics which fail to distinguish chaos from integrability.

In dissipative systems described by Lindblad dynamics, the bi-Lanczos tridiagonalizes the generator for efficient and interpretable computation of complexity and spectral features. In the strong dissipation regime, distinctions between chaotic and integrable behavior vanish at late times—operator complexity equilibrates to a universal saturated value, indicating the breakdown of late-time chaos as a meaningful concept in such open quantum systems (Bhattacharya et al., 2023).

5. Numerical Robustness and Exactness in Finite Precision

For specific model problems—when the operator is a signed permutation of a tridiagonal matrix and the initial vectors are aligned accordingly—the bi-Lanczos algorithm computes all vectors and tridiagonal coefficients \emph{exactly} in standard floating point arithmetic. There is no accumulation of rounding error, orthogonality loss, or spurious spectral pollution, provided the coefficients remain within representable range (Šimonová et al., 2021). This property is invaluable for algorithm validation and theoretical investigations, as it equates practical numerics with idealized mathematics on such instances.

Algorithm	Structure for Exactness	Floating Point Robustness
Symmetric Lanczos	$A = PTP^\top$, $v = P e_1$	Yes
Bi-Lanczos (nonsymmetric)	$A = PTP^\top$, $v, w = P e_1$	Yes

For general inputs, the algorithm is subject to the usual vulnerabilities to rounding error and instability; the exactness result holds specifically for structured instances.

6. Applications: Linear Solvers, Eigenproblems, Tensor Algorithms, and Beyond

The bi-Lanczos framework underpins a rich hierarchy of methods:

• Linear systems: BiCG, QMR, BiLQ (the latter being robust to singular tridiagonal projections and ill-conditioning, ensuring the existence of quasi-minimum error approximations in an implicit biorthogonal norm) (Montoison et al., 2019).
• Eigenvalue problems: Extraction of eigenvalues and eigenvectors via the tridiagonal projection; repeated and restarted approaches (eigBiCG) for large-scale and multiple right-hand sides (Abdel-Rehim et al., 2013).
• Nonlinear eigenproblems: Infinite-dimensional bi-Lanczos provides short-recurrence, block-wise representations encoding nonlinear structure with efficient storage and computation (Gaaf et al., 2016).
• Tensor computations: Generalization to 4-mode tensors enables efficient approximation of matrix-valued ODEs and time-ordered exponentials, preserving moment-matching and short-recurrence properties, and improving tractability for high-dimensional and operator-valued problems (Cipolla et al., 2022).

Application Domain	Bi-Lanczos Role
Linear systems ($Ax=b$)	Basis for BiCG, QMR, BiLQ, robust error minimization
Non-Hermitian quantum dynamics	Chaos/integrability diagnostics, Krylov complexity
Time-ordered exponentials, bilinear forms	Tensor generalization, time-dependent ODEs
Nonlinear eigenvalue problems	Infinite-dimensional companion linearization

7. Universality, Extensions, and Ongoing Developments

Universal empirical relations in bi-Lanczos coefficients signal underlying random matrix universality in quantum chaos. Across models—non-Hermitian Sachdev-Ye-Kitaev (nHSYK) and GinUE ensembles—the same diagnostic behaviors and coefficient relations manifest for a broad class of dynamics, establishing the method’s generality and versatility (Baggioli et al., 19 Aug 2025).

Recent research explores:

• Look-ahead and block variants for more robust operation under breakdown conditions (Pozza et al., 2019).
• Formal orthogonal polynomial generalizations with higher-degree recurrence relations, broadening the landscape of feasible Lanczos-type algorithms and iterative solvers (Farooq et al., 2014).
• Reduced-order modeling in inverse problems, where bi-Lanczos-type Gram-Schmidt processes yield weakly medium-dependent bases critical for stabilized inversion, and extensions to Bi-Lanczos appear for non-Hermitian or non-self-adjoint scenarios (Baker et al., 2023).

A plausible implication is the increasing unification of theory and numerics in non-Hermitian problem domains, driven both by problems in quantum science and large-scale computational mathematics.

---

The bi-Lanczos algorithm thus forms the theoretical and algorithmic backbone for robust and universal Krylov subspace methods in nonsymmetric settings, with deep connections to orthogonal polynomials, quadrature, moment matching, operator complexity, system realization, and modern computational science.