PaperlayerView exact revision on arXiv

Source-linked AI summary

On the stability of basis construction in extended Krylov subspace methods

Kathryn Lund, Marcel Schweitzer

arXiv:2610.21603v1math.NAcs.NA

TL;DR

The space Em,m(A, b) has a unique shape relative to Km(A, b) or Qm(A, b), requiring unique implementation decisions. The paper classifies mixed- and split-channel EKSM variants and analyzes their loss of orthogonality using an easily computable metric. For symmetric positive definite A, mixed channels have better orthogonality than split channels, while experiments show no predictable ranking for indefinite or nonsymmetric matrices.

  • Problem

    The space Em,m(A, b) has a unique shape relative to Km(A, b) or Qm(A, b), requiring unique implementation decisions.

    Source: §1 · p. 2

  • Method

    The paper classifies mixed- and split-channel EKSM variants and analyzes their loss of orthogonality using an easily computable metric.

    Source: §3 · p. 10, §5 · pp. 20–21

  • Results

    For symmetric positive definite A, mixed channels have better orthogonality than split channels, while experiments show no predictable ranking for indefinite or nonsymmetric matrices.

    Source: §1 · p. 2, §5 · p. 21

  • Takeaways & Limitations

    When in doubt, fully orthogonalizing the basis reduces the strong effects of starting vector, T-method, and orthogonalization-backbone choices.

    Source: §5 · p. 21

    Mixed-channel Lanczos convergence on nonsymmetric problems was tested on only two matrices and two functions.

    Source: §5 · p. 21

Abstract

from arXiv · show

Extended Krylov subspace methods (EKSMs) are a well known approach for approximating the action of a matrix function on a vector $f(A)\mathbf{b}$, and their block versions are often used for solving Lyapunov and Sylvester matrix equations. EKSMs produce Laurent approximations to the underlying operator, which are generally more accurate per iteration than polynomial KSMs, without requiring a pole selection procedure as rational KSMs do. At the same time, due to the expansion of the search space by not only $A$ but also by its inverse, a number of implementation choices must be made that can affect the stability of the algorithm in floating-point arithmetic. We concentrate on two variants that often appear in the literature and in code, based on whether new basis vectors are built via what we denote as ``mixed channels" or ``split channels." We characterize their loss of orthogonality with an easily computable metric and prove that when $A$ is symmetric positive definite, mixed channels have better orthogonality than split channels; however, numerical experiments demonstrate that straightforward predictions are not possible for indefinite or nonsymmetric matrices. We also provide a modular MATLAB toolbox, \texttt{ektoolbox}, that contains new versions of EKSMs with low-synchronization Gram-Schmidt backbones.

1. Introduction.

EKSMs offer a balance between polynomial and rational Krylov methods, but expanding the search space with both A and A^-1 creates implementation choices that can affect numerical stability. This paper classifies those choices and analyzes their effects on basis orthogonality, with results that depend on the matrix setting.

1. Introduction.

  • EKSMs balance polynomial methods with rational methods, avoiding rational methods’ pole-selection procedure while offering more accurate approximations per iteration than polynomial methods.Source: §1 · p. 1
  • The analysis addresses operator ordering and batch versus individual vector generation, assumes access to accurate factorizations, and leaves linear-system solution choices for future work.The study assumes equal numbers of applications of A and A^-1.Source: §1 · p. 2
  • Mixed channels apply A or A^-1 to the most recently generated basis vector; split channels apply each operator to its own channel and can batch both directions.Source: §1 · p. 2
  • For symmetric positive definite A, mixed channels have better orthogonality than split channels; experiments do not yield predictable rankings for indefinite or nonsymmetric A.The comparison is based on an easily computable metric.Source: §1 · p. 2
  • For symmetric A, mixed-channel methods have 3-term recurrences, whereas split-channel methods have 5-term recurrences.The paper also shows that the two approaches span the same space and produce essentially the same basis up to signs.Source: §1 · p. 2

2. Implementation variants.

