From 48c43add5612a9e885edee5aad3c9e475d3a87ad Mon Sep 17 00:00:00 2001 From: Gabriel Santos Date: Sat, 1 Aug 2026 17:37:06 +0000 Subject: [PATCH 1/3] add jacobi eigen solver --- ROADMAP.md | 3 +- doc/solvers/JacobiEigenSolver.md | 133 ++++++++++++++ doc/solvers/README.md | 1 + numerical/math/GivensRotation.hpp | 15 ++ numerical/math/MatrixNorms.hpp | 11 ++ numerical/math/test/TestGivensRotation.cpp | 21 +++ numerical/math/test/TestMatrixNorms.cpp | 12 ++ numerical/solvers/CMakeLists.txt | 2 + numerical/solvers/JacobiEigenSolver.cpp | 8 + numerical/solvers/JacobiEigenSolver.hpp | 170 ++++++++++++++++++ numerical/solvers/test/CMakeLists.txt | 1 + .../solvers/test/TestJacobiEigenSolver.cpp | 120 +++++++++++++ 12 files changed, 495 insertions(+), 2 deletions(-) create mode 100644 doc/solvers/JacobiEigenSolver.md create mode 100644 numerical/solvers/JacobiEigenSolver.cpp create mode 100644 numerical/solvers/JacobiEigenSolver.hpp create mode 100644 numerical/solvers/test/TestJacobiEigenSolver.cpp diff --git a/ROADMAP.md b/ROADMAP.md index 5786f348..ff809977 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -28,7 +28,6 @@ Difficulty legend: | # | Component | Target module | Difficulty | |----|------------------------------------------------------|---------------------------|------------| | 39 | Square-root / Information Kalman filter | `filters/active` | ★★★★☆ | -| 40 | Feedback linearization | `nonlinear_control` (new) | ★★★★☆ | | 41 | Backstepping controller | `nonlinear_control` (new) | ★★★★☆ | | 42 | Symmetric eigenvalue solver (Jacobi) | `solvers` | ★★★★★ | | 43 | Singular Value Decomposition (Golub-Kahan) | `solvers` | ★★★★★ | @@ -237,7 +236,7 @@ the library does not yet expose. Detailed below under ## Tier 5 — Hard / research-grade ★★★★★ -### 42. Symmetric eigenvalue solver (Jacobi) *(float-first)* +### ~~42. Symmetric eigenvalue solver (Jacobi)~~ *(float-first)* ✓ Done - **What:** Cyclic Jacobi rotations for the eigenvalues/vectors of a symmetric matrix. - **Embedded value:** PCA/feature extraction, modal analysis, covariance conditioning, Gramian analysis. - **Algorithm / paper:** Golub & Van Loan, *Matrix Computations*, Ch. 8 (symmetric eigenproblem / cyclic Jacobi). diff --git a/doc/solvers/JacobiEigenSolver.md b/doc/solvers/JacobiEigenSolver.md new file mode 100644 index 00000000..cbbe63d5 --- /dev/null +++ b/doc/solvers/JacobiEigenSolver.md @@ -0,0 +1,133 @@ +# Jacobi Eigenvalue Solver + +## Overview & Motivation + +Many embedded estimation and analysis tasks reduce to the **symmetric eigenvalue problem**: given a real symmetric matrix $A$, find scalars $\lambda_k$ and orthonormal vectors $v_k$ such that $A v_k = \lambda_k v_k$. Covariance matrices, Gramians, inertia tensors, and modal-analysis stiffness matrices are all symmetric, and their eigen-decomposition exposes principal directions (PCA), energy distribution, conditioning, and vibration modes. + +The **cyclic Jacobi method** is the oldest and one of the most robust ways to solve this problem. It repeatedly applies planar rotations that annihilate one off-diagonal entry at a time, gradually driving $A$ toward a diagonal matrix whose entries are the eigenvalues. Its appeal for resource-constrained systems is threefold: it needs only a fixed-size working copy of the matrix (no heap), it produces the full set of eigenvalues **and** a fully orthonormal eigenvector basis, and it is exceptionally accurate for symmetric matrices — including the small relative eigenvalues that QR-based methods can lose. + +## Mathematical Theory + +### Jacobi Rotations + +A Jacobi rotation $G(p, q, \theta)$ is an orthogonal matrix equal to the identity except in the $2\times 2$ block on rows/columns $p$ and $q$: + +$$ +\begin{bmatrix} c & s \\ -s & c \end{bmatrix}, \qquad c = \cos\theta, \; s = \sin\theta . +$$ + +Applying the similarity transform $A' = G^\top A G$ leaves the spectrum unchanged (orthogonal similarity) but modifies only rows/columns $p$ and $q$. + +### Annihilating the Pivot + +The angle is chosen so that the off-diagonal entry $a_{pq}$ becomes zero. Requiring $a'_{pq} = 0$ gives + +$$ +(c^2 - s^2)\,a_{pq} + c\,s\,(a_{pp} - a_{qq}) = 0 . +$$ + +Defining + +$$ +\theta = \frac{a_{qq} - a_{pp}}{2\,a_{pq}}, \qquad +t = \frac{\operatorname{sign}(\theta)}{|\theta| + \sqrt{\theta^2 + 1}}, \qquad +c = \frac{1}{\sqrt{t^2 + 1}}, \qquad s = t\,c , +$$ + +selects the **smaller root** $t = \tan\theta$. Choosing the smaller rotation angle keeps the transformation close to the identity, which is what guarantees numerical stability and monotone convergence. + +### Convergence Measure + +Let $\operatorname{off}(A) = \sqrt{\sum_{i \neq j} a_{ij}^2}$ be the Frobenius norm of the off-diagonal part. Each rotation that zeroes $a_{pq}$ reduces $\operatorname{off}(A)^2$ by exactly $2\,a_{pq}^2$, since orthogonal similarity preserves the total Frobenius norm and the diagonal absorbs the removed mass. Thus $\operatorname{off}(A)$ decreases monotonically toward zero. + +### Cyclic Sweeps + +Rather than searching for the largest off-diagonal entry each step (classical Jacobi, $O(n^2)$ search per rotation), the **cyclic** variant sweeps every pair $(p, q)$ with $p < q$ in fixed row-major order. A full sweep touches all $n(n-1)/2$ pairs. Near convergence the method is **quadratically convergent**: the off-diagonal norm is roughly squared each sweep, so a handful of sweeps suffices even for ill-conditioned inputs. + +### Eigenvectors + +Accumulating the rotations $V \leftarrow V\,G$ starting from $V = I$ yields the orthogonal matrix whose columns are the eigenvectors. On termination $A$ is (numerically) diagonal with $\lambda_k = a_{kk}$, and $A \approx V \operatorname{diag}(\lambda) V^\top$. + +## Complexity Analysis + +| Case | Time | Space | Notes | +|---------|-------------|----------|----------------------------------------------------------------| +| Best | $O(n^3)$ | $O(n^2)$ | Already near-diagonal; one sweep to confirm convergence | +| Average | $O(n^3)$ | $O(n^2)$ | Typically 6–10 sweeps; each sweep costs $O(n^3)$ | +| Worst | $O(S n^3)$ | $O(n^2)$ | $S$ sweeps capped by a fixed maximum; each rotation is $O(n)$ | + +**Why $O(n^3)$ per sweep:** a sweep performs $n(n-1)/2 = O(n^2)$ rotations, and each rotation updates two rows and two columns at $O(n)$ cost. Space is a single $n\times n$ working copy plus the $n\times n$ eigenvector accumulator — both stack/static, no dynamic allocation. + +## Step-by-Step Walkthrough + +**Input:** + +$$ +A = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} +$$ + +**Step 1 — Pick pivot** $(p, q) = (0, 1)$, $a_{pq} = 1$. + +**Step 2 — Compute angle.** $\theta = (a_{qq} - a_{pp}) / (2 a_{pq}) = (2 - 2)/2 = 0$, so $t = 1$, $c = s = 1/\sqrt{2} \approx 0.7071$ (a $45^\circ$ rotation). + +**Step 3 — Rotate.** The updated diagonal is + +$$ +a'_{pp} = c^2 a_{pp} - 2 s c\, a_{pq} + s^2 a_{qq} = 1, \qquad +a'_{qq} = s^2 a_{pp} + 2 s c\, a_{pq} + c^2 a_{qq} = 3, +$$ + +and $a'_{pq} = 0$. The matrix is now diagonal. + +**Step 4 — Read results.** Eigenvalues $\{1, 3\}$; eigenvectors (columns of the accumulated rotation) $\tfrac{1}{\sqrt 2}(1, -1)^\top$ and $\tfrac{1}{\sqrt 2}(1, 1)^\top$. After sorting ascending, $\lambda = (1, 3)$. + +## Pitfalls & Edge Cases + +- **Symmetry is assumed.** Only the symmetric part is meaningful; a non-symmetric input silently has its lower triangle mirrored by the rotations. Callers must pass a genuinely symmetric matrix. +- **Zero pivot skip.** When $a_{pq}$ is already zero the rotation is skipped, avoiding a division by zero in the $\theta$ formula. +- **Degenerate / repeated eigenvalues.** Convergence is unaffected, but the eigenvectors within a degenerate subspace are only determined up to rotation — any orthonormal basis of that subspace is valid. +- **Termination threshold.** Convergence is declared when the off-diagonal norm falls below a small multiple of the matrix scale (its Frobenius/diagonal magnitude). A fixed maximum sweep count bounds worst-case runtime; failure to converge within it is reported to the caller rather than looping forever. +- **Float precision.** In single precision the achievable off-diagonal residual is limited by machine epsilon times the largest eigenvalue; tolerances on reconstruction should be set accordingly. + +## Variants & Generalizations + +| Variant | Key Difference | +|--------------------------------|--------------------------------------------------------------------------------------------------| +| **Classical Jacobi** | Zeroes the *largest* off-diagonal entry each step; fewer rotations but an $O(n^2)$ search each time | +| **Threshold Jacobi** | Skips pivots below a per-sweep threshold, cheaper early sweeps on sparse-ish matrices | +| **One-sided Jacobi** | Applies rotations to a factor only; the basis of Jacobi SVD for tall matrices | +| **QR / tridiagonal method** | Reduces to tridiagonal form then iterates — faster asymptotically but less accurate on tiny eigenvalues | + +## Applications + +- **Principal Component Analysis** — eigen-decompose a covariance matrix to obtain principal directions and variances. +- **Modal analysis** — natural frequencies and mode shapes from mass/stiffness matrices. +- **Covariance conditioning** — inspect or floor eigenvalues to keep estimator covariances positive-definite. +- **Inertia and orientation** — principal axes of an inertia tensor from its eigenvectors. +- **Gramian analysis** — controllability/observability energy directions in control systems. + +## Connections to Other Algorithms + +```mermaid +graph LR + JAC["Jacobi Eigen"] + QR["QR Decomposition"] + CHOL["Cholesky"] + SVD["SVD (Golub-Kahan)"] + QR -.->|"shares Givens rotations"| JAC + JAC -->|"symmetric building block"| SVD + JAC -.->|"eigenvalue floor keeps SPD"| CHOL +``` + +| Algorithm | Relationship | +|----------------------------------------------|------------------------------------------------------------------------------------------------| +| [QR Decomposition](QrDecomposition.md) | Both are built from orthogonal (Givens/Householder) transforms; QR underlies the alternative tridiagonal eigen-method | +| [Cholesky Decomposition](CholeskyDecomposition.md) | Requires symmetric positive-definite input; Jacobi eigenvalues certify or restore definiteness | +| [Spectral Radius](SpectralRadius.md) | Returns only the dominant eigenvalue magnitude; Jacobi returns the full spectrum and vectors | + +## References & Further Reading + +- Jacobi, C.G.J., "Über ein leichtes Verfahren, die in der Theorie der Säcularstörungen vorkommenden Gleichungen numerisch aufzulösen", *Crelle's Journal*, 30, 1846. +- Golub, G.H. & Van Loan, C.F., *Matrix Computations*, 4th ed., Johns Hopkins University Press, 2013 — Chapter 8 (symmetric eigenproblem, cyclic Jacobi). +- Press, W.H. et al., *Numerical Recipes*, 3rd ed., Cambridge University Press, 2007 — Section 11.1 (Jacobi transformations of a symmetric matrix). +- Demmel, J. & Veselić, K., "Jacobi's method is more accurate than QR", *SIAM J. Matrix Anal. Appl.*, 13(4), 1992. diff --git a/doc/solvers/README.md b/doc/solvers/README.md index e2ef999c..daebab6a 100644 --- a/doc/solvers/README.md +++ b/doc/solvers/README.md @@ -17,3 +17,4 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations. | [QR Decomposition](QrDecomposition.md) | Householder factorization and Givens streaming row update for least-squares solves | | [LU Decomposition](LuDecomposition.md) | PA = LU factorization with partial pivoting for general dense linear systems | | [Lyapunov / Sylvester Solvers](LyapunovSylvester.md) | Sylvester AX+XB=C and continuous/discrete Lyapunov solvers via Kronecker vectorisation | +| [Jacobi Eigenvalue Solver](JacobiEigenSolver.md) | Cyclic Jacobi rotations for the full symmetric eigenvalue/eigenvector problem | diff --git a/numerical/math/GivensRotation.hpp b/numerical/math/GivensRotation.hpp index 3a538a06..83f89462 100644 --- a/numerical/math/GivensRotation.hpp +++ b/numerical/math/GivensRotation.hpp @@ -29,6 +29,21 @@ namespace math return { a / r, b / r }; } + template + [[nodiscard]] OPTIMIZE_FOR_SPEED GivensRotation ComputeJacobiRotation(T app, T aqq, T apq) + { + static_assert(std::is_floating_point_v, "ComputeJacobiRotation supports floating-point types only"); + + if (apq == T{}) + return { T{ 1 }, T{} }; + + T theta = (aqq - app) / (T{ 2 } * apq); + T t = ((theta >= T{}) ? T{ 1 } : T{ -1 }) / (std::abs(theta) + std::sqrt(theta * theta + T{ 1 })); + T c = T{ 1 } / std::sqrt(t * t + T{ 1 }); + + return { c, t * c }; + } + template OPTIMIZE_FOR_SPEED void ApplyGivens(const GivensRotation& g, T& x, T& y) { diff --git a/numerical/math/MatrixNorms.hpp b/numerical/math/MatrixNorms.hpp index 02b53b7e..405902fd 100644 --- a/numerical/math/MatrixNorms.hpp +++ b/numerical/math/MatrixNorms.hpp @@ -22,6 +22,17 @@ namespace math return std::sqrt(sum); } + template + [[nodiscard]] OPTIMIZE_FOR_SPEED T OffDiagonalFrobeniusNorm(const Matrix& a) + { + static_assert(std::is_floating_point_v, "MatrixNorms supports floating-point types"); + T sum{}; + for (std::size_t i = 0; i < N; ++i) + for (std::size_t j = i + 1; j < N; ++j) + sum += a.at(i, j) * a.at(i, j) + a.at(j, i) * a.at(j, i); + return std::sqrt(sum); + } + template [[nodiscard]] OPTIMIZE_FOR_SPEED T OneNorm(const Matrix& a) { diff --git a/numerical/math/test/TestGivensRotation.cpp b/numerical/math/test/TestGivensRotation.cpp index b5d18634..b5623c31 100644 --- a/numerical/math/test/TestGivensRotation.cpp +++ b/numerical/math/test/TestGivensRotation.cpp @@ -41,3 +41,24 @@ TEST_F(GivensRotationTest, DegenerateInputIsIdentity) EXPECT_NEAR(g.c, 1.0f, math::Tolerance()); EXPECT_NEAR(g.s, 0.0f, math::Tolerance()); } + +TEST_F(GivensRotationTest, JacobiRotationAnnihilatesSymmetricOffDiagonal) +{ + float app = 4.0f; + float aqq = 1.0f; + float apq = 2.0f; + + auto g = math::ComputeJacobiRotation(app, aqq, apq); + float offDiagonal = (g.c * g.c - g.s * g.s) * apq + g.c * g.s * (app - aqq); + + EXPECT_NEAR(offDiagonal, 0.0f, math::Tolerance()); + EXPECT_NEAR(g.c * g.c + g.s * g.s, 1.0f, math::Tolerance()); +} + +TEST_F(GivensRotationTest, JacobiRotationZeroOffDiagonalIsIdentity) +{ + auto g = math::ComputeJacobiRotation(3.0f, 5.0f, 0.0f); + + EXPECT_NEAR(g.c, 1.0f, math::Tolerance()); + EXPECT_NEAR(g.s, 0.0f, math::Tolerance()); +} diff --git a/numerical/math/test/TestMatrixNorms.cpp b/numerical/math/test/TestMatrixNorms.cpp index c609f7a9..27b428ff 100644 --- a/numerical/math/test/TestMatrixNorms.cpp +++ b/numerical/math/test/TestMatrixNorms.cpp @@ -21,6 +21,18 @@ TEST_F(MatrixNormsTest, FrobeniusNorm) EXPECT_NEAR(result, 3.87298f, math::Tolerance()); } +TEST_F(MatrixNormsTest, OffDiagonalFrobeniusNorm) +{ + math::Matrix m{ + { 5.0f, 3.0f, 0.0f }, + { -3.0f, 7.0f, 4.0f }, + { 0.0f, -4.0f, 9.0f } + }; + + float result = math::OffDiagonalFrobeniusNorm(m); + EXPECT_NEAR(result, 7.0710678f, math::Tolerance()); +} + TEST_F(MatrixNormsTest, OneNorm) { float result = math::OneNorm(a); diff --git a/numerical/solvers/CMakeLists.txt b/numerical/solvers/CMakeLists.txt index 6abc38d9..0d6d60af 100644 --- a/numerical/solvers/CMakeLists.txt +++ b/numerical/solvers/CMakeLists.txt @@ -17,6 +17,7 @@ target_sources(numerical.solver PRIVATE DormandPrince45.hpp DurandKerner.hpp GaussianElimination.hpp + JacobiEigenSolver.hpp LevinsonDurbin.hpp LuDecomposition.hpp LyapunovSylvester.hpp @@ -31,6 +32,7 @@ numerical_add_coverage_sources(numerical.solver DiscreteAlgebraicRiccatiEquation.cpp DurandKerner.cpp GaussianElimination.cpp + JacobiEigenSolver.cpp LuDecomposition.cpp LyapunovSylvester.cpp QrDecomposition.cpp diff --git a/numerical/solvers/JacobiEigenSolver.cpp b/numerical/solvers/JacobiEigenSolver.cpp new file mode 100644 index 00000000..076ac337 --- /dev/null +++ b/numerical/solvers/JacobiEigenSolver.cpp @@ -0,0 +1,8 @@ +#include "numerical/solvers/JacobiEigenSolver.hpp" + +namespace solvers +{ + template class JacobiEigenSolver; + template class JacobiEigenSolver; + template class JacobiEigenSolver; +} diff --git a/numerical/solvers/JacobiEigenSolver.hpp b/numerical/solvers/JacobiEigenSolver.hpp new file mode 100644 index 00000000..ee9e9d05 --- /dev/null +++ b/numerical/solvers/JacobiEigenSolver.hpp @@ -0,0 +1,170 @@ +#pragma once + +#if defined(__GNUC__) || defined(__clang__) +#pragma GCC optimize("O3", "fast-math") +#endif + +#include "numerical/math/CompilerOptimizations.hpp" +#include "numerical/math/GivensRotation.hpp" +#include "numerical/math/Matrix.hpp" +#include "numerical/math/MatrixNorms.hpp" +#include +#include +#include + +namespace solvers +{ + template + class JacobiEigenSolver + { + static_assert(std::is_floating_point_v, "JacobiEigenSolver supports floating-point types only"); + static_assert(N >= 1, "JacobiEigenSolver requires N >= 1"); + + public: + JacobiEigenSolver() = default; + + OPTIMIZE_FOR_SPEED bool Solve(const math::Matrix& a); + + const math::Vector& Eigenvalues() const; + const math::Matrix& Eigenvectors() const; + std::size_t Sweeps() const; + + private: + static void ApplyPivotBlock(math::Matrix& a, std::size_t p, std::size_t q, const math::GivensRotation& rot); + void Rotate(math::Matrix& a, std::size_t p, std::size_t q); + void SortAscending(); + + math::Vector eigenvalues{}; + math::Matrix eigenvectors{}; + std::size_t sweeps{}; + bool solved{ false }; + + static constexpr std::size_t maxSweeps{ 50 }; + }; + + template + void JacobiEigenSolver::ApplyPivotBlock( + math::Matrix& a, std::size_t p, std::size_t q, const math::GivensRotation& rot) + { + T c = rot.c; + T s = rot.s; + T app = a.at(p, p); + T aqq = a.at(q, q); + T apq = a.at(p, q); + + a.at(p, p) = c * c * app - T{ 2 } * s * c * apq + s * s * aqq; + a.at(q, q) = s * s * app + T{ 2 } * s * c * apq + c * c * aqq; + a.at(p, q) = T{}; + a.at(q, p) = T{}; + } + + template + OPTIMIZE_FOR_SPEED void JacobiEigenSolver::Rotate(math::Matrix& a, std::size_t p, std::size_t q) + { + if (a.at(p, q) == T{}) + return; + + math::GivensRotation rot = math::ComputeJacobiRotation(a.at(p, p), a.at(q, q), a.at(p, q)); + ApplyPivotBlock(a, p, q, rot); + + math::GivensRotation transposed{ rot.c, -rot.s }; + + for (std::size_t i = 0; i < N; ++i) + { + if (i != p && i != q) + { + math::ApplyGivens(transposed, a.at(i, p), a.at(i, q)); + a.at(p, i) = a.at(i, p); + a.at(q, i) = a.at(i, q); + } + + math::ApplyGivens(transposed, eigenvectors.at(i, p), eigenvectors.at(i, q)); + } + } + + template + void JacobiEigenSolver::SortAscending() + { + for (std::size_t i = 0; i + 1 < N; ++i) + { + std::size_t minIndex = i; + for (std::size_t j = i + 1; j < N; ++j) + if (eigenvalues.at(j, 0) < eigenvalues.at(minIndex, 0)) + minIndex = j; + + if (minIndex != i) + { + T tmp = eigenvalues.at(i, 0); + eigenvalues.at(i, 0) = eigenvalues.at(minIndex, 0); + eigenvalues.at(minIndex, 0) = tmp; + + for (std::size_t r = 0; r < N; ++r) + { + T v = eigenvectors.at(r, i); + eigenvectors.at(r, i) = eigenvectors.at(r, minIndex); + eigenvectors.at(r, minIndex) = v; + } + } + } + } + + template + OPTIMIZE_FOR_SPEED bool JacobiEigenSolver::Solve(const math::Matrix& a) + { + math::Matrix work = a; + eigenvectors = math::Matrix::Identity(); + sweeps = 0; + solved = false; + + T scale = math::OffDiagonalFrobeniusNorm(a); + for (std::size_t i = 0; i < N; ++i) + scale += std::abs(a.at(i, i)); + + T threshold = ((scale > T{}) ? scale : T{ 1 }) * T{ 1e-7f }; + + while (sweeps < maxSweeps) + { + if (math::OffDiagonalFrobeniusNorm(work) <= threshold) + solved = true; + + if (solved) + break; + + for (std::size_t p = 0; p < N; ++p) + for (std::size_t q = p + 1; q < N; ++q) + Rotate(work, p, q); + + ++sweeps; + } + + for (std::size_t i = 0; i < N; ++i) + eigenvalues.at(i, 0) = work.at(i, i); + + SortAscending(); + return solved; + } + + template + const math::Vector& JacobiEigenSolver::Eigenvalues() const + { + return eigenvalues; + } + + template + const math::Matrix& JacobiEigenSolver::Eigenvectors() const + { + return eigenvectors; + } + + template + std::size_t JacobiEigenSolver::Sweeps() const + { + return sweeps; + } + +#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD + extern template class JacobiEigenSolver; + extern template class JacobiEigenSolver; + extern template class JacobiEigenSolver; +#endif +} diff --git a/numerical/solvers/test/CMakeLists.txt b/numerical/solvers/test/CMakeLists.txt index 2a2a029b..7271a521 100644 --- a/numerical/solvers/test/CMakeLists.txt +++ b/numerical/solvers/test/CMakeLists.txt @@ -12,6 +12,7 @@ target_sources(numerical.solvers_test PRIVATE TestDiscreteAlgebraicRiccatiEquation.cpp TestDurandKerner.cpp TestGaussianElimination.cpp + TestJacobiEigenSolver.cpp TestLevinsonDurbin.cpp TestLuDecomposition.cpp TestLyapunovSylvester.cpp diff --git a/numerical/solvers/test/TestJacobiEigenSolver.cpp b/numerical/solvers/test/TestJacobiEigenSolver.cpp new file mode 100644 index 00000000..3be1cc8a --- /dev/null +++ b/numerical/solvers/test/TestJacobiEigenSolver.cpp @@ -0,0 +1,120 @@ +#include "numerical/math/Tolerance.hpp" +#include "numerical/solvers/JacobiEigenSolver.hpp" +#include + +namespace +{ + class TestJacobiEigenSolver : public ::testing::Test + { + protected: + solvers::JacobiEigenSolver solver2{}; + solvers::JacobiEigenSolver solver3{}; + }; +} + +TEST_F(TestJacobiEigenSolver, diagonal_matrix_returns_sorted_diagonal) +{ + math::Matrix a{ + { 3.0f, 0.0f, 0.0f }, + { 0.0f, 1.0f, 0.0f }, + { 0.0f, 0.0f, 2.0f } + }; + + EXPECT_TRUE(solver3.Solve(a)); + + EXPECT_NEAR(solver3.Eigenvalues().at(0, 0), 1.0f, 1e-5f); + EXPECT_NEAR(solver3.Eigenvalues().at(1, 0), 2.0f, 1e-5f); + EXPECT_NEAR(solver3.Eigenvalues().at(2, 0), 3.0f, 1e-5f); +} + +TEST_F(TestJacobiEigenSolver, symmetric_2x2_matches_closed_form) +{ + math::Matrix a{ + { 2.0f, 1.0f }, + { 1.0f, 2.0f } + }; + + EXPECT_TRUE(solver2.Solve(a)); + + EXPECT_NEAR(solver2.Eigenvalues().at(0, 0), 1.0f, 1e-5f); + EXPECT_NEAR(solver2.Eigenvalues().at(1, 0), 3.0f, 1e-5f); +} + +TEST_F(TestJacobiEigenSolver, recovers_known_eigenvalues_of_dense_symmetric) +{ + math::Matrix a{ + { 4.0f, 1.0f, 1.0f }, + { 1.0f, 4.0f, 1.0f }, + { 1.0f, 1.0f, 4.0f } + }; + + EXPECT_TRUE(solver3.Solve(a)); + + EXPECT_NEAR(solver3.Eigenvalues().at(0, 0), 3.0f, 1e-5f); + EXPECT_NEAR(solver3.Eigenvalues().at(1, 0), 3.0f, 1e-5f); + EXPECT_NEAR(solver3.Eigenvalues().at(2, 0), 6.0f, 1e-5f); +} + +TEST_F(TestJacobiEigenSolver, eigenvectors_are_orthonormal) +{ + math::Matrix a{ + { 2.0f, -1.0f, 0.0f }, + { -1.0f, 2.0f, -1.0f }, + { 0.0f, -1.0f, 2.0f } + }; + + solver3.Solve(a); + auto v = solver3.Eigenvectors(); + auto vtv = v.Transpose() * v; + + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR(vtv.at(i, j), (i == j) ? 1.0f : 0.0f, 1e-5f); +} + +TEST_F(TestJacobiEigenSolver, satisfies_eigen_equation) +{ + math::Matrix a{ + { 2.0f, -1.0f, 0.0f }, + { -1.0f, 2.0f, -1.0f }, + { 0.0f, -1.0f, 2.0f } + }; + + solver3.Solve(a); + auto v = solver3.Eigenvectors(); + + for (std::size_t k = 0; k < 3; ++k) + { + float lambda = solver3.Eigenvalues().at(k, 0); + for (std::size_t i = 0; i < 3; ++i) + { + float av{ 0.0f }; + for (std::size_t j = 0; j < 3; ++j) + av += a.at(i, j) * v.at(j, k); + + EXPECT_NEAR(av, lambda * v.at(i, k), 1e-5f); + } + } +} + +TEST_F(TestJacobiEigenSolver, reconstructs_matrix_from_spectral_decomposition) +{ + math::Matrix a{ + { 6.0f, 2.0f, 1.0f }, + { 2.0f, 3.0f, 1.0f }, + { 1.0f, 1.0f, 1.0f } + }; + + solver3.Solve(a); + auto v = solver3.Eigenvectors(); + + math::Matrix d{}; + for (std::size_t i = 0; i < 3; ++i) + d.at(i, i) = solver3.Eigenvalues().at(i, 0); + + auto reconstructed = v * d * v.Transpose(); + + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR(reconstructed.at(i, j), a.at(i, j), 1e-4f); +} From 719fbe2d8ce259b8de80b2008c621d7d600a3fb4 Mon Sep 17 00:00:00 2001 From: gfs Date: Sat, 1 Aug 2026 19:41:39 +0200 Subject: [PATCH 2/3] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- doc/solvers/JacobiEigenSolver.md | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/doc/solvers/JacobiEigenSolver.md b/doc/solvers/JacobiEigenSolver.md index cbbe63d5..24da218d 100644 --- a/doc/solvers/JacobiEigenSolver.md +++ b/doc/solvers/JacobiEigenSolver.md @@ -50,11 +50,11 @@ Accumulating the rotations $V \leftarrow V\,G$ starting from $V = I$ yields the ## Complexity Analysis -| Case | Time | Space | Notes | -|---------|-------------|----------|----------------------------------------------------------------| -| Best | $O(n^3)$ | $O(n^2)$ | Already near-diagonal; one sweep to confirm convergence | -| Average | $O(n^3)$ | $O(n^2)$ | Typically 6–10 sweeps; each sweep costs $O(n^3)$ | -| Worst | $O(S n^3)$ | $O(n^2)$ | $S$ sweeps capped by a fixed maximum; each rotation is $O(n)$ | +| Case | Time | Space | Notes | +|---------|------------|----------|---------------------------------------------------------------| +| Best | $O(n^3)$ | $O(n^2)$ | Already near-diagonal; one sweep to confirm convergence | +| Average | $O(n^3)$ | $O(n^2)$ | Typically 6–10 sweeps; each sweep costs $O(n^3)$ | +| Worst | $O(S n^3)$ | $O(n^2)$ | $S$ sweeps capped by a fixed maximum; each rotation is $O(n)$ | **Why $O(n^3)$ per sweep:** a sweep performs $n(n-1)/2 = O(n^2)$ rotations, and each rotation updates two rows and two columns at $O(n)$ cost. Space is a single $n\times n$ working copy plus the $n\times n$ eigenvector accumulator — both stack/static, no dynamic allocation. @@ -91,12 +91,12 @@ and $a'_{pq} = 0$. The matrix is now diagonal. ## Variants & Generalizations -| Variant | Key Difference | -|--------------------------------|--------------------------------------------------------------------------------------------------| -| **Classical Jacobi** | Zeroes the *largest* off-diagonal entry each step; fewer rotations but an $O(n^2)$ search each time | -| **Threshold Jacobi** | Skips pivots below a per-sweep threshold, cheaper early sweeps on sparse-ish matrices | -| **One-sided Jacobi** | Applies rotations to a factor only; the basis of Jacobi SVD for tall matrices | -| **QR / tridiagonal method** | Reduces to tridiagonal form then iterates — faster asymptotically but less accurate on tiny eigenvalues | +| Variant | Key Difference | +|-----------------------------|---------------------------------------------------------------------------------------------------------| +| **Classical Jacobi** | Zeroes the *largest* off-diagonal entry each step; fewer rotations but an $O(n^2)$ search each time | +| **Threshold Jacobi** | Skips pivots below a per-sweep threshold, cheaper early sweeps on sparse-ish matrices | +| **One-sided Jacobi** | Applies rotations to a factor only; the basis of Jacobi SVD for tall matrices | +| **QR / tridiagonal method** | Reduces to tridiagonal form then iterates — faster asymptotically but less accurate on tiny eigenvalues | ## Applications From b6e54672a22308b1eece606306549717ee499cdb Mon Sep 17 00:00:00 2001 From: gfs Date: Sat, 1 Aug 2026 19:41:59 +0200 Subject: [PATCH 3/3] Update doc/solvers/JacobiEigenSolver.md Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- doc/solvers/JacobiEigenSolver.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/solvers/JacobiEigenSolver.md b/doc/solvers/JacobiEigenSolver.md index 24da218d..b53e4e95 100644 --- a/doc/solvers/JacobiEigenSolver.md +++ b/doc/solvers/JacobiEigenSolver.md @@ -119,11 +119,11 @@ graph LR JAC -.->|"eigenvalue floor keeps SPD"| CHOL ``` -| Algorithm | Relationship | -|----------------------------------------------|------------------------------------------------------------------------------------------------| -| [QR Decomposition](QrDecomposition.md) | Both are built from orthogonal (Givens/Householder) transforms; QR underlies the alternative tridiagonal eigen-method | -| [Cholesky Decomposition](CholeskyDecomposition.md) | Requires symmetric positive-definite input; Jacobi eigenvalues certify or restore definiteness | -| [Spectral Radius](SpectralRadius.md) | Returns only the dominant eigenvalue magnitude; Jacobi returns the full spectrum and vectors | +| Algorithm | Relationship | +|----------------------------------------------------|-----------------------------------------------------------------------------------------------------------------------| +| [QR Decomposition](QrDecomposition.md) | Both are built from orthogonal (Givens/Householder) transforms; QR underlies the alternative tridiagonal eigen-method | +| [Cholesky Decomposition](CholeskyDecomposition.md) | Requires symmetric positive-definite input; Jacobi eigenvalues certify or restore definiteness | +| [Spectral Radius](SpectralRadius.md) | Returns only the dominant eigenvalue magnitude; Jacobi returns the full spectrum and vectors | ## References & Further Reading