blob: 08e17fa6d9057353bce9651b676be94e3246ce18 (
plain)
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
|
#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 {
// 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).
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.
// 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 = {});
} // namespace linalg
|