This section presents a position-vector framework for categorizing EKSM implementation variants by operator ordering and whether A and A⁻¹ act in split or mixed channels, and compares their basis constructions. It also places prior algorithms in these categories, noting that most identified implementations use split channels.

2. Implementation variants.

  • Two choices define four EKSM variants: whether A or A^-1 leads, and whether operators act within separate channels or on the most recently computed vector; operator order does not strongly affect stability.The authors briefly examine both leading-operator choices in separate sections.Source: §2 · pp. 3–4

2.1. Power-leading variants.

  • The power-leading MCPL process applies A first and then A^-1, but the authors do not pursue power-leading variants further, except noting their possible use when relations in A^-1 are needed.MCPL normalizes b, applies A to the first basis vector, and then applies A^-1 to the next.Source: §2.1 · p. 4

2.2. Inverse-leading variants.

  • Assuming no breakdown, MCIL and SCIL generate bases equal up to sign and spanning the same extended Krylov subspace.The sign relating corresponding basis vectors is ±1.Source: §2.2 · p. 6
  • Despite generating essentially the same basis, MCIL and SCIL differ in finite precision; for symmetric A, MCIL has a three-term recurrence while SCIL has a five-term recurrence.The paper cautions that implementation variants are not equivalent once floating-point arithmetic is used.Source: §2.2 · pp. 7–8

2.3. Categorization of variants in the literature.

  • Prior algorithms map to these categories: Druskin–Knizhnerman’s scheme is SeqSCIL, K-PIK is SCIL, and Jagels and Reichel mix channels, unlike Schweitzer’s SCIL approach.The classification also includes a block SCIL method for differential Lyapunov equations [46].Source: §2.3 · p. 9
  • For shifted-inverse variants, is classified as SeqMCIL and as SCIL, while does not describe basis generation enough to determine its channel type.The cited works use (A − γI)^-1 or another shifted inverse in their EKSM formulations.Source: §2.3 · p. 10
  • Mertens and Vandebril’s position-vector framework can describe any Table 2.2 variant and could support runtime adaptation, but choosing an optimal ordering is difficult.Their default ordering appears to be MCPL.Source: §2.3 · p. 10
  • Most implementations identified in the literature follow a split-channel approach.Source: §2.3 · p. 10

3. On the loss of orthogonality.

The section studies how floating-point rounding affects loss of orthogonality (LOO) in mixed- and split-channel EKSM variants. Its metric predicts their relative local error amplification, with a rigorous advantage for mixed channels when A is definite, but no consistent ranking for indefinite matrices.

3. On the loss of orthogonality.

  • The metric estimates local round-off error amplification: normalization can amplify perturbations when the post-orthogonalization remainder is small relative to the candidate vector.It does not account for every source of LOO, including effects of the orthogonalization routine.Source: §3 · pp. 10–11
  • In the general case, LOO increases in MCIL coincide with γk ≪1, while γk > 1 occurs less often and only mildly in Example 3.1.This suggests MCIL may sometimes fare worse than SCIL, while the opposite behavior seems less likely in that example.Source: §3 · pp. 11–12, Figure 3.1 · p. 12
  • For the symmetric positive definite Example 3.2, γk > 1 for all k ≥3, and SCIL consistently has larger LOO than MCIL.The figure compares LOO and γk across the two variants.Source: §3 · p. 13, Figure 3.2 · p. 13
  • For symmetric positive definite A, γ2 = 1 and γk > 1 for all k ≥3, indicating larger expected LOO for SCIL from step 3 onward.The same comparison holds for negative definite A after replacing A with −A.Source: §3 · p. 13
  • For symmetric indefinite A, there is no clear MCIL-versus-SCIL LOO ranking, and the comparison can depend on the starting vector.In Example 3.5, the two starting vectors yield different outcomes: SCIL has slightly larger overall LOO for b(1), while MCIL has much larger LOO for b(2).Source: §3 · pp. 13–14, Figure 3.3 · p. 14
  • The instability ratio is a local indicator for round-off amplification at step k, not a measure of accumulated LOO up to that step.The authors conclude that a general comparison like the definite-case result appears difficult for indefinite or nonsymmetric matrices.Source: §3 · p. 14

