trueno/eigen/
methods.rs

1//! Public API methods and iterator for `SymmetricEigen`.
2
3use super::SymmetricEigen;
4use crate::{Backend, Matrix, TruenoError, Vector};
5
6impl SymmetricEigen {
7    /// Returns the eigenvalues in descending order
8    ///
9    /// # Example
10    ///
11    /// ```
12    /// use trueno::{Matrix, SymmetricEigen};
13    ///
14    /// let m = Matrix::from_vec(2, 2, vec![3.0, 1.0, 1.0, 3.0])?;
15    /// let eigen = SymmetricEigen::new(&m)?;
16    ///
17    /// let values = eigen.eigenvalues();
18    /// assert!((values[0] - 4.0).abs() < 1e-5);  // λ₁ = 4
19    /// assert!((values[1] - 2.0).abs() < 1e-5);  // λ₂ = 2
20    /// # Ok::<(), trueno::TruenoError>(())
21    /// ```
22    pub fn eigenvalues(&self) -> &[f32] {
23        &self.eigenvalues
24    }
25
26    /// Returns the eigenvector matrix
27    ///
28    /// Columns are eigenvectors, ordered to correspond with eigenvalues.
29    /// Column `i` is the eigenvector for `eigenvalues()[i]`.
30    ///
31    /// # Example
32    ///
33    /// ```
34    /// use trueno::{Matrix, SymmetricEigen};
35    ///
36    /// let m = Matrix::identity(3);
37    /// let eigen = SymmetricEigen::new(&m)?;
38    ///
39    /// // Identity matrix has eigenvectors that are the standard basis
40    /// let vectors = eigen.eigenvectors();
41    /// assert_eq!(vectors.rows(), 3);
42    /// assert_eq!(vectors.cols(), 3);
43    /// # Ok::<(), trueno::TruenoError>(())
44    /// ```
45    pub fn eigenvectors(&self) -> &Matrix<f32> {
46        &self.eigenvectors
47    }
48
49    /// Returns an iterator over (eigenvalue, eigenvector) pairs
50    ///
51    /// # Example
52    ///
53    /// ```
54    /// use trueno::{Matrix, SymmetricEigen};
55    ///
56    /// let m = Matrix::from_vec(2, 2, vec![2.0, 0.0, 0.0, 1.0])?;
57    /// let eigen = SymmetricEigen::new(&m)?;
58    ///
59    /// for (value, vector) in eigen.iter() {
60    ///     println!("λ = {}, v = {:?}", value, vector.as_slice());
61    /// }
62    /// # Ok::<(), trueno::TruenoError>(())
63    /// ```
64    pub fn iter(&self) -> EigenIterator<'_> {
65        EigenIterator { eigen: self, index: 0 }
66    }
67
68    /// Returns the number of eigenvalue/eigenvector pairs
69    pub fn len(&self) -> usize {
70        self.eigenvalues.len()
71    }
72
73    /// Returns true if there are no eigenvalues
74    pub fn is_empty(&self) -> bool {
75        self.eigenvalues.is_empty()
76    }
77
78    /// Returns the backend used for computation
79    pub fn backend(&self) -> Backend {
80        self.backend
81    }
82
83    /// Get a specific eigenvector by index
84    ///
85    /// # Arguments
86    ///
87    /// * `i` - Index of the eigenvector (0 = largest eigenvalue)
88    ///
89    /// # Returns
90    ///
91    /// The eigenvector as a Vector, or None if index out of bounds
92    pub fn eigenvector(&self, i: usize) -> Option<Vector<f32>> {
93        contract_pre_eigenvector!();
94        if i >= self.eigenvalues.len() {
95            return None;
96        }
97
98        let n = self.eigenvectors.rows();
99        let mut data = Vec::with_capacity(n);
100
101        for row in 0..n {
102            if let Some(&val) = self.eigenvectors.get(row, i) {
103                data.push(val);
104            }
105        }
106
107        Some(Vector::from_slice(&data))
108    }
109
110    /// Reconstruct the original matrix from eigendecomposition
111    ///
112    /// Computes `V × D × V^T` where D is the diagonal matrix of eigenvalues.
113    /// This is useful for verifying the decomposition accuracy.
114    ///
115    /// # Example
116    ///
117    /// ```
118    /// use trueno::{Matrix, SymmetricEigen};
119    ///
120    /// let m = Matrix::from_vec(2, 2, vec![4.0, 2.0, 2.0, 4.0])?;
121    /// let eigen = SymmetricEigen::new(&m)?;
122    /// let reconstructed = eigen.reconstruct()?;
123    ///
124    /// // Should be approximately equal to original
125    /// assert!((reconstructed.get(0, 0).unwrap() - 4.0).abs() < 1e-5);
126    /// # Ok::<(), trueno::TruenoError>(())
127    /// ```
128    pub fn reconstruct(&self) -> Result<Matrix<f32>, TruenoError> {
129        let n = self.eigenvalues.len();
130
131        // V × D (scale each column by its eigenvalue)
132        let mut vd_data = vec![0.0f32; n * n];
133        for i in 0..n {
134            let lambda = self.eigenvalues[i];
135            for j in 0..n {
136                if let Some(&v) = self.eigenvectors.get(j, i) {
137                    vd_data[j * n + i] = v * lambda;
138                }
139            }
140        }
141
142        let vd = Matrix::from_vec(n, n, vd_data)?;
143        let vt = self.eigenvectors.transpose();
144
145        vd.matmul(&vt)
146    }
147}
148
149/// Iterator over eigenvalue-eigenvector pairs
150pub struct EigenIterator<'a> {
151    pub(super) eigen: &'a SymmetricEigen,
152    pub(super) index: usize,
153}
154
155impl<'a> Iterator for EigenIterator<'a> {
156    type Item = (f32, Vector<f32>);
157
158    fn next(&mut self) -> Option<Self::Item> {
159        if self.index >= self.eigen.len() {
160            return None;
161        }
162
163        let value = self.eigen.eigenvalues[self.index];
164        let vector = self.eigen.eigenvector(self.index)?;
165        self.index += 1;
166
167        Some((value, vector))
168    }
169
170    fn size_hint(&self) -> (usize, Option<usize>) {
171        let remaining = self.eigen.len() - self.index;
172        (remaining, Some(remaining))
173    }
174}
175
176impl<'a> ExactSizeIterator for EigenIterator<'a> {}
trueno/eigen/methods.rs

trueno/eigen/
methods.rs
