PaperlayerView exact revision on arXiv

Source-linked AI summary

A Fast and Stable Square-Root-Free Unitary Core-Chasing Algorithm

Mónica Esquivel-Rosado, Jared L. Aurentz, Giovanni Barbarino, María C. Quintana

arXiv:2610.12228v1math.NAcs.NA

TL;DR

Unitary core chasing had an unnoticed symmetry, and its matrix-based formulation left that symmetry unexploited. This paper derives a square-root-free variant and analyzes its stability, reporting comparable accuracy with reduced runtime while showing that strict structured backward stability is impossible in general.

  • Problem

    The unitary core-chasing formulation contained a previously unnoticed symmetry.

    Source: §1 · p. 2

  • Method

    The paper identifies the symmetry in a single-shift core-chasing step and uses it to derive a square-root-free variant that preserves near-unitarity without additional square roots.

    Source: §1 · p. 2

  • Results

    The algorithms are stable for nearly-unitary input and achieve accuracy comparable to LAPACK and [3] while reducing overall runtime; strict structured backward stability is impossible in general.

    Source: §1 · p. 2, §7 · p. 25

  • Takeaways & Limitations

    The experiments show that exploiting the symmetry while removing square roots reduces runtime without degrading accuracy.

    Source: §7 · p. 25

    Strict structured backward error is impossible in general, and the square-root-free variant does not store the transformations needed to compute eigenvectors.

    Source: §4.1 · p. 17, §5 · p. 18

Abstract

from arXiv · show

Symmetric and unitary matrices are among the most important classes of structured matrices admitting efficient and stable eigenvalue algorithms. Notable examples include structure-preserving variants of the QR algorithm and their square-root-free counterparts. For unitary matrices, these algorithms were historically derived from recurrence relations for orthogonal polynomials on the unit circle. A matrix-based derivation using core chasing was later given in [3]; this family of algorithms is backward stable, but the underlying symmetries needed to cast out the square roots were not identified, and the question of structured backward stability was left open. In this paper we present a matrix-based derivation of the symmetric unitary core-chasing algorithm and use it to obtain a square-root-free variant. We give a matrix-based proof of backward stability and show that strict structured backward stability is impossible in general. Open-source Fortran code and numerical experiments demonstrate that the combination of exploiting the symmetry and removing the square roots reduces computation time without degrading accuracy.

1. Introduction.

The paper identifies a previously unnoticed symmetry in unitary core chasing, uses it to construct a square-root-free algorithm, and analyzes its floating-point stability. It also reports numerical experiments comparing the algorithms with LAPACK and prior unitary core chasing.

1. Introduction.

  • The core-chasing formulation contains a previously unnoticed symmetry, preserved in exact arithmetic by a symmetric turnover.The symmetry is generated by the first fusion in a single-shift step.Source: §1 · p. 2
  • A modified symmetric turnover is nearly exact for nearly-unitary data, where naïve symmetric-turnover formulas can be unstable.The modified formula is analogous to Gragg’s stabilized recurrence.Source: §1 · p. 2
  • The turnover symmetry yields a square-root-free unitary core-chasing algorithm that preserves near-unitarity without additional square roots.Source: §1 · p. 2
  • Both turnovers are forward stable for nearly-unitary floating-point input, but strict structured backward error is impossible in general.The forward error can be converted into a backward error result for the turnover.Source: §1 · p. 2
  • The experiments compare the symmetric and square-root-free algorithms with LAPACK and the unitary core-chasing algorithm in [3].Source: §1 · p. 2

2. The Unitary Core-Chasing Algorithm.

Unitary upper Hessenberg matrices can be represented by local core transformations and updated through a sequence of turnovers. A symmetry in the initial fusion is preserved through the chasing step, enabling an efficient symmetric turnover.

2.1. Core transformations.

  • A unitary upper Hessenberg matrix factors into core transformations, enabling O(n) storage and an O(n)-flop QR step.The factorization is the one used by core-chasing algorithms.Source: §2.1 · p. 3

2.2. Unitary core-chasing.

  • A core-chasing iteration starts with a shift-dependent transformation and uses turnovers and similarity transformations to restore upper Hessenberg form.The initial transformation is chosen from the first column of U − ρI.Source: §2.2 · pp. 3–4

