1 Francis QR algorithm basics

1.1 Eigenvalue problems and matrix setup

1.1.1 Standard forms (dense, sparse, and structured matrices)

The Francis QR algorithm targets eigenvalues of a square matrix \(A\), typically in real arithmetic. In practice, \(A\) may originate as a dense matrix, a sparse matrix, or a structured matrix (for example, banded or arising from discretizations). The algorithm’s core iterations are most efficient when the matrix has been transformed into a compact form—most notably Hessenberg form—so preprocessing often converts the original representation into an appropriate dense working format for the iteration.

Structured matrices can sometimes be exploited to reduce memory traffic and improve cache locality, even when the iteration itself proceeds on a dense banded profile. For sparse input, many implementations first apply dense-stage preprocessing to obtain a condensed form suitable for QR steps.

1.1.2 Similarity transformations and invariants

QR iterations are based on similarity transformations, which preserve the spectrum: if \(A_{k+1} = Q_k^{T} A_k Q_k\) with orthogonal \(Q_k\), then \(A_k\) and \(A_{k+1}\) have the same eigenvalues. More broadly, similarity preserves many matrix invariants, such as the characteristic polynomial and algebraic multiplicities of eigenvalues (over the complex numbers). Numerically, the orthogonality of the similarity transforms also helps limit the amplification of rounding errors during the process.

1.2 QR factorization in iterations

1.2.1 Q-R decomposition

A basic QR factorization expresses a matrix \(A_k\) as \[ A_k = Q_k R_k, \] where \(Q_k\) is orthogonal (or unitary in complex arithmetic) and \(R_k\) is upper triangular. The Francis QR algorithm modifies this factorization by incorporating shifts, but the underlying computation remains the repeated extraction of an orthogonal–triangular decomposition followed by a recombination step.

In dense implementations, the QR factorization is performed implicitly via orthogonal transformations (such as Householder reflectors) to avoid forming \(Q_k\) explicitly.

1.2.2 From A_k to A_{k+1} via similarity

Once \(A_k = Q_k R_k\) is available, the unshifted QR step sets \[ A_{k+1} = R_k Q_k. \] Since \(R_k Q_k\) is similar to \(Q_k^{T} (Q_k R_k) Q_k = Q_k^{T} A_k Q_k\), the step preserves eigenvalues. With shifts, the iteration is applied to a shifted matrix, and the recombination is arranged so that the net transformation remains a similarity of \(A_k\). This is the mechanism behind convergence: although eigenvalues are invariant, the matrix iterates toward a form where eigenvalues are read off.

1.3 Convergence and stopping criteria

1.3.1 Residual norms

A common practical measure is how small certain subdiagonal elements become. For Hessenberg (or quasi-triangular Schur) structures, convergence is often identified when the magnitude of entries below the first subdiagonal becomes negligible and the first subdiagonal entries indicate separation into nearly invariant subspaces. Some implementations track residual-like quantities related to eigenpairs; however, for large-scale QR iterations, monitoring subdiagonal norms is frequently the primary signal.

Residual norms can be interpreted as indicating how close the matrix is to an exact block triangular form, where eigenvalues of the blocks decouple. As the iteration progresses, the off-block coupling fades due to the shift strategy.

1.3.2 Deflation and extracting converged eigenvalues

Deflation refers to splitting off parts of the matrix that have effectively converged to an invariant block. When a subdiagonal entry \(a_{i+1,i}\) becomes sufficiently small, it is treated as zero, and the matrix is partitioned into top-left and bottom-right blocks. Each block then undergoes its own QR iteration.

This block separation both accelerates convergence and makes eigenvalue extraction straightforward: once a block reaches size 1, its diagonal entry is an eigenvalue; for a size-2 block in real arithmetic, the diagonal and superdiagonal entries form a quadratic whose roots give two complex-conjugate eigenvalues.

2 Hessenberg reduction and preparation

2.1 Why Hessenberg form

2.1.1 Computational efficiency

