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> Eig for ArrayRef<A, Ix2>
43where
44    A: Scalar + Lapack,
45{
46    type EigVal = Array1<A::Complex>;
47    type EigVec = Array2<A::Complex>;
48
49    fn eig(&self) -> Result<(Self::EigVal, Self::EigVec)> {
50        let mut a = self.to_owned();
51        let layout = a.square_layout()?;
52        let (s, t) = A::eig(true, layout, a.as_allocated_mut()?)?;
53        let n = layout.len() as usize;
54        Ok((
55            ArrayBase::from(s),
56            Array2::from_shape_vec((n, n).f(), t).unwrap(),
57        ))
58    }
59}
60
61/// Calculate eigenvalues without eigenvectors
62pub trait EigVals {
63    type EigVal;
64    fn eigvals(&self) -> Result<Self::EigVal>;
65}
66
67impl<A> EigVals for ArrayRef<A, Ix2>
68where
69    A: Scalar + Lapack,
70{
71    type EigVal = Array1<A::Complex>;
72
73    fn eigvals(&self) -> Result<Self::EigVal> {
74        let mut a = self.to_owned();
75        let (s, _) = A::eig(false, a.square_layout()?, a.as_allocated_mut()?)?;
76        Ok(ArrayBase::from(s))
77    }
78}
79
80#[cfg_attr(doc, katexit::katexit)]
81/// Eigenvalue decomposition of general matrix reference
82pub trait EigGeneralized {
83    /// EigVec is the right eivenvector
84    type EigVal;
85    type EigVec;
86    type Real;
87    /// Calculate eigenvalues with the right eigenvector
88    ///
89    /// $$ A u_i = \lambda_i B u_i $$
90    ///
91    /// ```
92    /// use ndarray::*;
93    /// use ndarray_linalg::*;
94    ///
95    /// let a: Array2<f64> = array![
96    ///     [-1.01,  0.86, -4.60,  3.31, -4.81],
97    ///     [ 3.98,  0.53, -7.04,  5.29,  3.55],
98    ///     [ 3.30,  8.26, -3.89,  8.20, -1.51],
99    ///     [ 4.43,  4.96, -7.66, -7.33,  6.18],
100    ///     [ 7.31, -6.43, -6.16,  2.47,  5.58],
101    /// ];
102    /// let b: Array2<f64> = array![
103    ///     [ 1.23, -4.56,  7.89,  0.12, -3.45],
104    ///     [ 6.78, -9.01,  2.34, -5.67,  8.90],
105    ///     [-1.11,  3.33, -6.66,  9.99, -2.22],
106    ///     [ 4.44, -7.77,  0.00,  1.11,  5.55],
107    ///     [-8.88,  6.66, -3.33,  2.22, -9.99],
108    /// ];
109    /// let (geneigs, vecs) = (a.clone(), b.clone()).eig_generalized(None).unwrap();
110    ///
111    /// let a = a.map(|v| v.as_c());
112    /// let b = b.map(|v| v.as_c());
113    /// for (ge, vec) in geneigs.iter().zip(vecs.axis_iter(Axis(1))) {
114    ///     if let GeneralizedEigenvalue::Finite(e, _) = ge {
115    ///         let ebv = b.dot(&vec).map(|v| v * e);
116    ///         let av = a.dot(&vec);
117    ///         assert_close_l2!(&av, &ebv, 1e-5);
118    ///     }
119    /// }
120    /// ```
121    ///
122    /// # Arguments
123    ///
124    /// * `thresh_opt` - An optional threshold for determining approximate zero |β| values when
125    /// computing the eigenvalues as α/β. If `None`, no approximate comparisons to zero will be
126    /// made.
127    fn eig_generalized(
128        self,
129        thresh_opt: Option<Self::Real>,
130    ) -> Result<(Self::EigVal, Self::EigVec)>;
131}
132
133/// Turn arrays, references to arrays, and [`ArrayRef`]s into owned arrays
134pub trait MaybeOwnedMatrix {
135    type Elem;
136
137    /// Convert into an owned array, cloning only when necessary.
138    fn into_owned(self) -> Array2<Self::Elem>;
139}
140
141impl<S> MaybeOwnedMatrix for ArrayBase<S, Ix2>
142where
143    S: Data,
144    S::Elem: Clone,
145{
146    type Elem = S::Elem;
147
148    fn into_owned(self) -> Array2<S::Elem> {
149        ArrayBase::into_owned(self)
150    }
151}
152
153impl<S> MaybeOwnedMatrix for &ArrayBase<S, Ix2>
154where
155    S: Data,
156    S::Elem: Clone,
157{
158    type Elem = S::Elem;
159
160    fn into_owned(self) -> Array2<S::Elem> {
161        self.to_owned()
162    }
163}
164
165impl<A> MaybeOwnedMatrix for &ArrayRef2<A>
166where
167    A: Clone,
168{
169    type Elem = A;
170
171    fn into_owned(self) -> Array2<A> {
172        self.to_owned()
173    }
174}
175
176impl<T1, T2> EigGeneralized for (T1, T2)
177where
178    T1: MaybeOwnedMatrix,
179    T1::Elem: Lapack + Scalar,
180    T2: MaybeOwnedMatrix<Elem = T1::Elem>,
181{
182    type EigVal = Array1<GeneralizedEigenvalue<<T1::Elem as Scalar>::Complex>>;
183    type EigVec = Array2<<T1::Elem as Scalar>::Complex>;
184    type Real = <T1::Elem as Scalar>::Real;
185
186    fn eig_generalized(
187        self,
188        thresh_opt: Option<Self::Real>,
189    ) -> Result<(Self::EigVal, Self::EigVec)> {
190        let (mut a, mut b) = (self.0.into_owned(), self.1.into_owned());
191        let layout = a.square_layout()?;
192        let (s, t) = T1::Elem::eig_generalized(
193            true,
194            layout,
195            a.as_allocated_mut()?,
196            b.as_allocated_mut()?,
197            thresh_opt,
198        )?;
199        let n = layout.len() as usize;
200        Ok((
201            ArrayBase::from(s),
202            Array2::from_shape_vec((n, n).f(), t).unwrap(),
203        ))
204    }
205}
206
207#[cfg(test)]
208mod tests {
209    use crate::MaybeOwnedMatrix;
210
211    #[test]
212    fn test_maybe_owned_matrix() {
213        let a = array![[1.0, 2.0], [3.0, 4.0]];
214        let a_ptr = a.as_ptr();
215        let a1 = MaybeOwnedMatrix::into_owned(a);
216        assert_eq!(a_ptr, a1.as_ptr());
217
218        let b = a1.clone();
219        let b1 = MaybeOwnedMatrix::into_owned(&b);
220        assert_eq!(b, b1);
221        assert_ne!(b.as_ptr(), b1.as_ptr());
222
223        let b2 = MaybeOwnedMatrix::into_owned(&*b);
224        assert_eq!(b, b2);
225        assert_ne!(b.as_ptr(), b2.as_ptr());
226    }
227}
ndarray_linalg/eig.rs

ndarray_linalg/
eig.rs
