aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/qr.cpp134
-rw-r--r--src/qr_iteration.cpp258
2 files changed, 392 insertions, 0 deletions
diff --git a/src/qr.cpp b/src/qr.cpp
index 7e32ce3..a907bab 100644
--- a/src/qr.cpp
+++ b/src/qr.cpp
@@ -37,6 +37,24 @@ QRResult qr_modified_gs(const Matrix& A, double zero_tolerance = 1e-14);
// Throws DimensionMismatchError if rows < cols.
QRResult qr_householder(const Matrix& A);
+struct QRColPivResult {
+ Matrix Q;
+ Matrix R;
+ std::vector<std::size_t> perm;
+ std::size_t rank;
+};
+
+// Householder QR with column pivoting (rank-revealing).
+// At each step, the column with largest remaining norm is selected as pivot.
+// The numerical rank is determined by comparing diagonal entries of R to
+// rank_tolerance * |R(0,0)|.
+//
+// Returns Q (m x n), R (n x n upper triangular), perm (column permutation),
+// and rank (numerical rank estimate).
+//
+// Throws DimensionMismatchError if rows < cols.
+QRColPivResult qr_colpiv(const Matrix& A, double rank_tolerance = 1e-12);
+
} // namespace linalgebra
namespace {
@@ -194,4 +212,120 @@ QRResult qr_householder(const Matrix& A) {
return QRResult{std::move(Q), std::move(R)};
}
+QRColPivResult qr_colpiv(const Matrix& A, double rank_tolerance) {
+ require_tall(A, "qr_colpiv");
+ const std::size_t m = A.rows();
+ const std::size_t n = A.cols();
+
+ Matrix work = A;
+ Matrix Q_full = Matrix::identity(m);
+
+ std::vector<std::size_t> perm(n);
+ std::iota(perm.begin(), perm.end(), std::size_t{0});
+
+ // Precompute column norms squared.
+ std::vector<double> col_norms_sq(n);
+ for (std::size_t j = 0; j < n; ++j) {
+ double s = 0.0;
+ for (std::size_t i = 0; i < m; ++i) s += work(i, j) * work(i, j);
+ col_norms_sq[j] = s;
+ }
+
+ std::size_t rank = n;
+
+ for (std::size_t k = 0; k < n; ++k) {
+ // Find pivot: column with largest remaining norm.
+ std::size_t pivot = k;
+ double max_norm = col_norms_sq[k];
+ for (std::size_t j = k + 1; j < n; ++j) {
+ if (col_norms_sq[j] > max_norm) {
+ max_norm = col_norms_sq[j];
+ pivot = j;
+ }
+ }
+
+ // Swap columns k and pivot.
+ if (pivot != k) {
+ for (std::size_t i = 0; i < m; ++i) {
+ std::swap(work(i, k), work(i, pivot));
+ }
+ std::swap(col_norms_sq[k], col_norms_sq[pivot]);
+ std::swap(perm[k], perm[pivot]);
+ }
+
+ // Householder reflector for column k.
+ const std::size_t p = m - k;
+
+ std::vector<double> u(p);
+ for (std::size_t i = 0; i < p; ++i) u[i] = work(k + i, k);
+
+ const double x_norm = [&] {
+ double s = 0.0;
+ for (double v : u) s += v * v;
+ return std::sqrt(s);
+ }();
+
+ if (x_norm == 0.0) {
+ // Remaining columns are zero — rank determined.
+ rank = k;
+ break;
+ }
+
+ // Check rank: if this pivot norm is small relative to R(0,0).
+ if (k > 0) {
+ const double r00 = std::abs(work(0, 0));
+ if (x_norm <= rank_tolerance * r00) {
+ rank = k;
+ break;
+ }
+ }
+
+ const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm;
+ u[0] += sigma;
+
+ const double utu = [&] {
+ double s = 0.0;
+ for (double v : u) s += v * v;
+ return s;
+ }();
+ const double tau = 2.0 / utu;
+
+ // Apply reflector to work columns k..n-1.
+ for (std::size_t j = k; j < n; ++j) {
+ double d = 0.0;
+ for (std::size_t i = 0; i < p; ++i) d += u[i] * work(k + i, j);
+ const double coeff = tau * d;
+ for (std::size_t i = 0; i < p; ++i) work(k + i, j) -= coeff * u[i];
+ }
+
+ // Apply reflector to Q_full.
+ for (std::size_t j = 0; j < m; ++j) {
+ double d = 0.0;
+ for (std::size_t i = 0; i < p; ++i) d += u[i] * Q_full(k + i, j);
+ const double coeff = tau * d;
+ for (std::size_t i = 0; i < p; ++i) Q_full(k + i, j) -= coeff * u[i];
+ }
+
+ // Update column norms (downdate).
+ for (std::size_t j = k + 1; j < n; ++j) {
+ const double val = work(k, j);
+ col_norms_sq[j] -= val * val;
+ if (col_norms_sq[j] < 0.0) col_norms_sq[j] = 0.0;
+ }
+ }
+
+ // Extract Q (m x n) and R (n x n).
+ Matrix Q(m, n);
+ for (std::size_t i = 0; i < m; ++i)
+ for (std::size_t j = 0; j < n; ++j)
+ Q(i, j) = Q_full(j, i);
+
+ Matrix R(n, n);
+ for (std::size_t i = 0; i < n; ++i)
+ for (std::size_t j = 0; j < n; ++j)
+ R(i, j) = work(i, j);
+
+ return QRColPivResult{std::move(Q), std::move(R), std::move(perm), rank};
+}
+
} // namespace linalgebra
diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp
index 64fd25f..f8a3be5 100644
--- a/src/qr_iteration.cpp
+++ b/src/qr_iteration.cpp
@@ -65,6 +65,13 @@ void hessenberg_qr_step(Matrix& H, double sigma);
[[nodiscard]] QRIterationResult eigenvalues_hessenberg(const Matrix& A,
QRIterationOptions opts = {});
+// Francis double-shift QR — implicit bulge chasing on upper Hessenberg form.
+// Handles real matrices with complex conjugate eigenvalue pairs without
+// complex arithmetic. Uses robust deflation (subdiagonal + 2×2 block).
+// Reference: GVL §7.5, T&B Lecture 29.
+[[nodiscard]] QRIterationResult eigenvalues_francis(const Matrix& A,
+ QRIterationOptions opts = {});
+
} // namespace linalgebra
namespace {
@@ -502,4 +509,255 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, QRIterationOptions opt
return result;
}
+// ---------------------------------------------------------------------------
+// Francis double-shift QR with implicit bulge chasing
+// ---------------------------------------------------------------------------
+
+QRIterationResult eigenvalues_francis(const Matrix& A, QRIterationOptions opts) {
+ require_square(A, "eigenvalues_francis");
+ 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;
+ }
+
+ if (n == 2) {
+ const double a = A(0, 0), b = A(0, 1), c = A(1, 0), d = A(1, 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);
+ result.eigenvalues_real[0] = 0.5 * (tr + sq);
+ result.eigenvalues_real[1] = 0.5 * (tr - sq);
+ } else {
+ result.eigenvalues_real[0] = 0.5 * tr;
+ result.eigenvalues_imag[0] = 0.5 * std::sqrt(-disc);
+ result.eigenvalues_real[1] = 0.5 * tr;
+ result.eigenvalues_imag[1] = -0.5 * std::sqrt(-disc);
+ }
+ return result;
+ }
+
+ HessenbergResult hr = hessenberg_reduction(A);
+ Matrix& H = hr.H;
+
+ 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;
+ };
+
+ // Robust deflation: checks both subdiagonal magnitude and 2x2 block.
+ auto deflation_tol = [&](std::size_t i) -> double {
+ const double scale = std::abs(H(i - 1, i - 1)) + std::abs(H(i, i));
+ return opts.tolerance * (scale > 0.0 ? scale : 1.0);
+ };
+
+ 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;
+ };
+
+ // Find the start of the active unreduced block (split from the top).
+ auto find_block_start = [&]() -> std::size_t {
+ for (std::size_t i = active - 1; i >= 1; --i) {
+ if (std::abs(H(i, i - 1)) < deflation_tol(i)) {
+ H(i, i - 1) = 0.0;
+ return i;
+ }
+ }
+ return 0;
+ };
+
+ int exceptional_shift_count = 0;
+
+ for (int k = 0; k < opts.max_iterations; ++k) {
+ // Deflate converged eigenvalues from bottom.
+ while (active >= 2) {
+ if (std::abs(H(active - 1, active - 2)) < deflation_tol(active - 1)) {
+ H(active - 1, active - 2) = 0.0;
+ store_real(H(active - 1, active - 1));
+ --active;
+ } else {
+ break;
+ }
+ }
+
+ if (active == 0) break;
+ if (active == 1) { store_real(H(0, 0)); active = 0; break; }
+ if (active == 2) { close_2x2(); break; }
+
+ // Check for 2x2 block deflation (complex pair at bottom).
+ if (active >= 3 && std::abs(H(active - 2, active - 3)) < deflation_tol(active - 2)) {
+ H(active - 2, active - 3) = 0.0;
+ // The bottom 2x2 has converged.
+ 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;
+ exceptional_shift_count = 0;
+ continue;
+ }
+
+ std::size_t block_start = find_block_start();
+
+ // Compute Francis double shift from bottom 2x2 of active block.
+ const double a11 = H(active - 2, active - 2);
+ const double a12 = H(active - 2, active - 1);
+ const double a21 = H(active - 1, active - 2);
+ const double a22 = H(active - 1, active - 1);
+ double s = a11 + a22; // trace of bottom 2x2
+ double t = a11 * a22 - a12 * a21; // determinant of bottom 2x2
+
+ // Exceptional shift (Wilkinson's ad hoc) every 10 iterations to break stalls.
+ if (exceptional_shift_count > 0 && exceptional_shift_count % 10 == 0) {
+ const double w = std::abs(H(active - 1, active - 2))
+ + std::abs(H(block_start + 1, block_start));
+ s = 1.5 * w;
+ t = w * w;
+ }
+
+ // First column of M = H^2 - sH + tI (implicit).
+ const double h00 = H(block_start, block_start);
+ const double h01 = H(block_start, block_start + 1);
+ const double h10 = H(block_start + 1, block_start);
+ const double h11 = H(block_start + 1, block_start + 1);
+ const double h21 = (block_start + 2 < active) ? H(block_start + 2, block_start + 1) : 0.0;
+
+ double x = h00 * h00 + h01 * h10 - s * h00 + t;
+ double y = h10 * (h00 + h11 - s);
+ double z = h10 * h21;
+
+ // Chase the bulge through the Hessenberg matrix.
+ for (std::size_t i = block_start; i + 2 < active; ++i) {
+ // Determine Householder reflector P such that P * [x; y; z]^T = [*; 0; 0]^T.
+ const std::size_t p = (i + 3 <= active) ? 3 : 2;
+
+ double norm_v = std::sqrt(x * x + y * y + (p == 3 ? z * z : 0.0));
+ if (norm_v == 0.0) break;
+
+ const double sign = (x >= 0.0) ? 1.0 : -1.0;
+ double v0 = x + sign * norm_v;
+ double v1 = y;
+ double v2 = (p == 3) ? z : 0.0;
+
+ const double vdot = v0 * v0 + v1 * v1 + v2 * v2;
+ const double tau = 2.0 / vdot;
+
+ // Apply P from left to H rows [i, i+p-1], columns [max(i-1,0), active-1].
+ const std::size_t col_start = (i > 0) ? i - 1 : 0;
+ for (std::size_t j = col_start; j < active; ++j) {
+ double d = v0 * H(i, j) + v1 * H(i + 1, j);
+ if (p == 3) d += v2 * H(i + 2, j);
+ const double coeff = tau * d;
+ H(i, j) -= coeff * v0;
+ H(i + 1, j) -= coeff * v1;
+ if (p == 3) H(i + 2, j) -= coeff * v2;
+ }
+
+ // Apply P from right to H rows [0, min(i+p, active-1)], columns [i, i+p-1].
+ const std::size_t row_end = std::min(i + p + 1, active);
+ for (std::size_t j = 0; j < row_end; ++j) {
+ double d = v0 * H(j, i) + v1 * H(j, i + 1);
+ if (p == 3) d += v2 * H(j, i + 2);
+ const double coeff = tau * d;
+ H(j, i) -= coeff * v0;
+ H(j, i + 1) -= coeff * v1;
+ if (p == 3) H(j, i + 2) -= coeff * v2;
+ }
+
+ // Prepare for next bulge step.
+ if (i + 3 < active) {
+ x = H(i + 1, i);
+ y = H(i + 2, i);
+ z = (i + 3 < active) ? H(i + 3, i) : 0.0;
+ }
+ }
+
+ // Final 2x2 reflector to restore Hessenberg form at bottom.
+ {
+ const std::size_t i = active - 2;
+ const double xi = H(i, i - 1);
+ const double yi = H(i + 1, i - 1);
+ const double r = std::hypot(xi, yi);
+ if (r > 0.0) {
+ const double c = xi / r;
+ const double s_val = yi / r;
+ // Apply Givens from left.
+ for (std::size_t j = i - 1; j < active; ++j) {
+ const double t0 = H(i, j);
+ const double t1 = H(i + 1, j);
+ H(i, j) = c * t0 + s_val * t1;
+ H(i + 1, j) = -s_val * t0 + c * t1;
+ }
+ // Apply Givens from right.
+ for (std::size_t j = 0; j < std::min(i + 3, active); ++j) {
+ const double t0 = H(j, i);
+ const double t1 = H(j, i + 1);
+ H(j, i) = c * t0 + s_val * t1;
+ H(j, i + 1) = -s_val * t0 + c * t1;
+ }
+ }
+ }
+
+ ++exceptional_shift_count;
+
+ if (opts.track_convergence) {
+ double s_norm = 0.0;
+ for (std::size_t ii = 1; ii < active; ++ii)
+ s_norm += H(ii, ii - 1) * H(ii, ii - 1);
+ result.convergence_history.push_back(std::sqrt(s_norm));
+ }
+ ++result.iterations;
+ }
+
+ if (n_found > 0) {
+ std::ostringstream oss;
+ oss << "eigenvalues_francis: did not converge in "
+ << opts.max_iterations << " iterations ("
+ << n_found << " eigenvalue(s) not yet deflated).";
+ throw NonConvergenceError(oss.str());
+ }
+ return result;
+}
+
} // namespace linalgebra
generated by cgit v1.2.3 (git 2.51.0) at 2026-09-01 21:17:29 +0000