Hessenberg form is the sparsest dense form preserved by orthogonal similarity transforms in the sense that it has zeros below the first subdiagonal. For an \(n \times n\) matrix, a general Hessenberg matrix has \(O(n^2)\) nonzeros rather than \(O(n^3)\) operations per iteration that would arise from dense QR steps. With Hessenberg structure, QR iterations can be implemented with \(O(n^2)\) cost per sweep, rather than the \(O(n^3)\) cost typical of generic QR factorization on fully dense matrices.

Additionally, bulge-chasing operations introduced by shifted QR steps naturally preserve the near-banded structure and keep computational work localized.

2.1.2 Preservation under QR iterations

When the current iterate \(A_k\) is in Hessenberg form and the iteration uses orthogonal transformations that act locally near the subdiagonal, the next iterate remains in (or extremely close to) Hessenberg form. This invariance is crucial: it ensures that sparsity is not lost and that subsequent iterations remain efficient.

In implementations, small numerical fill-in from rounding is tolerated, then cleaned via the next transformations so the matrix effectively stays Hessenberg.

2.2 Householder reduction

2.2.1 Building the orthogonal similarity transforms

Householder reduction transforms a general matrix \(A\) into an upper Hessenberg matrix \(H\) by successive orthogonal similarities. Each step uses a Householder reflector to annihilate elements below the first subdiagonal in one column. After applying the reflector \(P_k\) on the left and the corresponding similarity update on the right, the transformation maintains orthogonality and improves numerical robustness.

If eigenvectors are requested, the algorithm may also accumulate these reduction transforms so that eigenvectors of the Hessenberg matrix can be mapped back to eigenvectors of the original matrix.

2.2.2 Complexity considerations

The reduction cost is typically \(O(n^3)\) for dense matrices, but it is a one-time preprocessing expense. For large \(n\), this cost is comparable to other dense eigenvalue approaches; however, QR iterations afterward are cheaper due to \(O(n^2)\) per iteration cost, and the number of iterations is moderated by shift strategies.

For structured or sparse inputs, the reduction stage can be adapted or hybridized in advanced software, but the conceptual role of Hessenberg preparation remains central.

2.3 Tridiagonal case (special simplification)

2.3.1 Relation to symmetric eigenproblems

For real symmetric matrices, orthogonal reduction naturally yields tridiagonal form rather than general Hessenberg form. The symmetric tridiagonal QR algorithm is closely related in spirit: its iteration steps exploit the additional structure to reduce work and improve stability properties, often leading to superior performance and easier convergence behavior.

In the symmetric case, the underlying transformations can be taken to be symmetric-aware, and eigenvalues are guaranteed real (though finite precision still requires careful handling).

2.3.2 Efficient implementations for structured matrices

Tridiagonal or banded matrices reduce the number of operations per iteration and simplify the bulge-chasing pathway. Many practical software packages implement specialized kernels for tridiagonal QR, using careful rotation sequences (often based on Givens rotations) to update only the few affected bands.

For near-tridiagonal matrices (e.g., resulting from particular discretizations), exploiting that near-band structure can reduce runtime while maintaining accuracy.

3 Shifted QR iterations

3.1 The role of shifts

3.1.1 Accelerating convergence

The unshifted QR method often converges slowly because it does not directly steer the iteration toward eigenvalues. Shifts modify the QR step so that the iteration acts as if it is trying to make a particular scalar close to an eigenvalue. When shifts are well chosen, the convergence to eigenvalues becomes much faster, especially for well-separated eigenvalues.

Shifts are also used to reduce the tendency of eigenvalues with large magnitude to dominate numerical behavior, thereby improving the balance of convergence across the spectrum.

3.1.2 Choice of shift parameters

The shift parameters are typically extracted from local information in the current iterate, frequently near the bottom-right corner where deflation is expected next. A common goal is to approximate an eigenvalue of the trailing principal submatrix. Different choices lead to different convergence profiles: some steps produce rapid progress for clusters, while others are more reliable but potentially slower.

