diff options
Diffstat (limited to 'include/qr_iteration.hpp')
| -rw-r--r-- | include/qr_iteration.hpp | 165 |
1 files changed, 0 insertions, 165 deletions
diff --git a/include/qr_iteration.hpp b/include/qr_iteration.hpp deleted file mode 100644 index ccf48cd..0000000 --- a/include/qr_iteration.hpp +++ /dev/null @@ -1,165 +0,0 @@ -#pragma once - -// QR iteration for eigenvalue computation -// -// Refs: -// Trefethen & Bau, "Numerical Linear Algebra" (T&B) -// Lecture 25 — Eigenvalue algorithms -// Lecture 26 — Schur factorisation -// Lecture 28 — The QR algorithm (unshifted) -// Lecture 29 — The QR algorithm with shifts -// Golub & Van Loan, "Matrix Computations" 4th ed. (GVL) -// §7.3 — The Unshifted QR Algorithm -// §7.4 — The Shifted QR Algorithm -// §7.4.2 — Wilkinson shift - -#include <vector> - -#include "linalg_error.hpp" -#include "matrix.hpp" -#include "vector.hpp" - -namespace linalg { - -// --------------------------------------------------------------------------- -// Options -// --------------------------------------------------------------------------- - -// All defaults are consistent with the recommendations in T&B Lecture 28. -struct QRIterationOptions { - // Convergence threshold. Iteration halts once the Frobenius norm of the - // strict lower triangle of A_k falls below this value. - // Ref: T&B §28; GVL §7.3. - double tolerance = 1e-10; - - int max_iterations = 1000; - - // When true, the Frobenius norm of the strict lower triangle is recorded - // after every QR step and returned in QRIterationResult::convergence_history. - bool track_convergence = false; -}; - - -struct QRIterationResult { - // For symmetric inputs all imaginary parts are zero. - // Complex-conjugate pairs from 2×2 Schur blocks appear as +/-imag entries. - Vector eigenvalues_real; - Vector eigenvalues_imag; - - int iterations = 0; - - // Populated only when QRIterationOptions::track_convergence is true. - // Entry k is ||lower(A_k)||_F after the k-th QR step. - // std::vector is used here because linalg::Vector has no push_back; - // convergence_history is a plain time-series container, not a math object. - std::vector<double> convergence_history; -}; - -// Algorithm (T&B Algorithm 28.1): -// -// A_0 = A -// for k = 1, 2, ...: -// factor A_{k-1} = Q_k R_k (Householder QR) -// set A_k = R_k Q_k (orthogonal similarity: preserves eigenvalues) -// -// The iterates A_k converge to the real Schur form of A: a quasi-upper- -// triangular matrix with 1×1 blocks (real eigenvalue) and 2×2 blocks -// (complex-conjugate pair) on the diagonal. -// -// Convergence rate: linear. Per-step factor ~= |lambda_{j+1} / lambda_j| -// for the off-diagonal entries linking eigenvalue clusters j and j+1. -// (T&B Lecture 28, Theorem 28.2) -// -// Throws DimensionMismatchError if A is not square. -// Throws NonConvergenceError if convergence is not achieved within -// opts.max_iterations steps. -[[nodiscard]] QRIterationResult eigenvalues_unshifted(const Matrix& A, - QRIterationOptions opts = {}); - -// --------------------------------------------------------------------------- -// Wilkinson-shifted QR iteration -// --------------------------------------------------------------------------- -// -// Same outer loop as Stage 1, but each step applies a shift σ chosen as -// the eigenvalue of the bottom-right 2×2 block of A_{k-1} that is closest -// to the (n,n) entry, then unshifts after the QR step: -// -// factor (A_{k-1} - σI) = Q_k R_k -// set A_k = R_k Q_k + σI -// -// The Wilkinson shift (T&B Lecture 29; GVL §7.4.2): -// Given the bottom-right 2×2 block | a b | -// | b c | -// δ = (a - c) / 2 -// σ = c - sign(δ) * b² / (|δ| + sqrt(δ² + b²)) -// equivalently: the eigenvalue of the block closer to c. -// -// Convergence rate: typically cubic near a simple eigenvalue. -// (T&B Lecture 29; GVL §7.5.1) -// -// Same exceptions as eigenvalues_unshifted. -[[nodiscard]] QRIterationResult eigenvalues_shifted(const Matrix& A, - QRIterationOptions opts = {}); - -// --------------------------------------------------------------------------- -// Stage 3: Hessenberg reduction algorithm -// --------------------------------------------------------------------------- - -// Givens rotation G acting on rows/columns i and i+1: -// -// G = | c s | chosen so that G * [x; y]^T = [r; 0]^T -// | -s c | with c = x/r, s = y/r, r = hypot(x, y) -// -// T&B Lecture 10 (Givens rotations). -struct GivensRotation { - double c; // cos(theta) - double s; // sin(theta) - std::size_t i; // first row/column index (second is i+1) - - // Construct the rotation that maps [x, y]^T → [hypot(x,y), 0]^T. - // Returns the identity (c=1, s=0) when x == y == 0. - [[nodiscard]] static GivensRotation make(double x, double y, - std::size_t row_index); - - // Apply G from the left to rows i and i+1 of M, columns [col_start, n). - // M[i:i+2, col_start:] ← G * M[i:i+2, col_start:] - void apply_left(Matrix& M, std::size_t col_start = 0) const; - - // Apply G^T from the right to columns i and i+1 of M, rows [0, row_end). - // M[0:row_end, i:i+2] ← M[0:row_end, i:i+2] * G^T - void apply_right(Matrix& M, std::size_t row_end) const; -}; - -// Result of reducing A to upper Hessenberg form. -// H is upper Hessenberg: H(i,j) = 0 for all i > j+1. -// Q is orthogonal and A = Q H Q^T. -// Ref: GVL Algorithm 7.4.2; T&B Lecture 26. -struct HessenbergResult { - Matrix H; // upper Hessenberg similarity of A - Matrix Q; // accumulated orthogonal transformation -}; - -// Reduce A to upper Hessenberg form via Householder reflectors applied -// from both sides. Costs O(10n³/3) flops; done once before QR iteration. -// Ref: GVL §7.4.2 (Algorithm 7.4.2). -// -// Throws DimensionMismatchError if A is not square. -[[nodiscard]] HessenbergResult hessenberg_reduction(const Matrix& A); - -// Apply one shifted QR step to an upper Hessenberg matrix H in-place, -// using n-1 Givens rotations. Costs O(n²) vs O(n³) for Householder QR. -// H remains upper Hessenberg after the step. -// Ref: GVL §7.4.2; T&B Lecture 29. -void hessenberg_qr_step(Matrix& H, double sigma); - -// Full practical QR algorithm: -// 1. Reduce A to Hessenberg H = Q^T A Q (O(n³), done once). -// 2. Run Wilkinson-shifted QR on H using Givens steps (O(n²) each). -// Substantially faster than eigenvalues_shifted for n ≥ 50. -// Ref: T&B Lecture 29. -// -// Same exceptions as eigenvalues_unshifted. -[[nodiscard]] QRIterationResult eigenvalues_hessenberg(const Matrix& A, - QRIterationOptions opts = {}); - -} // namespace linalg |