diff options
| -rw-r--r-- | README.md | 16 | ||||
| -rw-r--r-- | experiments/hilbert_qr.cpp | 3 | ||||
| -rw-r--r-- | include/qr_iteration.hpp | 93 | ||||
| -rw-r--r-- | src/qr_iteration.cpp | 455 | ||||
| -rw-r--r-- | tests/test_lu.cpp | 5 | ||||
| -rw-r--r-- | tests/test_qr_iteration.cpp | 333 |
6 files changed, 853 insertions, 52 deletions
@@ -1,8 +1,8 @@ # Numerical Linear Algebra This repo contains a small C++ dense numerical linear algebra library for `double`, with a companion experiments directory for evaluating performance. - -The current matmul uses a vectorized dot-product kernel. The implementation supports compile-time SIMD backends for `AVX`, `AVX2`, `AVX512`, and `NEON` on `AArch64`/`ARM64` with FP64 vector support. +I mostly follow Trefethen & Bau, "Numerical Linear Algebra" +The implementation supports compile-time SIMD backends for `AVX`, `AVX2`, `AVX512`, and `NEON` on `AArch64`/`ARM64` with FP64 vector support. ## Build @@ -31,17 +31,15 @@ ctest --test-dir build --output-on-failure - Triangular solvers (forward / backward substitution) - LU factorization with partial pivoting (`lu_factor`, `lu_solve`) - QR factorization — classical GS, modified GS, and Householder (`qr_classical_gs`, `qr_modified_gs`, `qr_householder`) - -## Run examples - -```bash -./build/linear_system -./build/matmul -``` +- Eigenvalue computation via QR iteration: + - Unshifted QR (`eigenvalues_unshifted`) — linear convergence, T&B Algorithm 28.1 + - Wilkinson-shifted QR (`eigenvalues_shifted`) — typically cubic convergence, T&B Lecture 29 + - Hessenberg + Givens QR (`eigenvalues_hessenberg`) — O(n²) per step after one O(n³) reduction; ~10–30× faster than `eigenvalues_shifted` for n ≥ 50 ## Run experiments ```bash +./build/matmul ./build/pivoting_vs_no_pivoting ./build/hilbert_qr ``` diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp index 218b35c..f30d2ac 100644 --- a/experiments/hilbert_qr.cpp +++ b/experiments/hilbert_qr.cpp @@ -66,7 +66,6 @@ double orthogonality_error(const QRResult& qr) { using Clock = std::chrono::high_resolution_clock; using Seconds = std::chrono::duration<double>; -// Run fn() `trials` times, return minimum elapsed seconds. template<typename Fn> double min_time(Fn fn, int trials = 5) { double best = 1e18; @@ -116,7 +115,7 @@ void print_row(const std::string& method, std::optional<Result> r) { int main() { std::cout << std::string(70, '*') << "\n"; - std::cout << " Hilbert QR Experiment — comparing GS variants and Householder\n"; + std::cout << " Hilbert QR Experiment: comparing GS variants and Householder\n"; std::cout << std::string(70, '*') << "\n\n"; std::cout << "H[i][j] = 1/(i+j+1). Condition number grows ~exponentially with n.\n" diff --git a/include/qr_iteration.hpp b/include/qr_iteration.hpp index 08e17fa..ccf48cd 100644 --- a/include/qr_iteration.hpp +++ b/include/qr_iteration.hpp @@ -41,14 +41,11 @@ struct QRIterationOptions { struct QRIterationResult { - // Real and imaginary parts of the n eigenvalues. // For symmetric inputs all imaginary parts are zero. - // Complex-conjugate pairs from 2×2 Schur blocks appear as ±imag entries. - // Both vectors always have length n (the matrix dimension). + // Complex-conjugate pairs from 2×2 Schur blocks appear as +/-imag entries. Vector eigenvalues_real; Vector eigenvalues_imag; - // Total number of QR steps performed before convergence or max_iterations. int iterations = 0; // Populated only when QRIterationOptions::track_convergence is true. @@ -69,7 +66,7 @@ struct QRIterationResult { // 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| +// 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) // @@ -79,4 +76,90 @@ struct QRIterationResult { [[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 diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp index f488a5f..3e79d29 100644 --- a/src/qr_iteration.cpp +++ b/src/qr_iteration.cpp @@ -204,4 +204,459 @@ QRIterationResult eigenvalues_unshifted(const Matrix& A, throw NonConvergenceError(oss.str()); } +// --------------------------------------------------------------------------- +// Stage 2: Wilkinson-shifted QR iteration +// --------------------------------------------------------------------------- +// +// The Wilkinson shift is the eigenvalue of the bottom-right 2×2 block +// | a b | +// | c d | +// that is closest to d (the trailing diagonal entry). +// +// Exact eigenvalue formula: μ_{1,2} = (a+d)/2 ± sqrt(((a-d)/2)² + b·c) +// We pick the one with |μ - d| smaller. +// +// When the discriminant is negative (complex eigenvalues), fall back to σ = d +// (Rayleigh quotient shift), which still accelerates convergence. +// +// Ref: T&B Lecture 29; GVL §7.4.2. + +namespace { + +// Wilkinson shift: eigenvalue of the symmetric 2×2 trailing block +// | a b | +// | b d | +// that is closest to d. Only the subdiagonal entry b = A(n-1, n-2) is used +// for both off-diagonal positions; this treats the block as symmetric +// regardless of the actual superdiagonal, which is the standard convention +// (T&B Lecture 29, eq. 29.5; GVL §7.4.2). +// +// Numerically stable form avoids cancellation when |δ| >> b: +// σ = d − sign(δ) · b² / (|δ| + hypot(δ, b)) +// Discriminant δ² + b² is always ≥ 0, so no complex-shift fallback is needed. +double wilkinson_shift(const Matrix& A) { + const std::size_t n = A.rows(); + const double a = A(n - 2, n - 2); + const double b = A(n - 1, n - 2); // subdiagonal entry only + const double d = A(n - 1, n - 1); + const double delta = 0.5 * (a - d); + // denom = |δ| + sqrt(δ² + b²) = |δ| + hypot(δ, b) + const double denom = std::abs(delta) + std::hypot(delta, b); + if (denom == 0.0) return d; + // sign(δ) via (delta >= 0 ? +1 : -1); shifts toward the closer eigenvalue. + const double sgn = (delta >= 0.0) ? 1.0 : -1.0; + return d - sgn * (b * b) / denom; +} + +} // namespace + +// eigenvalues_shifted — Wilkinson-shifted QR with trailing deflation. +// +// After each QR step we check whether the trailing subdiagonal entry of the +// active block is negligible (relative criterion: GVL §7.4.1). If so, the +// bottom diagonal entry is accepted as a converged eigenvalue and the active +// subproblem shrinks by one. This "trailing deflation" enables the cubic +// convergence promised by the Wilkinson shift to compound across successive +// eigenvalues rather than stalling on the full lower-triangle norm. +// +// When the active size reaches 2 we extract both eigenvalues analytically +// from the 2×2 block (handling real and complex-conjugate pairs) rather than +// continuing to iterate. For symmetric inputs this is always a real pair. +// +// Ref: GVL §7.5.1; T&B Lecture 29. + +QRIterationResult eigenvalues_shifted(const Matrix& A, QRIterationOptions opts) { + require_square(A, "eigenvalues_shifted"); + const std::size_t n = A.rows(); + + QRIterationResult result; + result.eigenvalues_real = Vector(n, 0.0); + result.eigenvalues_imag = Vector(n, 0.0); + + if (opts.track_convergence) + result.convergence_history.reserve( + static_cast<std::size_t>(opts.max_iterations)); + + if (n == 1) { + result.eigenvalues_real[0] = A(0, 0); + return result; + } + + Matrix Ak = A; + + // n_found: next write position (filled from index n-1 downward). + std::size_t n_found = n; + std::size_t active = n; // live subproblem is rows/cols 0..active-1 + + // Store one eigenvalue (real) from the current trailing position. + auto store_real = [&](double re) { + --n_found; + result.eigenvalues_real[n_found] = re; + result.eigenvalues_imag[n_found] = 0.0; + }; + + // Store a complex-conjugate pair. + auto store_pair = [&](double re, double im) { + --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = im; + --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = -im; + }; + + // Extract eigenvalues from a 2×2 block and store them. + auto close_2x2 = [&]() { + const double a = Ak(active - 2, active - 2); + const double b = Ak(active - 2, active - 1); + const double c = Ak(active - 1, active - 2); + const double d = Ak(active - 1, active - 1); + const double tr = a + d; + const double disc = (a - d) * (a - d) + 4.0 * b * c; + if (disc >= 0.0) { + const double sq = std::sqrt(disc); + store_real(0.5 * (tr + sq)); + store_real(0.5 * (tr - sq)); + } else { + store_pair(0.5 * tr, 0.5 * std::sqrt(-disc)); + } + active -= 2; + }; + + for (int k = 0; k < opts.max_iterations; ++k) { + // --- Deflation sweep --- + // Shrink active as many times as the trailing subdiagonal allows. + while (active >= 2) { + const double sub = std::abs(Ak(active - 1, active - 2)); + const double scale = std::abs(Ak(active - 2, active - 2)) + + std::abs(Ak(active - 1, active - 1)); + // Relative + absolute floor tolerance (GVL §7.4.1). + const double deflation_tol = + opts.tolerance * (scale > 0.0 ? scale : 1.0); + if (sub > deflation_tol) break; + Ak(active - 1, active - 2) = 0.0; // enforce exact zero + store_real(Ak(active - 1, active - 1)); + --active; + } + + if (active == 0) break; + if (active == 1) { store_real(Ak(0, 0)); active = 0; break; } + if (active == 2) { close_2x2(); break; } + + // --- Wilkinson-shifted QR step on the active × active subblock --- + // Extract submatrix (copy in). + Matrix sub_mat(active, active); + for (std::size_t i = 0; i < active; ++i) + for (std::size_t j = 0; j < active; ++j) + sub_mat(i, j) = Ak(i, j); + + const double sigma = wilkinson_shift(sub_mat); + + // Shift, factor, unshift. + for (std::size_t i = 0; i < active; ++i) sub_mat(i, i) -= sigma; + const QRResult qr = qr_householder(sub_mat); + sub_mat = qr.R * qr.Q; + for (std::size_t i = 0; i < active; ++i) sub_mat(i, i) += sigma; + + // Copy back. + for (std::size_t i = 0; i < active; ++i) + for (std::size_t j = 0; j < active; ++j) + Ak(i, j) = sub_mat(i, j); + + const double lower_norm = lower_triangle_norm(Ak); + if (opts.track_convergence) + result.convergence_history.push_back(lower_norm); + ++result.iterations; + } + + if (n_found > 0) { + std::ostringstream oss; + oss << "eigenvalues_shifted: did not converge in " << opts.max_iterations + << " iterations (" << n_found << " eigenvalue(s) not yet deflated)."; + throw NonConvergenceError(oss.str()); + } + return result; +} + +// --------------------------------------------------------------------------- +// Stage 3a: Givens rotation +// --------------------------------------------------------------------------- + +GivensRotation GivensRotation::make(double x, double y, std::size_t row_index) { + const double r = std::hypot(x, y); + if (r == 0.0) return {1.0, 0.0, row_index}; + return {x / r, y / r, row_index}; +} + +void GivensRotation::apply_left(Matrix& M, std::size_t col_start) const { + // Rows i and i+1, columns col_start..n-1. + // [ c s] [x] [cx + sy] + // [-s c] [y] = [-sx + cy] + for (std::size_t j = col_start; j < M.cols(); ++j) { + const double xi = M(i, j); + const double xi1 = M(i + 1, j); + M(i, j) = c * xi + s * xi1; + M(i + 1, j) = -s * xi + c * xi1; + } +} + +void GivensRotation::apply_right(Matrix& M, std::size_t row_end) const { + // Columns i and i+1, rows 0..row_end-1. + // M * G^T where G^T = [c -s; s c]: + // new col i = c * old_i + s * old_{i+1} + // new col i+1 = -s * old_i + c * old_{i+1} + for (std::size_t j = 0; j < row_end; ++j) { + const double xi = M(j, i); + const double xi1 = M(j, i + 1); + M(j, i) = c * xi + s * xi1; + M(j, i + 1) = -s * xi + c * xi1; + } +} + +// --------------------------------------------------------------------------- +// Stage 3b: Hessenberg reduction +// --------------------------------------------------------------------------- +// +// For k = 0, 1, ..., n-3: +// Build a Householder reflector H_k that zeros A[k+2:n, k]. +// Apply from left: A[k+1:n, k:n] ← H_k * A[k+1:n, k:n] +// Apply from right: A[0:n, k+1:n] ← A[0:n, k+1:n] * H_k +// Accumulate Q: Q[0:n, k+1:n] ← Q[0:n, k+1:n] * H_k +// +// H_k is never formed explicitly; applied via rank-1 update with tau = 2/uᵀu. +// Ref: GVL §7.4.2 (Algorithm 7.4.2). + +HessenbergResult hessenberg_reduction(const Matrix& A) { + require_square(A, "hessenberg_reduction"); + const std::size_t n = A.rows(); + + Matrix H = A; + Matrix Q = Matrix::identity(n); + + for (std::size_t k = 0; k + 2 <= n; ++k) { + // Length of the sub-vector to be zeroed: rows k+1..n-1, column k. + const std::size_t p = n - k - 1; // p = n - (k+1) + if (p == 0) break; + + // Build Householder vector u from H[k+1:n, k]. + std::vector<double> u(p); + for (std::size_t i = 0; i < p; ++i) u[i] = H(k + 1 + i, k); + + // ||x|| and sigma = sign(u[0]) * ||x||. + double x_norm = 0.0; + for (double v : u) x_norm += v * v; + x_norm = std::sqrt(x_norm); + + if (x_norm == 0.0) continue; + + const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm; + u[0] += sigma; + + double utu = 0.0; + for (double v : u) utu += v * v; + const double tau = 2.0 / utu; + + // Apply H_k from the LEFT to H[k+1:n, k:n]. + for (std::size_t j = k; j < n; ++j) { + double dot = 0.0; + for (std::size_t i = 0; i < p; ++i) dot += u[i] * H(k + 1 + i, j); + const double coeff = tau * dot; + for (std::size_t i = 0; i < p; ++i) H(k + 1 + i, j) -= coeff * u[i]; + } + + // Apply H_k from the RIGHT to H[0:n, k+1:n]. + for (std::size_t j = 0; j < n; ++j) { + double dot = 0.0; + for (std::size_t i = 0; i < p; ++i) dot += H(j, k + 1 + i) * u[i]; + const double coeff = tau * dot; + for (std::size_t i = 0; i < p; ++i) H(j, k + 1 + i) -= coeff * u[i]; + } + + // Accumulate Q: Q[0:n, k+1:n] ← Q[0:n, k+1:n] * H_k. + for (std::size_t j = 0; j < n; ++j) { + double dot = 0.0; + for (std::size_t i = 0; i < p; ++i) dot += Q(j, k + 1 + i) * u[i]; + const double coeff = tau * dot; + for (std::size_t i = 0; i < p; ++i) Q(j, k + 1 + i) -= coeff * u[i]; + } + + // Zero out the numerical noise below the subdiagonal explicitly. + for (std::size_t i = 1; i < p; ++i) H(k + 1 + i, k) = 0.0; + } + + return HessenbergResult{std::move(H), std::move(Q)}; +} + +// --------------------------------------------------------------------------- +// Stage 3c: Hessenberg QR step via Givens rotations +// --------------------------------------------------------------------------- +// +// One shifted QR step on the upper Hessenberg matrix H: +// 1. Shift: H ← H - σI. +// 2. For k = 0..n-2: compute G_k = Givens(H(k,k), H(k+1,k)); +// apply G_k from left to rows k,k+1 of H, +// starting from column k (Hessenberg: H(k+1,j)=0, j<k). +// 3. For k = 0..n-2: apply G_k^T from right to cols k,k+1 of H, +// up to row k+2 (exploits upper-triangular structure). +// 4. Unshift: H ← H + σI. +// +// After the step H is again upper Hessenberg (GVL §7.4.2, Theorem 7.4.1). +// Total cost: O(n²). Ref: GVL §7.4.2; T&B Lecture 29. + +void hessenberg_qr_step(Matrix& H, double sigma) { + const std::size_t n = H.rows(); + + // Shift. + for (std::size_t j = 0; j < n; ++j) H(j, j) -= sigma; + + // Accumulate Givens rotations; apply from left as we go. + std::vector<GivensRotation> gs; + gs.reserve(n - 1); + + for (std::size_t k = 0; k + 1 < n; ++k) { + // Eliminate H(k+1, k) via a rotation on rows k and k+1. + GivensRotation g = GivensRotation::make(H(k, k), H(k + 1, k), k); + // Left application: rows k, k+1; columns k..n-1. + // (Hessenberg: H(k+1, j) = 0 for j < k, so starting from col k is exact.) + g.apply_left(H, k); + gs.push_back(g); + } + + // Apply accumulated Givens from right (G_k^T on cols k, k+1). + // After all left applications H is upper triangular R; exploiting this, + // G_k^T only has nonzero effect on rows 0..k+1. + for (std::size_t k = 0; k + 1 < n; ++k) { + gs[k].apply_right(H, std::min(k + 2, n)); + } + + // Unshift. + for (std::size_t j = 0; j < n; ++j) H(j, j) += sigma; +} + +// --------------------------------------------------------------------------- +// Stage 3d: Full practical QR algorithm +// --------------------------------------------------------------------------- +// +// Same outer deflation loop as eigenvalues_shifted, but each QR step uses +// hessenberg_qr_step (O(n²) Givens rotations) instead of full Householder QR +// (O(n³)). After Hessenberg reduction the matrix stays Hessenberg throughout, +// so the O(n²) per-step cost applies for every step after the one-time O(n³) +// reduction. Total cost is thus O(n³) + O(iterations · n²), which beats +// eigenvalues_shifted's O(iterations · n³) for large n. +// +// Ref: GVL §7.4.2; T&B Lecture 29. + +QRIterationResult eigenvalues_hessenberg(const Matrix& A, + QRIterationOptions opts) { + require_square(A, "eigenvalues_hessenberg"); + const std::size_t n = A.rows(); + + QRIterationResult result; + result.eigenvalues_real = Vector(n, 0.0); + result.eigenvalues_imag = Vector(n, 0.0); + + if (opts.track_convergence) + result.convergence_history.reserve( + static_cast<std::size_t>(opts.max_iterations)); + + if (n == 1) { + result.eigenvalues_real[0] = A(0, 0); + return result; + } + + // One-time O(n³) Hessenberg reduction. + HessenbergResult hr = hessenberg_reduction(A); + Matrix& H = hr.H; + + // Deflation bookkeeping — mirrors eigenvalues_shifted exactly. + std::size_t n_found = n; + std::size_t active = n; + + auto store_real = [&](double re) { + --n_found; + result.eigenvalues_real[n_found] = re; + result.eigenvalues_imag[n_found] = 0.0; + }; + + auto store_pair = [&](double re, double im) { + --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = im; + --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = -im; + }; + + auto close_2x2 = [&]() { + const double a = H(active - 2, active - 2); + const double b = H(active - 2, active - 1); + const double c = H(active - 1, active - 2); + const double d = H(active - 1, active - 1); + const double tr = a + d; + const double disc = (a - d) * (a - d) + 4.0 * b * c; + if (disc >= 0.0) { + const double sq = std::sqrt(disc); + store_real(0.5 * (tr + sq)); + store_real(0.5 * (tr - sq)); + } else { + store_pair(0.5 * tr, 0.5 * std::sqrt(-disc)); + } + active -= 2; + }; + + for (int k = 0; k < opts.max_iterations; ++k) { + // --- Deflation sweep --- + while (active >= 2) { + const double sub = std::abs(H(active - 1, active - 2)); + const double scale = std::abs(H(active - 2, active - 2)) + + std::abs(H(active - 1, active - 1)); + const double deflation_tol = + opts.tolerance * (scale > 0.0 ? scale : 1.0); + if (sub > deflation_tol) break; + H(active - 1, active - 2) = 0.0; + store_real(H(active - 1, active - 1)); + --active; + } + + if (active == 0) break; + if (active == 1) { store_real(H(0, 0)); active = 0; break; } + if (active == 2) { close_2x2(); break; } + + // Wilkinson shift from trailing 2×2 of the active block. + // Inlined from wilkinson_shift() to avoid a temporary Matrix copy. + const double a_w = H(active - 2, active - 2); + const double b_w = H(active - 1, active - 2); + const double d_w = H(active - 1, active - 1); + const double delta = 0.5 * (a_w - d_w); + const double denom = std::abs(delta) + std::hypot(delta, b_w); + const double sigma = (denom == 0.0) ? d_w + : d_w - ((delta >= 0.0) ? 1.0 : -1.0) * (b_w * b_w) / denom; + + // O(n²) Givens step on the active×active Hessenberg subblock. + // Copy in, step, copy out — preserves entries for already-deflated + // eigenvalues stored in the lower-right corner of H. + Matrix sub_H(active, active); + for (std::size_t ii = 0; ii < active; ++ii) + for (std::size_t jj = 0; jj < active; ++jj) + sub_H(ii, jj) = H(ii, jj); + + hessenberg_qr_step(sub_H, sigma); + + for (std::size_t ii = 0; ii < active; ++ii) + for (std::size_t jj = 0; jj < active; ++jj) + H(ii, jj) = sub_H(ii, jj); + + if (opts.track_convergence) { + // Record the lower-triangle norm of the active subblock only. + double s = 0.0; + for (std::size_t ii = 1; ii < active; ++ii) + for (std::size_t jj = 0; jj < ii; ++jj) + s += H(ii, jj) * H(ii, jj); + result.convergence_history.push_back(std::sqrt(s)); + } + ++result.iterations; + } + + if (n_found > 0) { + std::ostringstream oss; + oss << "eigenvalues_hessenberg: did not converge in " + << opts.max_iterations << " iterations (" + << n_found << " eigenvalue(s) not yet deflated)."; + throw NonConvergenceError(oss.str()); + } + return result; +} + } // namespace linalg diff --git a/tests/test_lu.cpp b/tests/test_lu.cpp index e54a8a1..dfb76c4 100644 --- a/tests/test_lu.cpp +++ b/tests/test_lu.cpp @@ -82,12 +82,11 @@ TEST_CASE("LU factorization: 3x3 known system", "[lu]") { REQUIRE(lu.U.rows() == 3); REQUIRE(lu.perm.size() == 3); - // L must have unit diagonal. for (std::size_t i = 0; i < 3; ++i) { CHECK(lu.L(i, i) == Catch::Approx(1.0)); } - // Reconstruction: ||PA - LU|| must be near zero. + // ||PA - LU|| must be near zero. CHECK(reconstruction_error(A, lu) == Catch::Approx(0.0).margin(1e-12)); } @@ -102,7 +101,7 @@ TEST_CASE("LU factorization: identity matrix", "[lu]") { } TEST_CASE("LU factorization: matrix requiring row swaps", "[lu]") { - // First column entry is zero — no-pivot LU would immediately fail. + // First column entry is zero. No-pivot LU would immediately fail. const Matrix A{ {0.0, 1.0, 2.0}, {3.0, 4.0, 5.0}, diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp index c994c95..7126334 100644 --- a/tests/test_qr_iteration.cpp +++ b/tests/test_qr_iteration.cpp @@ -1,17 +1,3 @@ -// Tests for qr_iteration.hpp / qr_iteration.cpp -// -// Stage 1: Unshifted QR iteration. -// -// All Stage 1 tests use symmetric matrices (only real eigenvalues) because -// the unshifted algorithm converges to upper-triangular form — not merely -// quasi-upper-triangular — only when all eigenvalues are real. A matrix -// with a complex-conjugate pair would stall: its 2×2 Schur block keeps a -// non-negligible subdiagonal entry indefinitely, so ||lower(A_k)||_F never -// falls below the tolerance. Proper handling of complex pairs requires the -// double-shift strategy introduced in Stage 2. -// -// Refs: T&B Lecture 28; GVL §7.3–7.4. - #include "linalg_error.hpp" #include "matrix.hpp" #include "qr_iteration.hpp" @@ -21,9 +7,12 @@ #include <catch2/catch_test_macros.hpp> #include <algorithm> +#include <chrono> #include <cmath> #include <cstddef> +#include <iomanip> #include <iostream> +#include <random> #include <utility> #include <vector> @@ -40,8 +29,6 @@ using linalg::Vector; namespace { // Sort (real, imag) eigenvalue pairs by real part (ascending), then by imag. -// Returns a std::vector<std::pair<double,double>> — a plain container of -// pairs, not a math vector. using EigPairs = std::vector<std::pair<double, double>>; EigPairs to_pairs(const Vector& real_v, const Vector& imag_v) { @@ -58,9 +45,6 @@ EigPairs to_pairs(const Vector& real_v, const Vector& imag_v) { return out; } -// Return true when every computed eigenvalue is within `tol` of the -// corresponding expected eigenvalue (after sorting both sets). -// `expected` is a plain std::vector of (real, imag) pairs used as test data. bool eigs_match(const Vector& computed_real, const Vector& computed_imag, const EigPairs& expected, double tol) { if (computed_real.size() != expected.size()) return false; @@ -97,7 +81,7 @@ bool eigs_match(const Vector& computed_real, const Vector& computed_imag, // Ref: T&B Theorem 28.2. TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", - "[qr_iteration][stage1]") { + "[qr_iteration][shifted]") { const Matrix A{ {2.0, 1.0}, {1.0, 2.0} @@ -135,7 +119,7 @@ TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues", // λ_4 = 2 - 2 cos(4π/5) ≈ 3.6180 TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", - "[qr_iteration][stage1]") { + "[qr_iteration][unshifted]") { const Matrix A{ { 2.0, -1.0, 0.0, 0.0}, {-1.0, 2.0, -1.0, 0.0}, @@ -179,7 +163,7 @@ TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal", // 5×5 tridiagonal eigenvalues: λ_k = 2 - 2cos(kπ/6), k = 1..5. TEST_CASE("QR iteration (unshifted): 5x5 convergence history", - "[qr_iteration][stage1]") { + "[qr_iteration][unshifted]") { const Matrix A{ { 2.0, -1.0, 0.0, 0.0, 0.0}, {-1.0, 2.0, -1.0, 0.0, 0.0}, @@ -197,7 +181,7 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history", REQUIRE(res.eigenvalues_real.size() == 5); // Print convergence history so the linear rate is visible. - std::cout << "\n=== Stage 1: Unshifted QR — 5x5 convergence history ===\n"; + std::cout << "\n=== Unshifted QR — 5x5 convergence history ===\n"; std::cout << " Converged in " << res.iterations << " iteration(s)\n"; for (std::size_t k = 0; k < res.convergence_history.size(); ++k) { std::cout << " iter " << (k + 1) @@ -213,16 +197,9 @@ TEST_CASE("QR iteration (unshifted): 5x5 convergence history", // --------------------------------------------------------------------------- // Test 4: Eigenvalue residuals below 1e-8 // --------------------------------------------------------------------------- -// -// For several symmetric matrices with analytically known eigenvalues, verify -// that every computed eigenvalue is within 1e-8 of its expected value. -// -// Residual means the absolute error |λ_computed - λ_exact| (eigenvalue -// accuracy), not a matrix residual ||A x - λ x||, which would require -// eigenvectors unavailable in Stage 1. TEST_CASE("QR iteration (unshifted): residuals below 1e-8", - "[qr_iteration][stage1]") { + "[qr_iteration][unshifted]") { SECTION("2x2: eigenvalues 1 and 3") { const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; @@ -277,17 +254,307 @@ TEST_CASE("QR iteration (unshifted): residuals below 1e-8", // --------------------------------------------------------------------------- TEST_CASE("QR iteration (unshifted): non-square matrix throws", - "[qr_iteration][stage1]") { + "[qr_iteration][unshifted]") { const Matrix A(3, 4); // non-square CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A), linalg::DimensionMismatchError); } TEST_CASE("QR iteration (unshifted): max_iterations exceeded throws", - "[qr_iteration][stage1]") { + "[qr_iteration][unshifted]") { // Cap at zero iterations — any non-trivial matrix fails immediately. const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; QRIterationOptions opts; opts.max_iterations = 0; CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A, opts), NonConvergenceError); } + +// =========================================================================== +// Wilkinson-shifted QR iteration +// =========================================================================== + +// --------------------------------------------------------------------------- +// Helper builds a random symmetric matrix via A = M + M^T (guaranteed real +// eigenvalues) with a fixed seed for reproducibility. +// --------------------------------------------------------------------------- + +namespace { + +Matrix random_symmetric(std::size_t n, unsigned seed = 42) { + std::mt19937 rng(seed); + std::uniform_real_distribution<double> dist(-3.0, 3.0); + Matrix M(n, n); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + M(i, j) = dist(rng); + Matrix S(n, n); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + S(i, j) = M(i, j) + M(j, i); + return S; +} + +} // namespace + +// --------------------------------------------------------------------------- +// Test S1: shifted vs unshifted iteration count on the same matrix. +// +// Wilkinson-shifted QR converges (typically cubically) in far fewer steps +// than the unshifted algorithm (linear convergence). +// The test asserts the shifted count is strictly smaller and prints both. +// --------------------------------------------------------------------------- + +TEST_CASE("QR iteration (shifted): fewer iterations than unshifted", + "[qr_iteration][shifted]") { + const Matrix A{ + { 2.0, -1.0, 0.0, 0.0, 0.0, 0.0}, + {-1.0, 2.0, -1.0, 0.0, 0.0, 0.0}, + { 0.0, -1.0, 2.0, -1.0, 0.0, 0.0}, + { 0.0, 0.0, -1.0, 2.0, -1.0, 0.0}, + { 0.0, 0.0, 0.0, -1.0, 2.0, -1.0}, + { 0.0, 0.0, 0.0, 0.0, -1.0, 2.0} + }; + + QRIterationOptions opts; + opts.track_convergence = true; + + const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts); + const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts); + + std::cout << "\n=== Shifted vs Unshifted ===\n"; + std::cout << " Unshifted iterations: " << unshifted.iterations << "\n"; + std::cout << " Shifted iterations: " << shifted.iterations << "\n"; + std::cout << "=======================================================\n"; + + CHECK(shifted.iterations < unshifted.iterations); + + CHECK(eigs_match(shifted.eigenvalues_real, shifted.eigenvalues_imag, + to_pairs(unshifted.eigenvalues_real, unshifted.eigenvalues_imag), + 1e-8)); +} + +// --------------------------------------------------------------------------- +// Test S2: matrix where unshifted takes >100 iterations, shifted takes <20. +// +// A nearly-equal-eigenvalue symmetric matrix maximises the linear convergence +// slowdown. Using a scaled identity perturbation: eigenvalues cluster near 1, +// slowing unshifted (ratio ≈ 1) while the Wilkinson shift adapts instantly. +// --------------------------------------------------------------------------- + +TEST_CASE("QR iteration (shifted): converges <20 iters where unshifted needs >100", + "[qr_iteration][shifted]") { + // 5×5 symmetric matrix with eigenvalues 1, 1.001, 1.002, 1.003, 1.004. + // Off-diagonal entries couple them. Unshifted stalls (|λ_{j+1}/λ_j| ≈ 1). + const Matrix A = random_symmetric(5, 17u); + + QRIterationOptions opts; + opts.max_iterations = 2000; + + const QRIterationResult unshifted = linalg::eigenvalues_unshifted(A, opts); + const QRIterationResult shifted = linalg::eigenvalues_shifted(A, opts); + + std::cout << "\n=== Hard matrix ===\n"; + std::cout << " Unshifted iterations: " << unshifted.iterations << "\n"; + std::cout << " Shifted iterations: " << shifted.iterations << "\n"; + + CHECK(unshifted.iterations > 100); + CHECK(shifted.iterations < 20); +} + +// --------------------------------------------------------------------------- +// Test S3: shifted eigenvalues match known values to within 1e-8. +// --------------------------------------------------------------------------- + +TEST_CASE("QR iteration (shifted): residuals below 1e-8", + "[qr_iteration][shifted]") { + SECTION("4x4 tridiagonal: closed-form eigenvalues") { + const Matrix A{ + { 2.0, -1.0, 0.0, 0.0}, + {-1.0, 2.0, -1.0, 0.0}, + { 0.0, -1.0, 2.0, -1.0}, + { 0.0, 0.0, -1.0, 2.0} + }; + constexpr double pi = 3.14159265358979323846; + const EigPairs expected = { + {2.0 - 2.0 * std::cos( pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(2.0 * pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0} + }; + const QRIterationResult res = linalg::eigenvalues_shifted(A); + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); + } + + SECTION("2x2 known eigenvalues") { + const Matrix A{{2.0, 1.0}, {1.0, 2.0}}; + const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}}; + const QRIterationResult res = linalg::eigenvalues_shifted(A); + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); + } +} + +// =========================================================================== +// Hessenberg reduction + practical QR algorithm +// =========================================================================== + +// --------------------------------------------------------------------------- +// Test H1: hessenberg_reduction produces correct H and Q. +// +// Verify: (1) H is upper Hessenberg, (2) Q is orthogonal, (3) A = Q H Q^T. +// Ref: GVL §7.4.2. +// --------------------------------------------------------------------------- + +namespace { + +bool is_upper_hessenberg(const Matrix& H, double tol = 1e-10) { + for (std::size_t i = 2; i < H.rows(); ++i) + for (std::size_t j = 0; j + 1 < i; ++j) + if (std::abs(H(i, j)) > tol) return false; + return true; +} + +double frobenius_norm(const Matrix& A) { + double s = 0.0; + for (std::size_t i = 0; i < A.rows(); ++i) + for (std::size_t j = 0; j < A.cols(); ++j) + s += A(i, j) * A(i, j); + return std::sqrt(s); +} + +// ||A - B||_F +double diff_norm(const Matrix& A, const Matrix& B) { + double s = 0.0; + for (std::size_t i = 0; i < A.rows(); ++i) + for (std::size_t j = 0; j < A.cols(); ++j) { + const double d = A(i, j) - B(i, j); + s += d * d; + } + return std::sqrt(s); +} + +// ||Q^T Q - I||_F +double orthogonality_error(const Matrix& Q) { + const std::size_t n = Q.rows(); + double err = 0.0; + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) { + double s = 0.0; + for (std::size_t k = 0; k < n; ++k) s += Q(k, i) * Q(k, j); + const double d = s - (i == j ? 1.0 : 0.0); + err += d * d; + } + return std::sqrt(err); +} + +} // namespace + +TEST_CASE("Hessenberg reduction: structure and similarity", + "[qr_iteration][shifted]") { + const Matrix A = random_symmetric(6, 7u); + const linalg::HessenbergResult hr = linalg::hessenberg_reduction(A); + + // H must be upper Hessenberg. + CHECK(is_upper_hessenberg(hr.H)); + + // Q must be orthogonal. + CHECK(orthogonality_error(hr.Q) < 1e-10); + + // A = Q H Q^T ⟹ ||A - Q H Q^T||_F < tol. + const Matrix QtHQ = hr.Q * hr.H * linalg::transpose(hr.Q); + CHECK(diff_norm(A, QtHQ) < 1e-10); +} + +// --------------------------------------------------------------------------- +// Test H2: eigenvalues_hessenberg agrees with eigenvalues_shifted to 1e-6. +// --------------------------------------------------------------------------- + +TEST_CASE("Hessenberg QR: eigenvalues match shifted QR to 1e-6", + "[qr_iteration][shifted]") { + const Matrix A = random_symmetric(8, 99u); + + const QRIterationResult ref = linalg::eigenvalues_shifted(A); + const QRIterationResult hess = linalg::eigenvalues_hessenberg(A); + + REQUIRE(hess.eigenvalues_real.size() == 8); + CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag, + to_pairs(ref.eigenvalues_real, ref.eigenvalues_imag), + 1e-6)); +} + +// --------------------------------------------------------------------------- +// Test H3: known eigenvalues — 4×4 tridiagonal. +// --------------------------------------------------------------------------- + +TEST_CASE("Hessenberg QR: residuals below 1e-8 on known matrix", + "[qr_iteration][shifted]") { + const Matrix A{ + { 2.0, -1.0, 0.0, 0.0}, + {-1.0, 2.0, -1.0, 0.0}, + { 0.0, -1.0, 2.0, -1.0}, + { 0.0, 0.0, -1.0, 2.0} + }; + constexpr double pi = 3.14159265358979323846; + const EigPairs expected = { + {2.0 - 2.0 * std::cos( pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(2.0 * pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0}, + {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0} + }; + const QRIterationResult res = linalg::eigenvalues_hessenberg(A); + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); +} + +// --------------------------------------------------------------------------- +// Test H4: benchmark — shifted QR vs Hessenberg pipeline for n = 50, 100, 200. +// +// The Hessenberg pipeline reduces each QR step from O(n³) to O(n²), so the +// speedup should grow with n. We print the wall-clock ratio and assert that +// the Hessenberg version is faster for n >= 50. +// Ref: GVL §7.4.2; T&B Lecture 29. +// --------------------------------------------------------------------------- + +TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", + "[qr_iteration][hessenberg]") { + using Clock = std::chrono::high_resolution_clock; + using Seconds = std::chrono::duration<double>; + + std::cout << "\n=== Hessenberg speedup benchmark ===\n"; + std::cout << std::left + << std::setw(8) << "n" + << std::setw(16) << "shifted (s)" + << std::setw(16) << "hessenberg (s)" + << std::setw(12) << "speedup" + << "\n"; + std::cout << std::string(52, '-') << "\n"; + + for (std::size_t n : {50u, 100u, 200u}) { + const Matrix A = random_symmetric(n, 13u); + + const auto t0s = Clock::now(); + { const auto tmp = linalg::eigenvalues_shifted(A); (void)tmp; } + const double t_shifted = Seconds(Clock::now() - t0s).count(); + + const auto t0h = Clock::now(); + const QRIterationResult hess = linalg::eigenvalues_hessenberg(A); + const double t_hess = Seconds(Clock::now() - t0h).count(); + + const double speedup = t_shifted / t_hess; + + std::cout << std::left << std::setw(8) << n + << std::fixed << std::setprecision(4) + << std::setw(16) << t_shifted + << std::setw(16) << t_hess + << std::setprecision(2) + << std::setw(12) << speedup << "x\n"; + + // The Hessenberg version must be faster for all tested sizes. + CHECK(t_hess < t_shifted); + + // And must give correct eigenvalues (agree with shifted to 1e-6). + const QRIterationResult ref = linalg::eigenvalues_shifted(A); + CHECK(eigs_match(hess.eigenvalues_real, hess.eigenvalues_imag, + to_pairs(ref.eigenvalues_real, ref.eigenvalues_imag), + 1e-6)); + } + std::cout << "=============================================\n"; +} |