Practical implementations often combine deterministic shift rules with safeguards based on monitoring whether the step produces sufficient reduction in subdiagonal entries.

3.2 Single-shift QR method

3.2.1 Implicit update idea

In the single-shift method, the QR factorization is effectively performed on \(A_k - \mu I\) for a chosen shift \(\mu\). The resulting update is arranged so that \[ A_{k+1} = Q^{T}(A_k - \mu I)Q + \mu I, \] where \(Q\) comes from the QR factorization of the shifted matrix. The key is that QR steps can be applied without explicitly forming \(Q\), using orthogonal transformations that act through the matrix.

The “implicit” aspect matters for performance: applying transformations directly to the Hessenberg matrix maintains sparsity and keeps arithmetic count moderate.

3.2.2 Practical numerical stability

Orthogonal transformations are numerically stable because they preserve norms up to rounding. However, shifts can introduce cancellation effects when \(A_k\) is close to \(\mu I\) in some entries. To mitigate this, computations are typically organized to avoid subtracting nearly equal numbers in ways that amplify relative error.

Implementations also use thresholds relative to matrix norms to decide when an entry is “small enough” for deflation, preventing premature splitting due to rounding noise.

3.3 Double-shift (bulge-chasing) approach

3.3.1 Motivation for faster progress

Double shifts are especially effective when the target eigenvalues form a complex-conjugate pair in real arithmetic, or when trailing submatrices produce better spectral approximations using a second-order shift model. The implicit double shift is a cornerstone technique associated with the “bulge-chasing” mechanism in shifted Hessenberg QR methods.

Compared with single shifts, double shifts can produce more consistent reduction of the bottom subdiagonal elements and yield faster deflation in many cases.

3.3.2 Bulge introduction and elimination

The implicit double shift introduces a “bulge” in the matrix profile: although the matrix starts in Hessenberg form, the shifted operations temporarily create nonzero elements outside the Hessenberg bandwidth. These elements are then chased downward through a sequence of local orthogonal rotations that restore the Hessenberg pattern after each chase step.

This procedure concentrates computation: each rotation updates only a small set of entries (a few adjacent rows and columns), so total work per iteration remains \(O(n^2)\) rather than growing to \(O(n^3)\).

3.3.3 Handling complex conjugate pairs

Over the reals, complex eigenvalues arise in conjugate pairs. The double-shift strategy can produce a real arithmetic iteration whose converged form is a quasi-triangular matrix (Schur form) containing a \(2\times 2\) block for each conjugate pair. The bulge-chasing rotations are designed so that this \(2\times 2\) block structure emerges naturally without leaving the real field.

As a result, eigenvalues can be extracted without switching to complex arithmetic during the iterative phase.

4 Implicit QR and the bulge-chasing mechanism

4.1 Implicitly shifted QR

4.1.1 Avoiding explicit formation of Q

Forming the orthogonal matrix \(Q\) at each step would be too costly and memory-intensive. Instead, the algorithm applies the same sequence of orthogonal transformations to the Hessenberg matrix directly. This preserves the Hessenberg profile and reduces the overhead of storing dense orthogonal factors.

The method effectively simulates the effect of the QR factorization and recombination while operating purely on the banded structure.

4.1.2 Using transformations through the matrix

Because the matrix is Hessenberg, each orthogonal transformation affects only a limited region. The implicit QR update can be seen as applying a chain of rotations or reflectors that “push” information along the subdiagonal. The transformations propagate local changes through the matrix while maintaining the invariant that previously cleared entries remain (numerically) small.

This organization also helps prevent spurious fill-in, as each step immediately restores the expected sparsity pattern.

4.2 Bulge-chasing steps

4.2.1 Local rotations and updates

The bulge-chasing phase uses a sequence of Givens rotations (or equivalent orthogonal transformations). Each rotation is selected to annihilate one element of the bulge, typically targeting a particular subdiagonal or near-subdiagonal entry. After applying the rotation, the bulge shifts one position further down, ready to be eliminated by the next rotation.

