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}