mathru/algebra/linear/matrix/uppertriangular/eigendec/
native.rs

1use crate::algebra::abstr::{AbsDiffEq, Field, Scalar};
2use crate::algebra::linear::matrix::Diagonal;
3use crate::algebra::linear::matrix::{EigenDec, EigenDecomposition, General, UpperTriangular};
4use crate::elementary::Power;
5
6impl<T> EigenDecomposition<T> for UpperTriangular<T>
7where
8    T: Field + Scalar + Power + AbsDiffEq<Epsilon = T>,
9{
10    /// Computes the eigen decomposition
11    ///
12    /// # Arguments
13    ///
14    /// # Example
15    ///
16    /// ```
17    /// use mathru::algebra::linear::matrix::{EigenDec, EigenDecomposition, General, UpperTriangular};
18    /// use mathru::matrix;
19    ///
20    /// let a: UpperTriangular<f64> = matrix![  -3.0, 3.0, 6.0;
21    ///                                         0.0, -5.0, -6.0;
22    ///                                         0.0, 0.0, 4.0].into();
23    ///
24    /// let eigen: EigenDec<f64> = a.dec_eigen().unwrap();
25    /// ```
26    fn dec_eigen(&self) -> Result<EigenDec<T>, String> {
27        let (m, _): (usize, usize) = self.dim();
28
29        let mut eigen_values = Vec::with_capacity(m);
30        for i in 0..m {
31            eigen_values.push(self[[i, i]]);
32        }
33        let values: Diagonal<T> = Diagonal::new(&eigen_values);
34        let vectors: General<T> = self.calc_eigenvector(&eigen_values);
35
36        Ok(EigenDec::new(values, vectors))
37    }
38}
39
40impl<T> UpperTriangular<T>
41where
42    T: Field + Scalar + Power + AbsDiffEq<Epsilon = T>,
43{
44    pub fn calc_eigenvector(&self, eigen_values: &[T]) -> General<T> {
45        let mut x: General<T> = General::zero(eigen_values.len(), eigen_values.len());
46
47        for (idx, lambda) in eigen_values.iter().enumerate().rev() {
48            x[[idx, idx]] = T::one();
49
50            for n in (0..idx).rev() {
51                let mut sum = T::zero();
52
53                for i in n + 1..eigen_values.len() {
54                    sum += self.matrix[[n, i]] * x[[i, idx]];
55                }
56
57                let divisor = *lambda - self.matrix[[n, n]];
58
59                x[[n, idx]] = sum / divisor;
60            }
61        }
62
63        UpperTriangular::norm(x)
64    }
65
66    fn norm(mut m: General<T>) -> General<T> {
67        let (rows, cols) = m.dim();
68        for c in 0..cols {
69            let norm = m.get_column(c).p_norm(&T::from_f32(2.0));
70
71            for r in 0..rows {
72                m[[r, c]] /= norm;
73            }
74        }
75
76        m
77    }
78}
mathru/algebra/linear/matrix/uppertriangular/eigendec/native.rs

mathru/algebra/linear/matrix/uppertriangular/eigendec/
native.rs