Only a few rows and columns are updated per rotation, leading to efficient computation and reduced numerical noise compared with broader operations.

4.2.2 Maintaining Hessenberg structure

The goal is that after each rotation update, the matrix returns to Hessenberg form (except for the moving bulge). When the bulge reaches the bottom of the active submatrix, the result is again an upper Hessenberg matrix for the next iteration.

This repeated maintenance of structure is central to the Francis QR algorithm’s efficiency: the cost per step depends on the matrix staying narrow-banded.

4.3 Computational cost per iteration

4.3.1 Operation counts

For an \(n\times n\) Hessenberg matrix, each QR sweep with bulge chasing costs on the order of \(O(n^2)\) floating-point operations. More precisely, constant factors depend on whether the algorithm uses single or double shifts, the use of rotations versus reflectors, and whether eigenvectors are being accumulated.

The cost is typically much smaller than the one-time Hessenberg reduction, which is \(O(n^3)\).

4.3.2 Practical performance considerations

Actual runtime depends on memory access patterns, the ability to use Level-2/Level-3 BLAS operations for dense subblocks, and the overhead of accumulating transformations for eigenvectors. For eigenvalues-only problems, implementations can omit much bookkeeping and run faster.

Deflation reduces the effective matrix size over time, shrinking the work required for later iterations and improving overall efficiency.

5 Extracting eigenvalues and deflation

5.1 Deflation strategy

5.1.1 Detecting negligible subdiagonal entries

Deflation relies on deciding when a subdiagonal element is “effectively zero.” Practical rules compare \(a_{i+1,i}\) against a tolerance derived from local or global norms and machine precision. This prevents misclassification when entries are small merely due to scaling rather than actual convergence.

Once an entry is judged negligible, the matrix is treated as block upper triangular, and the corresponding split is made explicit in the iteration logic.

5.1.2 Splitting the matrix into blocks

After deflation, the matrix is divided into smaller blocks whose eigenvalues can be computed independently. Each block proceeds through its own QR iterations, and the algorithm focuses computational effort on unresolved trailing blocks.

This block decomposition is not only an efficiency improvement; it also provides a robust path to handle clusters by isolating subspaces as soon as they decouple sufficiently.

5.2 Recovering eigenvectors (when requested)

5.2.1 Backward transformation framework

If eigenvectors of the original matrix are required, the algorithm must map eigenvectors from the final structured form back through all similarity transformations used during reduction and iteration. In the QR framework, this typically involves applying accumulated orthogonal transformations to the eigenvectors obtained from the quasi-triangular (Schur) form.

When deflation creates blocks, eigenvectors can often be computed progressively, aligning with the same splitting used for eigenvalues.

5.2.2 Accuracy and conditioning considerations

Eigenvector computation is more sensitive than eigenvalue extraction, particularly when eigenvalues are clustered. The conditioning of eigenvectors depends on the separation between eigenvalues and on the non-normality of the matrix. Even though orthogonal transformations are stable, the final step that solves for eigenvectors can magnify errors when eigenvalues are close.

To address this, implementations may use careful normalization, backward error checks, and additional steps based on the computed Schur form to improve the quality of returned eigenvectors.

5.3 Complex eigenvalues handling

5.3.1 Real arithmetic formulations

In real arithmetic QR iterations, complex eigenvalues appear as converged \(2\times2\) blocks. Extracting eigenvalues from such blocks avoids explicit complex iterations, allowing the algorithm to remain efficient on real-valued data structures.

The eigenvalues are determined by the characteristic polynomial of each \(2\times2\) block, yielding a pair of complex conjugates when the discriminant is negative.

5.3.2 2×2 blocks in the quasi-triangular Schur form

