rustyqlib/core/linalg/decomp/
eigen.rs

1//! Eigendecomposition of symmetric matrices by the cyclic Jacobi
2//! method — robust and accurate for the small matrices of correlation
3//! and covariance work (the PSD projection inside Higham's
4//! nearest-correlation algorithm runs on this).
5
6/// Eigendecomposition of a symmetric matrix: returns
7/// `(eigenvalues, eigenvectors)` with the eigenvectors in the columns
8/// (`A = V diag(vals) V^T`). Order is unspecified.
9pub fn symmetric_eigen(a: &[Vec<f64>]) -> (Vec<f64>, Vec<Vec<f64>>) {
10    let n = a.len();
11    let mut m = a.to_vec();
12    let mut v: Vec<Vec<f64>> = (0..n)
13        .map(|i| (0..n).map(|j| if i == j { 1.0 } else { 0.0 }).collect())
14        .collect();
15    for _sweep in 0..100 {
16        let off: f64 = (0..n)
17            .flat_map(|i| ((i + 1)..n).map(move |j| (i, j)))
18            .map(|(i, j)| m[i][j] * m[i][j])
19            .sum();
20        if off < 1e-24 {
21            break;
22        }
23        for p in 0..n {
24            for q in (p + 1)..n {
25                if m[p][q].abs() < 1e-300 {
26                    continue;
27                }
28                let theta = (m[q][q] - m[p][p]) / (2.0 * m[p][q]);
29                let t = theta.signum() / (theta.abs() + (theta * theta + 1.0).sqrt());
30                let c = 1.0 / (t * t + 1.0).sqrt();
31                let s = t * c;
32                for k in 0..n {
33                    let (mkp, mkq) = (m[k][p], m[k][q]);
34                    m[k][p] = c * mkp - s * mkq;
35                    m[k][q] = s * mkp + c * mkq;
36                }
37                for k in 0..n {
38                    let (mpk, mqk) = (m[p][k], m[q][k]);
39                    m[p][k] = c * mpk - s * mqk;
40                    m[q][k] = s * mpk + c * mqk;
41                }
42                for row in v.iter_mut() {
43                    let (vp, vq) = (row[p], row[q]);
44                    row[p] = c * vp - s * vq;
45                    row[q] = s * vp + c * vq;
46                }
47            }
48        }
49    }
50    ((0..n).map(|i| m[i][i]).collect(), v)
51}
52
53#[cfg(test)]
54mod tests {
55    use super::*;
56
57    #[test]
58    fn reproduces_known_eigenvalues() {
59        // [[2,1],[1,2]] has eigenvalues 1 and 3
60        let (mut vals, _) = symmetric_eigen(&[vec![2.0, 1.0], vec![1.0, 2.0]]);
61        vals.sort_by(f64::total_cmp);
62        assert!((vals[0] - 1.0).abs() < 1e-10 && (vals[1] - 3.0).abs() < 1e-10);
63    }
64
65    #[test]
66    fn reconstructs_the_matrix() {
67        let a = vec![
68            vec![3.0, 1.0, 0.5],
69            vec![1.0, 2.0, -0.4],
70            vec![0.5, -0.4, 1.5],
71        ];
72        let (vals, vecs) = symmetric_eigen(&a);
73        for i in 0..3 {
74            for j in 0..3 {
75                let recon: f64 = (0..3).map(|k| vecs[i][k] * vals[k] * vecs[j][k]).sum();
76                assert!((recon - a[i][j]).abs() < 1e-10, "[{i}][{j}]");
77            }
78        }
79    }
80}
rustyqlib/core/linalg/decomp/eigen.rs

rustyqlib/core/linalg/decomp/
eigen.rs
