1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
|
#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
|