The output of shifted QR iterations is often described in terms of the (real) Schur decomposition: the matrix is transformed into a quasi-triangular form with \(1\times1\) blocks for real eigenvalues and \(2\times2\) blocks for complex-conjugate pairs. The double-shift mechanism naturally drives the iteration toward this structure.

Once the Schur form is obtained (implicitly or explicitly), eigenvalue extraction is systematic: diagonal entries for \(1\times1\) blocks and quadratic roots for \(2\times2\) blocks.

6.1 Symmetric vs nonsymmetric variants

6.1.1 Symmetric tridiagonal QR (Lanczos-style context)

For symmetric matrices, the algorithm simplifies substantially when reduced to tridiagonal form. Iterations exploit symmetry to preserve a strictly tridiagonal or near-tridiagonal structure and can use shift rules tailored to symmetric problems. In such settings, convergence behavior is typically favorable, and eigenvalues are real.

The symmetric QR approach is conceptually connected to other Krylov methods (often discussed in the broader context of Lanczos-type procedures) but differs in that QR iterations directly target the full spectral decomposition of the structured matrix.

6.1.2 Nonsymmetric Hessenberg QR generalization

For general nonsymmetric matrices, Hessenberg QR methods apply with shifts and deflation but require handling potentially complex eigenvalues. The resulting quasi-triangular Schur form may contain \(2\times2\) blocks even when the input is real.

The nonsymmetric case is more delicate in eigenvector accuracy because non-normality can cause large transient growth and sensitivity. Nonetheless, the orthogonal similarity framework remains the basis for numerical stability.

6.2 Connection to the Schur decomposition

6.2.1 QR iterations as a route to quasi-triangular form

QR iterations with shifts are commonly viewed as a process that transforms \(A\) toward real or complex Schur form. The quasi-triangular output provides a reliable representation of eigenvalues and, when requested, a platform for computing eigenvectors.

In this perspective, deflation and block splitting correspond to identifying invariant subspaces that have nearly converged to Schur blocks.

6.3 Alternative shift strategies

6.3.1 Wilkinson-type shifts

A widely used shift strategy is the Wilkinson shift, which selects \(\mu\) based on eigenvalue approximations of a small trailing submatrix (commonly \(2\times2\) in symmetric settings). Its purpose is to maximize the local convergence rate and reduce the number of iterations needed for deflation.

For nonsymmetric problems, analogous local rational approximations can be employed to guide shifts effectively without requiring full spectral information.

6.3.2 Heuristics used in implementations

Real implementations often blend theoretically motivated shift rules with safeguards. Examples include preventing stagnation when shifts produce insufficient reduction, switching between shift types, or adjusting tolerances in response to scaling. Heuristics may also address cases where the matrix has nearly equal eigenvalues (clusters), which can otherwise slow deflation.

These choices are typically tuned for robustness across diverse matrix types and conditioning regimes.

7 Numerical properties

7.1 Stability considerations

7.1.1 Orthogonality and rounding error control

Orthogonal similarity transforms limit growth in rounding error because they preserve norms. In the QR algorithm, each step uses orthogonal transformations that are norm-preserving (up to machine precision), which makes the iteration backward stable in many practical senses.

Bulge-chasing preserves this property by using localized orthogonal rotations rather than large dense updates, keeping numerical behavior controlled throughout the iteration.

7.1.2 Conditioning of eigenvalue computation

Eigenvalues can be relatively well-conditioned compared with eigenvectors, but conditioning depends on eigenvalue separation and on matrix non-normality. QR iterations still produce accurate eigenvalues in typical settings, yet clusters can lead to larger relative errors in eigenvector directions and, in some cases, in the separation of nearby eigenvalues.

Overall, the method’s use of orthogonal transformations helps ensure that the computed spectrum corresponds closely to the exact spectrum of a nearby matrix.

7.2 Accuracy of eigenvalues

7.2.1 Error estimation concepts