2.3. Symmetry in the unitary core-chasing algorithm.

  • The initial fusion creates a symmetry that reduces turnover arithmetic and is maintained throughout the iteration.The symmetry follows from the special form of the initial transforming matrix.Source: §2.2 · p. 4

2.3.1. Symmetry in the turnover.

  • The symmetric turnover refactors the structured product while preserving the initial symmetry as the transformations move through the core-chasing step.The shift requires an extra parameter, which allows the symmetry to persist.Source: §2.3.1 · pp. 5, 7
  • The turnover formulas expose a structural resemblance between unitary core chasing and the QR algorithm for symmetric tridiagonal matrices.The paper uses a 3 × 3 symmetric tridiagonal matrix as an illustrative example.Source: §2.4 · p. 8

2.4. Connection to symmetric tridiagonal QR.

  • The symmetric tridiagonal QR iteration begins with a real shift and a similarity transformation that creates a bulge.Source: §2.4 · p. 8

1 TQ1 where Q1 is the following core transformation,

The symmetric tridiagonal QR step creates a bulge, then removes it with a core transformation to restore tridiagonal form. The exact-arithmetic formulas are not stable in floating point.

1 TQ1 where Q1 is the following core transformation,

  • A real-shift similarity transformation creates a bulge in the symmetric tridiagonal matrix, which a core transformation then removes.Source: §1 · p. 8
  • The formulas compute the next core transformation and the first row and column of the updated tridiagonal matrix in exact arithmetic.The formulas use normalized parameters for the transformation.Source: §1 · p. 8
  • The exact-arithmetic turnover formulas are unstable in floating point when unit-norm parameter relations drift or cancellation occurs.These formulas therefore do not permit a square-root-free symmetric turnover.Source: §3 · p. 9

3. Nearly-Unitary Matrices.

The section stabilizes the symmetric turnover by changing variables, then proves it remains nearly exact for nearly-unitary inputs. It also explains why an alternative formula for z is needed in that setting.

3. Nearly-Unitary Matrices.

  • The section’s changed-variable formulation stabilizes the symmetric turnover on nearly-unitary floating-point inputs and supports a square-root-free implementation.The square-root-free derivation is developed in the next section.Source: §3 · p. 9, §3.1 · p. 9, §4 · p. 15

3.1. A change of variables.

  • The original symmetric turnover can produce an O(1) error even when its input parameters are unitary to working precision.The instability arises in floating-point evaluation; the section then develops revised formulas to avoid it.Source: The corresponding computed output… · p. 10, §3.2 · p. 12
  • Writing c in polar form and introducing t and z avoids cancellation in key formulas and lets the update compute |bc| and bσ directly.The formula for z is evaluated using a branch chosen to avoid cancellation, following an idea from.Source: The corresponding computed output… · pp. 10–11

3.2. Nearly-unitary matrices.

  • For sufficiently small unitary defects, the symmetric turnover maps nearly-unitary inputs to nearly-unitary outputs, with nearby exactly unitary inputs producing close outputs.The theorem establishes this mixed-stability property for the turnover.Source: §3.2 · pp. 13–14

3.3. Two formulas for z.

  • The alternative formula for z is necessary because, for nearly-unitary inputs, z = 1 + t can be zero even when n > 0.The alternative formula guarantees that n ≠ 0 implies z ≠ 0; the two formulas agree for unitary core transformations.Source: §3.3 · p. 15

4. A Square-Root-Free Symmetric Turnover.

The symmetric turnover can be computed without square roots by working with squared variables. The section also describes renormalization to control unitary defects and notes limits of the square-root-free formulation.

4.1. Casting out the square roots.

  • The turnover can be computed without square roots by working directly with squared variables and combining γ and σ into ω = γσ^2.This removes two square-root computations, while the phase-tracking variable γ is absorbed into ω.Source: §4.1 · pp. 15–16
  • The square-root-free turnover computes the same outputs as the turnover with square roots for corresponding unsquared inputs.The equivalence is stated for inputs satisfying the corollary’s conditions.Source: §4.1 · pp. 16–17
  • The square-root-free formulation does not store the core transformations Qi explicitly, so it cannot compute eigenvectors.The passage describes the method in the context of eigenvalue computation.Source: §4.1 · p. 17

