Skip to main content

symplex/domains/
decompositions.rs

1//! Named results of matrix decompositions and normal forms.
2//!
3//! Every factorisation here used to come back as a tuple of two or three
4//! matrices — `(Q, R)`, `(P, D)`, `(H, U)`, `(S, U, V)` — where nothing but
5//! memory said which position was which, and a transposition compiled
6//! silently.  These structs name each factor and document the identity it
7//! satisfies, at zero runtime cost.  See CONTRIBUTING.md, "Tuples versus
8//! structs".
9//!
10//! The structs are generic over the matrix type so that the symbolic
11//! [`Matrix`](crate::matrix::Matrix) and the exact
12//! [`ZMatrix`](crate::matrix::ZMatrix) / [`QMatrix`](crate::matrix::QMatrix)
13//! share one vocabulary: `Qr<Matrix>`, `HermiteNormalForm<ZMatrix>`, ….
14//!
15//! # Examples
16//!
17//! ```
18//! use symplex::prelude::*;
19//!
20//! let ctx = Context::new();
21//! let a = matrix![ctx, [4, 2], [2, 3]];
22//! let ldl = a.ldl().unwrap();
23//! assert_eq!((&(&ldl.l * &ldl.d) * &ldl.l.transpose()).eval(), a);
24//!
25//! let Qr { q, r } = matrix![ctx, [1, 1], [0, 1]].qr().unwrap();
26//! assert_eq!((&q * &r).simplify(), matrix![ctx, [1, 1], [0, 1]]);
27//! ```
28
29/// QR decomposition `A = Q·R`: `Q` has orthonormal columns, `R` is upper
30/// triangular.
31#[derive(Clone, Debug, PartialEq)]
32pub struct Qr<M> {
33    /// `Q` (`m × n`), orthonormal columns.
34    pub q: M,
35    /// `R` (`n × n`), upper triangular.
36    pub r: M,
37}
38
39/// LDLᵀ decomposition `A = L·D·Lᵀ` of a symmetric matrix: `L` unit lower
40/// triangular, `D` diagonal.
41#[derive(Clone, Debug, PartialEq)]
42pub struct Ldl<M> {
43    /// `L`, unit lower triangular.
44    pub l: M,
45    /// `D`, diagonal.
46    pub d: M,
47}
48
49/// LU decomposition with partial pivoting `P·A = L·U`: `L` unit lower
50/// triangular, `U` upper triangular, and `perm` the row permutation —
51/// row `i` of `P·A` is row `perm[i]` of `A`.
52#[derive(Clone, Debug, PartialEq)]
53pub struct Lu<M> {
54    /// `L`, unit lower triangular.
55    pub l: M,
56    /// `U`, upper triangular.
57    pub u: M,
58    /// Row permutation: `perm[i]` is the original index of row `i` of `P·A`.
59    pub perm: Vec<usize>,
60}
61
62/// Eigendecomposition `A = P·D·P⁻¹`: the columns of `P` are eigenvectors,
63/// `D` is diagonal with the eigenvalues in the same order.
64#[derive(Clone, Debug, PartialEq)]
65pub struct Diagonalization<M> {
66    /// `P`, invertible, eigenvectors as columns.
67    pub p: M,
68    /// `D`, diagonal of eigenvalues.
69    pub d: M,
70}
71
72/// Jordan normal form `A = P·J·P⁻¹`: `J` is block diagonal with Jordan
73/// blocks `J_k(λ)` (eigenvalue on the diagonal, ones on the superdiagonal),
74/// `P` holds the (generalized) eigenvectors.
75#[derive(Clone, Debug, PartialEq)]
76pub struct JordanForm<M> {
77    /// `P`, invertible, (generalized) eigenvectors as columns.
78    pub p: M,
79    /// `J`, block diagonal of Jordan blocks.
80    pub j: M,
81}
82
83/// Upper Hessenberg form by similarity: `H = P⁻¹·A·P` (equivalently
84/// `A·P = P·H`) with `h_ij = 0` for `i > j + 1`.
85#[derive(Clone, Debug, PartialEq)]
86pub struct Hessenberg<M> {
87    /// `H`, upper Hessenberg.
88    pub h: M,
89    /// `P`, the invertible similarity transform.
90    pub p: M,
91}
92
93/// Full-rank factorisation `A = C·F` with `r = rank A`: `C` (`m × r`) is
94/// made of the pivot columns of `A`, `F` (`r × n`) of the nonzero rows of
95/// `rref(A)`.
96#[derive(Clone, Debug, PartialEq)]
97pub struct RankDecomposition<M> {
98    /// `C` (`m × r`), the pivot columns of `A`.
99    pub c: M,
100    /// `F` (`r × n`), the nonzero rows of `rref(A)`.
101    pub f: M,
102}
103
104/// Row-style Hermite normal form `H = U·A` with `U` unimodular
105/// (`det U = ±1`).
106#[derive(Clone, Debug, PartialEq)]
107pub struct HermiteNormalForm<M> {
108    /// `H`, the Hermite normal form of `A`.
109    pub h: M,
110    /// `U`, unimodular, with `H = U·A`.
111    pub u: M,
112}
113
114/// Smith normal form `S = U·A·V` with `U`, `V` unimodular
115/// (`det U = det V = ±1`) and `S = diag(d₁, …, dᵣ, 0, …)`, `dᵢ | dᵢ₊₁`.
116#[derive(Clone, Debug, PartialEq)]
117pub struct SmithNormalForm<M> {
118    /// `S`, the diagonal Smith normal form of `A`.
119    pub s: M,
120    /// `U`, unimodular row transform.
121    pub u: M,
122    /// `V`, unimodular column transform.
123    pub v: M,
124}
125
126/// LLL reduction of the lattice basis formed by the rows of `A`:
127/// `reduced = transform·A` with `transform` unimodular (`det = ±1`), so
128/// both span the same lattice.
129#[derive(Clone, Debug, PartialEq)]
130pub struct LllReduction<M> {
131    /// The LLL-reduced basis, one lattice vector per row.
132    pub reduced: M,
133    /// `T`, unimodular, with `reduced = T·A`.
134    pub transform: M,
135}