Skip to main content

tenferro_linalg/
rank_revealing_qr.rs

1use tenferro_tensor::Error;
2
3use crate::QrGauge;
4
5/// Four fixed-shape outputs of a rank-revealing QR factorization.
6///
7/// `Q` and `R` use `T`; the column permutation and numerical rank use `M`.
8/// For dynamic, eager, and traced tensors, `T` and `M` are the same tensor
9/// type. Typed tensors use [`crate::TypedRankRevealingQrResult`], whose metadata
10/// tensors have scalar type `i64`.
11///
12/// # Examples
13///
14/// ```rust
15/// use tenferro_linalg::RankRevealingQrResult;
16///
17/// let result = RankRevealingQrResult {
18///     q: "q",
19///     r: "r",
20///     column_permutation: vec![1_i64, 0],
21///     rank: vec![2_i64],
22/// };
23/// assert_eq!(result.column_permutation, [1, 0]);
24/// assert_eq!(result.rank, [2]);
25/// ```
26#[derive(Clone, Debug)]
27pub struct RankRevealingQrResult<T, M = T> {
28    /// Thin orthonormal factor with shape `[m, min(m, n), batch...]`.
29    pub q: T,
30    /// Upper-trapezoidal factor with shape `[min(m, n), n, batch...]`.
31    pub r: T,
32    /// Zero-based original column at each factor column, shaped `[n, batch...]`.
33    pub column_permutation: M,
34    /// Leading-prefix numerical rank, shaped `[batch...]`.
35    pub rank: M,
36}
37
38/// Options for column-pivoted rank-revealing QR.
39///
40/// Rank is the length of the leading diagonal prefix satisfying
41/// `abs(R[i,i]) > max(atol, rtol * abs(R[0,0]))`. Both tolerances must be
42/// finite and non-negative. The default uses zero tolerances so callers choose
43/// an application-specific numerical threshold explicitly.
44///
45/// # Examples
46///
47/// ```rust
48/// use tenferro_linalg::{QrGauge, RankRevealingQrOptions};
49///
50/// let options = RankRevealingQrOptions::default()
51///     .gauge(QrGauge::PositiveDiagonal)
52///     .rtol(1.0e-10)
53///     .atol(1.0e-14);
54/// assert_eq!(options.rtol, 1.0e-10);
55/// assert_eq!(options.atol, 1.0e-14);
56/// ```
57#[derive(Clone, Copy, Debug, PartialEq)]
58pub struct RankRevealingQrOptions {
59    /// QR sign or phase convention. The default is [`QrGauge::Raw`].
60    pub gauge: QrGauge,
61    /// Relative diagonal threshold. The default is `0.0`.
62    pub rtol: f64,
63    /// Absolute diagonal threshold. The default is `0.0`.
64    pub atol: f64,
65}
66
67impl Default for RankRevealingQrOptions {
68    fn default() -> Self {
69        Self {
70            gauge: QrGauge::Raw,
71            rtol: 0.0,
72            atol: 0.0,
73        }
74    }
75}
76
77impl RankRevealingQrOptions {
78    /// Return options with the requested QR gauge.
79    ///
80    /// # Examples
81    ///
82    /// ```rust
83    /// use tenferro_linalg::{QrGauge, RankRevealingQrOptions};
84    ///
85    /// let options = RankRevealingQrOptions::default().gauge(QrGauge::PositiveDiagonal);
86    /// assert_eq!(options.gauge, QrGauge::PositiveDiagonal);
87    /// ```
88    pub fn gauge(mut self, gauge: QrGauge) -> Self {
89        self.gauge = gauge;
90        self
91    }
92
93    /// Return options with the requested relative rank tolerance.
94    ///
95    /// # Examples
96    ///
97    /// ```rust
98    /// use tenferro_linalg::RankRevealingQrOptions;
99    ///
100    /// let options = RankRevealingQrOptions::default().rtol(1.0e-8);
101    /// assert_eq!(options.rtol, 1.0e-8);
102    /// ```
103    pub fn rtol(mut self, rtol: f64) -> Self {
104        self.rtol = rtol;
105        self
106    }
107
108    /// Return options with the requested absolute rank tolerance.
109    ///
110    /// # Examples
111    ///
112    /// ```rust
113    /// use tenferro_linalg::RankRevealingQrOptions;
114    ///
115    /// let options = RankRevealingQrOptions::default().atol(1.0e-12);
116    /// assert_eq!(options.atol, 1.0e-12);
117    /// ```
118    pub fn atol(mut self, atol: f64) -> Self {
119        self.atol = atol;
120        self
121    }
122}
123
124pub(crate) fn validate_rank_revealing_qr_options(
125    op: &'static str,
126    options: RankRevealingQrOptions,
127) -> tenferro_tensor::Result<()> {
128    for (name, value) in [("rtol", options.rtol), ("atol", options.atol)] {
129        if !value.is_finite() || value < 0.0 {
130            return Err(Error::invalid_argument(
131                op,
132                name,
133                format!("must be finite and non-negative, got {value}"),
134            ));
135        }
136    }
137    Ok(())
138}