Many error estimates in numerical linear algebra relate computed eigenvalues to perturbations in \(A\). Since the QR iteration is driven by similarity transforms with orthogonal matrices, the effect resembles small structured perturbations over time. Practical checks often use residuals (for eigenpairs) or the smallness of off-diagonal entries (for Schur convergence).

For eigenvalues alone, convergence diagnostics via deflation and block structure provide an indirect but meaningful indicator of accuracy.

7.2.2 Effect of scaling and normalization

The scaling of input matrices can affect the numerical magnitude of subdiagonal elements and the choice of deflation thresholds. Proper scaling can improve the interpretability of tolerances and reduce the chance of underflow or overflow in intermediate computations.

Normalization of computed eigenvectors also influences reported errors, since relative accuracy metrics depend on the chosen scale.

7.3 Failure modes and mitigations

7.3.1 Slow convergence cases

Slow convergence can occur when eigenvalues are tightly clustered, when shifts do not align well with the local spectral structure, or when deflation tolerances are too strict relative to floating-point limitations. The algorithm mitigates this through better shift rules, adaptive tolerances, and fallback strategies that ensure progress continues.

Even with these measures, some pathological cases require more iterations, particularly for nonsymmetric matrices with complex dynamics.

7.3.2 Breakdown scenarios (and typical remedies)

Numerical “breakdown” in QR iteration typically refers to stagnation—no further reduction in key subdiagonal entries. Remedies include changing shift parameters, switching from double to single shifts (or vice versa), or performing alternative transformations that restart the bulge-chasing process.

Robust implementations also rely on careful handling of tiny pivots and on conservative decisions about when to declare deflation.

8 Practical implementation notes

8.1 Data structures and efficiency

8.1.1 Working with Hessenberg matrices

Implementations store only the relevant bands of the Hessenberg matrix or, in many cases, the full dense matrix while restricting updates to the Hessenberg profile. The former saves memory and bandwidth; the latter simplifies indexing but may waste work updating entries that remain structurally zero.

Efficient kernels assume the Hessenberg pattern when applying rotations and reflectors, reducing arithmetic and improving cache efficiency.

8.1.2 In-place updates

In-place updates reduce memory overhead by overwriting matrix entries during each QR sweep. Bulge-chasing is well suited to in-place operation because each rotation affects a local region that can be updated immediately.

Care must be taken with indexing order so that later computations do not rely on overwritten data that should still be available.

8.2 Interface with linear algebra libraries

8.2.1 Common API expectations (eigenvalues only vs eigenvectors)

Software libraries often provide options for “eigenvalues only” and “eigenvalues and eigenvectors.” Eigenvalues-only mode skips accumulation of reduction and iteration transformations, thereby lowering runtime and memory usage. Eigenvectors mode includes additional bookkeeping and postprocessing to map computed Schur vectors back to eigenvectors of the original matrix.

Return formats may differ depending on whether a Schur form is computed explicitly or whether eigenvectors are computed directly.

8.2.2 Output formats (real vs complex, Schur forms)

For real input, output can be presented as:

  • real eigenvalues plus complex conjugate pairs, or
  • a real Schur form containing \(1\times1\) and \(2\times2\) blocks, from which eigenvalues are derived.

Some interfaces can return a Schur form explicitly, which is useful for subsequent tasks requiring invariant subspaces, while others return only the eigenvalues (and optionally eigenvectors).

8.3 Example workflow

8.3.1 Reduction → QR iterations → deflation → result extraction

A typical end-to-end workflow begins by transforming \(A\) to Hessenberg form using orthogonal Householder reflections. Next, the algorithm performs shifted QR iterations on the Hessenberg matrix, repeatedly applying implicit single or double shifts and bulge-chasing rotations. As the iteration proceeds, deflation detects nearly zero subdiagonal entries and splits the matrix into smaller blocks, reducing the active problem size.

Finally, the algorithm extracts eigenvalues from the converged blocks (diagonal entries for real eigenvalues, \(2\times2\) blocks for complex pairs) and, when requested, computes eigenvectors through accumulated transformations back to the original basis.