diff --git a/src/rsrs/rsrs_factors/commutative_factors.rs b/src/rsrs/rsrs_factors/commutative_factors.rs index 595460e..a020196 100644 --- a/src/rsrs/rsrs_factors/commutative_factors.rs +++ b/src/rsrs/rsrs_factors/commutative_factors.rs @@ -19,8 +19,8 @@ use crate::utils::{ row_perm, row_subs, }, linear_algebra::{ - block_extraction_into, streaming_chunk_rows, BlockExtractionMethod, - NormalEquationAccumulator, + block_extraction_into, fixed_rank_relative_tol, streaming_chunk_rows, + BlockExtractionMethod, NormalEquationAccumulator, NullMethod, }, memory::{matrix_bytes, trace_memory_event, trace_memory_growth}, }; @@ -536,6 +536,7 @@ where where StandardNormal: Distribution, Standard: Distribution, + Real: num::NumCast, LuDecomposition, 2>>: MatrixLuDecomposition, QrDecomposition, 2>>: @@ -582,6 +583,12 @@ where }; return (None, Times::Id(id_times)); } + let fixed_rank_accuracy = |rank: usize| match id_options.null_method { + NullMethod::Svd | NullMethod::Qr => { + Accuracy::MinRank(fixed_rank_relative_tol::(), rank) + } + NullMethod::Projection => Accuracy::FixedRank(rank), + }; let id_sketch = match rank_par { BoxType::Full(tol) => { // for a box that hasn't been merged yet it @@ -605,7 +612,7 @@ where .r_mut() .into_subview([0, 0], null_shape) .into_id_alloc_no_skel( - Accuracy::FixedRank(loc_rank), + fixed_rank_accuracy(loc_rank), id_options.qr_method.clone(), TransMode::Trans, ) @@ -619,7 +626,7 @@ where .r_mut() .into_subview([0, 0], null_shape) .into_id_alloc_no_skel( - Accuracy::FixedRank((*rank).min(max_rank)), + fixed_rank_accuracy((*rank).min(max_rank)), id_options.qr_method.clone(), TransMode::Trans, ) diff --git a/src/utils/linear_algebra.rs b/src/utils/linear_algebra.rs index e8c0e84..42a4528 100644 --- a/src/utils/linear_algebra.rs +++ b/src/utils/linear_algebra.rs @@ -1,4 +1,4 @@ -use num::{Float, One, Zero}; +use num::{Float, NumCast, One, Zero}; use rlst::dense::linalg::{lu::MatrixLu, null_space::Method}; pub use rlst::prelude::*; use serde::Deserialize; @@ -93,6 +93,17 @@ pub fn add_diagonal( } } +pub fn fixed_rank_relative_tol() -> Real +where + Real: NumCast, +{ + if std::mem::size_of::>() <= std::mem::size_of::() { + NumCast::from(1e-5f64).unwrap() + } else { + NumCast::from(1e-10f64).unwrap() + } +} + fn normal_equation_scale(normal: &DynamicArray) -> Real { let shape = normal.shape(); let view = normal.r(); @@ -421,4 +432,10 @@ mod tests { .iter() .all(|entry| Float::abs(*entry) <= f64::EPSILON)); } + + #[test] + fn fixed_rank_relative_tol_depends_on_precision() { + assert_eq!(fixed_rank_relative_tol::(), 1e-5f32); + assert_eq!(fixed_rank_relative_tol::(), 1e-10f64); + } }