eig.cpp File Reference

Full symmetric eigendecomposition – cyclic Jacobi algorithm. More...

#include "linalg/eigen/jacobi_eig.hpp"
#include <algorithm>
#include <cmath>
#include <stdexcept>

Go to the source code of this file.

Namespaces
namespace	num

Functions
EigenResult	num::eig_sym (const Matrix &A_in, real tol, idx max_sweeps)

Detailed Description

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)

Definition in file eig.cpp.