diff --git a/PATCHSET.md b/PATCHSET.md index c798cecb..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 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 | 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 d6d23408..b4ae02f3 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,96 @@ 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] { + 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 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}" + ); + } + } + + /// 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), + ); + assert!( + relative_error <= 1e-13, + "3x{n}: relative reconstruction error {relative_error:e}" + ); + } + } + #[test] fn test_rank_deficient() { let i = c64::new(0.0, 1.0);