eig.cpp

Go to the documentation of this file.

/// @file eigen/eig.cpp
/// @brief Full symmetric eigendecomposition  -- cyclic Jacobi algorithm
///
/// Reference: Golub & Van Loan "Matrix Computations" §8.4
///
/// Each sweep visits every off-diagonal pair (p,q) and applies a 2x2
/// similarity (Givens) rotation that zeros A[p,q].  The rotation angle theta
/// satisfies:
///
///   tan(2theta) = 2*A[p,q] / (A[q,q] - A[p,p])
///
/// which is solved as:
///   tau = (A[q,q] - A[p,p]) / (2*A[p,q])
///   t = sign(tau) / (|tau| + sqrt(1 + tau^2))   <- smaller root of t^2 + 2taut - 1 = 0
///   c = 1/sqrt(1+t^2),  s = c*t
///
/// After the rotation:
///   A'[p,p] = c^2*A[p,p] - 2cs*A[p,q] + s^2*A[q,q]
///   A'[q,q] = s^2*A[p,p] + 2cs*A[p,q] + c^2*A[q,q]
///   A'[p,q] = A'[q,p] = 0                          <- zeroed by construction
///   A'[r,p] = c*A[r,p] - s*A[r,q]   for r != p,q
///   A'[r,q] = s*A[r,p] + c*A[r,q]   for r != p,q
///
/// Eigenvectors accumulate: V' = V * G  (G is the Givens matrix in (p,q) plane)
 
#include "linalg/eigen/jacobi_eig.hpp"
#include <algorithm>
#include <cmath>
#include <stdexcept>
 
namespace num {
 
EigenResult eig_sym(const Matrix& A_in, real tol, idx max_sweeps) {
    if (A_in.rows() != A_in.cols())
        throw std::invalid_argument("eig_sym: matrix must be square");
 
    idx n = A_in.rows();
    Matrix A = A_in;         // working copy  -- will be diagonalised
    Matrix V(n, n, 0.0);    // eigenvectors  -- initialised to identity
    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) {
        // Off-diagonal Frobenius norm (convergence check)
        real off = 0;
        for (idx p = 0; p < n; ++p)
            for (idx q = p+1; q < n; ++q)
                off += A(p,q) * A(p,q);
 
        if (std::sqrt(2.0 * off) < tol) {
            converged = true;
            break;
        }
 
        // Cyclic sweep over all (p,q) pairs
        for (idx p = 0; p < n-1; ++p) {
            for (idx q = p+1; q < n; ++q) {
                real apq = A(p,q);
                if (std::abs(apq) < 1e-15) continue;  // already negligible
 
                real app = A(p,p), aqq = A(q,q);
 
                // Jacobi rotation angle
                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;
 
                // Update diagonal and zero (p,q)
                A(p,p) = c*c*app - 2.0*c*s*apq + s*s*aqq;
                A(q,q) = s*s*app + 2.0*c*s*apq + c*c*aqq;
                A(p,q) = A(q,p) = 0.0;
 
                // Update off-diagonal rows/cols
                for (idx r = 0; r < n; ++r) {
                    if (r == p || r == q) continue;
                    real arp = A(r,p), arq = A(r,q);
                    A(r,p) = A(p,r) = c*arp - s*arq;
                    A(r,q) = A(q,r) = s*arp + c*arq;
                }
 
                // Accumulate eigenvectors: V' = V * G
                // G[:,p] = c*e_p - s*e_q,  G[:,q] = s*e_p + c*e_q
                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;
    }
 
    // Extract eigenvalues from diagonal
    Vector values(n);
    for (idx i = 0; i < n; ++i) values[i] = A(i,i);
 
    // Sort eigenvalues ascending, permute eigenvector columns to match
    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