From e0f266bbb335808d0331dfeb34162de00fb4b347 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jannik=20Schu=CC=88rg?= Date: Sat, 9 May 2020 16:36:57 +0200 Subject: [PATCH] Change the QR implementation, inspired by LAPACK --- CHANGELOG.md | 17 +++ src/linalg/decomposition.rs | 3 +- src/linalg/qr.rs | 236 ++++++++++++++++++++++-------------- src/linalg/svd3.rs | 20 ++- tests/linalg/qr.rs | 97 +++++++++++++-- 5 files changed, 272 insertions(+), 101 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8d64160d3..513f7151d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,23 @@ This project adheres to [Semantic Versioning](https://semver.org/). **nalgebra-lapack change log** is found [here](https://github.com/dimforge/nalgebra/blob/main/nalgebra-lapack/CHANGELOG.md) starting with `nalgebra-lapack` version `0.27.0`. +## Unreleased + +### Added + +- Add `QR::q_columns` which computes the first `k` columns of `Q`, up to the full `Q` of a tall matrix. The + existing `QR::q` still returns the first `min(nrows, ncols)` columns only. + +### Changed + +- The `QR` decomposition now scales the Householder vectors like LAPACK's `?GEQR2` and `?ORG2R`, which is more + accurate. This is a breaking change: + - The internal representation changed. `QR::qr_internal` now holds `R` in its upper trapezoidal part and the + scaled Householder vectors below it. The `diag` field became LAPACK's `tau` vector. + - The diagonal of `R` is no longer forced to be positive. Some columns of `Q` and rows of `R` can have a + different sign than before. + - `Matrix::qr` and `QR::new` need the additional `DefaultAllocator: Allocator` bound. + ## [0.35.0] (24 May 2026) ### Added diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index c22ffceae..07486a21a 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -57,7 +57,8 @@ impl> Matrix { pub fn qr(self) -> QR where R: DimMin, - DefaultAllocator: Allocator + Allocator + Allocator>, + DefaultAllocator: + Allocator + Allocator + Allocator + Allocator>, { QR::new(self.into_owned()) } diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs index 64983e890..8feb8848c 100644 --- a/src/linalg/qr.rs +++ b/src/linalg/qr.rs @@ -3,17 +3,20 @@ use num::Zero; use serde::{Deserialize, Serialize}; use crate::allocator::{Allocator, Reallocator}; -use crate::base::{DefaultAllocator, Matrix, OMatrix, OVector, Unit}; +use crate::base::{DefaultAllocator, Matrix, OMatrix, OVector}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Const, Dim, DimMin, DimMinimum}; use crate::storage::{Storage, StorageMut}; use simba::scalar::ComplexField; -use crate::geometry::Reflection; -use crate::linalg::householder; use std::mem::MaybeUninit; /// The QR decomposition of a general matrix. +/// +/// The decomposition is stored like in LAPACK's `?GEQR2`: the upper trapezoidal part of `qr` is +/// `R`, and the columns below its diagonal hold the Householder vectors `v` (their first component +/// is `1` and is not stored). `Q` is the product `H(0) * H(1) * ... ` of the reflections +/// `H(i) = I - tau(i) * v * v.adjoint()`. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize-no-std", @@ -36,7 +39,7 @@ where DefaultAllocator: Allocator + Allocator>, { qr: OMatrix, - diag: OVector>, + tau: OVector>, } impl, C: Dim> Copy for QR @@ -49,7 +52,7 @@ where impl, C: Dim> QR where - DefaultAllocator: Allocator + Allocator + Allocator>, + DefaultAllocator: Allocator + Allocator + Allocator + Allocator>, { /// Computes the QR decomposition using householder reflections. pub fn new(mut matrix: OMatrix) -> Self { @@ -59,20 +62,67 @@ where if min_nrows_ncols.value() == 0 { return QR { qr: matrix, - diag: Matrix::zeros_generic(min_nrows_ncols, Const::<1>), + tau: Matrix::zeros_generic(min_nrows_ncols, Const::<1>), }; } - let mut diag = Matrix::uninit(min_nrows_ncols, Const::<1>); + let mut tau = Matrix::uninit(min_nrows_ncols, Const::<1>); + let mut work = Matrix::zeros_generic(ncols, Const::<1>); for i in 0..min_nrows_ncols.value() { - diag[i] = - MaybeUninit::new(householder::clear_column_unchecked(&mut matrix, i, 0, None)); + let (mut left, mut right) = matrix.columns_range_pair_mut(i, i + 1..); + let mut axis = left.rows_range_mut(i..); + + // Compute the scaled Householder vector, cf. LAPACK's `?LARFG`. + let (beta, tau_i) = { + let alpha = unsafe { axis.vget_unchecked(0).clone() }; + let xnorm = axis.rows_range(1..).norm(); + + if xnorm.is_zero() && alpha.clone().imaginary().is_zero() { + // The column is already in the wanted form. + (alpha, T::zero()) + } else { + let a_r = alpha.clone().real(); + let a_i = alpha.clone().imaginary(); + // TODO: use LAPACK's `?LAPY3` once `RealField` has a `max` method. + let reflection_norm = + (a_r.clone() * a_r.clone() + a_i.clone() * a_i + xnorm.clone() * xnorm) + .sqrt(); + // TODO: use `reflection_norm.copysign(a_r)`. + let beta = -reflection_norm.abs() * a_r.signum(); + // TODO: rescale if `beta` is close to underflow, cf. LAPACK's `?LARFG`. + let tau_i = (T::from_real(beta.clone()) - alpha.clone()).unscale(beta.clone()); + // Scale the Householder vector such that its first component is `1`. + let tmp = alpha - T::from_real(beta.clone()); + axis.rows_range_mut(1..).apply(|x| *x /= tmp.clone()); + + (T::from_real(beta), tau_i) + } + }; + + tau[i] = MaybeUninit::new(tau_i.clone()); + + if !tau_i.is_zero() { + // Apply the Householder reflection to the remaining columns. + unsafe { + *axis.vget_unchecked_mut(0) = T::one(); + } + + let mut work = work.rows_range_mut(i + 1..); + work.gemv_ad(T::one(), &right.rows_range(i..), &axis, T::zero()); + right + .rows_range_mut(i..) + .gerc(-tau_i.conjugate(), &axis, &work, T::one()); + } + + unsafe { + *axis.vget_unchecked_mut(0) = beta; + } } - // Safety: diag is now fully initialized. - let diag = unsafe { diag.assume_init() }; - QR { qr: matrix, diag } + // Safety: tau is now fully initialized. + let tau = unsafe { tau.assume_init() }; + QR { qr: matrix, tau } } /// Retrieves the upper trapezoidal submatrix `R` of this decomposition. @@ -83,9 +133,7 @@ where DefaultAllocator: Allocator, C>, { let (nrows, ncols) = self.qr.shape_generic(); - let mut res = self.qr.rows_generic(0, nrows.min(ncols)).upper_triangle(); - res.set_partial_diagonal(self.diag.iter().map(|e| T::from_real(e.clone().modulus()))); - res + self.qr.rows_generic(0, nrows.min(ncols)).upper_triangle() } /// Retrieves the upper trapezoidal submatrix `R` of this decomposition. @@ -99,33 +147,76 @@ where let (nrows, ncols) = self.qr.shape_generic(); let mut res = self.qr.resize_generic(nrows.min(ncols), ncols, T::zero()); res.fill_lower_triangle(T::zero(), 1); - res.set_partial_diagonal(self.diag.iter().map(|e| T::from_real(e.clone().modulus()))); res } - /// Computes the orthogonal matrix `Q` of this decomposition. + /// Computes the first `ncols` columns of the orthogonal matrix `Q` of this decomposition. + /// + /// Use this to get the full `Q` of a tall matrix: `q_columns` accepts any `ncols` up to the + /// number of rows of the decomposed matrix, while [`QR::q`] returns the first + /// `min(nrows, ncols)` columns only. + /// + /// # Panics + /// Panics if `ncols` is bigger than the number of rows of the decomposed matrix. #[must_use] - pub fn q(&self) -> OMatrix> + pub fn q_columns(&self, ncols: K) -> OMatrix where - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator + Allocator, { - let (nrows, ncols) = self.qr.shape_generic(); + // This is LAPACK's `?ORG2R`. + let (q_nrows, q_ncols) = self.qr.shape_generic(); + assert!( + ncols.value() <= q_nrows.value(), + "The number of columns of Q cannot be bigger than the number of rows of the decomposed matrix." + ); + + let mut a = OMatrix::::identity_generic(q_nrows, ncols); + let mut work = Matrix::zeros_generic(ncols, Const::<1>); + // The reflections after the first `k` ones do not change the first `ncols` columns. + let k = q_nrows.value().min(q_ncols.value()).min(ncols.value()); + + a.view_range_mut(.., ..k) + .copy_from(&self.qr.view_range(.., ..k)); + + for i in (0..k).rev() { + let tau_i = unsafe { self.tau.vget_unchecked(i).clone() }; - // NOTE: we could build the identity matrix and call q_mul on it. - // Instead we don't so that we take in account the matrix sparseness. - let mut res = Matrix::identity_generic(nrows, nrows.min(ncols)); - let dim = self.diag.len(); + if i + 1 < ncols.value() { + // Apply the reflection to the columns computed so far. + unsafe { + *a.get_unchecked_mut((i, i)) = T::one(); + } + + let (left, mut right) = a.columns_range_pair_mut(i, i + 1..); + let axis = left.rows_range(i..); + let mut work = work.rows_range_mut(i + 1..); + work.gemv_ad(T::one(), &right.rows_range(i..), &axis, T::zero()); + right + .rows_range_mut(i..) + .gerc(-tau_i.clone(), &axis, &work, T::one()); + } - for i in (0..dim).rev() { - let axis = self.qr.view_range(i.., i); - // TODO: sometimes, the axis might have a zero magnitude. - let refl = Reflection::new(Unit::new_unchecked(axis), T::zero()); + if i + 1 < q_nrows.value() { + a.view_range_mut(i + 1.., i).apply(|x| *x *= -tau_i.clone()); + } - let mut res_rows = res.view_range_mut(i.., i..); - refl.reflect_with_sign(&mut res_rows, self.diag[i].clone().signum()); + unsafe { + *a.get_unchecked_mut((i, i)) = T::one() - tau_i; + } + a.view_range_mut(..i, i).fill(T::zero()); } - res + a + } + + /// Computes the orthogonal matrix `Q` of this decomposition. + #[must_use] + pub fn q(&self) -> OMatrix> + where + DefaultAllocator: Allocator> + Allocator>, + { + let (nrows, ncols) = self.qr.shape_generic(); + self.q_columns(nrows.min(ncols)) } /// Unpacks this decomposition into its two matrix factors. @@ -137,8 +228,9 @@ where ) where DimMinimum: DimMin>, - DefaultAllocator: - Allocator> + Reallocator, C>, + DefaultAllocator: Allocator> + + Allocator> + + Reallocator, C>, { (self.q(), self.unpack_r()) } @@ -148,25 +240,33 @@ where &self.qr } - #[must_use] - pub(crate) const fn diag_internal(&self) -> &OVector> { - &self.diag - } - /// Multiplies the provided matrix by the transpose of the `Q` matrix of this decomposition. pub fn q_tr_mul(&self, rhs: &mut Matrix) - // TODO: do we need a static constraint on the number of rows of rhs? where S2: StorageMut, + ShapeConstraint: SameNumberOfRows, { - let dim = self.diag.len(); + for i in 0..self.tau.len() { + let tau_i = unsafe { self.tau.vget_unchecked(i).clone() }; + + if tau_i.is_zero() { + continue; + } + + // The first component of the Householder vector is `1` and is not stored. + let axis = self.qr.view_range(i + 1.., i); - for i in 0..dim { - let axis = self.qr.view_range(i.., i); - let refl = Reflection::new(Unit::new_unchecked(axis), T::zero()); + for j in 0..rhs.ncols() { + let mut col = rhs.column_mut(j); + let dot = + unsafe { col.vget_unchecked(i).clone() } + axis.dotc(&col.rows_range(i + 1..)); + let factor = -(tau_i.clone().conjugate() * dot); - let mut rhs_rows = rhs.rows_range_mut(i..); - refl.reflect_with_sign(&mut rhs_rows, self.diag[i].clone().signum().conjugate()); + unsafe { + *col.vget_unchecked_mut(i) += factor.clone(); + } + col.rows_range_mut(i + 1..).axpy(factor, &axis, T::one()); + } } } } @@ -217,42 +317,7 @@ where ); self.q_tr_mul(b); - self.solve_upper_triangular_mut(b) - } - - // TODO: duplicate code from the `solve` module. - fn solve_upper_triangular_mut( - &self, - b: &mut Matrix, - ) -> bool - where - S2: StorageMut, - ShapeConstraint: SameNumberOfRows, - { - let dim = self.qr.nrows(); - - for k in 0..b.ncols() { - let mut b = b.column_mut(k); - for i in (0..dim).rev() { - let coeff; - - unsafe { - let diag = self.diag.vget_unchecked(i).clone().modulus(); - - if diag.is_zero() { - return false; - } - - coeff = b.vget_unchecked(i).clone().unscale(diag); - *b.vget_unchecked_mut(i) = coeff.clone(); - } - - b.rows_range_mut(..i) - .axpy(-coeff, &self.qr.view_range(..i, i), T::one()); - } - } - - true + self.qr.solve_upper_triangular_mut(b) } /// Computes the inverse of the decomposed matrix. @@ -283,14 +348,7 @@ where self.qr.is_square(), "QR: unable to test the invertibility of a non-square matrix." ); - - for i in 0..self.diag.len() { - if self.diag[i].is_zero() { - return false; - } - } - - true + (0..self.qr.ncols()).all(|i| unsafe { !self.qr.get_unchecked((i, i)).is_zero() }) } // /// Computes the determinant of the decomposed matrix. diff --git a/src/linalg/svd3.rs b/src/linalg/svd3.rs index a8c39d282..72bb49c64 100644 --- a/src/linalg/svd3.rs +++ b/src/linalg/svd3.rs @@ -46,10 +46,26 @@ pub fn svd_ordered3( } let qr = b.qr(); + // The QR decomposition does not force the diagonal of `R` to be positive, so the singular + // values are the absolute values of that diagonal, and the columns of `U` must get the + // sign of the corresponding diagonal element. + let r_diagonal = qr.qr_internal().diagonal(); Some(SVD { - u: if compute_u { Some(qr.q()) } else { None }, - singular_values: qr.diag_internal().map(|e| e.abs()), + u: if compute_u { + let mut u = qr.q(); + + for i in 0..3 { + if r_diagonal[i] < T::zero() { + u.column_mut(i).neg_mut(); + } + } + + Some(u) + } else { + None + }, + singular_values: r_diagonal.map(|e| e.abs()), v_t: if compute_v { Some(v.transpose()) } else { None }, }) } diff --git a/tests/linalg/qr.rs b/tests/linalg/qr.rs index f499b030d..2760a05c7 100644 --- a/tests/linalg/qr.rs +++ b/tests/linalg/qr.rs @@ -1,9 +1,67 @@ #![cfg(feature = "proptest-support")] +use na::{Matrix2, Matrix4x2, U3, U4}; + +#[test] +fn simple_qr() { + #[rustfmt::skip] + let a = Matrix4x2::new( + -0.8943285241224914 , 0.12787800716234649, + -0.37320804072796987, 0.21338804264385058, + 0. , -0.2456767687354977 , + 0.2456767687354977 , 0. ); + let qr = a.qr(); + // The reference values were generated by converting the input + // to the form `m * 2 ^ e` for integers m and e. This was then used to + // obtain the QR decomposition without rounding errors. The result was + // converted back to f64. + #[rustfmt::skip] + let r_ref = Matrix2::new( + 0.99973237689865724, -0.19405501632841561, + 0. , -0.2908383860381578); + assert_relative_eq!(qr.r(), r_ref); + + #[rustfmt::skip] + let q_ref = Matrix4x2::new( + -0.89456793116659196, 0.15719172406996297, + -0.3733079465583837 , -0.48461884587835711, + 0. , 0.8447191998351451, + 0.24574253511487697, -0.1639658791740342); + assert_relative_eq!(qr.q(), q_ref); +} + +#[test] +fn q_columns() { + let a = Matrix4x2::new(0., 1., 3., 3., 1., 1., 2., 1.); + let qr = a.qr(); + + // The full `Q` is orthogonal, and its first columns are the ones of `q()`. + let q_full = qr.q_columns(U4); + assert!(q_full.is_orthogonal(1.0e-15)); + assert_relative_eq!(q_full.fixed_columns::<2>(0).into_owned(), qr.q()); + assert_relative_eq!(qr.q_columns(U3), q_full.fixed_columns::<3>(0).into_owned()); +} + +#[test] +#[should_panic] +fn q_columns_panic() { + let _ = Matrix2::::zeros().qr().q_columns(U3); +} + +#[test] +fn qr_zero_matrix() { + // All the reflections are the identity, so `Q` is the identity too. + let qr = Matrix4x2::::zeros().qr(); + + assert_eq!(qr.r(), Matrix2::zeros()); + assert!(qr.q_columns(U4).is_identity(1.0e-15)); + assert!(qr.q().is_orthogonal(1.0e-15)); +} + macro_rules! gen_tests( ($module: ident, $scalar: expr, $scalar_type: ty) => { mod $module { - use na::{DMatrix, DVector, Matrix4x3, Vector4}; + use na::{DMatrix, DVector, Dyn, Matrix4x3, Vector4}; use std::cmp; #[allow(unused_imports)] use crate::core::helper::{RandScalar, RandComplex}; @@ -17,8 +75,8 @@ macro_rules! gen_tests( let q = qr.q(); let r = qr.r(); - prop_assert!(relative_eq!(m, &q * r, epsilon = 1.0e-7)); - prop_assert!(q.is_orthogonal(1.0e-7)); + prop_assert!(relative_eq!(m, &q * r, epsilon = 1.0e-9)); + prop_assert!(q.is_orthogonal(1.0e-13)); } #[test] @@ -27,8 +85,8 @@ macro_rules! gen_tests( let q = qr.q(); let r = qr.r(); - prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-7)); - prop_assert!(q.is_orthogonal(1.0e-7)); + prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-9)); + prop_assert!(q.is_orthogonal(1.0e-13)); } #[test] @@ -37,8 +95,8 @@ macro_rules! gen_tests( let q = qr.q(); let r = qr.r(); - prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-7)); - prop_assert!(q.is_orthogonal(1.0e-7)); + prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-9)); + prop_assert!(q.is_orthogonal(1.0e-13)); } #[test] @@ -47,8 +105,29 @@ macro_rules! gen_tests( let q = qr.q(); let r = qr.r(); - prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-7)); - prop_assert!(q.is_orthogonal(1.0e-7)); + prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-9)); + prop_assert!(q.is_orthogonal(1.0e-13)); + } + + #[test] + fn qr_q_columns(m in dmatrix_($scalar)) { + let nrows = m.nrows(); + let q = m.qr().q_columns(Dyn(nrows)); + + prop_assert!(q.is_orthogonal(1.0e-13)); + } + + #[test] + fn qr_q_tr_mul(m in dmatrix_($scalar)) { + let nrows = m.nrows(); + let qr = m.qr(); + let b = DMatrix::<$scalar_type>::new_random(nrows, 3).map(|e| e.0); + + let mut q_tr_b = b.clone(); + qr.q_tr_mul(&mut q_tr_b); + + let q = qr.q_columns(Dyn(nrows)); + prop_assert!(relative_eq!(q_tr_b, q.adjoint() * b, epsilon = 1.0e-9)); } #[test]