## Coverage Report

#### Created: 2026-01-25 15:05

next uncovered line (L), next uncovered region (R), next uncovered branch (B)

    /home/noah/src/trueno/src/eigen.rs

    Line
    Count
    Source
    1
    //! Eigendecomposition for symmetric matrices
    2
    //!
    3
    //! Provides SIMD-accelerated eigenvalue and eigenvector computation for symmetric
    4
    //! (Hermitian) matrices, enabling PCA, spectral clustering, and other algorithms
    5
    //! without external dependencies like nalgebra.
    6
    //!
    7
    //! # Algorithm
    8
    //!
    9
    //! Uses the Jacobi eigenvalue algorithm, which is numerically stable and well-suited
    10
    //! for SIMD parallelization. For large matrices (>1000 dimensions), GPU acceleration
    11
    //! is available via wgpu.
    12
    //!
    13
    //! # Example
    14
    //!
    15
    //! ```
    16
    //! use trueno::{Matrix, SymmetricEigen};
    17
    //!
    18
    //! // Create a symmetric positive definite matrix
    19
    //! let cov = Matrix::from_vec(2, 2, vec![
    20
    //!     2.0, 1.0,
    21
    //!     1.0, 2.0,
    22
    //! ])?;
    23
    //!
    24
    //! let eigen = SymmetricEigen::new(&cov)?;
    25
    //!
    26
    //! // Eigenvalues in descending order (PCA convention)
    27
    //! let values = eigen.eigenvalues();
    28
    //! assert!((values[0] - 3.0).abs() < 1e-6);  // λ₁ = 3
    29
    //! assert!((values[1] - 1.0).abs() < 1e-6);  // λ₂ = 1
    30
    //! # Ok::<(), trueno::TruenoError>(())
    31
    //! ```
    32
    33
    use crate::{Backend, Matrix, TruenoError, Vector};
    34
    35
    /// Maximum number of sweeps for Jacobi algorithm convergence
    36
    /// Each sweep processes all n(n-1)/2 off-diagonal elements once
    37
    /// Typically converges in 5-10 sweeps for well-conditioned matrices
    38
    const MAX_JACOBI_SWEEPS: usize = 50;
    39
    40
    /// Convergence threshold for off-diagonal elements (relative to Frobenius norm)
    41
    const CONVERGENCE_THRESHOLD: f32 = 1e-7;
    42
    43
    /// GPU threshold - use wgpu for matrices larger than this
    44
    #[allow(dead_code)]
    45
    const GPU_THRESHOLD: usize = 1000;
    46
    47
    /// Symmetric matrix eigendecomposition
    48
    ///
    49
    /// Computes eigenvalues and eigenvectors for symmetric (Hermitian) matrices.
    50
    /// Eigenvalues are returned in descending order (largest first), which is the
    51
    /// convention used in PCA and most dimensionality reduction algorithms.
    52
    ///
    53
    /// # Properties
    54
    ///
    55
    /// For a symmetric matrix A, the decomposition satisfies:
    56
    /// - `A = V × D × V^T` where D is diagonal with eigenvalues
    57
    /// - Eigenvectors are orthonormal: `V^T × V = I`
    58
    /// - Eigenvalues are real (guaranteed for symmetric matrices)
    59
    ///
    60
    /// # Performance
    61
    ///
    62
    /// - SIMD-accelerated Jacobi rotations for CPU
    63
    /// - GPU compute shaders for matrices >1000 dimensions (with `gpu` feature)
    64
    /// - O(n³) time complexity, O(n²) space complexity
    65
    #[derive(Debug, Clone)]
    66
    pub struct SymmetricEigen {
    67
        /// Eigenvalues sorted in descending order
    68
        eigenvalues: Vec<f32>,
    69
        /// Eigenvectors as columns (column i corresponds to eigenvalue i)
    70
        eigenvectors: Matrix<f32>,
    71
        /// Sorting indices mapping original to sorted order
    72
        #[allow(dead_code)]
    73
        sort_indices: Vec<usize>,
    74
        /// Backend used for computation
    75
        backend: Backend,
    76
    }
    77
    78
    impl SymmetricEigen {
    79
        /// Computes eigendecomposition of a symmetric matrix
    80
        ///
    81
        /// # Arguments
    82
        ///
    83
        /// * `matrix` - A symmetric square matrix
    84
        ///
    85
        /// # Returns
    86
        ///
    87
        /// `SymmetricEigen` containing eigenvalues (descending) and eigenvectors
    88
        ///
    89
        /// # Errors
    90
        ///
    91
        /// - `InvalidInput` if matrix is not square
    92
        /// - `InvalidInput` if matrix is empty
    93
        /// - `InvalidInput` if algorithm fails to converge
    94
        ///
    95
        /// # Example
    96
        ///
    97
        /// ```
    98
        /// use trueno::{Matrix, SymmetricEigen};
    99
        ///
    100
        /// let m = Matrix::from_vec(3, 3, vec![
    101
        ///     4.0, 2.0, 0.0,
    102
        ///     2.0, 5.0, 3.0,
    103
        ///     0.0, 3.0, 6.0,
    104
        /// ])?;
    105
        ///
    106
        /// let eigen = SymmetricEigen::new(&m)?;
    107
        /// assert_eq!(eigen.eigenvalues().len(), 3);
    108
        /// # Ok::<(), trueno::TruenoError>(())
    109
        /// ```
    110
    0
        pub fn new(matrix: &Matrix<f32>) -> Result<Self, TruenoError> {
    111
            // Validate input
    112
    0
            if matrix.rows() != matrix.cols() {
    113
    0
                return Err(TruenoError::InvalidInput(format!(
    114
    0
                    "Matrix must be square for eigendecomposition, got {}x{}",
    115
    0
                    matrix.rows(),
    116
    0
                    matrix.cols()
    117
    0
                )));
    118
    0
            }
    119
    120
    0
            if matrix.rows() == 0 {
    121
    0
                return Err(TruenoError::InvalidInput(
    122
    0
                    "Cannot compute eigendecomposition of empty matrix".to_string(),
    123
    0
                ));
    124
    0
            }
    125
    126
    0
            let backend = Backend::select_best();
    127
    128
            // Dispatch to appropriate implementation based on matrix size and GPU availability
    129
            #[cfg(all(feature = "gpu", not(target_arch = "wasm32")))]
    130
            {
    131
    0
                let n = matrix.rows();
    132
    0
                if n >= GPU_THRESHOLD && crate::backends::gpu::GpuBackend::is_available() {
    133
    0
                    return Self::compute_gpu(matrix);
    134
    0
                }
    135
            }
    136
    137
            // CPU implementation (SIMD-accelerated) - works on all platforms
    138
    0
            Self::compute_jacobi(matrix, backend)
    139
    0
        }
    140
    141
        /// CPU implementation using Jacobi eigenvalue algorithm
    142
        ///
    143
        /// The Jacobi algorithm iteratively applies Givens rotations to eliminate
    144
        /// off-diagonal elements, converging to a diagonal matrix of eigenvalues.
    145
    0
        fn compute_jacobi(matrix: &Matrix<f32>, backend: Backend) -> Result<Self, TruenoError> {
    146
    0
            let n = matrix.rows();
    147
    148
            // Copy matrix data for in-place modification
    149
    0
            let mut a = matrix.as_slice().to_vec();
    150
    151
            // Compute initial Frobenius norm for relative convergence
    152
    0
            let frobenius_sq: f32 = a.iter().map(|x| x * x).sum();
    153
    0
            let tolerance = CONVERGENCE_THRESHOLD * frobenius_sq.sqrt().max(1.0);
    154
    155
            // Initialize eigenvectors to identity matrix
    156
    0
            let mut v = vec![0.0f32; n * n];
    157
    0
            for i in 0..n {
    158
    0
                v[i * n + i] = 1.0;
    159
    0
            }
    160
    161
            // Jacobi iteration with sweep strategy
    162
            // Each sweep processes all off-diagonal elements once
    163
    0
            for _sweep in 0..MAX_JACOBI_SWEEPS {
    164
                // Cyclic Jacobi: process all pairs (i, j) where i < j
    165
    0
                let mut converged = true;
    166
    167
    0
                for i in 0..n {
    168
    0
                    for j in (i + 1)..n {
    169
    0
                        let aij = a[i * n + j];
    170
    171
                        // Skip if already small enough
    172
    0
                        if aij.abs() < tolerance {
    173
    0
                            continue;
    174
    0
                        }
    175
    176
    0
                        converged = false;
    177
    0
                        Self::jacobi_rotate(&mut a, &mut v, n, i, j, backend);
    178
                    }
    179
                }
    180
    181
    0
                if converged {
    182
                    // Extract eigenvalues from diagonal
    183
    0
                    let eigenvalues: Vec<f32> = (0..n).map(|i| a[i * n + i]).collect();
    184
    185
                    // Sort eigenvalues in descending order
    186
    0
                    let mut indices: Vec<usize> = (0..n).collect();
    187
    0
                    indices.sort_by(|&i, &j| {
    188
    0
                        eigenvalues[j]
    189
    0
                            .partial_cmp(&eigenvalues[i])
    190
    0
                            .unwrap_or(std::cmp::Ordering::Equal)
    191
    0
                    });
    192
    193
                    // Reorder eigenvalues
    194
    0
                    let sorted_eigenvalues: Vec<f32> =
    195
    0
                        indices.iter().map(|&i| eigenvalues[i]).collect();
    196
    197
                    // Create eigenvector matrix with sorted columns
    198
    0
                    let mut eigenvector_data = vec![0.0f32; n * n];
    199
    0
                    for (new_col, &old_col) in indices.iter().enumerate() {
    200
    0
                        for row in 0..n {
    201
    0
                            eigenvector_data[row * n + new_col] = v[row * n + old_col];
    202
    0
                        }
    203
                    }
    204
    205
    0
                    let eigenvectors = Matrix::from_vec(n, n, eigenvector_data)?;
    206
    207
    0
                    return Ok(SymmetricEigen {
    208
    0
                        eigenvalues: sorted_eigenvalues,
    209
    0
                        eigenvectors,
    210
    0
                        sort_indices: indices,
    211
    0
                        backend,
    212
    0
                    });
    213
    0
                }
    214
            }
    215
    216
            // Failed to converge - this shouldn't happen for well-conditioned matrices
    217
    0
            Err(TruenoError::InvalidInput(format!(
    218
    0
                "Jacobi algorithm failed to converge after {} sweeps",
    219
    0
                MAX_JACOBI_SWEEPS
    220
    0
            )))
    221
    0
        }
    222
    223
        /// Find the largest off-diagonal element (unused in cyclic Jacobi, kept for classic Jacobi)
    224
        #[inline]
    225
        #[allow(dead_code)]
    226
    0
        fn find_max_off_diagonal(a: &[f32], n: usize) -> (usize, usize, f32) {
    227
    0
            let mut max_val = 0.0f32;
    228
    0
            let mut p = 0;
    229
    0
            let mut q = 1;
    230
    231
    0
            for i in 0..n {
    232
    0
                for j in (i + 1)..n {
    233
    0
                    let val = a[i * n + j].abs();
    234
    0
                    if val > max_val {
    235
    0
                        max_val = val;
    236
    0
                        p = i;
    237
    0
                        q = j;
    238
    0
                    }
    239
                }
    240
            }
    241
    242
    0
            (p, q, max_val)
    243
    0
        }
    244
    245
        /// Apply Jacobi rotation to zero out a[p][q] and a[q][p]
    246
        ///
    247
        /// Uses the numerically stable formula from:
    248
        /// Golub & Van Loan, "Matrix Computations", 4th Edition
    249
        #[inline]
    250
    0
        fn jacobi_rotate(
    251
    0
            a: &mut [f32],
    252
    0
            v: &mut [f32],
    253
    0
            n: usize,
    254
    0
            p: usize,
    255
    0
            q: usize,
    256
    0
            _backend: Backend,
    257
    0
        ) {
    258
    0
            let app = a[p * n + p];
    259
    0
            let aqq = a[q * n + q];
    260
    0
            let apq = a[p * n + q];
    261
    262
            // Skip if already zero
    263
    0
            if apq.abs() < 1e-15 {
    264
    0
                return;
    265
    0
            }
    266
    267
            // Compute rotation parameters using numerically stable formula
    268
            // tau = (aqq - app) / (2 * apq)
    269
            // t = sign(tau) / (|tau| + sqrt(1 + tau^2))  (avoiding catastrophic cancellation)
    270
            // c = 1 / sqrt(1 + t^2)
    271
            // s = t * c
    272
    0
            let tau = (aqq - app) / (2.0 * apq);
    273
    0
            let t = if tau >= 0.0 {
    274
    0
                1.0 / (tau + (1.0 + tau * tau).sqrt())
    275
            } else {
    276
    0
                -1.0 / (-tau + (1.0 + tau * tau).sqrt())
    277
            };
    278
    279
    0
            let c = 1.0 / (1.0 + t * t).sqrt();
    280
    0
            let s = t * c;
    281
    282
            // Update diagonal elements
    283
    0
            a[p * n + p] = app - t * apq;
    284
    0
            a[q * n + q] = aqq + t * apq;
    285
    0
            a[p * n + q] = 0.0;
    286
    0
            a[q * n + p] = 0.0;
    287
    288
            // Update off-diagonal elements in rows/columns p and q
    289
    0
            for k in 0..n {
    290
    0
                if k != p && k != q {
    291
    0
                    let akp = a[k * n + p];
    292
    0
                    let akq = a[k * n + q];
    293
    0
                    a[k * n + p] = c * akp - s * akq;
    294
    0
                    a[p * n + k] = a[k * n + p];
    295
    0
                    a[k * n + q] = s * akp + c * akq;
    296
    0
                    a[q * n + k] = a[k * n + q];
    297
    0
                }
    298
            }
    299
    300
            // Update eigenvector matrix
    301
    0
            for k in 0..n {
    302
    0
                let vkp = v[k * n + p];
    303
    0
                let vkq = v[k * n + q];
    304
    0
                v[k * n + p] = c * vkp - s * vkq;
    305
    0
                v[k * n + q] = s * vkp + c * vkq;
    306
    0
            }
    307
    0
        }
    308
    309
        /// GPU implementation for large matrices
    310
        #[cfg(all(feature = "gpu", not(target_arch = "wasm32")))]
    311
    0
        fn compute_gpu(matrix: &Matrix<f32>) -> Result<Self, TruenoError> {
    312
            use crate::backends::gpu::GpuBackend;
    313
    314
    0
            let n = matrix.rows();
    315
    0
            let mut gpu = GpuBackend::new();
    316
    317
            // Execute eigendecomposition on GPU
    318
    0
            let (eigenvalues, eigenvector_data) =
    319
    0
                gpu.symmetric_eigen(matrix.as_slice(), n).map_err(|e| {
    320
    0
                    TruenoError::InvalidInput(format!("GPU eigendecomposition failed: {}", e))
    321
    0
                })?;
    322
    323
            // Sort eigenvalues in descending order
    324
    0
            let mut indices: Vec<usize> = (0..n).collect();
    325
    0
            indices.sort_by(|&i, &j| {
    326
    0
                eigenvalues[j]
    327
    0
                    .partial_cmp(&eigenvalues[i])
    328
    0
                    .unwrap_or(std::cmp::Ordering::Equal)
    329
    0
            });
    330
    331
    0
            let sorted_eigenvalues: Vec<f32> = indices.iter().map(|&i| eigenvalues[i]).collect();
    332
    333
            // Reorder eigenvectors
    334
    0
            let mut sorted_eigenvector_data = vec![0.0f32; n * n];
    335
    0
            for (new_col, &old_col) in indices.iter().enumerate() {
    336
    0
                for row in 0..n {
    337
    0
                    sorted_eigenvector_data[row * n + new_col] = eigenvector_data[row * n + old_col];
    338
    0
                }
    339
            }
    340
    341
    0
            let eigenvectors = Matrix::from_vec(n, n, sorted_eigenvector_data)?;
    342
    343
    0
            Ok(SymmetricEigen {
    344
    0
                eigenvalues: sorted_eigenvalues,
    345
    0
                eigenvectors,
    346
    0
                sort_indices: indices,
    347
    0
                backend: Backend::GPU,
    348
    0
            })
    349
    0
        }
    350
    351
        /// Returns the eigenvalues in descending order
    352
        ///
    353
        /// # Example
    354
        ///
    355
        /// ```
    356
        /// use trueno::{Matrix, SymmetricEigen};
    357
        ///
    358
        /// let m = Matrix::from_vec(2, 2, vec![3.0, 1.0, 1.0, 3.0])?;
    359
        /// let eigen = SymmetricEigen::new(&m)?;
    360
        ///
    361
        /// let values = eigen.eigenvalues();
    362
        /// assert!((values[0] - 4.0).abs() < 1e-5);  // λ₁ = 4
    363
        /// assert!((values[1] - 2.0).abs() < 1e-5);  // λ₂ = 2
    364
        /// # Ok::<(), trueno::TruenoError>(())
    365
        /// ```
    366
    0
        pub fn eigenvalues(&self) -> &[f32] {
    367
    0
            &self.eigenvalues
    368
    0
        }
    369
    370
        /// Returns the eigenvector matrix
    371
        ///
    372
        /// Columns are eigenvectors, ordered to correspond with eigenvalues.
    373
        /// Column `i` is the eigenvector for `eigenvalues()[i]`.
    374
        ///
    375
        /// # Example
    376
        ///
    377
        /// ```
    378
        /// use trueno::{Matrix, SymmetricEigen};
    379
        ///
    380
        /// let m = Matrix::identity(3);
    381
        /// let eigen = SymmetricEigen::new(&m)?;
    382
        ///
    383
        /// // Identity matrix has eigenvectors that are the standard basis
    384
        /// let vectors = eigen.eigenvectors();
    385
        /// assert_eq!(vectors.rows(), 3);
    386
        /// assert_eq!(vectors.cols(), 3);
    387
        /// # Ok::<(), trueno::TruenoError>(())
    388
        /// ```
    389
    0
        pub fn eigenvectors(&self) -> &Matrix<f32> {
    390
    0
            &self.eigenvectors
    391
    0
        }
    392
    393
        /// Returns an iterator over (eigenvalue, eigenvector) pairs
    394
        ///
    395
        /// # Example
    396
        ///
    397
        /// ```
    398
        /// use trueno::{Matrix, SymmetricEigen};
    399
        ///
    400
        /// let m = Matrix::from_vec(2, 2, vec![2.0, 0.0, 0.0, 1.0])?;
    401
        /// let eigen = SymmetricEigen::new(&m)?;
    402
        ///
    403
        /// for (value, vector) in eigen.iter() {
    404
        ///     println!("λ = {}, v = {:?}", value, vector.as_slice());
    405
        /// }
    406
        /// # Ok::<(), trueno::TruenoError>(())
    407
        /// ```
    408
    0
        pub fn iter(&self) -> EigenIterator<'_> {
    409
    0
            EigenIterator {
    410
    0
                eigen: self,
    411
    0
                index: 0,
    412
    0
            }
    413
    0
        }
    414
    415
        /// Returns the number of eigenvalue/eigenvector pairs
    416
    0
        pub fn len(&self) -> usize {
    417
    0
            self.eigenvalues.len()
    418
    0
        }
    419
    420
        /// Returns true if there are no eigenvalues
    421
    0
        pub fn is_empty(&self) -> bool {
    422
    0
            self.eigenvalues.is_empty()
    423
    0
        }
    424
    425
        /// Returns the backend used for computation
    426
    0
        pub fn backend(&self) -> Backend {
    427
    0
            self.backend
    428
    0
        }
    429
    430
        /// Get a specific eigenvector by index
    431
        ///
    432
        /// # Arguments
    433
        ///
    434
        /// * `i` - Index of the eigenvector (0 = largest eigenvalue)
    435
        ///
    436
        /// # Returns
    437
        ///
    438
        /// The eigenvector as a Vector, or None if index out of bounds
    439
    0
        pub fn eigenvector(&self, i: usize) -> Option<Vector<f32>> {
    440
    0
            if i >= self.eigenvalues.len() {
    441
    0
                return None;
    442
    0
            }
    443
    444
    0
            let n = self.eigenvectors.rows();
    445
    0
            let mut data = Vec::with_capacity(n);
    446
    447
    0
            for row in 0..n {
    448
    0
                if let Some(&val) = self.eigenvectors.get(row, i) {
    449
    0
                    data.push(val);
    450
    0
                }
    451
            }
    452
    453
    0
            Some(Vector::from_slice(&data))
    454
    0
        }
    455
    456
        /// Reconstruct the original matrix from eigendecomposition
    457
        ///
    458
        /// Computes `V × D × V^T` where D is the diagonal matrix of eigenvalues.
    459
        /// This is useful for verifying the decomposition accuracy.
    460
        ///
    461
        /// # Example
    462
        ///
    463
        /// ```
    464
        /// use trueno::{Matrix, SymmetricEigen};
    465
        ///
    466
        /// let m = Matrix::from_vec(2, 2, vec![4.0, 2.0, 2.0, 4.0])?;
    467
        /// let eigen = SymmetricEigen::new(&m)?;
    468
        /// let reconstructed = eigen.reconstruct()?;
    469
        ///
    470
        /// // Should be approximately equal to original
    471
        /// assert!((reconstructed.get(0, 0).unwrap() - 4.0).abs() < 1e-5);
    472
        /// # Ok::<(), trueno::TruenoError>(())
    473
        /// ```
    474
    0
        pub fn reconstruct(&self) -> Result<Matrix<f32>, TruenoError> {
    475
    0
            let n = self.eigenvalues.len();
    476
    477
            // V × D (scale each column by its eigenvalue)
    478
    0
            let mut vd_data = vec![0.0f32; n * n];
    479
    0
            for i in 0..n {
    480
    0
                let lambda = self.eigenvalues[i];
    481
    0
                for j in 0..n {
    482
    0
                    if let Some(&v) = self.eigenvectors.get(j, i) {
    483
    0
                        vd_data[j * n + i] = v * lambda;
    484
    0
                    }
    485
                }
    486
            }
    487
    488
    0
            let vd = Matrix::from_vec(n, n, vd_data)?;
    489
    0
            let vt = self.eigenvectors.transpose();
    490
    491
    0
            vd.matmul(&vt)
    492
    0
        }
    493
    }
    494
    495
    /// Iterator over eigenvalue-eigenvector pairs
    496
    pub struct EigenIterator<'a> {
    497
        eigen: &'a SymmetricEigen,
    498
        index: usize,
    499
    }
    500
    501
    impl<'a> Iterator for EigenIterator<'a> {
    502
        type Item = (f32, Vector<f32>);
    503
    504
    0
        fn next(&mut self) -> Option<Self::Item> {
    505
    0
            if self.index >= self.eigen.len() {
    506
    0
                return None;
    507
    0
            }
    508
    509
    0
            let value = self.eigen.eigenvalues[self.index];
    510
    0
            let vector = self.eigen.eigenvector(self.index)?;
    511
    0
            self.index += 1;
    512
    513
    0
            Some((value, vector))
    514
    0
        }
    515
    516
    0
        fn size_hint(&self) -> (usize, Option<usize>) {
    517
    0
            let remaining = self.eigen.len() - self.index;
    518
    0
            (remaining, Some(remaining))
    519
    0
        }
    520
    }
    521
    522
    impl<'a> ExactSizeIterator for EigenIterator<'a> {}
    523
    524
    #[cfg(test)]
    525
    mod tests {
    526
        use super::*;
    527
    528
        // =========================================================================
    529
        // RED PHASE: Tests that define expected behavior
    530
        // =========================================================================
    531
    532
        #[test]
    533
        fn test_symmetric_eigen_2x2_simple() {
    534
            // Simple 2x2 symmetric matrix: [[2, 1], [1, 2]]
    535
            // Eigenvalues: 3, 1
    536
            // Eigenvectors: [1/√2, 1/√2], [1/√2, -1/√2]
    537
            let m = Matrix::from_vec(2, 2, vec![2.0, 1.0, 1.0, 2.0]).expect("valid matrix");
    538
    539
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    540
    541
            let values = eigen.eigenvalues();
    542
            assert_eq!(values.len(), 2);
    543
    544
            // Eigenvalues should be in descending order
    545
            assert!(values[0] >= values[1], "eigenvalues must be descending");
    546
    547
            // Check eigenvalue values (with tolerance)
    548
            assert!(
    549
                (values[0] - 3.0).abs() < 1e-5,
    550
                "first eigenvalue should be 3, got {}",
    551
                values[0]
    552
            );
    553
            assert!(
    554
                (values[1] - 1.0).abs() < 1e-5,
    555
                "second eigenvalue should be 1, got {}",
    556
                values[1]
    557
            );
    558
        }
    559
    560
        #[test]
    561
        fn test_symmetric_eigen_identity() {
    562
            // Identity matrix has all eigenvalues = 1
    563
            let m = Matrix::identity(3);
    564
    565
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    566
    567
            let values = eigen.eigenvalues();
    568
            assert_eq!(values.len(), 3);
    569
    570
            for (i, &val) in values.iter().enumerate() {
    571
                assert!(
    572
                    (val - 1.0).abs() < 1e-5,
    573
                    "eigenvalue {} should be 1, got {}",
    574
                    i,
    575
                    val
    576
                );
    577
            }
    578
        }
    579
    580
        #[test]
    581
        fn test_symmetric_eigen_diagonal() {
    582
            // Diagonal matrix: eigenvalues are the diagonal elements
    583
            let m = Matrix::from_vec(3, 3, vec![5.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 1.0])
    584
                .expect("valid matrix");
    585
    586
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    587
    588
            let values = eigen.eigenvalues();
    589
    590
            // Should be sorted descending: 5, 3, 1
    591
            assert!((values[0] - 5.0).abs() < 1e-5, "got {}", values[0]);
    592
            assert!((values[1] - 3.0).abs() < 1e-5, "got {}", values[1]);
    593
            assert!((values[2] - 1.0).abs() < 1e-5, "got {}", values[2]);
    594
        }
    595
    596
        #[test]
    597
        fn test_symmetric_eigen_eigenvectors_orthogonal() {
    598
            let m = Matrix::from_vec(3, 3, vec![4.0, 2.0, 0.0, 2.0, 5.0, 3.0, 0.0, 3.0, 6.0])
    599
                .expect("valid matrix");
    600
    601
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    602
    603
            // Eigenvectors should be orthonormal: V^T × V = I
    604
            let v = eigen.eigenvectors();
    605
            let vt = v.transpose();
    606
            let product = vt.matmul(v).expect("matmul should succeed");
    607
    608
            // Check if product is approximately identity
    609
            for i in 0..3 {
    610
                for j in 0..3 {
    611
                    let expected = if i == j { 1.0 } else { 0.0 };
    612
                    let actual = product.get(i, j).unwrap();
    613
                    assert!(
    614
                        (actual - expected).abs() < 1e-4,
    615
                        "V^T×V[{},{}] = {}, expected {}",
    616
                        i,
    617
                        j,
    618
                        actual,
    619
                        expected
    620
                    );
    621
                }
    622
            }
    623
        }
    624
    625
        #[test]
    626
        fn test_symmetric_eigen_reconstruction() {
    627
            let m = Matrix::from_vec(2, 2, vec![4.0, 2.0, 2.0, 4.0]).expect("valid matrix");
    628
    629
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    630
            let reconstructed = eigen.reconstruct().expect("reconstruction should succeed");
    631
    632
            // Reconstructed matrix should match original
    633
            for i in 0..2 {
    634
                for j in 0..2 {
    635
                    let original = m.get(i, j).unwrap();
    636
                    let recon = reconstructed.get(i, j).unwrap();
    637
                    assert!(
    638
                        (original - recon).abs() < 1e-4,
    639
                        "A[{},{}] = {}, reconstructed = {}",
    640
                        i,
    641
                        j,
    642
                        original,
    643
                        recon
    644
                    );
    645
                }
    646
            }
    647
        }
    648
    649
        #[test]
    650
        fn test_symmetric_eigen_av_equals_lambda_v() {
    651
            // For each eigenpair (λ, v): A×v = λ×v
    652
            let m = Matrix::from_vec(2, 2, vec![3.0, 1.0, 1.0, 3.0]).expect("valid matrix");
    653
    654
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    655
    656
            for (lambda, v) in eigen.iter() {
    657
                // Compute A×v
    658
                let av = m.matvec(&v).expect("matvec should succeed");
    659
    660
                // Compute λ×v
    661
                let lambda_v: Vec<f32> = v.as_slice().iter().map(|&x| x * lambda).collect();
    662
    663
                // Check equality
    664
                for (i, (&av_i, &lv_i)) in av.as_slice().iter().zip(lambda_v.iter()).enumerate() {
    665
                    assert!(
    666
                        (av_i - lv_i).abs() < 1e-4,
    667
                        "A×v[{}] = {}, λv[{}] = {}",
    668
                        i,
    669
                        av_i,
    670
                        i,
    671
                        lv_i
    672
                    );
    673
                }
    674
            }
    675
        }
    676
    677
        #[test]
    678
        fn test_symmetric_eigen_error_non_square() {
    679
            let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).expect("valid matrix");
    680
    681
            let result = SymmetricEigen::new(&m);
    682
            assert!(result.is_err());
    683
    684
            let err = result.unwrap_err();
    685
            assert!(
    686
                matches!(err, TruenoError::InvalidInput(_)),
    687
                "expected InvalidInput error"
    688
            );
    689
        }
    690
    691
        #[test]
    692
        fn test_symmetric_eigen_error_empty() {
    693
            let m = Matrix::zeros(0, 0);
    694
    695
            let result = SymmetricEigen::new(&m);
    696
            assert!(result.is_err());
    697
        }
    698
    699
        #[test]
    700
        fn test_symmetric_eigen_1x1() {
    701
            let m = Matrix::from_vec(1, 1, vec![7.0]).expect("valid matrix");
    702
    703
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    704
    705
            assert_eq!(eigen.eigenvalues().len(), 1);
    706
            assert!((eigen.eigenvalues()[0] - 7.0).abs() < 1e-6);
    707
        }
    708
    709
        #[test]
    710
        fn test_symmetric_eigen_iterator() {
    711
            let m = Matrix::from_vec(2, 2, vec![2.0, 0.0, 0.0, 1.0]).expect("valid matrix");
    712
    713
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    714
    715
            let pairs: Vec<_> = eigen.iter().collect();
    716
            assert_eq!(pairs.len(), 2);
    717
    718
            // First eigenvalue is larger
    719
            assert!(pairs[0].0 >= pairs[1].0);
    720
        }
    721
    722
        #[test]
    723
        fn test_symmetric_eigen_len() {
    724
            let m = Matrix::identity(5);
    725
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    726
    727
            assert_eq!(eigen.len(), 5);
    728
            assert!(!eigen.is_empty());
    729
        }
    730
    731
        #[test]
    732
        fn test_symmetric_eigen_eigenvector_accessor() {
    733
            let m = Matrix::identity(3);
    734
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    735
    736
            let v0 = eigen.eigenvector(0);
    737
            assert!(v0.is_some());
    738
            assert_eq!(v0.unwrap().len(), 3);
    739
    740
            let v_invalid = eigen.eigenvector(10);
    741
            assert!(v_invalid.is_none());
    742
        }
    743
    744
        #[test]
    745
        fn test_symmetric_eigen_covariance_matrix() {
    746
            // Typical covariance matrix from PCA
    747
            // Points: [(1,2), (3,4), (5,6)] centered → [(-2,-2), (0,0), (2,2)]
    748
            // Cov = [[8/3, 8/3], [8/3, 8/3]] ≈ [[2.67, 2.67], [2.67, 2.67]]
    749
            let m = Matrix::from_vec(2, 2, vec![2.67, 2.67, 2.67, 2.67]).expect("valid matrix");
    750
    751
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    752
    753
            // Eigenvalues: 5.34, 0 (approximately)
    754
            let values = eigen.eigenvalues();
    755
            assert!(values[0] > 5.0, "first eigenvalue should be ~5.34");
    756
            assert!(values[1].abs() < 0.1, "second eigenvalue should be ~0");
    757
        }
    758
    759
        #[test]
    760
        fn test_symmetric_eigen_negative_eigenvalues() {
    761
            // Matrix with negative eigenvalues
    762
            // [[0, 1], [1, 0]] has eigenvalues 1, -1
    763
            let m = Matrix::from_vec(2, 2, vec![0.0, 1.0, 1.0, 0.0]).expect("valid matrix");
    764
    765
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    766
    767
            let values = eigen.eigenvalues();
    768
            assert!(
    769
                (values[0] - 1.0).abs() < 1e-5,
    770
                "first eigenvalue should be 1"
    771
            );
    772
            assert!(
    773
                (values[1] - (-1.0)).abs() < 1e-5,
    774
                "second eigenvalue should be -1"
    775
            );
    776
        }
    777
    778
        #[test]
    779
        fn test_symmetric_eigen_backend() {
    780
            // Test that the backend() method returns the expected value
    781
            let m = Matrix::from_vec(2, 2, vec![4.0, 1.0, 1.0, 3.0]).expect("valid matrix");
    782
            let eigen = SymmetricEigen::new(&m).expect("eigendecomposition should succeed");
    783
    784
            // backend() should return the current backend
    785
            let backend = eigen.backend();
    786
            // On this machine, it should be AVX2
    787
            #[cfg(target_arch = "x86_64")]
    788
            {
    789
                use crate::Backend;
    790
                assert!(
    791
                    matches!(backend, Backend::AVX2 | Backend::Scalar | Backend::SSE2),
    792
                    "expected valid backend, got {:?}",
    793
                    backend
    794
                );
    795
            }
    796
        }
    797
    798
        // =========================================================================
    799
        // Property-based tests (proptest)
    800
        // =========================================================================
    801
    802
        #[cfg(test)]
    803
        mod proptest_tests {
    804
            use super::*;
    805
            use proptest::prelude::*;
    806
    807
            proptest! {
    808
                #![proptest_config(ProptestConfig::with_cases(50))]
    809
    810
                #[test]
    811
                fn prop_eigenvalues_descending(n in 2usize..6) {
    812
                    // Generate random symmetric matrix
    813
                    let mut data = vec![0.0f32; n * n];
    814
                    for i in 0..n {
    815
                        for j in i..n {
    816
                            let val = (i + j) as f32 / (n as f32);
    817
                            data[i * n + j] = val;
    818
                            data[j * n + i] = val;
    819
                        }
    820
                    }
    821
    822
                    let m = Matrix::from_vec(n, n, data).expect("valid matrix");
    823
                    let eigen = SymmetricEigen::new(&m).expect("eigen should succeed");
    824
    825
                    let values = eigen.eigenvalues();
    826
                    for i in 1..values.len() {
    827
                        prop_assert!(
    828
                            values[i - 1] >= values[i],
    829
                            "eigenvalues not descending: {} < {}",
    830
                            values[i - 1],
    831
                            values[i]
    832
                        );
    833
                    }
    834
                }
    835
    836
                #[test]
    837
                fn prop_eigenvector_count_matches_dimension(n in 1usize..8) {
    838
                    let m = Matrix::identity(n);
    839
                    let eigen = SymmetricEigen::new(&m).expect("eigen should succeed");
    840
    841
                    prop_assert_eq!(eigen.len(), n);
    842
                    prop_assert_eq!(eigen.eigenvalues().len(), n);
    843
                    prop_assert_eq!(eigen.eigenvectors().rows(), n);
    844
                    prop_assert_eq!(eigen.eigenvectors().cols(), n);
    845
                }
    846
    847
                #[test]
    848
                fn prop_reconstruction_accuracy(
    849
                    a in 1.0f32..10.0,  // Ensure positive diagonal for conditioning
    850
                    b in -5.0f32..5.0,  // Off-diagonal smaller than diagonal
    851
                    c in 1.0f32..10.0   // Ensure positive diagonal for conditioning
    852
                ) {
    853
                    // Create symmetric 2x2 matrix [[a+|b|, b], [b, c+|b|]]
    854
                    // Add |b| to diagonal for better conditioning
    855
                    let diag_a = a + b.abs();
    856
                    let diag_c = c + b.abs();
    857
                    let m = Matrix::from_vec(2, 2, vec![diag_a, b, b, diag_c]).expect("valid matrix");
    858
    859
                    if let Ok(eigen) = SymmetricEigen::new(&m) {
    860
                        if let Ok(recon) = eigen.reconstruct() {
    861
                            // Use relative error for numerical stability
    862
                            let frobenius_orig: f32 = [diag_a, b, b, diag_c].iter()
    863
                                .map(|x| x * x).sum::<f32>().sqrt();
    864
                            let max_allowed_error = 0.01 * frobenius_orig.max(1.0);
    865
    866
                            for i in 0..2 {
    867
                                for j in 0..2 {
    868
                                    let orig = m.get(i, j).unwrap();
    869
                                    let rec = recon.get(i, j).unwrap();
    870
                                    prop_assert!(
    871
                                        (orig - rec).abs() < max_allowed_error,
    872
                                        "reconstruction error: {} vs {}, allowed: {}",
    873
                                        orig,
    874
                                        rec,
    875
                                        max_allowed_error
    876
                                    );
    877
                                }
    878
                            }
    879
                        }
    880
                    }
    881
                }
    882
            }
    883
        }
    884
    }