4.2. Square-root-free renormalization.

  • A degree-one Taylor correction renormalizes small defects without square roots, and is sufficient when defects are on the order of machine precision.The correction controls defect growth; floating-point representation still limits how small defects can become.Source: §4.2 · pp. 17–18
  • The turnovers are backward stable, but strict structured backward error is impossible in general for floating-point core transformations.The paper instead proves normwise backward stability for the input core transformations.Source: §5 · p. 18

5. Numerical Stability.

The analysis establishes forward and backward stability for the symmetric and square-root-free turnovers when inputs are nearly unitary. It also shows that strict structured backward stability is impossible in general.

5. Numerical Stability.

  • For floating-point inputs, naïve symmetric turnover formulas can be unstable, so the analysis uses a modified formula analogous to Gragg’s stabilized recurrence.The resulting symmetric turnover is proved nearly exact for nearly-unitary data.Source: §5.1 · p. 18

5.1. Forward stability.

  • With nearly-unitary floating-point inputs having O(µ) defects, both turnovers produce factorizations within O(µ) of the input product.The result covers the symmetric turnover in Algorithm 4.1 and the square-root-free turnover in Algorithm 4.2.Source: §5.1 · pp. 19–20

5.2. Backward stability.

  • For nearly-unitary inputs with O(µ) defects, Algorithms 4.1 and 4.2 are backward stable, producing factorizations exact for a nearby input matrix.The nearby input matrix need not itself be a core transformation, so the backward errors are not confined to its nontrivial block.Source: §5.2 · p. 21

5.3. Structured backward error.

  • Strict structured backward stability is impossible in general: nearly-unitary outputs need not be exact turnovers of nearby nearly-unitary input transformations.The obstruction applies even when the turnover is exact for inputs that are only nearly unitary.Source: §5.3 · pp. 21–22

6. Numerical Experiments.

The experiments compare the cost of basic operations and four eigenvalue methods, finding that exploiting structure, symmetry, and square-root removal reduces runtime while retaining comparable accuracy.

6.1. Square roots are expensive.

  • Division and square-root operations take more time on average than addition and multiplication in the sequential-operation experiment.The comparison executes 108 real, double-precision operations of each type on an Intel Core i7-4790 processor.Source: §6.1 · p. 22, Figure 6.1 · p. 23

6.2. Square-root-free unitary core-chasing.

  • The four methods show comparable accuracy, with computed roots within a modest multiple of machine precision of the true roots.The runtime plot contrasts O(n^3) dense QR with O(n^2) algorithms exploiting unitary upper Hessenberg structure.Source: §6.2 · pp. 22–23
  • At n = 4096, exploiting unitary upper Hessenberg structure reduces runtime versus LAPACK, while symmetry and square-root removal further reduce it versus structure-exploiting methods.Table 6.1 compares runtime and accuracy for the four methods at this size.Source: §6.2 · p. 23, Table 6.1 · p. 23
  • A change of variables stabilizes the symmetric turnover for nearly-unitary data and also enables a square-root-free version of the algorithm.The derivation uses a symmetry in the turnover operation to reduce arithmetic in the core-chasing step.Source: §7 · p. 24

7. Conclusions.

The paper derives a square-root-free unitary core-chasing algorithm and analyzes its floating-point stability. Experiments show comparable accuracy with reduced runtime, while strict structured backward error is impossible in general.

7. Conclusions.

  • The matrix-based derivation uses turnover symmetry to obtain a square-root-free algorithm, with stable turnovers for nearly-unitary inputs and a local backward error result.The resulting forward error can be converted into a local backward error for the turnover.Source: §7 · p. 25
  • Strict structured backward error is impossible in general.Source: §7 · p. 25
  • The experiments report accuracy comparable to LAPACK and the unitary core-chasing algorithm in [3], while reducing overall runtime.Source: §7 · p. 25
Loading 2610.12228v1…