rlst/dense/linalg/lapack/
lq.rs

1//! LQ decomposition and full-row-rank minimum-norm solves.
2
3use crate::UnsafeRandom1DAccessByValue;
4use crate::base_types::{RlstError, RlstResult};
5use crate::dense::array::{Array, DynArray};
6use crate::dense::linalg::lapack::interface::trsm::{TrsmDiag, TrsmSide, TrsmTransA, TrsmUplo};
7use crate::traits::base_operations::Shape;
8use crate::traits::linalg::{decompositions::Lq, lapack::Lapack};
9
10pub use super::qr::QMode;
11pub use crate::dense::linalg::lapack::interface::mqr::{MqrSide, MqrTransMode};
12
13/// Stores the packed reflectors and scalar factors of an LQ decomposition.
14///
15/// For an `m` by `n` input and `k = min(m, n)`, the decomposition satisfies
16/// `A = L Q`. [`Self::l_mat`] has shape `m` by `k`. Compact Q has shape
17/// `k` by `n`, while full Q has shape `n` by `n`. Thus compact factors
18/// reconstruct A directly; with full Q, use its first `k` rows or pad L with
19/// zero columns to shape `m` by `n`.
20pub struct LqDecomposition<Item> {
21    a: DynArray<Item, 2>,
22    tau: Vec<Item>,
23}
24
25impl<Item, ArrayImpl> Lq for Array<ArrayImpl, 2>
26where
27    Item: Lapack,
28    ArrayImpl: UnsafeRandom1DAccessByValue<Item = Item> + Shape<2>,
29{
30    type Item = Item;
31
32    fn lq(&self) -> RlstResult<LqDecomposition<Item>> {
33        let mut a = DynArray::new_from(self);
34        let [m, n] = a.shape();
35        let mut tau = vec![Item::zero(); m.min(n)];
36        if m != 0 && n != 0 {
37            Item::gelqf(m, n, a.data_mut().unwrap(), m, &mut tau)?;
38        }
39        Ok(LqDecomposition { a, tau })
40    }
41}
42
43impl<Item> LqDecomposition<Item>
44where
45    Item: Lapack,
46{
47    /// Form Q explicitly.
48    ///
49    /// Compact Q has shape `min(m, n)` by `n`; full Q has shape `n` by `n`.
50    /// For an empty reflector set, full Q is the identity matrix.
51    pub fn q_mat(&self, mode: QMode) -> RlstResult<DynArray<Item, 2>> {
52        let [m, n] = self.a.shape();
53        let k = m.min(n);
54        let q_rows = match mode {
55            QMode::Compact => k,
56            QMode::Full => n,
57        };
58        let mut q = DynArray::<Item, 2>::from_shape([q_rows, n]);
59        for row in 0..k {
60            for col in row + 1..n {
61                q[[row, col]] = self.a[[row, col]];
62            }
63        }
64        if k == 0 {
65            if matches!(mode, QMode::Full) {
66                q.set_identity();
67            }
68            return Ok(q);
69        }
70        Item::orglq(
71            q_rows,
72            n,
73            k,
74            q.data_mut().unwrap(),
75            q_rows.max(1),
76            &self.tau,
77        )?;
78        Ok(q)
79    }
80
81    /// Extract the `m` by `min(m, n)` lower trapezoidal factor L.
82    pub fn l_mat(&self) -> RlstResult<DynArray<Item, 2>> {
83        let [m, n] = self.a.shape();
84        let k = m.min(n);
85        let mut l = DynArray::<Item, 2>::from_shape([m, k]);
86        for col in 0..k {
87            for row in col..m {
88                l[[row, col]] = self.a[[row, col]];
89            }
90        }
91        Ok(l)
92    }
93
94    /// Apply full Q or Q^H implicitly to a two-dimensional input.
95    ///
96    /// [`MqrSide::Left`] computes `Q rhs` and
97    /// [`MqrSide::Right`] computes `rhs Q`; selecting
98    /// [`MqrTransMode::ConjugateTranspose`] substitutes `Q^H` (which is `Q^T`
99    /// for real scalars). The affected input dimension must equal the number
100    /// of columns of the factorized matrix. The input and packed factorization
101    /// are not modified.
102    pub fn apply_q<ArrayImpl>(
103        &self,
104        side: MqrSide,
105        trans: MqrTransMode,
106        rhs: &Array<ArrayImpl, 2>,
107    ) -> RlstResult<DynArray<Item, 2>>
108    where
109        ArrayImpl: UnsafeRandom1DAccessByValue<Item = Item> + Shape<2>,
110    {
111        let [rows, cols] = rhs.shape();
112        let n = self.a.shape()[1];
113        let affected = match side {
114            MqrSide::Left => rows,
115            MqrSide::Right => cols,
116        };
117        if affected != n {
118            return Err(RlstError::GeneralError(format!(
119                "LQ Q application requires affected dimension {n}, got {affected}."
120            )));
121        }
122        let mut result = DynArray::new_from(rhs);
123        if rows == 0 || cols == 0 || self.tau.is_empty() {
124            return Ok(result);
125        }
126        Item::mlq(
127            side,
128            trans,
129            rows,
130            cols,
131            self.tau.len(),
132            self.a.data().unwrap(),
133            self.a.shape()[0].max(1),
134            &self.tau,
135            result.data_mut().unwrap(),
136            rows.max(1),
137        )?;
138        Ok(result)
139    }
140
141    /// Solve `A X = B` for the minimum Euclidean/Frobenius-norm solution.
142    ///
143    /// The factorized matrix must satisfy `m <= n` and have full row rank.
144    /// A one-dimensional right-hand side has length `m` and returns length
145    /// `n`; a two-dimensional right-hand side has shape `m` by `nrhs` and
146    /// returns shape `n` by `nrhs`. The result has minimum Euclidean norm for
147    /// a vector and minimum Frobenius norm for a matrix.
148    ///
149    /// This factorization is not rank revealing: near-rank-deficient inputs
150    /// have no reliability guarantee. Only exact zero diagonal entries of the
151    /// triangular factor are rejected.
152    pub fn solve<ArrayImpl, const NDIM: usize>(
153        &self,
154        rhs: &Array<ArrayImpl, NDIM>,
155    ) -> RlstResult<DynArray<Item, NDIM>>
156    where
157        ArrayImpl: UnsafeRandom1DAccessByValue<Item = Item> + Shape<NDIM>,
158    {
159        if NDIM != 1 && NDIM != 2 {
160            return Err(RlstError::GeneralError(
161                "LQ solve is only implemented for 1D and 2D arrays.".to_string(),
162            ));
163        }
164        let [m, n] = self.a.shape();
165        if m > n {
166            return Err(RlstError::GeneralError(
167                "LQ minimum-norm solve requires m <= n.".to_string(),
168            ));
169        }
170        if rhs.shape()[0] != m {
171            return Err(RlstError::GeneralError(format!(
172                "LQ solve requires {m} right-hand-side rows, got {}.",
173                rhs.shape()[0]
174            )));
175        }
176        for diagonal in 0..m {
177            if self.a[[diagonal, diagonal]] == Item::zero() {
178                return Err(RlstError::GeneralError(
179                    "LQ minimum-norm solve requires full row rank.".to_string(),
180                ));
181            }
182        }
183
184        let nrhs = if NDIM == 1 { 1 } else { rhs.shape()[1] };
185        let mut shape = rhs.shape();
186        shape[0] = n;
187        let mut solution = DynArray::<Item, NDIM>::from_shape(shape);
188        let rhs_copy = DynArray::new_from(rhs);
189        for col in 0..nrhs {
190            for row in 0..m {
191                solution.data_mut().unwrap()[row + n * col] =
192                    rhs_copy.data().unwrap()[row + m * col];
193            }
194        }
195        if m == 0 || nrhs == 0 {
196            return Ok(solution);
197        }
198
199        let mut l11 = DynArray::<Item, 2>::from_shape([m, m]);
200        for col in 0..m {
201            for row in col..m {
202                l11[[row, col]] = self.a[[row, col]];
203            }
204        }
205        Item::trsm(
206            TrsmSide::Left,
207            TrsmUplo::Lower,
208            TrsmTransA::NoTrans,
209            TrsmDiag::NonUnit,
210            m,
211            nrhs,
212            Item::one(),
213            l11.data_mut().unwrap(),
214            m,
215            solution.data_mut().unwrap(),
216            n,
217        )?;
218        Item::mlq(
219            MqrSide::Left,
220            MqrTransMode::ConjugateTranspose,
221            n,
222            nrhs,
223            m,
224            self.a.data().unwrap(),
225            m.max(1),
226            &self.tau,
227            solution.data_mut().unwrap(),
228            n,
229        )?;
230        Ok(solution)
231    }
232}
233
234#[cfg(test)]
235mod tests {
236    use super::{MqrSide, MqrTransMode, QMode};
237    use crate::base_types::{c32, c64};
238    use crate::dense::array::DynArray;
239    use crate::{Lq, MultIntoResize, RlstScalar, empty_array};
240    use paste::paste;
241
242    macro_rules! lq_scalar_tests {
243        ($scalar:ty, $tol:expr) => {
244            paste! {
245                #[test]
246                fn [<lq_reconstruct_apply_and_solve_ $scalar>]() {
247                    let mut a = DynArray::<$scalar, 2>::from_shape([2, 3]);
248                    a.fill_from_seed_equally_distributed(17);
249                    a[[0, 0]] += <$scalar>::from_real(3.0);
250                    a[[1, 1]] += <$scalar>::from_real(3.0);
251                    let lq = a.lq().unwrap();
252                    let l = lq.l_mat().unwrap();
253                    let qc = lq.q_mat(QMode::Compact).unwrap();
254                    let reconstructed = empty_array().simple_mult_into_resize(l.r(), qc.r());
255                    crate::assert_array_relative_eq!(reconstructed, a, $tol);
256                    let qch = DynArray::new_from(&qc.r().conj().transpose());
257                    let qcqch = empty_array().simple_mult_into_resize(qc.r(), qch.r());
258                    let mut compact_identity = DynArray::<$scalar, 2>::from_shape([2, 2]);
259                    compact_identity.set_identity();
260                    crate::assert_array_relative_eq!(qcqch, compact_identity, $tol);
261
262                    for shape in [[2, 2], [3, 2]] {
263                        let mut other = DynArray::<$scalar, 2>::from_shape(shape);
264                        other.fill_from_seed_equally_distributed(31);
265                        other[[0, 0]] += <$scalar>::from_real(3.0);
266                        other[[1, 1]] += <$scalar>::from_real(3.0);
267                        let other_lq = other.lq().unwrap();
268                        let other_reconstructed = empty_array().simple_mult_into_resize(
269                            other_lq.l_mat().unwrap().r(),
270                            other_lq.q_mat(QMode::Compact).unwrap().r(),
271                        );
272                        crate::assert_array_relative_eq!(other_reconstructed, other, $tol);
273                    }
274
275                    let q = lq.q_mat(QMode::Full).unwrap();
276                    let qh = DynArray::new_from(&q.r().conj().transpose());
277                    let mut identity = DynArray::<$scalar, 2>::from_shape([3, 3]);
278                    identity.set_identity();
279                    let qqh = empty_array().simple_mult_into_resize(q.r(), qh.r());
280                    crate::assert_array_relative_eq!(qqh, identity, $tol);
281
282                    let mut left_rhs = DynArray::<$scalar, 2>::from_shape([3, 2]);
283                    left_rhs.fill_from_seed_equally_distributed(41);
284                    let left = lq.apply_q(MqrSide::Left, MqrTransMode::NoTranspose, &left_rhs).unwrap();
285                    let explicit_left = empty_array().simple_mult_into_resize(q.r(), left_rhs.r());
286                    crate::assert_array_relative_eq!(left, explicit_left, $tol);
287                    let left_adjoint = lq.apply_q(MqrSide::Left, MqrTransMode::ConjugateTranspose, &left_rhs).unwrap();
288                    let explicit_left_adjoint = empty_array().simple_mult_into_resize(qh.r(), left_rhs.r());
289                    crate::assert_array_relative_eq!(left_adjoint, explicit_left_adjoint, $tol);
290                    let mut right_rhs = DynArray::<$scalar, 2>::from_shape([2, 3]);
291                    right_rhs.fill_from_seed_equally_distributed(43);
292                    let right = lq.apply_q(MqrSide::Right, MqrTransMode::ConjugateTranspose, &right_rhs).unwrap();
293                    let explicit_right = empty_array().simple_mult_into_resize(right_rhs.r(), qh.r());
294                    crate::assert_array_relative_eq!(right, explicit_right, $tol);
295                    let right_plain = lq.apply_q(MqrSide::Right, MqrTransMode::NoTranspose, &right_rhs).unwrap();
296                    let explicit_right_plain = empty_array().simple_mult_into_resize(right_rhs.r(), q.r());
297                    crate::assert_array_relative_eq!(right_plain, explicit_right_plain, $tol);
298                    let q_again = lq.q_mat(QMode::Full).unwrap();
299                    crate::assert_array_relative_eq!(q_again, q, $tol);
300
301                    let mut expected = DynArray::<$scalar, 2>::from_shape([3, 2]);
302                    expected.fill_from_seed_equally_distributed(23);
303                    let rhs = empty_array().simple_mult_into_resize(a.r(), expected.r());
304                    let x = lq.solve(&rhs).unwrap();
305                    let residual = empty_array().simple_mult_into_resize(a.r(), x.r());
306                    crate::assert_array_relative_eq!(residual, rhs, $tol);
307                    let qx = empty_array().simple_mult_into_resize(q.r(), x.r());
308                    for col in 0..2 {
309                        assert!(qx[[2, col]].abs() < $tol * 10.0);
310                    }
311
312                    let mut rhs_vector = DynArray::<$scalar, 1>::from_shape([2]);
313                    rhs_vector[[0]] = rhs[[0, 0]];
314                    rhs_vector[[1]] = rhs[[1, 0]];
315                    let vector_x = lq.solve(&rhs_vector).unwrap();
316                    assert_eq!(vector_x.shape(), [3]);
317                    for row in 0..3 {
318                        assert!((vector_x[[row]] - x[[row, 0]]).abs() < $tol);
319                    }
320                }
321            }
322        };
323    }
324
325    lq_scalar_tests!(f32, 2e-4);
326    lq_scalar_tests!(f64, 2e-12);
327    lq_scalar_tests!(c32, 3e-4);
328    lq_scalar_tests!(c64, 3e-12);
329
330    #[test]
331    fn lq_shapes_views_errors_and_empty_contracts() {
332        let mut source = DynArray::<f64, 2>::from_shape([3, 2]);
333        source.fill_from_seed_equally_distributed(5);
334        let view = source.transpose();
335        let lq = view.lq().unwrap();
336        assert_eq!(lq.l_mat().unwrap().shape(), [2, 2]);
337        assert_eq!(lq.q_mat(QMode::Compact).unwrap().shape(), [2, 3]);
338        assert_eq!(lq.q_mat(QMode::Full).unwrap().shape(), [3, 3]);
339        assert!(
340            lq.apply_q(
341                MqrSide::Left,
342                MqrTransMode::NoTranspose,
343                &DynArray::from_shape([2, 1])
344            )
345            .is_err()
346        );
347        assert!(lq.solve(&DynArray::<f64, 1>::from_shape([3])).is_err());
348        let rhs_storage = DynArray::<f64, 2>::from_shape([1, 2]);
349        assert_eq!(lq.solve(&rhs_storage.transpose()).unwrap().shape(), [3, 1]);
350        let apply_storage = DynArray::<f64, 2>::from_shape([1, 3]);
351        assert_eq!(
352            lq.apply_q(
353                MqrSide::Left,
354                MqrTransMode::NoTranspose,
355                &apply_storage.transpose(),
356            )
357            .unwrap()
358            .shape(),
359            [3, 1]
360        );
361
362        let tall = DynArray::<f64, 2>::from_shape([3, 2]).lq().unwrap();
363        assert!(tall.solve(&DynArray::<f64, 1>::from_shape([3])).is_err());
364        let singular = DynArray::<f64, 2>::from_shape([2, 3]).lq().unwrap();
365        assert!(
366            singular
367                .solve(&DynArray::<f64, 1>::from_shape([2]))
368                .is_err()
369        );
370
371        let empty = DynArray::<f64, 2>::from_shape([0, 3]).lq().unwrap();
372        let q = empty.q_mat(QMode::Full).unwrap();
373        let mut identity = DynArray::<f64, 2>::from_shape([3, 3]);
374        identity.set_identity();
375        crate::assert_array_abs_diff_eq!(q, identity, 0.0);
376        let rhs = DynArray::<f64, 2>::from_shape([0, 2]);
377        let x = empty.solve(&rhs).unwrap();
378        assert_eq!(x.shape(), [3, 2]);
379        assert!(x.iter_value().all(|value| value == 0.0));
380        let zero_columns = lq.solve(&DynArray::<f64, 2>::from_shape([2, 0])).unwrap();
381        assert_eq!(zero_columns.shape(), [3, 0]);
382        let applied = lq
383            .apply_q(
384                MqrSide::Left,
385                MqrTransMode::NoTranspose,
386                &DynArray::<f64, 2>::from_shape([3, 0]),
387            )
388            .unwrap();
389        assert_eq!(applied.shape(), [3, 0]);
390        assert!(
391            lq.solve(&DynArray::<f64, 3>::from_shape([2, 1, 1]))
392                .is_err()
393        );
394    }
395}
rlst/dense/linalg/lapack/lq.rs

rlst/dense/linalg/lapack/
lq.rs
