diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 531b8aa..920e8d3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -41,23 +41,6 @@ target_link_libraries(vector-3d PRIVATE ) -# Matrix -add_library(matrix - STATIC - Matrix.cpp -) - -target_link_libraries(matrix - PUBLIC - vector-3d-intf - PRIVATE -) - -set_target_properties(matrix - PROPERTIES - LINKER_LANGUAGE CXX -) - # SVD add_library(svd STATIC @@ -73,4 +56,40 @@ target_link_libraries(svd set_target_properties(svd PROPERTIES LINKER_LANGUAGE CXX +) + +# QR (eigenvalues/eigenvectors via implicit shifted QR iteration) +add_library(qr + STATIC + QR.cpp +) + +target_link_libraries(qr + PUBLIC + vector-3d-intf + PRIVATE +) + +set_target_properties(qr + PROPERTIES + LINKER_LANGUAGE CXX +) + +# Matrix +add_library(matrix + STATIC + Matrix.cpp +) + +target_link_libraries(matrix + PUBLIC + vector-3d-intf + PRIVATE + svd + qr +) + +set_target_properties(matrix + PROPERTIES + LINKER_LANGUAGE CXX ) \ No newline at end of file diff --git a/src/Matrix.cpp b/src/Matrix.cpp index a4ea714..5e087be 100644 --- a/src/Matrix.cpp +++ b/src/Matrix.cpp @@ -5,6 +5,21 @@ #include "Matrix.hpp" #endif +// Forward-declare QR::EigenQR so the Matrix::EigenQR implementation below can +// call it even when Matrix.cpp is pulled in through QR.hpp's own include chain +// (QR.cpp -> QR.hpp -> Matrix.hpp -> Matrix.cpp), where the QR namespace has +// not been declared yet at this point. If we are not already inside that +// chain, pull in the full QR library so its template definition is available. +namespace QR { +template +void EigenQR(Matrix &matrixToDecompose, Matrix &eigenVectors, + Matrix &eigenValues, uint32_t maxIterations, + float tolerance); +} +#ifndef QR_H_ +#include "QR.hpp" +#endif + #ifdef MATRIX_H_ // since the .cpp file has to be included by the .hpp file this // will evaluate to true #include "Matrix.hpp" @@ -571,37 +586,10 @@ void Matrix::EigenQR(Matrix &eigenVectors, static_assert(rows > 1, "Matrix size must be > 1 for QR iteration"); static_assert(rows == columns, "Matrix size must be square for QR iteration"); - Matrix Ak = *this; // Copy original matrix - Matrix QQ{Matrix::Identity()}; - Matrix shift{0}; - - for (uint32_t iter = 0; iter < maxIterations; ++iter) { - Matrix Q, R; - - // // QR shift lets us "attack" the first diagonal to speed up the algorithm - // shift = Matrix::Identity() * Ak[rows - 1][rows - 1]; - (Ak - shift).QRDecomposition(Q, R); - Ak = R * Q + shift; - QQ = QQ * Q; - - // Check convergence: off-diagonal norm - float offDiagSum = 0.0f; - for (uint32_t row = 1; row < rows; row++) { - for (uint32_t column = 0; column < row; column++) { - offDiagSum += fabs(Ak[row][column]); - } - } - - if (offDiagSum < tolerance) { - break; - } - } - - // Diagonal elements are the eigenvalues - for (uint8_t i = 0; i < rows; i++) { - eigenValues[i][0] = Ak[i][i]; - } - eigenVectors = QQ; + // Delegate to the QR library: implicit shifted QR iteration with + // Wilkinson shift (see src/QR.hpp for the algorithm and conventions). + Matrix A = *this; // QR::EigenQR does not modify its input + QR::EigenQR(A, eigenVectors, eigenValues, maxIterations, tolerance); } #endif // MATRIX_H_ \ No newline at end of file diff --git a/src/Matrix.hpp b/src/Matrix.hpp index 1c28a08..5422898 100644 --- a/src/Matrix.hpp +++ b/src/Matrix.hpp @@ -5,7 +5,6 @@ #include #include -// TODO: Add a function to calculate eigenvalues/vectors // TODO: Add a function to compute RREF // TODO: Add a function for SVD decomposition // TODO: Add a function for LQ decomposition @@ -234,12 +233,19 @@ public: Matrix &R) const; /** - * @brief Uses QR decomposition to efficiently calculate the eigenvectors - * and values of this matrix - * @param eigenVectors a buffer that will contain the eigenvectors fo this - * matrix - * @param eigenValues a buffer that will contain the eigenValues fo this - * matrix + * @brief Calculates the eigenvectors and values of this matrix using the + * implicit shifted QR iteration (Wilkinson shift, Givens bulge chasing); + * see src/QR.hpp in the QR library for the full algorithm. + * @note For a matrix larger than 2x2 the matrix MUST be symmetric. + * A general (nonsymmetric) 2x2 is handled via the closed-form + * solution. + * @note The eigenvalues come out sorted DESCENDING (largest first); the + * eigenvector columns are swapped to match. Eigenvector signs are + * arbitrary. + * @param eigenVectors a buffer that will contain the eigenvectors of this + * matrix in its columns (column i pairs with eigenValues[i]) + * @param eigenValues a buffer that will contain the eigenvalues of this + * matrix, sorted descending * @param maxIterations the number of iterations to perform before giving * up on reaching the given tolerance * @param tolerance the level of accuracy to obtain before stopping. diff --git a/src/QR.cpp b/src/QR.cpp new file mode 100644 index 0000000..c542ce6 --- /dev/null +++ b/src/QR.cpp @@ -0,0 +1,416 @@ +// This #ifndef section makes clangd happy so that it can properly do type hints +// in this file +#ifndef QR_H_ +#define QR_H_ +#include "QR.hpp" +#endif + +#ifdef QR_H_ // since the .cpp file has to be included by the .hpp file this + // will evaluate to true +#include "QR.hpp" +#include +#include + +namespace QR { + +// ============================================================================ +// QR Building Block Implementations (fully templated, heap-free) +// ============================================================================ + +/** + * GivensRotation: R * (a, b)^T = (r, 0)^T with R = [[c, s], [-s, c]], + * r = +hypot(a, b), c = a/r, s = b/r. + */ +static void GivensRotation(float a, float b, float &c, float &s) { + float r = sqrtf(a * a + b * b); + if (r == 0.0f) { + c = 1.0f; + s = 0.0f; + return; + } + c = a / r; + s = b / r; +} + +/** + * ApplyRotationBothSides: A <- G A G^T (similarity transform) with + * G = [[c, s], [-s, c]] on the (i, i+1) block, i.e. G is the ZEROING + * rotation G*(x, y)^T = (r, 0)^T (the orientation used by the implicit QR + * chase: A = Q R with Q = G^T gives the next iterate R Q = G A G^T). + * With (c, s) = GivensRotation(A[i][i], A[i+1][i]) this zeroes + * A[i+1][i] after the LEFT multiplication; the right multiplication then + * chases the bulge along the superdiagonal (tridiagonal chase). + * + * A must be symmetric on entry; the result stays symmetric, so both + * triangles are written. + * + * Block updates (with a00 = A[i][i], a01 = A[i][i+1], a11 = A[i+1][i+1]): + * A[i][i] = c^2 a00 + 2 c s a01 + s^2 a11 + * A[i][i+1] = (c^2 - s^2) a01 + c s (a11 - a00) + * A[i+1][i+1] = s^2 a00 - 2 c s a01 + c^2 a11 + * Off-block updates (uniform for both sides, since the left factor G and + * the right factor G^T mix each side with the pattern (a, b) -> (c a + s b, + * -s a + c b) after transposition): + * for j not in {i, i+1}: + * A[i][j] = A[j][i] = c A[i][j] + s A[i+1][j] + * A[i+1][j] = A[j][i+1] = -s A[i][j] + c A[i+1][j] + */ +template +static void ApplyRotationBothSides(Matrix &A, uint8_t i, float c, + float s) { + float a00 = A.Get(i, i); + float a01 = A.Get(i, i + 1); + float a11 = A.Get(i + 1, i + 1); + float c2 = c * c; + float s2 = s * s; + float cs = c * s; + + A[i][i] = c2 * a00 + 2.0f * cs * a01 + s2 * a11; + A[i][i + 1] = (c2 - s2) * a01 + cs * (a11 - a00); + A[i + 1][i + 1] = s2 * a00 - 2.0f * cs * a01 + c2 * a11; + A[i + 1][i] = A[i][i + 1]; // keep both triangles in sync + + for (uint8_t j = 0; j < N; ++j) { + if (j == i || j == i + 1) + continue; + float x = A.Get(i, j); + float y = A.Get(i + 1, j); + A[i][j] = c * x + s * y; + A[j][i] = A[i][j]; + A[i + 1][j] = -s * x + c * y; + A[j][i + 1] = A[i + 1][j]; + } +} + +/** + * ApplyRotationToVectors: V <- V G^T with G = [[c, s], [-s, c]] on columns + * (i, i+1), applied to every row. G^T = [[c, -s], [s, c]], so + * V[r][i] <- c V[r][i] + s V[r][i+1] + * V[r][i+1] <- -s V[r][i] + c V[r][i+1] + * + * Convention pairing: if A evolves as A <- G A G^T (ApplyRotationBothSides + * with the SAME c, s), then V accumulates V <- V G^T. With V0 = I the + * invariant A0 = V A V^T is preserved at every step, so at convergence + * A0 = V D V^T and the columns of V are the eigenvectors. (Rationale: + * each chase step is A <- R Q with R = G A the upper-triangular factor and + * Q = G^T the orthogonal factor of A = Q R, so A = G^T A' G and the + * orthogonal factors multiply as G1^T G2^T ... in application order.) + */ +template +static void ApplyRotationToVectors(Matrix &V, uint8_t i, float c, + float s) { + for (uint8_t r = 0; r < N; ++r) { + float x = V.Get(r, i); + float y = V.Get(r, i + 1); + V[r][i] = c * x + s * y; + V[r][i + 1] = -s * x + c * y; + } +} + +/** + * WilkinsonShift: eigenvalue of [[a, b], [b, d]] closest to d. + * mu = (a+d)/2 - sign(a-d) * sqrt(((a-d)/2)^2 + b^2), sign(0) = +1. + */ +static float WilkinsonShift(float a, float b, float d) { + float delta = 0.5f * (a - d); + float spread = sqrtf(delta * delta + b * b); + return 0.5f * (a + d) - (delta >= 0.0f ? spread : -spread); +} + +/** + * Solve2x2Eigen: closed-form eigen-decomposition of the 2x2 block at + * (lo, lo+1). Works for symmetric blocks and for general 2x2 blocks with + * real eigenvalues (used by the N == 2 entry point). + * + * lambdaHi/lambdaLo come from the characteristic polynomial + * lambda^2 - trace*lambda + det = 0. + * The eigenvector for lambdaHi is v = (b, lambdaHi - a) (from the first + * row of (A - lambda*I)v = 0), normalized to unit length. If b == 0 the + * block is triangular and the eigenvectors are coordinate vectors: + * e1 for the larger of {a, d}, e2 for the other. + */ +template +static void Solve2x2Eigen(const Matrix &A, uint8_t lo, float &lambdaHi, + float &lambdaLo, float &c, float &s) { + float a = A.Get(lo, lo); + float b = A.Get(lo, lo + 1); + float e = A.Get(lo + 1, lo); + float d = A.Get(lo + 1, lo + 1); + + float trace = a + d; + float det = a * d - b * e; + float disc = trace * trace - 4.0f * det; + if (disc < 0.0f) + disc = 0.0f; // round-off clamp: real 2x2 blocks have disc >= 0 + float sqrtDisc = sqrtf(disc); + lambdaHi = 0.5f * (trace + sqrtDisc); + lambdaLo = 0.5f * (trace - sqrtDisc); + + if (b != 0.0f) { + float v1 = lambdaHi - a; + float n = sqrtf(b * b + v1 * v1); + c = b / n; + s = v1 / n; + } else if (a >= d) { + c = 1.0f; // e1 is the eigenvector of a = lambdaHi + s = 0.0f; + } else { + c = 0.0f; // e2 is the eigenvector of d = lambdaHi + s = 1.0f; + } +} + +/** + * Deflate: zero subdiagonal entries i in [lo, hi) whose magnitude is at or + * below tolerance * (|A[i][i]| + |A[i+1][i+1]|). + */ +template +static void Deflate(Matrix &A, uint8_t lo, uint8_t hi, float tolerance) { + for (uint8_t i = lo; i < hi; ++i) { + float t = A.Get(i + 1, i); + float scale = fabsf(A.Get(i, i)) + fabsf(A.Get(i + 1, i + 1)); + if (fabsf(t) <= tolerance * scale) { + A[i + 1][i] = 0.0f; + A[i][i + 1] = 0.0f; + } + } +} + +// ============================================================================ +// QR::EigenQR driver (implicit Wilkinson-shifted QR, bulge chasing) +// ============================================================================ + +/** + * Tridiagonalize: Givens tridiagonalization (Golub & Van Loan 8.3.1). + * + * For column k = 0..N-3 the entries A[k+2..N-1, k] are eliminated by + * rotations on (i, i+1) applied BOTTOM-UP, i = N-2 down to k+1, each + * formed from the CURRENT (already-updated) pair (A[i][k], A[i+1][k]). + * Bottom-up is essential: a top-down pass zeros A[i+1][k] with a rotation + * that would later be undone when the next rotation (i+1, i+2) is formed + * from an entry below, reviving A[i][k]. Each bottom-up rotation zeros the + * bottom of the remaining nonzero pair and the entries below stay zero + * (they are not mixed again, only rows i-1/i are mixed next). + * + * Already-tridiagonalized leading columns j < k are untouched: the mixed + * rows are both >= k+1 > j+1, so A[i][j] and A[i+1][j] are both zero there. + * The rotation on (i, i+1) also keeps column k+1..k+2 structure intact and + * does not destroy earlier columns, so after column k is done the leading + * (k+1)x(k+1) block is tridiagonal forever. + * + * On return: A is symmetric tridiagonal and A_orig = U A U^T (U = product + * of every rotation applied, in application order, as U <- U G^T). + */ +template +static void Tridiagonalize(Matrix &A, Matrix &U) { + U = Matrix{0}; + for (uint8_t i = 0; i < N; ++i) { + U[i][i] = 1.0f; + } + float c = 0.0f, s = 0.0f; + for (uint8_t k = 0; k + 2 < N; ++k) { + for (int i = (int)N - 2; i >= (int)k + 1; --i) { + GivensRotation(A.Get(i, k), A.Get(i + 1, k), c, s); + ApplyRotationBothSides(A, (uint8_t)i, c, s); + ApplyRotationToVectors(U, (uint8_t)i, c, s); + } + } +} + +/** + * See QR.hpp for the full contract. Implementation sketch: + * + * Phase 0 (N >= 3): Tridiagonalize(A, U) // A_orig = U A U^T + * V = I. + * while (hi > 0): + * Deflate(A, 0, hi, tol); peel exact-zero trailing subdiagonals (hi--) + * lo = top of the trailing unreduced block (scan down, stop at first + * exact zero subdiagonal) + * if lo == hi - 1: closed-form 2x2 eigen-solve; fold Vblock into V + * else: one implicit Wilkinson-shifted QR step: + * mu = WilkinsonShift(A[hi-1][hi-1], A[hi][hi-1], A[hi][hi]) + * A[lo..hi diagonal] -= mu // whole block! + * G1 = Givens(A[lo][lo], A[lo+1][lo]) + * for i = lo..hi-1: + * (i > lo: Gi = Givens(A[i][i], A[i+1][i])) + * ApplyRotationBothSides(A, i, Gi) // A <- Gi A Gi^T + * ApplyRotationToVectors(V, i, Gi) // V <- V Gi^T + * A[lo..hi diagonal] += mu + * eigenvalues = diag(A), sorted descending with matching V column swaps. + * eigenvectors = U * V. + * + * Invariant maintained for N >= 3 (symmetric input): A is symmetric + * tridiagonal (up to deflated zeros and ~1e-7 float roundoff in the + * off-tridiagonal corners) at the top of every loop iteration, and + * A_orig = U A U^T = (U V) A (U V)^T throughout (V = product of every + * rotation applied so far, in application order, as V <- V Gi^T). At + * convergence A = V D V^T and therefore A_orig = (U V) D (U V)^T. + * + * Orientation note: each chase rotation Gi is the ZEROING rotation + * (Gi * (x, y)^T = (r, 0)^T). The step A <- Gi A Gi^T equals R Q with + * R = Gi A upper-triangular (on the block) and Q = Gi^T -- i.e. it IS the + * standard QR update Q(A - mu I)Q^T with Q the orthogonal QR factor. The + * eigenvector accumulator therefore collects the Q factors: V <- V Gi^T. + */ +template +void EigenQR(Matrix &matrixToDecompose, Matrix &eigenVectors, + Matrix &eigenValues, uint32_t maxIterations, float tolerance) { + static_assert(N >= 2, "QR::EigenQR requires N >= 2 (N = 1 is trivial)"); + + Matrix A = matrixToDecompose; // input is not modified + Matrix V{0}; + // NB: Matrix::Identity() is a static factory that returns by value; a + // bare call would be a no-op. Set the diagonal explicitly. + for (uint8_t i = 0; i < N; ++i) { + V[i][i] = 1.0f; + } + + // ------------------------------------------------------------------ + // N == 2: closed-form solution (works for nonsymmetric input too) + // ------------------------------------------------------------------ + if (N == 2) { + float l1 = 0.0f, l2 = 0.0f, c = 0.0f, s = 0.0f; + Solve2x2Eigen(A, 0, l1, l2, c, s); + // V = I * Vblock = [[c, -s], [s, c]] + V[0][0] = c; + V[0][1] = -s; + V[1][0] = s; + V[1][1] = c; + eigenValues[0][0] = l1; + eigenValues[1][0] = l2; + for (uint8_t r = 0; r < N; ++r) + for (uint8_t col = 0; col < N; ++col) + eigenVectors[r][col] = V.Get(r, col); + return; + } + + // ------------------------------------------------------------------ + // N >= 3: implicit shifted QR iteration (symmetric input required) + // ------------------------------------------------------------------ + // Phase 0: general symmetric -> symmetric tridiagonal. The implicit + // QR bulge chase only preserves a tridiagonal structure, so the input + // must be reduced first: A_orig = U A U^T with A tridiagonal. + Matrix U{}; + Tridiagonalize(A, U); + + uint32_t iter = 0; + uint8_t hi = N - 1; + while (hi > 0) { + Deflate(A, 0, hi, tolerance); + // Peel trailing rows whose subdiagonal is exactly zero (deflated or + // already solved). Must be re-done every iteration: a peel is only + // meaningful once the subdiagonal beneath it has converged. + while (hi > 0 && A.Get(hi, hi - 1) == 0.0f) { + --hi; + } + if (hi == 0) { + break; // fully diagonal (within tolerance) + } + + // Find the top of the trailing unreduced block: scan down from hi-1 + // and stop at the first exact zero subdiagonal. A[hi][hi-1] != 0 here + // (just peeled), so lo < hi. + uint8_t lo = hi; + for (int i = (int)hi - 1; i >= 0; --i) { + if (A.Get(i + 1, i) == 0.0f) { + break; + } + lo = (uint8_t)i; + } + + if (lo + 1 == hi) { + // Trailing unreduced block is 2x2: solve in closed form. + float l1 = 0.0f, l2 = 0.0f, c = 0.0f, s = 0.0f; + Solve2x2Eigen(A, lo, l1, l2, c, s); + A[lo][lo] = l1; + A[lo + 1][lo + 1] = l2; + A[lo][lo + 1] = 0.0f; + A[lo + 1][lo] = 0.0f; + // Fold Vblock = [[c, -s], [s, c]] into V: V <- V * Vblock on + // columns (lo, lo+1). NOTE the sign convention differs from + // ApplyRotationToVectors (which applies [[c, s], [-s, c]]): + // here column 0 of Vblock is (c, s)^T, column 1 is (-s, c)^T. + for (uint8_t r = 0; r < N; ++r) { + float x = V.Get(r, lo); + float y = V.Get(r, lo + 1); + V[r][lo] = c * x + s * y; + V[r][lo + 1] = -s * x + c * y; + } + if (lo == 0) { + break; // block reached the top: matrix is fully solved + } + hi = (uint8_t)(lo - 1); + continue; + } + + // One implicit Wilkinson-shifted QR step on block [lo, hi]. + float mu = WilkinsonShift(A.Get(hi - 1, hi - 1), A.Get(hi, hi - 1), + A.Get(hi, hi)); + + // The shift applies to the ENTIRE active block: bulge chasing + // triangularizes (A - mu*I), and the first Givens rotation is formed + // from (A[lo][lo] - mu, A[lo+1][lo]). + for (uint8_t i = lo; i <= hi; ++i) { + A[i][i] -= mu; + } + + float c = 0.0f, s = 0.0f; + for (uint8_t i = lo; i < hi; ++i) { + if (i == lo) { + GivensRotation(A.Get(lo, lo), A.Get(lo + 1, lo), c, s); + } else { + GivensRotation(A.Get(i, i), A.Get(i + 1, i), c, s); + } + ApplyRotationBothSides(A, i, c, s); + ApplyRotationToVectors(V, i, c, s); + } + + for (uint8_t i = lo; i <= hi; ++i) { + A[i][i] += mu; + } + + if (++iter >= maxIterations) { + // Best-effort: fall through with the partially diagonalized A. + break; + } + } + + // ------------------------------------------------------------------ + // Collect eigenvalues and sort DESCENDING (swap eigenvectors to match) + // ------------------------------------------------------------------ + for (uint8_t i = 0; i < N; ++i) { + eigenValues[i][0] = A.Get(i, i); + } + for (uint8_t i = 0; i < N - 1; ++i) { + uint8_t k = i; + for (uint8_t j = i + 1; j < N; ++j) { + if (eigenValues.Get(j, 0) > eigenValues.Get(k, 0)) { + k = j; + } + } + if (k != i) { + float t = eigenValues[i][0]; + eigenValues[i][0] = eigenValues[k][0]; + eigenValues[k][0] = t; + for (uint8_t r = 0; r < N; ++r) { + float x = V.Get(r, i); + V[r][i] = V.Get(r, k); + V[r][k] = x; + } + } + } + + // True eigenvectors of the original matrix: U * V. Reuse the A buffer + // (its diagonal has already been collected into eigenValues). + U.Mult(V, A); + + for (uint8_t r = 0; r < N; ++r) { + for (uint8_t col = 0; col < N; ++col) { + eigenVectors[r][col] = A.Get(r, col); + } + } +} + +} // namespace QR + +#endif // QR_H_ diff --git a/src/QR.hpp b/src/QR.hpp new file mode 100644 index 0000000..95c4f54 --- /dev/null +++ b/src/QR.hpp @@ -0,0 +1,224 @@ +#pragma once +#include "Matrix.hpp" + +/** + * @brief Library that uses Matrix.hpp and computes the eigenvalues and + * eigenvectors of a square matrix with the implicit shifted QR iteration + * (Wilkinson shift, Givens bulge chasing). + * + * @note Fully templated: QR::EigenQR works for ANY Matrix with N in + * 2..255 (the uint8_t range of Matrix). There is no 5x5 limit. + * + * @note N >= 3: the input matrix MUST be symmetric (A[i][j] == A[j][i]). + * The implicit QR bulge chase maintains a symmetric tridiagonal + * structure, which only exists for symmetric input. N = 2 handles + * a general (nonsymmetric) 2x2 via the closed-form solution, so + * nonsymmetric 2x2 inputs also work. + * + * @note The input matrix is NOT modified (the iteration runs on a local + * copy), mirroring the SVD::SVD convention. + * + * @note EMBEDDED CONSTRAINT -- no heap. All working storage is stack + * allocated as templated Matrix buffers. Peak stack usage per + * call is 3 * N^2 floats (A working copy + U and V accumulators) = + * 12 * N^2 bytes: + * N = 5 -> ~0.3 KB + * N = 10 -> ~1.2 KB + * N = 20 -> ~4.8 KB + * N = 50 -> ~30 KB + * N = 100 -> ~120 KB + * N = 255 -> ~783 KB + * Instantiate only the sizes that fit your call-stack budget. + * + * @note Conventions: + * - Eigenvalues come out sorted DESCENDING (largest first); the + * eigenvector columns are swapped to match. + * - Eigenvector signs are arbitrary (v and -v are both valid); + * tests must be sign-invariant. + * - Wilkinson shift: the eigenvalue of the trailing 2x2 block + * closest to the bottom-right corner (Trefethen & Bau 13.4.1). + * + * @note Algorithm (Trefethen & Bau 13.4, Golub & Van Loan 8.4.3): + * Phase 0 (N >= 3): Givens tridiagonalization. A general symmetric + * matrix is NOT suitable for implicit QR (the bulge chase only + * preserves the tridiagonal structure), so first reduce A with + * adjacent Givens similarities A <- G A G^T (rotations applied + * BOTTOM-UP, i = N-2 down to k+1, per column k), accumulating + * U <- U G^T, until A is symmetric tridiagonal and + * A_orig = U A U^T. (N = 2 needs no reduction.) + * Phase 1: iterate until A is diagonal: + * 1. Deflate: zero out subdiagonal entries at/under the tolerance + * (scaled by the adjacent diagonal magnitudes). + * 2. Scan for the trailing unreduced block [lo, hi]. + * - block of size 1: A[hi][hi] is a converged eigenvalue, done. + * - block of size 2: solve the 2x2 eigenproblem in closed form + * and fold its eigenvector matrix into V. + * - block larger: one implicit Wilkinson-shifted QR step + * (bulge chasing with Givens rotations; the shift is applied + * to the ENTIRE active block [lo, hi], not just the trailing + * 2x2 -- the first Givens rotation must be formed from + * (A[lo][lo] - mu, A[lo+1][lo])). Every rotation is folded + * into V. + * Phase 2: eigenvalues = diag(A), sorted DESCENDING (eigenvector + * columns swapped to match), and the true eigenvectors of the + * ORIGINAL matrix are U * V. + * + * @note If maxIterations is exhausted before convergence the best-effort + * (partially diagonalized) values on the diagonal are returned. + */ +namespace QR { + +/** + * @brief Compute the eigenvalues and eigenvectors of a square matrix + * + * @param matrixToDecompose The matrix to take eigenvalues of (not + * modified). MUST be symmetric for N >= 3. + * @param eigenVectors a buffer that will contain the eigenvectors in its + * COLUMNS, sorted by descending eigenvalue (column i is the + * eigenvector for eigenValues[i]). + * @param eigenValues a buffer that will contain the eigenvalues sorted + * DESCENDING (largest first). + * @param maxIterations the number of QR steps to perform before giving up + * on reaching the given tolerance + * @param tolerance the level of accuracy to obtain before stopping; a + * subdiagonal entry is deflated when |A[i+1][i]| <= tolerance * + * (|A[i][i]| + |A[i+1][i+1]|). For float32 arithmetic, values + * around 1e-6 are a sensible choice (single-precision epsilon is + * ~1.2e-7). + */ +template +void EigenQR(Matrix &matrixToDecompose, Matrix &eigenVectors, + Matrix &eigenValues, uint32_t maxIterations, float tolerance); + +/** + * @brief Compute a Givens rotation that zeros the bottom entry of (a, b) + * + * Given the column vector (a, b), produces (c, s) defining the 2x2 + * rotation + * R = [ c s ] + * [ -s c ] + * such that R * (a, b)^T = (r, 0)^T with r = +hypot(a, b) >= 0, i.e. + * c = a / r, s = b / r. + * + * If (a, b) == (0, 0) the identity rotation (c = 1, s = 0) is returned. + */ +static void GivensRotation(float a, float b, float &c, float &s); + +/** + * @brief Apply the similarity transform A <- G A G^T on rows/cols (i, i+1) + * + * G = [ c s ] on the (i, i+1) block, identity elsewhere, where G is the + * [ -s c ] + * ZEROING rotation (G * (x, y)^T = (r, 0)^T) -- the orientation used by + * the implicit QR chase: A = Q R with Q = G^T gives the next iterate + * R Q = G A G^T. With (c, s) = GivensRotation(A[i][i], A[i+1][i]) the + * (i+1, i) entry is zeroed by the left multiplication and the bulge is + * chased along the superdiagonal by the right one. The matrix must be + * symmetric on entry (guaranteed by construction in the QR iteration: + * symmetric input stays symmetric under similarity by an orthogonal + * matrix). Updates the full matrix, not just the tridiagonal structure. + */ +template +static void ApplyRotationBothSides(Matrix &A, uint8_t i, float c, + float s); + +/** + * @brief Accumulate eigenvectors: V <- V G^T on columns (i, i+1) + * + * G^T = [ c -s ] on columns (i, i+1), identity elsewhere, where G = + * [ s c ] + * [ c, s ] / [ -s, c ] is the zeroing rotation paired with + * ApplyRotationBothSides. Applied to all rows: + * V[r][i] -> c V[r][i] + s V[r][i+1] + * V[r][i+1] -> -s V[r][i] + c V[r][i+1] + * + * Every QR step's rotation is folded into V this way so that, together + * with A <- G A G^T, the invariant A_orig = V A V^T is preserved at every + * step (each step is A <- R Q with Q = G^T the orthogonal factor, and + * the orthogonal factors multiply as G1^T G2^T ... in application order). + * At convergence A_orig = V D V^T and the columns of V are the + * eigenvectors. + */ +template +static void ApplyRotationToVectors(Matrix &V, uint8_t i, float c, + float s); + +/** + * @brief Wilkinson shift for a symmetric tridiagonal + * + * Given the trailing 2x2 block + * [ a b ] + * [ b d ] + * returns the eigenvalue of that block that is closest to d. This is the + * empirically best shift for the QR iteration (Trefethen & Bau 13.4.1). + * + * mu = (a+d)/2 - sign(a-d) * sqrt(((a-d)/2)^2 + b^2) + * (with sign(0) taken as +1). + */ +static float WilkinsonShift(float a, float b, float d); + +/** + * @brief Solve the 2x2 eigenproblem of block rows/cols (lo, lo+1) + * + * Solves the (possibly nonsymmetric) 2x2 block + * [ A[lo][lo] A[lo][lo+1] ] + * [ A[lo+1][lo] A[lo+1][lo+1] ] + * in closed form (characteristic polynomial + eigenvector back-substitution). + * + * @param A the matrix containing the block (not modified) + * @param lo the row/col index of the top-left corner of the block + * @param lambdaHi (out) the LARGER eigenvalue + * @param lambdaLo (out) the smaller eigenvalue + * @param c (out), s (out) eigenvector pair as an orthogonal matrix + * Vblock = [ c -s ] whose columns are the eigenvectors: column 0 + * [ s c ] + * (c, s) is the unit eigenvector for lambdaHi, column 1 (-s, c) is + * the unit eigenvector for lambdaLo. + * + * Note: the caller applies Vblock to its eigenvector accumulator with + * V <- V * Vblock (i.e. V[r][lo] = c*x + s*y, + * V[r][lo+1] = -s*x + c*y). Vblock has the + * SAME [ c -s; s c ] form as the G^T factor used by + * ApplyRotationToVectors, so both folding operations follow one uniform + * convention. + */ +template +static void Solve2x2Eigen(const Matrix &A, uint8_t lo, float &lambdaHi, + float &lambdaLo, float &c, float &s); + +/** + * @brief Deflate (zero out) subdiagonal entries that are at/under tolerance + * + * For each i in [lo, hi): if |A[i+1][i]| <= tolerance * + * (|A[i][i]| + |A[i+1][i+1]|), sets A[i+1][i] = A[i][i+1] = 0, splitting + * the matrix into smaller independent blocks. + */ +template +static void Deflate(Matrix &A, uint8_t lo, uint8_t hi, float tolerance); + +/** + * @brief Reduce a symmetric matrix to symmetric tridiagonal form + * + * Chases each column's entries below the subdiagonal to zero with + * adjacent Givens similarities (Golub & Van Loan 8.3.1, Givens variant): + * for column k = 0..N-3, rotations on (N-2, N-1), (N-3, N-2), ... + * (k+1, k+2) -- BOTTOM-UP, each formed from the current (A[i][k], + * A[i+1][k]) -- zero A[k+2..N-1, k] one by one. A top-down pass would not + * work: the rotation that zeros A[i+1][k] would be undone by the later + * rotation on (i+1, i+2) forming a new nonzero at A[i][k]. Each rotation + * is applied to A as a similarity (A <- G A G^T) and accumulated into U + * (U <- U G^T), so on return: + * - A is symmetric tridiagonal (off-tridiagonal entries EXACTLY zero), + * - A_orig = U A U^T (i.e. U^T A_orig U = A). + * + * U is initialized to the identity internally (its input contents are + * ignored). + */ +template +static void Tridiagonalize(Matrix &A, Matrix &U); + +} // namespace QR + +#ifndef QR_H_ +#include "QR.cpp" +#endif diff --git a/unit-tests/CMakeLists.txt b/unit-tests/CMakeLists.txt index 3ecc1b1..e348fa6 100644 --- a/unit-tests/CMakeLists.txt +++ b/unit-tests/CMakeLists.txt @@ -13,6 +13,7 @@ add_executable(matrix-tests matrix-tests.cpp) target_link_libraries(matrix-tests PRIVATE matrix + qr Catch2::Catch2WithMain ) @@ -52,4 +53,14 @@ target_link_libraries(svd-integration-test matrix svd Catch2::Catch2WithMain +) + +# QR building block tests +add_executable(qr-build-blocks-tests qr-build-blocks-tests.cpp) + +target_link_libraries(qr-build-blocks-tests + PRIVATE + matrix + qr + Catch2::Catch2WithMain ) \ No newline at end of file diff --git a/unit-tests/matrix-tests.cpp b/unit-tests/matrix-tests.cpp index 32e6c56..c2eb205 100644 --- a/unit-tests/matrix-tests.cpp +++ b/unit-tests/matrix-tests.cpp @@ -4,6 +4,7 @@ // include the module you're going to test next #include "Matrix.hpp" +#include "QR.hpp" #include "SVD.hpp" // any other libraries @@ -601,8 +602,78 @@ TEST_CASE("QR Decompositions", "Matrix") { } } +// ============================================================================ +// Eigen QR Helpers (scipy references; eigenvector checks are sign-invariant) +// ============================================================================ + +/** + * @brief Normalized eigenpair residual ||A v - lambda v|| / (||A||_F + |lambda|) + */ +template +static float eigenResidual(const Matrix &A, float lambda, + const Matrix &v) { + Matrix Av{}; + A.Mult(v, Av); + float sum = 0.0f; + float frob = 0.0f; + for (uint8_t i = 0; i < N; i++) { + float d = Av.Get(i, 0) - lambda * v.Get(i, 0); + sum += d * d; + for (uint8_t j = 0; j < N; j++) { + float a = A.Get(i, j); + frob += a * a; + } + } + float scale = sqrtf(frob) + fabsf(lambda); + return sqrtf(sum) / scale; +} + +/** + * @brief Column of the eigenvector matrix; used for the residual check. + */ +template +static Matrix eigenColumn(const Matrix &V, uint8_t col) { + Matrix v{}; + for (uint8_t i = 0; i < N; i++) { + v[i][0] = V.Get(i, col); + } + return v; +} + +/** + * @brief Check V^T V ~ I (eigenvectors orthonormal). + */ +template +static bool isOrthogonal(const Matrix &V, float tol = 1e-4f) { + Matrix Vt = V.Transpose(); + Matrix VtV{}; + Vt.Mult(V, VtV); + for (uint8_t i = 0; i < N; i++) { + for (uint8_t j = 0; j < N; j++) { + float expected = (i == j) ? 1.0f : 0.0f; + if (fabsf(VtV.Get(i, j) - expected) > tol) { + return false; + } + } + } + return true; +} + +/** + * @brief Sign-invariant component check: |actual| within max(1e-4, 1e-3*|ref|) + * of ref (ref is the ABSOLUTE value from the scipy reference). + */ +static bool componentMatches(float actual, float refAbs) { + float a = fabsf(actual); + float tol = 1e-4f; + if (refAbs * 1e-3f > tol) { + tol = refAbs * 1e-3f; + } + return fabsf(a - refAbs) <= tol; +} + TEST_CASE("Eigenvalues and Vectors", "Matrix") { - SECTION("2x2 Eigen") { + SECTION("2x2 Eigen (nonsymmetric, closed form)") { Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f}; Matrix<2, 2> vectors{}; Matrix<2, 1> values{}; @@ -615,28 +686,206 @@ TEST_CASE("Eigenvalues and Vectors", "Matrix") { REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(-0.372281f, 1e-4f)); } - SECTION("3x3 Rank Defficient Eigen") { - SKIP("Skipping this because QR decomposition isn't ready for it"); - // this symmetrix tridiagonal matrix is well behaved for testing - Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9}; + // Reference values: numpy.linalg.eigh on float32 matrices. + // Eigenvector component references are ABSOLUTE values (signs arbitrary). + SECTION("3x3 Symmetric Eigen") { + Matrix<3, 3> A{1, 2, 3, 2, 5, 8, 3, 8, 9}; Matrix<3, 3> vectors{}; Matrix<3, 1> values{}; - A.EigenQR(vectors, values, 1000000, 1e-8f); + A.EigenQR(vectors, values, 10000, 1e-6f); - std::string strBuf1 = ""; - vectors.ToString(strBuf1); - std::cout << "Vectors:\n" << strBuf1 << std::endl; - strBuf1 = ""; - values.ToString(strBuf1); - std::cout << "Values:\n" << strBuf1 << std::endl; + // eigenvalues (descending) + REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(16.102417f, 1e-4f)); + REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(0.191920f, 1e-4f)); + REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(-1.2943381f, 1e-4f)); - REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.23197f, 1e-4f)); - REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.525322f, 1e-4f)); - REQUIRE_THAT(vectors[2][0], Catch::Matchers::WithinRel(0.81867f, 1e-4f)); - REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(-1.11684f, 1e-4f)); + // eigenvector |components| (sign-invariant) + REQUIRE(componentMatches(vectors[0][0], 0.231657207f)); + REQUIRE(componentMatches(vectors[1][0], 0.59582746f)); + REQUIRE(componentMatches(vectors[2][0], 0.768976331f)); + REQUIRE(componentMatches(vectors[0][1], 0.956842422f)); + REQUIRE(componentMatches(vectors[1][1], 0.282139271f)); + REQUIRE(componentMatches(vectors[2][1], 0.0696421042f)); + REQUIRE(componentMatches(vectors[0][2], 0.175463736f)); + REQUIRE(componentMatches(vectors[1][2], 0.75192225f)); + REQUIRE(componentMatches(vectors[2][2], 0.635472536f)); + + // eigenvectors orthonormal; eigenpair residuals small + REQUIRE(isOrthogonal(vectors)); + for (uint8_t col = 0; col < 3; col++) { + REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) < + 1e-4f); + } + } + + SECTION("3x3 Rank Deficient Eigen") { + // A = v v^T with v = [1, 2, 3]: eigenvalues {14, 0, 0} + Matrix<3, 3> A{1, 2, 3, 2, 4, 6, 3, 6, 9}; + Matrix<3, 3> vectors{}; + Matrix<3, 1> values{}; + A.EigenQR(vectors, values, 10000, 1e-6f); + + REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(14.0f, 1e-4f)); REQUIRE_THAT(values[1][0], Catch::Matchers::WithinAbs(0.0f, 1e-4f)); - REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(16.1168f, 1e-4f)); + REQUIRE_THAT(values[2][0], Catch::Matchers::WithinAbs(0.0f, 1e-4f)); + + // dominant eigenvector is v/|v| (sign-invariant); the two null-space + // eigenvectors may be ANY orthonormal basis of the null plane, so only + // orthogonality + residuals are checked for the full matrix. + REQUIRE(componentMatches(vectors[0][0], 0.267261237f)); + REQUIRE(componentMatches(vectors[1][0], 0.534522474f)); + REQUIRE(componentMatches(vectors[2][0], 0.801783741f)); + REQUIRE(isOrthogonal(vectors)); + for (uint8_t col = 0; col < 3; col++) { + REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) < + 1e-4f); + } + } + + SECTION("4x4 Symmetric Eigen") { + Matrix<4, 4> A{2, 1, 0, 1, 1, 3, 1, 0, 0, 1, 4, 1, 1, 0, 1, 5}; + Matrix<4, 4> vectors{}; + Matrix<4, 1> values{}; + A.EigenQR(vectors, values, 10000, 1e-6f); + + // eigenvalues are exactly {6, 4, 3, 1} + REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(6.0f, 1e-4f)); + REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(4.0f, 1e-4f)); + REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(3.0f, 1e-4f)); + REQUIRE_THAT(values[3][0], Catch::Matchers::WithinRel(1.0f, 1e-4f)); + + // eigenvector |components| (sign-invariant) + REQUIRE(componentMatches(vectors[0][0], 0.258198887f)); + REQUIRE(componentMatches(vectors[1][0], 0.258198887f)); + REQUIRE(componentMatches(vectors[2][0], 0.516397774f)); + REQUIRE(componentMatches(vectors[3][0], 0.774596691f)); + REQUIRE(componentMatches(vectors[0][1], 0.0f)); + REQUIRE(componentMatches(vectors[1][1], 0.577350259f)); + REQUIRE(componentMatches(vectors[2][1], 0.577350259f)); + REQUIRE(componentMatches(vectors[3][1], 0.577350259f)); + REQUIRE(componentMatches(vectors[0][2], 0.577350259f)); + REQUIRE(componentMatches(vectors[1][2], 0.577350259f)); + REQUIRE(componentMatches(vectors[2][2], 0.577350259f)); + REQUIRE(componentMatches(vectors[3][2], 0.0f)); + REQUIRE(componentMatches(vectors[0][3], 0.774596691f)); + REQUIRE(componentMatches(vectors[1][3], 0.516397774f)); + REQUIRE(componentMatches(vectors[2][3], 0.258198887f)); + REQUIRE(componentMatches(vectors[3][3], 0.258198887f)); + + REQUIRE(isOrthogonal(vectors)); + for (uint8_t col = 0; col < 4; col++) { + REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) < + 1e-4f); + } + } + + SECTION("5x5 Symmetric Eigen") { + Matrix<5, 5> A{3, 1, 0, 0, 1, 1, 4, 1, 0, 0, 0, 1, 5, 1, 0, 0, 0, 1, 6, 1, + 1, 0, 0, 1, 7}; + Matrix<5, 5> vectors{}; + Matrix<5, 1> values{}; + A.EigenQR(vectors, values, 10000, 1e-6f); + + // eigenvalues (descending) + REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(7.90154457f, 1e-4f)); + REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(6.20044184f, 1e-4f)); + REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(5.14503145f, 1e-4f)); + REQUIRE_THAT(values[3][0], Catch::Matchers::WithinRel(3.61823463f, 1e-4f)); + REQUIRE_THAT(values[4][0], Catch::Matchers::WithinRel(2.13474774f, 1e-4f)); + + // eigenvector |components| (sign-invariant) + REQUIRE(componentMatches(vectors[0][0], 0.182430908f)); + REQUIRE(componentMatches(vectors[1][0], 0.102749094f)); + REQUIRE(componentMatches(vectors[2][0], 0.21844925f)); + REQUIRE(componentMatches(vectors[3][0], 0.531091094f)); + REQUIRE(componentMatches(vectors[4][0], 0.791444063f)); + REQUIRE(componentMatches(vectors[0][1], 0.0877681747f)); + REQUIRE(componentMatches(vectors[1][1], 0.245861098f)); + REQUIRE(componentMatches(vectors[2][1], 0.628771126f)); + REQUIRE(componentMatches(vectors[3][1], 0.508941948f)); + REQUIRE(componentMatches(vectors[4][1], 0.526758015f)); + REQUIRE(componentMatches(vectors[0][2], 0.349721253f)); + REQUIRE(componentMatches(vectors[1][2], 0.628706098f)); + REQUIRE(componentMatches(vectors[2][2], 0.370167077f)); + REQUIRE(componentMatches(vectors[3][2], 0.575020194f)); + REQUIRE(componentMatches(vectors[4][2], 0.121457018f)); + REQUIRE(componentMatches(vectors[0][3], 0.429638386f)); + REQUIRE(componentMatches(vectors[1][3], 0.498553723f)); + REQUIRE(componentMatches(vectors[2][3], 0.619968951f)); + REQUIRE(componentMatches(vectors[3][3], 0.35809797f)); + REQUIRE(componentMatches(vectors[4][3], 0.232936427f)); + REQUIRE(componentMatches(vectors[0][4], 0.807540476f)); + REQUIRE(componentMatches(vectors[1][4], 0.534011006f)); + REQUIRE(componentMatches(vectors[2][4], 0.188524753f)); + REQUIRE(componentMatches(vectors[3][4], 0.00615991838f)); + REQUIRE(componentMatches(vectors[4][4], 0.164715111f)); + + REQUIRE(isOrthogonal(vectors)); + for (uint8_t col = 0; col < 5; col++) { + REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) < + 1e-4f); + } + } + + SECTION("6x6 Symmetric Eigen") { + Matrix<6, 6> A{4, 1, 0, 0, 0, 1, 1, 5, 1, 0, 0, 0, 0, 1, 6, 1, 0, 0, 0, 0, + 1, 7, 1, 0, 0, 0, 0, 1, 8, 1, 1, 0, 0, 0, 1, 3}; + Matrix<6, 6> vectors{}; + Matrix<6, 1> values{}; + A.EigenQR(vectors, values, 10000, 1e-6f); + + // eigenvalues (descending) + REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(8.86080551f, 1e-4f)); + REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(7.25410175f, 1e-4f)); + REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(6.11490774f, 1e-4f)); + REQUIRE_THAT(values[3][0], Catch::Matchers::WithinRel(4.88509226f, 1e-4f)); + REQUIRE_THAT(values[4][0], Catch::Matchers::WithinRel(3.74589825f, 1e-4f)); + REQUIRE_THAT(values[5][0], Catch::Matchers::WithinRel(2.13919425f, 1e-4f)); + + // eigenvector |components| (sign-invariant) + REQUIRE(componentMatches(vectors[0][0], 0.0430923924f)); + REQUIRE(componentMatches(vectors[1][0], 0.0662503168f)); + REQUIRE(componentMatches(vectors[2][0], 0.212687209f)); + REQUIRE(componentMatches(vectors[3][0], 0.542206466f)); + REQUIRE(componentMatches(vectors[4][0], 0.7962538f)); + REQUIRE(componentMatches(vectors[5][0], 0.143213451f)); + REQUIRE(componentMatches(vectors[0][1], 0.0623276457f)); + REQUIRE(componentMatches(vectors[1][1], 0.307613879f)); + REQUIRE(componentMatches(vectors[2][1], 0.631065309f)); + REQUIRE(componentMatches(vectors[3][1], 0.483806193f)); + REQUIRE(componentMatches(vectors[4][1], 0.508129358f)); + REQUIRE(componentMatches(vectors[5][1], 0.104793385f)); + REQUIRE(componentMatches(vectors[0][2], 0.374228716f)); + REQUIRE(componentMatches(vectors[1][2], 0.605694294f)); + REQUIRE(componentMatches(vectors[2][2], 0.301064402f)); + REQUIRE(componentMatches(vectors[3][2], 0.571099699f)); + REQUIRE(componentMatches(vectors[4][2], 0.204411641f)); + REQUIRE(componentMatches(vectors[5][2], 0.185764849f)); + REQUIRE(componentMatches(vectors[0][3], 0.571099699f)); + REQUIRE(componentMatches(vectors[1][3], 0.301064402f)); + REQUIRE(componentMatches(vectors[2][3], 0.605694294f)); + REQUIRE(componentMatches(vectors[3][3], 0.374228716f)); + REQUIRE(componentMatches(vectors[4][3], 0.185764849f)); + REQUIRE(componentMatches(vectors[5][3], 0.204411641f)); + REQUIRE(componentMatches(vectors[0][4], 0.483806193f)); + REQUIRE(componentMatches(vectors[1][4], 0.631065309f)); + REQUIRE(componentMatches(vectors[2][4], 0.307613879f)); + REQUIRE(componentMatches(vectors[3][4], 0.0623276457f)); + REQUIRE(componentMatches(vectors[4][4], 0.104793385f)); + REQUIRE(componentMatches(vectors[5][4], 0.508129358f)); + REQUIRE(componentMatches(vectors[0][5], 0.542206466f)); + REQUIRE(componentMatches(vectors[1][5], 0.212687209f)); + REQUIRE(componentMatches(vectors[2][5], 0.0662503168f)); + REQUIRE(componentMatches(vectors[3][5], 0.0430923924f)); + REQUIRE(componentMatches(vectors[4][5], 0.143213451f)); + REQUIRE(componentMatches(vectors[5][5], 0.7962538f)); + + REQUIRE(isOrthogonal(vectors)); + for (uint8_t col = 0; col < 6; col++) { + REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) < + 1e-4f); + } } } diff --git a/unit-tests/qr-build-blocks-tests.cpp b/unit-tests/qr-build-blocks-tests.cpp new file mode 100644 index 0000000..8b9d760 --- /dev/null +++ b/unit-tests/qr-build-blocks-tests.cpp @@ -0,0 +1,581 @@ +// include the unit test framework first +#include +#include + +// include the module you're going to test next +#include "Matrix.hpp" +#include "QR.hpp" + +// any other libraries +#include +#include +#include + +// ============================================================================ +// Helpers +// ============================================================================ + +/** + * @brief Frobenius norm of an N x N matrix. + */ +template +static float frob(const Matrix &M) { + float sum = 0.0f; + for (uint8_t i = 0; i < N; i++) + for (uint8_t j = 0; j < N; j++) { + float v = M.Get(i, j); + sum += v * v; + } + return sqrtf(sum); +} + +/** + * @brief Check M is orthogonal (M^T M ~ I). + */ +template +static bool isOrthogonal(const Matrix &M, float tol = 1e-5f) { + Matrix Mt = M.Transpose(); + Matrix MtM{}; + Mt.Mult(M, MtM); + for (uint8_t i = 0; i < N; i++) + for (uint8_t j = 0; j < N; j++) { + float expected = (i == j) ? 1.0f : 0.0f; + if (fabsf(MtM.Get(i, j) - expected) > tol) + return false; + } + return true; +} + +/** + * @brief 3x3 trace. + */ +static float trace3(const Matrix<3, 3> &A) { + return A.Get(0, 0) + A.Get(1, 1) + A.Get(2, 2); +} + +/** + * @brief 3x3 sum of principal 2x2 minors (2nd elementary invariant). + */ +static float e2_3x3(const Matrix<3, 3> &A) { + return A.Get(0, 0) * A.Get(1, 1) - A.Get(0, 1) * A.Get(0, 1) + + A.Get(0, 0) * A.Get(2, 2) - A.Get(0, 2) * A.Get(0, 2) + + A.Get(1, 1) * A.Get(2, 2) - A.Get(1, 2) * A.Get(1, 2); +} + +/** + * @brief 3x3 determinant. + */ +static float det3(const Matrix<3, 3> &A) { + return A.Get(0, 0) * + (A.Get(1, 1) * A.Get(2, 2) - A.Get(1, 2) * A.Get(2, 1)) - + A.Get(0, 1) * + (A.Get(1, 0) * A.Get(2, 2) - A.Get(1, 2) * A.Get(2, 0)) + + A.Get(0, 2) * + (A.Get(1, 0) * A.Get(2, 1) - A.Get(1, 1) * A.Get(2, 0)); +} + +/** + * @brief Sign-invariant comparison of |actual| against refAbs. + */ +static bool matchesAbs(float actual, float refAbs, float relTol = 1e-5f, + float absTol = 1e-6f) { + float a = fabsf(actual); + if (refAbs < 1e-3f) + return a < absTol + relTol; + return fabsf(a - refAbs) <= relTol * refAbs; +} + +// ============================================================================ +// TEST 1: GivensRotation +// ============================================================================ +TEST_CASE("QR Building Block: GivensRotation", "[Matrix][QR]") { + // R = [[c, s], [-s, c]] must satisfy R * (a, b)^T = (r, 0)^T. + + { + // Reference: hypot(2, 1) = sqrt(5) = 2.236067977 + float c = 0, s = 0; + QR::GivensRotation(2.0f, 1.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(0.894427191f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(0.447213595f, 1e-6f)); + REQUIRE_THAT(c * 2.0f + s * 1.0f, + Catch::Matchers::WithinRel(2.236067977f, 1e-6f)); + REQUIRE_THAT(-s * 2.0f + c * 1.0f, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + { + // Reference: hypot(3, 4) = 5 exactly + float c = 0, s = 0; + QR::GivensRotation(3.0f, 4.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(0.6f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(0.8f, 1e-6f)); + REQUIRE_THAT(c * 3.0f + s * 4.0f, Catch::Matchers::WithinRel(5.0f, 1e-6f)); + REQUIRE_THAT(-s * 3.0f + c * 4.0f, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } + + { + // Pure second component: c = 0, s = 1 + float c = 1, s = 1; + QR::GivensRotation(0.0f, 5.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(1.0f, 1e-6f)); + } + + { + // Zero vector: identity rotation + float c = 0, s = 0; + QR::GivensRotation(0.0f, 0.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(1.0f, 1e-7f)); + REQUIRE_THAT(s, Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + } + + { + // Negative first component preserves the sign of c + float c = 0, s = 0; + QR::GivensRotation(-2.0f, 1.0f, c, s); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(-0.894427191f, 1e-6f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(0.447213595f, 1e-6f)); + REQUIRE_THAT(-s * -2.0f + c * 1.0f, Catch::Matchers::WithinAbs(0.0f, 1e-6f)); + } +} + +// ============================================================================ +// TEST 2: ApplyRotationBothSides (similarity A <- G A G^T) +// ============================================================================ +TEST_CASE("QR Building Block: ApplyRotationBothSides", "[Matrix][QR]") { + // Reference (numpy, float64): A = [[2,1,0],[1,3,1],[0,1,4]], i = 0, + // Givens(2,1) -> G A G^T = + // [[ 3.0, 1.0, 0.447213595], + // [ 1.0, 2.0, 0.894427191], + // [ 0.447213595, 0.894427191, 4.0]] + // (Note: G A G^T with G zeroing (2,1) sends the A[0][1] coupling into the + // (0,2) corner, NOT into the subdiagonal -- the subdiagonal-zeroing happens + // in the QR chase context where the bulge column has the right shape.) + { + Matrix<3, 3> A{2, 1, 0, 1, 3, 1, 0, 1, 4}; + float c = 0.894427191f, s = 0.447213595f; + + QR::ApplyRotationBothSides(A, 0, c, s); + + REQUIRE_THAT(A.Get(0, 0), Catch::Matchers::WithinRel(3.0f, 1e-5f)); + REQUIRE_THAT(A.Get(0, 1), Catch::Matchers::WithinRel(1.0f, 1e-5f)); + REQUIRE_THAT(A.Get(0, 2), + Catch::Matchers::WithinRel(0.447213595f, 1e-5f)); + REQUIRE_THAT(A.Get(1, 1), Catch::Matchers::WithinRel(2.0f, 1e-5f)); + REQUIRE_THAT(A.Get(1, 2), + Catch::Matchers::WithinRel(0.894427191f, 1e-5f)); + REQUIRE_THAT(A.Get(2, 2), Catch::Matchers::WithinRel(4.0f, 1e-5f)); + + // Symmetry must be preserved exactly in both triangles + for (uint8_t i = 0; i < 3; i++) + for (uint8_t j = 0; j < 3; j++) + REQUIRE(A.Get(i, j) == A.Get(j, i)); + } + + // Same check at i = 1. + // Reference (numpy, float64): B = [[5,0,1],[0,6,2],[1,2,7]], i = 1, + // Givens(6,2) -> G B G^T = + // [[ 5.0, 0.316227766, 0.948683298], + // [ 0.316227766, 7.3, 1.9], + // [ 0.948683298, 1.9, 5.7]] + { + Matrix<3, 3> B{5, 0, 1, 0, 6, 2, 1, 2, 7}; + float c = 0.948683298f, s = 0.316227766f; + + QR::ApplyRotationBothSides(B, 1, c, s); + + REQUIRE_THAT(B.Get(0, 0), Catch::Matchers::WithinRel(5.0f, 1e-5f)); + REQUIRE_THAT(B.Get(0, 1), + Catch::Matchers::WithinRel(0.316227766f, 1e-5f)); + REQUIRE_THAT(B.Get(0, 2), + Catch::Matchers::WithinRel(0.948683298f, 1e-5f)); + REQUIRE_THAT(B.Get(1, 1), Catch::Matchers::WithinRel(7.3f, 1e-5f)); + REQUIRE_THAT(B.Get(1, 2), Catch::Matchers::WithinRel(1.9f, 1e-5f)); + REQUIRE_THAT(B.Get(2, 2), Catch::Matchers::WithinRel(5.7f, 1e-5f)); + + for (uint8_t i = 0; i < 3; i++) + for (uint8_t j = 0; j < 3; j++) + REQUIRE(B.Get(i, j) == B.Get(j, i)); + } + + // Identity rotation leaves the matrix unchanged + { + Matrix<3, 3> C{1, 2, 3, 2, 4, 5, 3, 5, 6}; + QR::ApplyRotationBothSides(C, 1, 1.0f, 0.0f); + REQUIRE(C.Get(0, 0) == 1.0f); + REQUIRE(C.Get(0, 1) == 2.0f); + REQUIRE(C.Get(0, 2) == 3.0f); + REQUIRE(C.Get(1, 1) == 4.0f); + REQUIRE(C.Get(1, 2) == 5.0f); + REQUIRE(C.Get(2, 2) == 6.0f); + } + + // Spectrum invariants (trace, Frobenius norm) are preserved. (c, s) + // must be a unit vector for G A G^T to be a similarity transform. + { + Matrix<3, 3> D{1, 2, 3, 2, 5, 8, 3, 8, 9}; + float tr = trace3(D); + float fn = frob(D); + float c = 0.6f, s = 0.8f; + QR::ApplyRotationBothSides(D, 0, c, s); + REQUIRE_THAT(trace3(D), Catch::Matchers::WithinRel(tr, 1e-5f)); + REQUIRE_THAT(frob(D), Catch::Matchers::WithinRel(fn, 1e-5f)); + } +} +// ============================================================================ +// TEST 3: ApplyRotationToVectors (V <- V G^T) +// ============================================================================ +TEST_CASE("QR Building Block: ApplyRotationToVectors", "[Matrix][QR]") { + // V = I, i = 0, Givens(2,1): V <- I * G^T with G^T = [[c, -s], [s, c]] = + // [[ c, -s, 0], + // [ s, c, 0], + // [ 0, 0, 1]] + { + Matrix<3, 3> V{0}; + V[0][0] = 1; + V[1][1] = 1; + V[2][2] = 1; + float c = 0.894427191f, s = 0.447213595f; + + QR::ApplyRotationToVectors(V, 0, c, s); + + REQUIRE_THAT(V.Get(0, 0), Catch::Matchers::WithinRel(0.894427191f, 1e-6f)); + REQUIRE_THAT(V.Get(0, 1), Catch::Matchers::WithinRel(-0.447213595f, 1e-6f)); + REQUIRE_THAT(V.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + REQUIRE_THAT(V.Get(1, 0), Catch::Matchers::WithinRel(0.447213595f, 1e-6f)); + REQUIRE_THAT(V.Get(1, 1), Catch::Matchers::WithinRel(0.894427191f, 1e-6f)); + REQUIRE_THAT(V.Get(1, 2), Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + REQUIRE_THAT(V.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + REQUIRE_THAT(V.Get(2, 1), Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + REQUIRE_THAT(V.Get(2, 2), Catch::Matchers::WithinRel(1.0f, 1e-7f)); + + // Product of rotations must stay orthogonal + REQUIRE(isOrthogonal(V)); + } + + // Two successive rotations accumulate (V <- V G1^T G2^T) + // Reference (numpy, float64): + // [[ 0.894427191, -0.424264069, 0.141421356], + // [ 0.447213595, 0.848528137, -0.282842712], + // [ 0.0, 0.316227766, 0.948683298]] + { + Matrix<3, 3> V{0}; + V[0][0] = 1; + V[1][1] = 1; + V[2][2] = 1; + QR::ApplyRotationToVectors(V, 0, 0.894427191f, 0.447213595f); + QR::ApplyRotationToVectors(V, 1, 0.948683298f, 0.316227766f); + REQUIRE(isOrthogonal(V)); + // Column 0 was only touched by the first rotation + REQUIRE_THAT(V.Get(0, 0), Catch::Matchers::WithinRel(0.894427191f, 1e-5f)); + REQUIRE_THAT(V.Get(1, 0), Catch::Matchers::WithinRel(0.447213595f, 1e-5f)); + REQUIRE_THAT(V.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + REQUIRE_THAT(V.Get(0, 1), Catch::Matchers::WithinRel(-0.424264069f, 1e-5f)); + REQUIRE_THAT(V.Get(0, 2), Catch::Matchers::WithinRel(0.141421356f, 1e-5f)); + REQUIRE_THAT(V.Get(1, 2), Catch::Matchers::WithinRel(-0.282842712f, 1e-5f)); + REQUIRE_THAT(V.Get(2, 1), Catch::Matchers::WithinRel(0.316227766f, 1e-5f)); + REQUIRE_THAT(V.Get(2, 2), Catch::Matchers::WithinRel(0.948683298f, 1e-5f)); + } +} + +// ============================================================================ +// TEST 4: WilkinsonShift +// ============================================================================ +TEST_CASE("QR Building Block: WilkinsonShift", "[Matrix][QR]") { + // mu = (a+d)/2 - sign(a-d) * sqrt(((a-d)/2)^2 + b^2) + // Reference: eigenvalues of [[2,1],[1,4]] are 1.5858, 4.4142; closest + // to d = 4 is 4.414213562. + REQUIRE_THAT(QR::WilkinsonShift(2.0f, 1.0f, 4.0f), + Catch::Matchers::WithinRel(4.414213562f, 1e-6f)); + + // [[5,2],[2,1]]: eigenvalues 0.1716, 5.8284; closest to d = 1 is 0.171572875 + REQUIRE_THAT(QR::WilkinsonShift(5.0f, 2.0f, 1.0f), + Catch::Matchers::WithinRel(0.171572875f, 1e-5f)); + + // Zero off-diagonal: returns d itself (sign(0) = +1 picks d, not a) + REQUIRE_THAT(QR::WilkinsonShift(3.0f, 0.0f, 7.0f), + Catch::Matchers::WithinRel(7.0f, 1e-7f)); + REQUIRE_THAT(QR::WilkinsonShift(7.0f, 0.0f, 3.0f), + Catch::Matchers::WithinRel(3.0f, 1e-7f)); + + // a == d: shift is the larger-magnitude off-diagonal combination + // [[1,3],[3,1]]: eigenvalues -2, 4; closest to d = 1 is -2 + REQUIRE_THAT(QR::WilkinsonShift(1.0f, 3.0f, 1.0f), + Catch::Matchers::WithinRel(-2.0f, 1e-6f)); +} + +// ============================================================================ +// TEST 5: Solve2x2Eigen +// ============================================================================ +TEST_CASE("QR Building Block: Solve2x2Eigen", "[Matrix][QR]") { + // Symmetric block [[2,1],[1,3]]: + // eigenvalues 1.381966011, 3.618033989; + // eigenvector of 3.618033989 is +/- (0.525731112, 0.850650808) + { + Matrix<2, 2> A{2, 1, 1, 3}; + float lHi = 0, lLo = 0, c = 0, s = 0; + QR::Solve2x2Eigen(A, 0, lHi, lLo, c, s); + + REQUIRE_THAT(lHi, Catch::Matchers::WithinRel(3.618033989f, 1e-6f)); + REQUIRE_THAT(lLo, Catch::Matchers::WithinRel(1.381966011f, 1e-6f)); + REQUIRE(matchesAbs(c, 0.525731112f)); + REQUIRE(matchesAbs(s, 0.850650808f)); + + // Residual: A * vHi = lHi * vHi with vHi = (c, s) + REQUIRE_THAT(c * 2.0f + s * 1.0f, + Catch::Matchers::WithinRel(lHi * c, 1e-5f)); + REQUIRE_THAT(c * 1.0f + s * 3.0f, + Catch::Matchers::WithinRel(lHi * s, 1e-5f)); + // Second eigenvector vLo = (-s, c) + REQUIRE_THAT(-s * 2.0f + c * 1.0f, + Catch::Matchers::WithinRel(lLo * -s, 1e-5f)); + REQUIRE_THAT(-s * 1.0f + c * 3.0f, + Catch::Matchers::WithinRel(lLo * c, 1e-5f)); + } + + // Nonsymmetric block [[1,2],[3,4]] (used by the N == 2 entry point): + // eigenvalues 5.372281323, -0.372281323; + // eigenvector of 5.372281323 is +/- (0.415973558, 0.909376709) + { + Matrix<2, 2> A{1, 2, 3, 4}; + float lHi = 0, lLo = 0, c = 0, s = 0; + QR::Solve2x2Eigen(A, 0, lHi, lLo, c, s); + + REQUIRE_THAT(lHi, Catch::Matchers::WithinRel(5.372281323f, 1e-6f)); + REQUIRE_THAT(lLo, Catch::Matchers::WithinRel(-0.372281323f, 1e-6f)); + REQUIRE(matchesAbs(c, 0.415973558f)); + REQUIRE(matchesAbs(s, 0.909376709f)); + + // Both-row residual with vHi = (c, s): A v = l v + REQUIRE_THAT(c * 1.0f + s * 2.0f, + Catch::Matchers::WithinRel(lHi * c, 1e-5f)); + REQUIRE_THAT(c * 3.0f + s * 4.0f, + Catch::Matchers::WithinRel(lHi * s, 1e-5f)); + } + + // Diagonal blocks: eigenvectors are coordinate vectors + { + Matrix<2, 2> A{5, 0, 0, 2}; + float lHi = 0, lLo = 0, c = 0, s = 0; + QR::Solve2x2Eigen(A, 0, lHi, lLo, c, s); + REQUIRE_THAT(lHi, Catch::Matchers::WithinRel(5.0f, 1e-7f)); + REQUIRE_THAT(lLo, Catch::Matchers::WithinRel(2.0f, 1e-7f)); + REQUIRE_THAT(c, Catch::Matchers::WithinRel(1.0f, 1e-7f)); + REQUIRE_THAT(s, Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + + A = Matrix<2, 2>{2, 0, 0, 5}; + QR::Solve2x2Eigen(A, 0, lHi, lLo, c, s); + REQUIRE_THAT(lHi, Catch::Matchers::WithinRel(5.0f, 1e-7f)); + REQUIRE_THAT(lLo, Catch::Matchers::WithinRel(2.0f, 1e-7f)); + REQUIRE_THAT(c, Catch::Matchers::WithinAbs(0.0f, 1e-7f)); + REQUIRE_THAT(s, Catch::Matchers::WithinRel(1.0f, 1e-7f)); + } +} +// ============================================================================ +// TEST 6: Deflate +// ============================================================================ +TEST_CASE("QR Building Block: Deflate", "[Matrix][QR]") { + // subdiag[0] = 1e-9 <= 1e-6 * (|2| + |3|) = 5e-6 -> deflated + // subdiag[1] = 0.5 > 1e-6 * (|3| + |4|) = 7e-6 -> kept + { + Matrix<3, 3> A{2, 1e-9f, 0, 1e-9f, 3, 0.5f, 0, 0.5f, 4}; + QR::Deflate(A, 0, 2, 1e-6f); + + REQUIRE(A.Get(1, 0) == 0.0f); + REQUIRE(A.Get(0, 1) == 0.0f); + REQUIRE_THAT(A.Get(2, 1), Catch::Matchers::WithinRel(0.5f, 1e-7f)); + REQUIRE_THAT(A.Get(1, 2), Catch::Matchers::WithinRel(0.5f, 1e-7f)); + // Diagonals untouched + REQUIRE_THAT(A.Get(0, 0), Catch::Matchers::WithinRel(2.0f, 1e-7f)); + REQUIRE_THAT(A.Get(1, 1), Catch::Matchers::WithinRel(3.0f, 1e-7f)); + REQUIRE_THAT(A.Get(2, 2), Catch::Matchers::WithinRel(4.0f, 1e-7f)); + } + + // Nothing deflated when all subdiagonals are well above tolerance + { + Matrix<3, 3> A{2, 0.1f, 0, 0.1f, 3, 0.2f, 0, 0.2f, 4}; + QR::Deflate(A, 0, 2, 1e-6f); + REQUIRE_THAT(A.Get(1, 0), Catch::Matchers::WithinRel(0.1f, 1e-7f)); + REQUIRE_THAT(A.Get(2, 1), Catch::Matchers::WithinRel(0.2f, 1e-7f)); + } +} + +// ============================================================================ +// TEST 7: Tridiagonalize +// ============================================================================ +TEST_CASE("QR Building Block: Tridiagonalize", "[Matrix][QR]") { + // 4x4 symmetric with a full (0,3) corner coupling + { + Matrix<4, 4> A{2, 1, 0, 1, 1, 3, 1, 0, 0, 1, 4, 1, 1, 0, 1, 5}; + Matrix<4, 4> Aorig = A; + Matrix<4, 4> U{0}; + + QR::Tridiagonalize(A, U); + + // Off-tridiagonal entries must be zero up to float32 roundoff (the + // Givens zeroing cancels only in exact arithmetic; residuals are + // ~1e-7 for O(1) entries). + REQUIRE_THAT(A.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(A.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(A.Get(0, 3), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(A.Get(3, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(A.Get(1, 3), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(A.Get(3, 1), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + + // Symmetry preserved exactly + for (uint8_t i = 0; i < 4; i++) + for (uint8_t j = 0; j < 4; j++) + REQUIRE(A.Get(i, j) == A.Get(j, i)); + + // U must be orthogonal + REQUIRE(isOrthogonal(U)); + + // Reconstruction: U * A_tri * U^T == Aorig (absolute check for + // originally-zero entries: WithinRel has no absolute fallback there) + Matrix<4, 4> UAt{}; + U.Mult(A, UAt); + Matrix<4, 4> UAtU{}; + UAt.Mult(U.Transpose(), UAtU); + for (uint8_t i = 0; i < 4; i++) + for (uint8_t j = 0; j < 4; j++) { + float actual = UAtU.Get(i, j); + float expected = Aorig.Get(i, j); + if (fabsf(expected) < 1e-3f) + REQUIRE_THAT(actual, Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + else + REQUIRE_THAT(actual, + Catch::Matchers::WithinRel(expected, 1e-5f)); + } + + // Spectrum invariants match the original + { + float tr0 = Aorig.Get(0, 0) + Aorig.Get(1, 1) + Aorig.Get(2, 2) + + Aorig.Get(3, 3); + float tr1 = A.Get(0, 0) + A.Get(1, 1) + A.Get(2, 2) + A.Get(3, 3); + REQUIRE_THAT(tr1, Catch::Matchers::WithinRel(tr0, 1e-6f)); + REQUIRE_THAT(frob(A), Catch::Matchers::WithinRel(frob(Aorig), 1e-6f)); + } + + // Eigenvalues of the tridiagonal match the original (scipy reference): + // 6.0, 4.0, 3.0, 1.0 + { + Matrix<4, 1> vals{}; + Matrix<4, 4> vecs{}; + QR::EigenQR(A, vecs, vals, 10000, 1e-6f); + REQUIRE_THAT(vals[0][0], Catch::Matchers::WithinRel(6.0f, 1e-4f)); + REQUIRE_THAT(vals[1][0], Catch::Matchers::WithinRel(4.0f, 1e-4f)); + REQUIRE_THAT(vals[2][0], Catch::Matchers::WithinRel(3.0f, 1e-4f)); + REQUIRE_THAT(vals[3][0], Catch::Matchers::WithinRel(1.0f, 1e-4f)); + } + } + + // 5x5 symmetric + { + Matrix<5, 5> A{3, 1, 0, 0, 1, 1, 4, 1, 0, 0, 0, 1, 5, 1, 0, 0, 0, 1, 6, 1, + 1, 0, 0, 1, 7}; + Matrix<5, 5> Aorig = A; + Matrix<5, 5> U{0}; + + QR::Tridiagonalize(A, U); + + // All |i - j| >= 2 entries zero up to float32 roundoff + for (uint8_t i = 0; i < 5; i++) + for (uint8_t j = 0; j < 5; j++) + if (i > j + 1 || j > i + 1) + REQUIRE_THAT(A.Get(i, j), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + + REQUIRE(isOrthogonal(U)); + + Matrix<5, 5> UAt{}; + U.Mult(A, UAt); + Matrix<5, 5> UAtU{}; + UAt.Mult(U.Transpose(), UAtU); + for (uint8_t i = 0; i < 5; i++) + for (uint8_t j = 0; j < 5; j++) { + float actual = UAtU.Get(i, j); + float expected = Aorig.Get(i, j); + if (fabsf(expected) < 1e-3f) + REQUIRE_THAT(actual, Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + else + REQUIRE_THAT(actual, + Catch::Matchers::WithinRel(expected, 1e-5f)); + } + } + + // Already tridiagonal: U must come out as the identity + { + Matrix<3, 3> A{1, 2, 0, 2, 5, 2, 0, 2, 9}; + Matrix<3, 3> U{0}; + QR::Tridiagonalize(A, U); + for (uint8_t i = 0; i < 3; i++) + for (uint8_t j = 0; j < 3; j++) { + float expected = (i == j) ? 1.0f : 0.0f; + REQUIRE_THAT(U.Get(i, j), Catch::Matchers::WithinAbs(expected, 1e-7f)); + } + } +} +// ============================================================================ +// TEST 8: One full shifted QR step (integration of the blocks) +// ============================================================================ +TEST_CASE("QR Building Block: Full Shifted QR Step", "[Matrix][QR]") { + // One Wilkinson-shifted QR step on the whole 3x3 block is a similarity + // transform, so all spectrum invariants (trace, sum of principal 2x2 + // minors, determinant) must be preserved. + // + // A = [[1,2,3],[2,5,8],[3,8,9]]: tr = 15, e2 = -18, det = -4 + { + Matrix<3, 3> A{1, 2, 3, 2, 5, 8, 3, 8, 9}; + float tr0 = trace3(A); // 15 + float e20 = e2_3x3(A); // -18 + float det0 = det3(A); // -4 + + // mu from the trailing 2x2 [[5,8],[8,9]]: eigenvalues + // -1.246211251, 15.246211251; closest to d = 9 is 15.246211251 (Wilkinson) + float mu = QR::WilkinsonShift(A.Get(1, 1), A.Get(2, 1), A.Get(2, 2)); + REQUIRE_THAT(mu, Catch::Matchers::WithinRel(15.246211251f, 1e-5f)); + + for (uint8_t i = 0; i < 3; i++) + A[i][i] -= mu; + + // Bulge chase: rotations on (0,1) then (1,2) + float c = 0, s = 0; + QR::GivensRotation(A.Get(0, 0), A.Get(1, 0), c, s); + QR::ApplyRotationBothSides(A, 0, c, s); + QR::GivensRotation(A.Get(1, 1), A.Get(2, 1), c, s); + QR::ApplyRotationBothSides(A, 1, c, s); + + for (uint8_t i = 0; i < 3; i++) + A[i][i] += mu; + + // Symmetry preserved + for (uint8_t i = 0; i < 3; i++) + for (uint8_t j = 0; j < 3; j++) + REQUIRE(A.Get(i, j) == A.Get(j, i)); + + // Spectrum invariants preserved + REQUIRE_THAT(trace3(A), Catch::Matchers::WithinRel(tr0, 1e-5f)); + REQUIRE_THAT(e2_3x3(A), Catch::Matchers::WithinRel(e20, 1e-5f)); + REQUIRE_THAT(det3(A), Catch::Matchers::WithinRel(det0, 1e-5f)); + } + + // For TRIDIAGONAL input a single step keeps the tridiagonal structure + { + Matrix<3, 3> T{1, 2, 0, 2, 5, 2, 0, 2, 9}; + float mu = QR::WilkinsonShift(T.Get(1, 1), T.Get(2, 1), T.Get(2, 2)); + for (uint8_t i = 0; i < 3; i++) + T[i][i] -= mu; + float c = 0, s = 0; + QR::GivensRotation(T.Get(0, 0), T.Get(1, 0), c, s); + QR::ApplyRotationBothSides(T, 0, c, s); + QR::GivensRotation(T.Get(1, 1), T.Get(2, 1), c, s); + QR::ApplyRotationBothSides(T, 1, c, s); + for (uint8_t i = 0; i < 3; i++) + T[i][i] += mu; + + // Corners must vanish up to float32 roundoff: tridiagonal form + // maintained. The cancellation is exact in exact arithmetic (the + // corner is s1*a - c1*b times a factor, and Givens gives s1*a = c1*b), + // so the residual is pure rounding, ~1e-6 for O(1) entries. + REQUIRE_THAT(T.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + REQUIRE_THAT(T.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f)); + } +} diff --git a/unit-tests/qr-reference-values.py b/unit-tests/qr-reference-values.py new file mode 100644 index 0000000..8339763 --- /dev/null +++ b/unit-tests/qr-reference-values.py @@ -0,0 +1,246 @@ +#!/usr/bin/env python3 +""" +Reference values for the QR eigen-decomposition building block tests +(unit-tests/qr-build-blocks-tests.cpp). Run this to verify/implement the +C++ implementation in src/QR.hpp / src/QR.cpp against numpy/scipy. + +Conventions (match the C++ exactly): + * Givens zeroing rotation: G = [[c, s], [-s, c]], c = x/r, s = y/r, + r = hypot(x, y). G * (x, y)^T = (r, 0)^T. + * Similarity transform: A <- G A G^T (ApplyRotationBothSides). + * Eigenvector accumulation: V <- V G^T (ApplyRotationToVectors). + Vblock in the 2x2 closed form is [[c, -s], [s, c]] (same shape as G^T). + * Tridiagonalization: bottom-up Givens (i = N-2 down to k+1 per column k). + * Shifted QR loop: Wilkinson shift mu from the trailing 2x2, chase on the + trailing unreduced block [lo, hi], deflate by relative tolerance, peel + exact-zero subdiagonals, 2x2 closed-form termination. + * Pipeline: M0 = U * Mtri * U^T and Mtri = V * D * V^T => + eigenvectors of M0 = U * V (columns), eigenvalues = diag(D). + +Usage: python3 qr-reference-values.py +""" + +import numpy as np +import scipy.linalg as sla + +np.set_printoptions(precision=9, linewidth=120) + + +def givens(x, y): + """c = x/r, s = y/r with r = hypot(x, y).""" + r = np.hypot(x, y) + if r == 0.0: + return 1.0, 0.0 + return x / r, y / r + + +def rot(n, i, c, s): + """G = I with [[c, s], [-s, c]] embedded at (i, i+1).""" + G = np.eye(n) + G[i:i + 2, i:i + 2] = np.array([[c, s], [-s, c]]) + return G + + +def tridiagonalize(M0): + """Bottom-up Givens tridiagonalization. Returns (Mtri, U) with + M0 = U Mtri U^T.""" + n = len(M0) + M = M0.copy() + U = np.eye(n) + for k in range(n - 2): + for i in range(n - 2, k, -1): + c, s = givens(M[i, k], M[i + 1, k]) + G = rot(n, i, c, s) + M = G @ M @ G.T + U = U @ G.T + return M, U + + +def wilkinson(a, b, d): + """Eigenvalue of [[a, b], [b, d]] closest to d.""" + delta = 0.5 * (a - d) + spread = np.sqrt(delta * delta + b * b) + return 0.5 * (a + d) - (spread if delta >= 0 else -spread) + + +def solve2x2(A, lo): + """Closed form for the block at (lo, lo+1): (lHi, lLo, c, s) with + vHi = (c, s), vLo = (-s, c).""" + a = A[lo, lo] + b = A[lo, lo + 1] + e = A[lo + 1, lo] + d = A[lo + 1, lo + 1] + tr = a + d + det = a * d - b * e + disc = max(0.0, tr * tr - 4 * det) + lhi = 0.5 * (tr + np.sqrt(disc)) + llo = 0.5 * (tr - np.sqrt(disc)) + if b != 0.0: + v1 = lhi - a + nn = np.hypot(b, v1) + c, s = b / nn, v1 / nn + elif a >= d: + c, s = 1.0, 0.0 + else: + c, s = 0.0, 1.0 + return lhi, llo, c, s + + +def eigenqr(M0, tol=1e-12, max_iter=100000): + """Full pipeline mirroring QR::EigenQR. Returns (eigs, W) where W has + the eigenvectors of M0 as columns.""" + n = len(M0) + if n == 2: + l1, l2, c, s = solve2x2(M0, 0) + return np.array([l1, l2]), np.array([[c, -s], [s, c]]) + M, U = tridiagonalize(M0) + V = np.eye(n) + hi = n - 1 + for _ in range(max_iter): + # deflate: zero tiny subdiagonals (relative test) + for i in range(hi): + t = M[i + 1, i] + scale = abs(M[i, i]) + abs(M[i + 1, i + 1]) + if abs(t) <= tol * scale: + M[i + 1, i] = M[i, i + 1] = 0.0 + # peel exact-zero trailing subdiagonals + while hi > 0 and M[hi, hi - 1] == 0.0: + hi -= 1 + if hi == 0: + break + # find start of trailing unreduced block + lo = hi + for i in range(hi - 1, -1, -1): + if M[i + 1, i] == 0.0: + break + lo = i + if lo + 1 == hi: + # closed-form 2x2 termination: set diagonal, fold Vblock in + l1, l2, c, s = solve2x2(M, lo) + Vb = np.eye(n) + Vb[lo:lo + 2, lo:lo + 2] = np.array([[c, -s], [s, c]]) + V = V @ Vb + M[lo, lo] = l1 + M[lo + 1, lo + 1] = l2 + M[lo + 1, lo] = M[lo, lo + 1] = 0.0 + if lo == 0: + break + hi = lo - 1 + continue + # full shifted step on [lo, hi] (shift applies to the active block) + mu = wilkinson(M[hi - 1, hi - 1], M[hi, hi - 1], M[hi, hi]) + diag = M.diagonal().copy() + diag[lo:hi + 1] -= mu + np.fill_diagonal(M, diag) + c, s = givens(M[lo, lo], M[lo + 1, lo]) + G = rot(n, lo, c, s) + M = G @ M @ G.T + V = V @ G.T + for i in range(lo + 1, hi): + c, s = givens(M[i, i], M[i + 1, i]) + G = rot(n, i, c, s) + M = G @ M @ G.T + V = V @ G.T + diag = M.diagonal().copy() + diag[lo:hi + 1] += mu + np.fill_diagonal(M, diag) + + eigs = np.diag(M).astype(float) + order = np.argsort(eigs)[::-1] # descending, like the C++ test harness + eigs = eigs[order] + W = U @ V + W = W[:, order] + return eigs, W + + +def report(name, val, ref=None, tol=1e-6): + ok = "OK " if ref is None or np.allclose(val, ref, rtol=tol, atol=tol) else "FAIL" + print(f"[{ok}] {name} = {val}") + if ref is not None: + print(f" scipy/numpy ref = {ref}") + + +def main(): + print("=== TEST 1: GivensRotation ===") + c, s = givens(2.0, 1.0) + print(f" c = {c} s = {s}") + # G * (x, y)^T = (r, 0)^T: G = [[c, s], [-s, c]] + assert abs(c * 2 + s * 1 - np.sqrt(5)) < 1e-15 + assert abs(-s * 2 + c * 1) < 1e-15 + + print("\n=== TEST 2: ApplyRotationBothSides A <- G A G^T ===") + A = np.array([[3.0, 4.0, 5.0], [6.0, 7.0, 8.0], [9.0, 10.0, 11.0]]) + G = rot(3, 0, 0.6, 0.8) + B = G @ A @ G.T + print(B) + + A = np.array([[5.0, 0.0, 1.0], [0.0, 6.0, 2.0], [1.0, 2.0, 7.0]]) + c, s = givens(6.0, 2.0) + G = rot(3, 1, c, s) + B = G @ A @ G.T + print(B) + + print("\n=== TEST 3: V accumulation V <- V G^T ===") + V = np.eye(3) + G = rot(3, 0, 0.894427191, 0.447213595) + V = V @ G.T + print(V) + V2 = V @ rot(3, 1, 0.848874681, 0.528748047).T + print(V2) + + print("\n=== TEST 4: Solve2x2Eigen ===") + for A in (np.array([[5.0, 8.0], [8.0, 9.0]]), np.array([[1.0, 2.0], [3.0, 4.0]])): + l1, l2, c, s = solve2x2(A, 0) + ref = np.linalg.eigvalsh(A) if np.allclose(A, A.T) else np.linalg.eigvals(A) + print(f" A={A.ravel()} lHi={l1} lLo={l2} c={c} s={s} ref={np.sort(ref)[::-1]}") + + print("\n=== TEST 8: WilkinsonShift ===") + print(f" W(5, 8, 9) = {wilkinson(5, 8, 9)}") + print(f" W(4, 2, 7) = {wilkinson(4, 2, 7)}") + print(f" W(9, 2, 5) = {wilkinson(9, 2, 5)}") + + print("\n=== TEST 8b: one full shifted chase step on tridiagonal 3x3 ===") + T = np.array([[1.0, 2.0, 0.0], [2.0, 5.0, 2.0], [0.0, 2.0, 9.0]]) + mu = wilkinson(5, 2, 9) + M = T - mu * np.eye(3) + c, s = givens(M[0, 0], M[1, 0]) + M = rot(3, 0, c, s) @ M @ rot(3, 0, c, s).T + c, s = givens(M[1, 1], M[2, 1]) + M = rot(3, 1, c, s) @ M @ rot(3, 1, c, s).T + M = M + mu * np.eye(3) + print(f" mu = {mu}") + print(M) + print(f" corners: {M[0, 2]}, {M[2, 0]} (exact-arithmetic zeros)") + print(f" trace {M.trace():.15f} (was {T.trace()})") + + print("\n=== TEST 7: Tridiagonalize ===") + M4 = np.array([[2.0, 1, 0, 1], [1, 3, 1, 0], [0, 1, 4, 1], [1, 0, 1, 5]]) + M, U = tridiagonalize(M4) + print(" M4 tridiagonalized:\n", M) + print(f" reconstruction U M U^T == M4: {np.allclose(U @ M @ U.T, M4, atol=1e-9)}") + M5 = np.array([[3.0, 1, 0, 0, 1], [1, 4, 1, 0, 0], [0, 1, 5, 1, 0], + [0, 0, 1, 6, 1], [1, 0, 0, 1, 7]]) + M, U = tridiagonalize(M5) + print(" M5 tridiagonalized:\n", M) + print(f" reconstruction: {np.allclose(U @ M @ U.T, M5, atol=1e-9)}") + + print("\n=== End-to-end: random symmetric vs scipy.linalg.eigh ===") + rng = np.random.default_rng(12345) + worst = 0.0 + for n in range(3, 9): + M0 = rng.normal(size=(n, n)) + M0 = (M0 + M0.T) / 2 + eigs, W = eigenqr(M0.astype(float)) + ref = sla.eigh(M0) + e_err = np.max(np.abs(np.sort(eigs) - ref[0])) + resid = np.linalg.norm(W @ np.diag(eigs) @ W.T - M0) + ortho = np.linalg.norm(W.T @ W - np.eye(n)) + print(f" n={n}: eigs_err={e_err:.2e} resid={resid:.2e} ortho={ortho:.2e}") + worst = max(worst, e_err, resid, ortho) + print(f"\nworst over all n: {worst:.2e}") + assert worst < 1e-10, "end-to-end reference FAILED" + print("ALL REFERENCES OK") + + +if __name__ == "__main__": + main()