From 2e6ea52fe72b3c20b14ed68f977f7e16d387bf42 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 8 Oct 2026 20:38:54 +0900 Subject: [PATCH 1/2] fix: retain nonzero columns in unpivoted QR --- PATCHSET.md | 2 +- faer/src/linalg/qr/no_pivoting/factor.rs | 64 +++++++++++++++++++++--- 2 files changed, 58 insertions(+), 8 deletions(-) diff --git a/PATCHSET.md b/PATCHSET.md index c798cecb..61b4f92b 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -25,7 +25,7 @@ upstream `faer`. | Experiment | What | Status | Issues | Notes | | --- | --- | --- | --- | --- | -| Unpivoted QR column skip | Householder QR drops a column whose part orthogonal to the previous ones is below `16·(m−k)·ε`. The tall/wide `svd` path goes through it too. LAPACK has no such skip. | TODO | tlinalg-rs#28, tensor4all-rs#836 | Likely an upstream bug. Any upstream issue or PR is opened only after the fork maintainer confirms. | +| Unpivoted QR column skip | Householder QR now applies the reflector and advances the row for every nonzero column, matching LAPACK; the tall/wide `svd` path goes through it too. | KEEP | tlinalg-rs#28, tensor4all-rs#836 | Adopted. Test: `linalg::qr::no_pivoting::factor::tests::test_near_dependent_columns_are_not_skipped` at m=100, 1000, and 8192. Existing columns are unchanged; formerly skipped nonzero columns now pay reflector/update flops (the threshold norm bookkeeping is removed). No existing benchmark suite was run. Release smoke, 3-run mean, before→after: 100000×8 2.984→2.986 ms; 1000×1000 19.793→19.921 ms. | | Typed uninitialized overwrite destination | A `MatUninitMut`-style destination that never reads the prior output, accepted by `matmul_with_conj`. | TODO | strided-rs#198 | strided-rs currently returns a Faer-only `Unsupported`. | | Native SVD / eigh direct output | At nbatch=1, stop writing U and the eigenvectors into an intermediate `Mat` and copying them out; write into the caller's output region instead. | TODO | tlinalg-rs#25 | Preserve output shape, initialization, and failure cleanup. | | Borrowed reflector columns | Stop copying reflectors out of the compact QR buffer into `basis` scratch; borrow them through a column split. | TODO | tlinalg-rs#17 | Numerics, the tau convention, and the restored diagonal beta stay identical. | diff --git a/faer/src/linalg/qr/no_pivoting/factor.rs b/faer/src/linalg/qr/no_pivoting/factor.rs index d6d23408..236b8a3e 100644 --- a/faer/src/linalg/qr/no_pivoting/factor.rs +++ b/faer/src/linalg/qr/no_pivoting/factor.rs @@ -22,7 +22,6 @@ fn qr_in_place_unblocked( let mut row = row_start; while row < Ord::min(size, m) && col < n { - let norm = A.rb().col(col).get(..row).norm_l2(); let (mut A00, A01, A10, A11) = A.rb_mut().split_at_mut(row + 1, col + 1); @@ -49,11 +48,6 @@ fn qr_in_place_unblocked( (info, A10l.rb().col(row)) }; - let norm = &info.norm.hypot(norm); - let eps = &eps::(); - let leeway = &from_f64::((m - row) as f64 * 16.0); - let threshold = eps * leeway * norm; - let tau_inv = &info.tau.recip(); H[row] = info.tau.to_cplx(); @@ -61,7 +55,7 @@ fn qr_in_place_unblocked( if info.norm > zero() { row += 1; } - } else if info.norm > threshold { + } else if info.norm > zero() { for (head, tail) in core::iter::zip(A01.iter_mut(), A11.col_iter_mut()) { @@ -536,6 +530,62 @@ mod tests { } } } + #[test] + fn test_near_dependent_columns_are_not_skipped() { + for m in [100, 1000, 8192] { + let A = Mat::from_fn(m, 3, |i, j| match j { + 0 => 1.0, + 1 => 1.0 + if i == 0 { 1e-10 } else { 0.0 }, + 2 => 1.0 + if i == 1 { 1e-6 } else { 0.0 }, + _ => unreachable!(), + }); + let block_size = recommended_block_size::(m, 3); + let mut QR = A.cloned(); + let mut H = Mat::zeros(block_size, 3); + qr_in_place( + QR.as_mut(), + H.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new(qr_in_place_scratch::( + m, + 3, + block_size, + Par::Seq, + default(), + ))), + default(), + ); + + let mut Q = Mat::::zeros(m, m); + let mut R = QR.as_ref().cloned(); + for j in 0..m { + Q[(j, j)] = 1.0; + } + for j in 0..3 { + for i in j + 1..m { + R[(i, j)] = 0.0; + } + } + householder::apply_block_householder_sequence_on_the_left_in_place_with_conj( + QR.as_ref(), + H.as_ref(), + Conj::No, + Q.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new( + householder::apply_block_householder_sequence_on_the_left_in_place_scratch::( + m, block_size, m, + ), + )), + ); + let relative_error = (&A - &Q * &R).norm_l2() / A.norm_l2(); + assert!( + relative_error <= 1e-13, + "m={m}: relative reconstruction error {relative_error:e}" + ); + } + } + #[test] fn test_rank_deficient() { let i = c64::new(0.0, 1.0); From d16929be4582c6b3679f5ec17223227bbfc289c9 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 8 Oct 2026 20:53:20 +0900 Subject: [PATCH 2/2] test(qr): cover the wide fixture; the skip has no second site The unpivoted QR column skip lived only in `qr_in_place_unblocked`. `qr_in_place_blocked` has no threshold of its own and recurses into that core, so removing the threshold there already covers the wide path; there is no wide-specific branch. Add the transposed fixture as a wide `3 x n` regression test (n = 100/1000/8192). That wide fixture reconstructs to ~4e-16 both with and without the change: at m = 3 the removed threshold `16 (m-k) eps ||column||` is ~1e-14 and never tripped. The tall m = 1000 case goes 2.58e-12 -> 1.98e-15. tlinalg's reported wide `3.33e-7` for this fixture is a reconstruction- indexing bug in that test: it writes entry (i, j) to flat `j + 3i` while the column-major factor uses `i + 3j`. faer's factors for the matrix are accurate to ~4e-16 under the correct indexing. --- PATCHSET.md | 2 +- faer/src/linalg/qr/no_pivoting/factor.rs | 110 +++++++++++++++-------- 2 files changed, 73 insertions(+), 39 deletions(-) diff --git a/PATCHSET.md b/PATCHSET.md index 61b4f92b..ee6f76d9 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -25,7 +25,7 @@ upstream `faer`. | Experiment | What | Status | Issues | Notes | | --- | --- | --- | --- | --- | -| Unpivoted QR column skip | Householder QR now applies the reflector and advances the row for every nonzero column, matching LAPACK; the tall/wide `svd` path goes through it too. | KEEP | tlinalg-rs#28, tensor4all-rs#836 | Adopted. Test: `linalg::qr::no_pivoting::factor::tests::test_near_dependent_columns_are_not_skipped` at m=100, 1000, and 8192. Existing columns are unchanged; formerly skipped nonzero columns now pay reflector/update flops (the threshold norm bookkeeping is removed). No existing benchmark suite was run. Release smoke, 3-run mean, before→after: 100000×8 2.984→2.986 ms; 1000×1000 19.793→19.921 ms. | +| Unpivoted QR column skip | The skip lived only in `qr_in_place_unblocked`; `qr_in_place_blocked` has no threshold of its own and recurses into that core, so removing it there covers the tall and wide paths at once. The reflector is applied and the row advances for every nonzero column, matching LAPACK; the tall/wide `svd` path goes through it too. | KEEP | tlinalg-rs#28, tensor4all-rs#836 | Adopted. Tests: `test_near_dependent_columns_are_not_skipped` (tall `m x 3`, m=100/1000/8192) and `test_near_dependent_rows_are_not_skipped` (its transpose, wide `3 x n`, n=100/1000/8192), both at `<=1e-13`; measured 1e-16 to 5e-15. The wide fixture was already accurate before the change: at m=3 the removed threshold `16 (m-k) eps ||column||` is ~1e-14 and never tripped it (a wide skip's relative reconstruction error is ~`16 eps sqrt(m) sqrt(m/n)`), so the wide path needed no separate treatment — the same core change covers it. tlinalg's reported wide `3.33e-7` is a reconstruction-indexing bug in that test (entry `(i,j)` written to flat `j+3i` while the factor is column-major `i+3j`); faer's factors for that matrix are accurate to ~4e-16. Existing columns are unchanged; formerly skipped nonzero columns now pay reflector/update flops (the threshold norm bookkeeping is removed). No existing benchmark suite was run. Release smoke, 3-run mean, before→after: 100000×8 2.984→2.986 ms; 1000×1000 19.793→19.921 ms. | | Typed uninitialized overwrite destination | A `MatUninitMut`-style destination that never reads the prior output, accepted by `matmul_with_conj`. | TODO | strided-rs#198 | strided-rs currently returns a Faer-only `Unsupported`. | | Native SVD / eigh direct output | At nbatch=1, stop writing U and the eigenvectors into an intermediate `Mat` and copying them out; write into the caller's output region instead. | TODO | tlinalg-rs#25 | Preserve output shape, initialization, and failure cleanup. | | Borrowed reflector columns | Stop copying reflectors out of the compact QR buffer into `basis` scratch; borrow them through a column split. | TODO | tlinalg-rs#17 | Numerics, the tau convention, and the restored diagonal beta stay identical. | diff --git a/faer/src/linalg/qr/no_pivoting/factor.rs b/faer/src/linalg/qr/no_pivoting/factor.rs index 236b8a3e..b4ae02f3 100644 --- a/faer/src/linalg/qr/no_pivoting/factor.rs +++ b/faer/src/linalg/qr/no_pivoting/factor.rs @@ -530,6 +530,52 @@ mod tests { } } } + /// Factor `A`, rebuild `Q` and `R`, and return `||A - QR|| / ||A||`. + fn qr_reconstruction_error(A: MatRef<'_, f64>, block_size: usize) -> f64 { + let (m, n) = (A.nrows(), A.ncols()); + let k = Ord::min(m, n); + let mut QR = A.cloned(); + let mut H = Mat::zeros(block_size, k); + qr_in_place( + QR.as_mut(), + H.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new(qr_in_place_scratch::( + m, + n, + block_size, + Par::Seq, + default(), + ))), + default(), + ); + + let mut Q = Mat::::zeros(m, m); + let mut R = QR.as_ref().cloned(); + for j in 0..m { + Q[(j, j)] = 1.0; + } + for j in 0..k { + for i in j + 1..m { + R[(i, j)] = 0.0; + } + } + householder::apply_block_householder_sequence_on_the_left_in_place_with_conj( + QR.as_ref().subcols(0, k), + H.as_ref(), + Conj::No, + Q.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new( + householder::apply_block_householder_sequence_on_the_left_in_place_scratch::( + m, block_size, m, + ), + )), + ); + (&A - &Q * &R).norm_l2() / A.norm_l2() + } + + /// `A = [a, a + 1e-10 e_0, a + 1e-6 e_1]`, `a = (1, ..., 1)`, tall `m x 3`. #[test] fn test_near_dependent_columns_are_not_skipped() { for m in [100, 1000, 8192] { @@ -539,49 +585,37 @@ mod tests { 2 => 1.0 + if i == 1 { 1e-6 } else { 0.0 }, _ => unreachable!(), }); - let block_size = recommended_block_size::(m, 3); - let mut QR = A.cloned(); - let mut H = Mat::zeros(block_size, 3); - qr_in_place( - QR.as_mut(), - H.as_mut(), - Par::Seq, - MemStack::new(&mut MemBuffer::new(qr_in_place_scratch::( - m, - 3, - block_size, - Par::Seq, - default(), - ))), - default(), + let relative_error = qr_reconstruction_error( + A.as_ref(), + recommended_block_size::(m, 3), ); + assert!( + relative_error <= 1e-13, + "m={m}: relative reconstruction error {relative_error:e}" + ); + } + } - let mut Q = Mat::::zeros(m, m); - let mut R = QR.as_ref().cloned(); - for j in 0..m { - Q[(j, j)] = 1.0; - } - for j in 0..3 { - for i in j + 1..m { - R[(i, j)] = 0.0; - } - } - householder::apply_block_householder_sequence_on_the_left_in_place_with_conj( - QR.as_ref(), - H.as_ref(), - Conj::No, - Q.as_mut(), - Par::Seq, - MemStack::new(&mut MemBuffer::new( - householder::apply_block_householder_sequence_on_the_left_in_place_scratch::( - m, block_size, m, - ), - )), + /// The same fixture transposed, so the nearly dependent columns are the + /// rows: wide `3 x n`. This is the shape that a caller sees as "wide", and + /// it exercises the blocked path's recursion into the unblocked core. + #[test] + fn test_near_dependent_rows_are_not_skipped() { + for n in [100, 1000, 8192] { + // W = A^T, i.e. W[(i, j)] = A[(j, i)] for the fixture above. + let W = Mat::from_fn(3, n, |i, j| match i { + 0 => 1.0, + 1 => 1.0 + if j == 0 { 1e-10 } else { 0.0 }, + 2 => 1.0 + if j == 1 { 1e-6 } else { 0.0 }, + _ => unreachable!(), + }); + let relative_error = qr_reconstruction_error( + W.as_ref(), + recommended_block_size::(3, n), ); - let relative_error = (&A - &Q * &R).norm_l2() / A.norm_l2(); assert!( relative_error <= 1e-13, - "m={m}: relative reconstruction error {relative_error:e}" + "3x{n}: relative reconstruction error {relative_error:e}" ); } }