ndarray_linalg/
eig.rs

1//! Eigenvalue decomposition for non-symmetric square matrices
2
3use crate::error::*;
4use crate::layout::*;
5use crate::types::*;
6pub use lax::GeneralizedEigenvalue;
7use ndarray::*;
8
9#[cfg_attr(doc, katexit::katexit)]
10/// Eigenvalue decomposition of general matrix reference
11pub trait Eig {
12    /// EigVec is the right eivenvector
13    type EigVal;
14    type EigVec;
15    /// Calculate eigenvalues with the right eigenvector
16    ///
17    /// $$ A u_i = \lambda_i u_i $$
18    ///
19    /// ```
20    /// use ndarray::*;
21    /// use ndarray_linalg::*;
22    ///
23    /// let a: Array2<f64> = array![
24    ///     [-1.01,  0.86, -4.60,  3.31, -4.81],
25    ///     [ 3.98,  0.53, -7.04,  5.29,  3.55],
26    ///     [ 3.30,  8.26, -3.89,  8.20, -1.51],
27    ///     [ 4.43,  4.96, -7.66, -7.33,  6.18],
28    ///     [ 7.31, -6.43, -6.16,  2.47,  5.58],
29    /// ];
30    /// let (eigs, vecs) = a.eig().unwrap();
31    ///
32    /// let a = a.map(|v| v.as_c());
33    /// for (&e, vec) in eigs.iter().zip(vecs.axis_iter(Axis(1))) {
34    ///     let ev = vec.map(|v| v * e);
35    ///     let av = a.dot(&vec);
36    ///     assert_close_l2!(&av, &ev, 1e-5);
37    /// }
38    /// ```
39    fn eig(&self) -> Result<(Self::EigVal, Self::EigVec)>;
40}
41
42impl<A, S> Eig for ArrayBase<S, Ix2>
43where
44    A: Scalar + Lapack,
45    S: Data<Elem = A>,
46{
47    type EigVal = Array1<A::Complex>;
48    type EigVec = Array2<A::Complex>;
49
50    fn eig(&self) -> Result<(Self::EigVal, Self::EigVec)> {
51        let mut a = self.to_owned();
52        let layout = a.square_layout()?;
53        let (s, t) = A::eig(true, layout, a.as_allocated_mut()?)?;
54        let n = layout.len() as usize;
55        Ok((
56            ArrayBase::from(s),
57            Array2::from_shape_vec((n, n).f(), t).unwrap(),
58        ))
59    }
60}
61
62/// Calculate eigenvalues without eigenvectors
63pub trait EigVals {
64    type EigVal;
65    fn eigvals(&self) -> Result<Self::EigVal>;
66}
67
68impl<A, S> EigVals for ArrayBase<S, Ix2>
69where
70    A: Scalar + Lapack,
71    S: Data<Elem = A>,
72{
73    type EigVal = Array1<A::Complex>;
74
75    fn eigvals(&self) -> Result<Self::EigVal> {
76        let mut a = self.to_owned();
77        let (s, _) = A::eig(false, a.square_layout()?, a.as_allocated_mut()?)?;
78        Ok(ArrayBase::from(s))
79    }
80}
81
82#[cfg_attr(doc, katexit::katexit)]
83/// Eigenvalue decomposition of general matrix reference
84pub trait EigGeneralized {
85    /// EigVec is the right eivenvector
86    type EigVal;
87    type EigVec;
88    type Real;
89    /// Calculate eigenvalues with the right eigenvector
90    ///
91    /// $$ A u_i = \lambda_i B u_i $$
92    ///
93    /// ```
94    /// use ndarray::*;
95    /// use ndarray_linalg::*;
96    ///
97    /// let a: Array2<f64> = array![
98    ///     [-1.01,  0.86, -4.60,  3.31, -4.81],
99    ///     [ 3.98,  0.53, -7.04,  5.29,  3.55],
100    ///     [ 3.30,  8.26, -3.89,  8.20, -1.51],
101    ///     [ 4.43,  4.96, -7.66, -7.33,  6.18],
102    ///     [ 7.31, -6.43, -6.16,  2.47,  5.58],
103    /// ];
104    /// let b: Array2<f64> = array![
105    ///     [ 1.23, -4.56,  7.89,  0.12, -3.45],
106    ///     [ 6.78, -9.01,  2.34, -5.67,  8.90],
107    ///     [-1.11,  3.33, -6.66,  9.99, -2.22],
108    ///     [ 4.44, -7.77,  0.00,  1.11,  5.55],
109    ///     [-8.88,  6.66, -3.33,  2.22, -9.99],
110    /// ];
111    /// let (geneigs, vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
112    ///
113    /// let a = a.map(|v| v.as_c());
114    /// let b = b.map(|v| v.as_c());
115    /// for (ge, vec) in geneigs.iter().zip(vecs.axis_iter(Axis(1))) {
116    ///     if let GeneralizedEigenvalue::Finite(e, _) = ge {
117    ///         let ebv = b.dot(&vec).map(|v| v * e);
118    ///         let av = a.dot(&vec);
119    ///         assert_close_l2!(&av, &ebv, 1e-5);
120    ///     }
121    /// }
122    /// ```
123    ///
124    /// # Arguments
125    ///
126    /// * `thresh_opt` - An optional threshold for determining approximate zero |β| values when
127    /// computing the eigenvalues as α/β. If `None`, no approximate comparisons to zero will be
128    /// made.
129    fn eig_generalized(
130        &self,
131        thresh_opt: Option<Self::Real>,
132    ) -> Result<(Self::EigVal, Self::EigVec)>;
133}
134
135impl<A, S> EigGeneralized for (ArrayBase<S, Ix2>, ArrayBase<S, Ix2>)
136where
137    A: Scalar + Lapack,
138    S: Data<Elem = A>,
139{
140    type EigVal = Array1<GeneralizedEigenvalue<A::Complex>>;
141    type EigVec = Array2<A::Complex>;
142    type Real = A::Real;
143
144    fn eig_generalized(
145        &self,
146        thresh_opt: Option<Self::Real>,
147    ) -> Result<(Self::EigVal, Self::EigVec)> {
148        let (mut a, mut b) = (self.0.to_owned(), self.1.to_owned());
149        let layout = a.square_layout()?;
150        let (s, t) = A::eig_generalized(
151            true,
152            layout,
153            a.as_allocated_mut()?,
154            b.as_allocated_mut()?,
155            thresh_opt,
156        )?;
157        let n = layout.len() as usize;
158        Ok((
159            ArrayBase::from(s),
160            Array2::from_shape_vec((n, n).f(), t).unwrap(),
161        ))
162    }
163}
ndarray_linalg/eig.rs

ndarray_linalg/
eig.rs
