jacobi_eig.hpp

Go to the documentation of this file.

/// @file eigen/jacobi_eig.hpp
/// @brief Full symmetric eigendecomposition via cyclic Jacobi sweeps
///
/// Computes ALL eigenvalues and eigenvectors of a real symmetric matrix.
///
/// Algorithm  -- cyclic Jacobi:
///   Repeatedly apply plane (Givens) rotations in the (p,q) plane to zero
///   the off-diagonal element A[p,q].  Each rotation is a similarity transform
///   so eigenvalues are preserved.  Convergence is quadratic: after each full
///   sweep the sum of squared off-diagonal elements decreases by at least a
///   constant factor.
///
/// Why Jacobi instead of the implicit QR algorithm?
///   Jacobi is conceptually simpler and each rotation is the exact same
///   2x2 eigendecomposition students encounter first.  For dense n < 1000 the
///   performance is acceptable.  LAPACK uses the Francis QR algorithm (dsyev)
///   which is O(n^2) per iteration vs O(n^2) per sweep here, but harder to
///   explain.
///
/// Complexity: O(n^2) per sweep, O(n) sweeps typical -> O(n^3) total.
#pragma once
 
#include "core/matrix.hpp"
#include "core/vector.hpp"
#include "core/policy.hpp"
#include <cmath>
#include <algorithm>
#include <stdexcept>
 
namespace num {
 
/// @brief Full eigendecomposition result: A = V * diag(values) * V^T
struct EigenResult {
    Vector values;    ///< Eigenvalues in ascending order
    Matrix vectors;   ///< Eigenvectors as columns of an nxn orthogonal matrix
    idx    sweeps;    ///< Number of Jacobi sweeps performed
    bool   converged; ///< Whether off-diagonal norm fell below tol
};
 
/// @brief Full eigendecomposition of a real symmetric matrix.
///
/// The rotation accumulation loop is parallelised when backend == Backend::omp.
///
/// @param A          Real symmetric nxn matrix (upper/lower triangle used)
/// @param backend    Execution backend (Backend::omp parallelises the rotation loop)
/// @param tol        Convergence: stop when Frobenius norm of off-diagonal < tol
/// @param max_sweeps Safety cap on Jacobi sweeps (typically < 20 for n < 500)
/// @return EigenResult with eigenvalues in ascending order and corresponding
///         eigenvectors as columns of the V matrix (A = V diag(lambda) Vᵀ)
inline EigenResult eig_sym(const Matrix& A,
                            real tol = 1e-12, idx max_sweeps = 100,
                            Backend backend = Backend::seq) {
    if (A.rows() != A.cols())
        throw std::invalid_argument("eig_sym: matrix must be square");
 
    constexpr real rotation_tol = 1e-15;
    idx n = A.rows();
    Matrix Aw = A;
    Matrix V(n, n, 0.0);
    for (idx i = 0; i < n; ++i) V(i,i) = 1.0;
 
    idx sweeps = 0;
    bool converged = false;
 
    for (idx sweep = 0; sweep < max_sweeps; ++sweep) {
        real off = 0;
        for (idx p = 0; p < n; ++p)
            for (idx q = p+1; q < n; ++q)
                off += Aw(p,q) * Aw(p,q);
 
        if (std::sqrt(2.0 * off) < tol) { converged = true; break; }
 
        for (idx p = 0; p < n-1; ++p) {
            for (idx q = p+1; q < n; ++q) {
                real apq = Aw(p,q);
                if (std::abs(apq) < rotation_tol) continue;
 
                real app = Aw(p,p), aqq = Aw(q,q);
                real tau = (aqq - app) / (2.0 * apq);
                real t   = std::copysign(1.0, tau)
                           / (std::abs(tau) + std::sqrt(1.0 + tau*tau));
                real c   = 1.0 / std::sqrt(1.0 + t*t);
                real s   = c * t;
 
                Aw(p,p) = c*c*app - 2.0*c*s*apq + s*s*aqq;
                Aw(q,q) = s*s*app + 2.0*c*s*apq + c*c*aqq;
                Aw(p,q) = Aw(q,p) = 0.0;
 
                // Update off-diagonal rows/cols; parallel when omp backend
#ifdef NUMERICS_HAS_OMP
                #pragma omp parallel for schedule(static) if(backend == Backend::omp)
#endif
                for (idx r = 0; r < n; ++r) {
                    if (r == p || r == q) continue;
                    real arp = Aw(r,p), arq = Aw(r,q);
                    Aw(r,p) = Aw(p,r) = c*arp - s*arq;
                    Aw(r,q) = Aw(q,r) = s*arp + c*arq;
                }
#ifdef NUMERICS_HAS_OMP
                #pragma omp parallel for schedule(static) if(backend == Backend::omp)
#endif
                for (idx r = 0; r < n; ++r) {
                    real vrp = V(r,p), vrq = V(r,q);
                    V(r,p) = c*vrp - s*vrq;
                    V(r,q) = s*vrp + c*vrq;
                }
            }
        }
        ++sweeps;
    }
 
    Vector values(n);
    for (idx i = 0; i < n; ++i) values[i] = Aw(i,i);
 
    for (idx i = 0; i < n-1; ++i) {
        idx min_j = i;
        for (idx j = i+1; j < n; ++j)
            if (values[j] < values[min_j]) min_j = j;
        if (min_j != i) {
            std::swap(values[i], values[min_j]);
            for (idx r = 0; r < n; ++r)
                std::swap(V(r,i), V(r,min_j));
        }
    }
 
    return {values, V, sweeps, converged};
}
 
} // namespace num