eigenvalues/
utils.rs

1/*!
2
3## Auxiliar functions to manipulate arrays
4
5 */
6
7use nalgebra as na;
8
9use approx::relative_eq;
10use na::linalg::SymmetricEigen;
11use na::Dynamic;
12use na::{DMatrix, DVector};
13
14/// Generate a random highly diagonal symmetric matrix
15pub fn generate_diagonal_dominant(dim: usize, sparsity: f64) -> DMatrix<f64> {
16    let xs = 1..=dim;
17    let it = xs.map(|x: usize| x as f64);
18    let vs = DVector::<f64>::from_iterator(dim, it);
19    let mut arr = DMatrix::<f64>::new_random(dim, dim);
20    arr += &arr.transpose();
21    arr *= sparsity;
22    arr.set_diagonal(&vs);
23    arr
24}
25
26/// Random symmetric matrix
27pub fn generate_random_symmetric(dim: usize, magnitude: f64) -> DMatrix<f64> {
28    let arr = DMatrix::<f64>::new_random(dim, dim) * magnitude;
29    &arr * arr.transpose()
30}
31
32/// Random Sparse matrix
33pub fn generate_random_sparse_symmetric(dim: usize, lim: usize, sparsity: f64) -> DMatrix<f64> {
34    let arr = generate_diagonal_dominant(dim, sparsity);
35    let lambda = |i, j| {
36        if j > i + lim && i > j + lim {
37            0.0
38        } else {
39            arr[(i, j)]
40        }
41    };
42    DMatrix::<f64>::from_fn(dim, dim, lambda)
43}
44
45/// Sort the eigenvalues and their corresponding eigenvectors in ascending order
46pub fn sort_eigenpairs(
47    eig: SymmetricEigen<f64, Dynamic>,
48    ascending: bool,
49) -> SymmetricEigen<f64, Dynamic> {
50    // Sort the eigenvalues
51    let mut vs: Vec<(f64, usize)> = eig
52        .eigenvalues
53        .iter()
54        .enumerate()
55        .map(|(idx, &x)| (x, idx))
56        .collect();
57    sort_vector(&mut vs, ascending);
58
59    // Sorted eigenvalues
60    let eigenvalues = DVector::<f64>::from_iterator(vs.len(), vs.iter().map(|t| t.0));
61
62    // Indices of the sorted eigenvalues
63    let indices: Vec<_> = vs.iter().map(|t| t.1).collect();
64
65    // Create sorted eigenvectors
66    let dim_rows = eig.eigenvectors.nrows();
67    let dim_cols = eig.eigenvectors.ncols();
68    let mut eigenvectors = DMatrix::<f64>::zeros(dim_rows, dim_cols);
69
70    for (k, i) in indices.iter().enumerate() {
71        eigenvectors.set_column(k, &eig.eigenvectors.column(*i));
72    }
73    SymmetricEigen {
74        eigenvalues,
75        eigenvectors,
76    }
77}
78
79pub fn sort_vector<T: PartialOrd>(vs: &mut Vec<T>, ascending: bool) {
80    if ascending {
81        vs.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
82    } else {
83        vs.sort_unstable_by(|a, b| b.partial_cmp(a).unwrap());
84    }
85}
86
87pub fn test_eigenpairs(
88    reference: &na::linalg::SymmetricEigen<f64, na::Dynamic>,
89    eigenpair: (na::DVector<f64>, na::DMatrix<f64>),
90    number: usize,
91) {
92    let (dav_eigenvalues, dav_eigenvectors) = eigenpair;
93    for i in 0..number {
94        // Test Eigenvalues
95        assert!(relative_eq!(
96            reference.eigenvalues[i],
97            dav_eigenvalues[i],
98            epsilon = 1e-6
99        ));
100        // Test Eigenvectors
101        let x = reference.eigenvectors.column(i);
102        let y = dav_eigenvectors.column(i);
103        // The autovectors may different in their sign
104        // They should be either parallel or antiparallel
105        let dot = x.dot(&y).abs();
106        assert!(relative_eq!(dot, 1.0, epsilon = 1e-6));
107    }
108}
109
110#[cfg(test)]
111mod test {
112    use nalgebra as na;
113    use std::f64;
114
115    #[test]
116    fn test_random_symmetric() {
117        let matrix = super::generate_random_symmetric(10, 2.5);
118        test_symmetric(matrix);
119    }
120    #[test]
121    fn test_diagonal_dominant() {
122        let matrix = super::generate_diagonal_dominant(10, 0.005);
123        test_symmetric(matrix);
124    }
125
126    fn test_symmetric(matrix: na::DMatrix<f64>) {
127        let rs = &matrix - &matrix.transpose();
128        assert!(rs.sum() < f64::EPSILON);
129    }
130}
eigenvalues/utils.rs

eigenvalues/
utils.rs
