Skip to main content

mdarray_linalg/
svd.rs

1//! Singular Value Decomposition (SVD)
2//!
3//! The matrix A is decomposed as `A = U * S * V^T` where:
4//! - `s` contains the singular values (1D vector)
5//! - `u` contains the left singular vectors (matrix U)
6//! - `vt` contains the transposed right singular vectors (matrix V^T)
7//!
8//! Singular values are mathematically real. Backends choose the scalar type
9//! used to represent them through [`SVD::SingularValue`].
10//!```rust,ignore
11//!// ----- Singular Value Decomposition (SVD) -----
12//!use mdarray_linalg::svd::SVDDecomp;
13//!use mdarray_linalg::prelude::*; // Import traits anonymously
14//!use mdarray_linalg_backend::Backend; // Use the real backend here, Lapack, Faer, ...
15//!
16//!let bd = Backend::default();
17//!let SVDDecomp { s, u, vt } = bd.svd(&mut a.clone()).expect("SVD failed");
18//!// Or the shorter ...
19//!let SVDDecomp { s, u, vt } = bd.svd(&mut a.clone()).expect("SVD failed");
20//!```
21use mdarray::{Array, Dim, Layout, Slice};
22use thiserror::Error;
23
24/// Error types related to singular value decomposition
25#[derive(Debug, Error)]
26pub enum SVDError {
27    #[error("Backend error code: {0}")]
28    BackendError(i32),
29
30    #[error("Inconsistent U and VT: must be both Some or both None")]
31    InconsistentUV,
32
33    #[error("Backend failed to converge: {superdiagonals} superdiagonals did not converge to zero")]
34    BackendDidNotConverge { superdiagonals: i32 },
35}
36
37/// Holds the results of a singular value decomposition, including
38/// singular values and the left and right singular vectors.
39///
40/// `T` is the matrix scalar type, `S` is the singular-value scalar type,
41/// and `D` is the matrix dimension type.
42pub struct SVDDecomp<T, S, D: Dim> {
43    pub s: Array<S, (D,)>,
44    pub u: Array<T, (D, D)>,
45    pub vt: Array<T, (D, D)>,
46}
47
48/// Singular value decomposition for matrix factorization and analysis
49pub trait SVD<T, D: Dim> {
50    /// Scalar type used for singular values.
51    ///
52    /// Singular values are mathematically real. Backends choose their
53    /// representation; some current backends use the matrix scalar type `T`.
54    type SingularValue;
55
56    /// Compute full SVD with new allocated matrices
57    fn svd<L: Layout>(
58        &self,
59        a: &mut Slice<T, (D, D), L>,
60    ) -> Result<SVDDecomp<T, Self::SingularValue, D>, SVDError>;
61
62    /// Compute thin SVD with new allocated matrices
63    fn svd_thin<L: Layout>(
64        &self,
65        a: &mut Slice<T, (D, D), L>,
66    ) -> Result<SVDDecomp<T, Self::SingularValue, D>, SVDError>;
67
68    /// Compute only singular values with new allocated matrix
69    fn svd_s<L: Layout>(
70        &self,
71        a: &mut Slice<T, (D, D), L>,
72    ) -> Result<Array<Self::SingularValue, (D,)>, SVDError>;
73
74    /// Compute SVD, overwriting existing matrices
75    fn svd_write<L: Layout, Ls: Layout, Lu: Layout, Lvt: Layout>(
76        &self,
77        a: &mut Slice<T, (D, D), L>,
78        s: &mut Slice<Self::SingularValue, (D,), Ls>,
79        u: &mut Slice<T, (D, D), Lu>,
80        vt: &mut Slice<T, (D, D), Lvt>,
81    ) -> Result<(), SVDError>;
82
83    /// Compute only singular values, overwriting existing matrix
84    fn svd_write_s<L: Layout, Ls: Layout>(
85        &self,
86        a: &mut Slice<T, (D, D), L>,
87        s: &mut Slice<Self::SingularValue, (D,), Ls>,
88    ) -> Result<(), SVDError>;
89}