arr_rs/linalg/operations/
eigen.rs

1use crate::{
2    core::prelude::*,
3    errors::prelude::*,
4    linalg::prelude::*,
5    numeric::prelude::*,
6    validators::prelude::*,
7};
8
9/// `ArrayTrait` - Array Linalg Eigen functions
10pub trait ArrayLinalgEigen<N: NumericOps> where Self: Sized + Clone {
11
12    /// Compute the eigenvalues of a square array
13    ///
14    /// # Examples
15    ///
16    /// ```
17    /// use arr_rs::prelude::*;
18    ///
19    /// let array = Array::new(vec![12., -51., 4., 6., 167., -68., -4., 24., -41.], vec![3, 3]).unwrap();
20    /// let vals = array.eigvals().unwrap()[0].clone();
21    /// assert_eq!(Array::flat(vec![156.20350762022625, -35.473244608879945, 17.26973698865371]).unwrap(), vals);
22    /// ```
23    ///
24    /// # Errors
25    ///
26    /// may returns `ArrayError`
27    fn eigvals(&self) -> Result<Vec<Array<N>>, ArrayError>;
28
29    /// Compute the eigenvalues and right eigenvectors of a square array
30    ///
31    /// # Examples
32    ///
33    /// ```
34    /// use arr_rs::prelude::*;
35    ///
36    /// let array = Array::new(vec![12., -51., 4., 6., 167., -68., -4., 24., -41.], vec![3, 3]).unwrap();
37    /// let (vals, vecs) = array.eig().unwrap()[0].clone();
38    /// assert_eq!(Array::flat(vec![156.20350762022625, -35.473244608879945, 17.26973698865371]).unwrap(), vals);
39    /// assert_eq!(Array::new(vec![0.3274744043621118, 0.2741043942532785, -0.9929229419852038, -0.9370666300603048, 0.3040145084128465, 0.08536232642252124, -0.12110592601150529, 0.9123825731158716, 0.08256696983166067], vec![3, 3]).unwrap(), vecs);
40    /// ```
41    ///
42    /// # Errors
43    ///
44    /// may returns `ArrayError`
45    fn eig(&self) -> LinalgResult<N>;
46}
47
48impl <N: NumericOps> ArrayLinalgEigen<N> for Array<N> {
49
50    fn eigvals(&self) -> Result<Vec<Self>, ArrayError> {
51        self.is_square()?;
52        if self.ndim()? == 2 {
53            let mut h = self.hessenberg_reduction()?;
54
55            let mut max_iter = 10000;
56            // fix condition - check convergence
57            while max_iter > 0 {
58                let (q, r) = h.qr()?[0].clone();
59                h = r.dot(&q)?.hessenberg_reduction()?;
60                if h.is_convergent()? { break }
61                max_iter -= 1;
62            }
63
64            let mut eigenvalues = vec![];
65            for i in 0..self.get_shape()?[0] {
66                eigenvalues.push(h.at(&[i, i])?);
67            }
68
69            Ok(vec![Self::flat(eigenvalues)?])
70        } else {
71            let shape = self.get_shape()?;
72            let sub_shape = shape[self.ndim()? - 2 ..].to_vec();
73            let vals = self
74                .ravel()?
75                .split(self.len()? / sub_shape.iter().product::<usize>(), None)?
76                .iter()
77                .map(|arr| arr.reshape(&sub_shape).eigvals())
78                .collect::<Vec<Result<Vec<Self>, _>>>()
79                .has_error()?.into_iter()
80                .flat_map(Result::unwrap)
81                .collect::<Vec<Self>>();
82            Ok(vals)
83        }
84    }
85
86    fn eig(&self) -> LinalgResult<N> {
87        let mut results = vec![];
88        for eigenvalues in self.eigvals()? {
89            let mut vectors = vec![];
90            for value in eigenvalues.get_elements()? {
91                let mut vector = (self.clone() - (Self::eye(self.get_shape()?[0], None, None)? * value)?)
92                    .solve(&Self::ones(vec![self.get_shape()?[0]])?)?;
93                vector /= vector
94                    .norm(None::<NormOrd>, None, None)
95                    .broadcast_to(vector.get_shape()?)?;
96                vectors.push(vector);
97            }
98            let eigenvectors = vectors.iter()
99                .flatten()
100                .copied()
101                .collect::<Self>()
102                .reshape(&[eigenvalues.len()?; 2])
103                .transpose(None)?;
104            results.push((eigenvalues, eigenvectors));
105        }
106
107        Ok(results)
108    }
109}
110
111impl <N: NumericOps> ArrayLinalgEigen<N> for Result<Array<N>, ArrayError> {
112
113    fn eigvals(&self) -> Result<Vec<Array<N>>, ArrayError> {
114        self.clone()?.eigvals()
115    }
116
117    fn eig(&self) -> LinalgResult<N> {
118        self.clone()?.eig()
119    }
120}
arr_rs/linalg/operations/eigen.rs

arr_rs/linalg/operations/
eigen.rs