4. Numerical Experiments.

The experiments compare EKSM variants on small SPD, symmetric-indefinite, and nonsymmetric matrices, tracking basis orthogonality and matrix-function error. Results vary by matrix type: mixed channels fare better on the SPD case, while indefinite and nonsymmetric cases resist straightforward prediction.

Toolbox and methods

  • The modular MATLAB toolbox implements EKSM variants with low-synchronization block Gram-Schmidt backbones; SCIL supports block routines even for a single starting vector, unlike MCIL.The experiments include Lanczos, modified Gram-Schmidt, and classical Gram-Schmidt with inner reorthogonalization.Source: §4 · p. 15

Experimental setup

  • The experiments use n=1000 SPD, symmetric-indefinite, or nonsymmetric matrices and compare LOO with relative error for exp(z) and z^-1/2 in double precision.Reference solutions use MATLAB expm or sqrtm with backslash; the matrices are built from a shared banded matrix.Source: §4 · p. 15, Table 4.1 · p. 15

Example 4.1 (exp(A1)b)

  • For the SPD example, SCIL has worse orthogonality than MCIL across T-computation strategies and orthogonalization backbones; explicit T can diverge unless orthogonalization is unconditionally stable.With bcgs iro, all methods attain machine-precision (O(10^-16)) LOO; MCIL-Lanczos is identified as the preferred approach for this example because of its shorter recurrence.Source: §4 · pp. 15–16, Figure 4.1 · p. 16, Figure 4.2 · p. 16

Example 4.2 (exp(A2)b)

  • For the symmetric-indefinite example, LOO is high for both Lanczos variants regardless of block size, while only MCIL appears to struggle with bmgs; convergence behavior is attributed to how T is computed.The authors identify MCIL-Lanczos as the best approach in this example because of its lower computational costs.Source: §4 · p. 17, Figure 4.3 · p. 17, Figure 4.4 · p. 18

Example 4.3 (exp(A3)b)

  • For nonsymmetric A3, all Lanczos methods quickly lose orthogonality, yet MCIL-blanczos converges while the other Lanczos methods fail or stagnate.The authors recommend SCIL as more reliable in this context and suggest low-synchronization Gram-Schmidt to avoid high orthogonalization costs.Source: §4 · p. 17, Figure 4.5 · p. 18, Figure 4.6 · p. 19

Example 4.4

  • For nonsymmetric A4, SCIL loses more orthogonality than MCIL with bmgs, but MCIL with blanczos can still converge; other non-Lanczos methods converge.A4 is tailored to avoid branch cuts in the matrix square root, and LOO appears to affect convergence only for explicit T.Source: §4 · pp. 18–20, Figure 4.8 · p. 20, Figure 4.9 · p. 20

5. Conclusions.

The paper analyzes how EKSM basis-building choices affect floating-point orthogonality, finding a mixed-channel advantage for symmetric positive definite matrices but less predictable behavior beyond that setting. Its experiments also show that implementation choices matter and that full orthogonalization reduces their effects.

5. Conclusions.

  • The investigation examines how basis-building choices affect loss of orthogonality and attainable accuracy in floating-point arithmetic.The methods are classified by whether new basis vectors are added through split or mixed channels.Source: §5 · pp. 20–21
  • For symmetric positive definite A, mixed channels have better orthogonality than split channels, while experiments show less predictable behavior for indefinite or nonsymmetric matrices.The comparison uses an easily computable ratio; starting vector, choice of T method, and orthogonalization backbone also strongly affect experimental behavior.Source: §5 · p. 21
  • Full orthogonalization is recommended when implementation choices produce strong effects on the basis.The reported effects include those of the starting vector, choice of T method, and orthogonalization backbone.Source: §5 · p. 21
  • Whether mixed channels with a Lanczos backbone converge for some nonsymmetric problems remains open; the authors tested only two matrices and two functions.They expect the modular ektoolbox to support further exploration.Source: §5 · p. 21
Loading 2610.21603v1…