Cap fixed-rank ID after QR/SVD nullification

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<Item::Real>,
         Standard: Distribution<Item::Real>,
+        Real<Item>: num::NumCast,
         LuDecomposition<Item, BaseArray<Item, VectorContainer<Item>, 2>>:
             MatrixLuDecomposition<Item = Item>,
         QrDecomposition<Item, BaseArray<Item, VectorContainer<Item>, 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::<Item>(), 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,
                 )

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<Item: RlstScalar>(
     }
 }
 
+pub fn fixed_rank_relative_tol<Item: RlstScalar>() -> Real<Item>
+where
+    Real<Item>: NumCast,
+{
+    if std::mem::size_of::<Real<Item>>() <= std::mem::size_of::<f32>() {
+        NumCast::from(1e-5f64).unwrap()
+    } else {
+        NumCast::from(1e-10f64).unwrap()
+    }
+}
+
 fn normal_equation_scale<Item: RlstScalar>(normal: &DynamicArray<Item, 2>) -> Real<Item> {
     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::<f32>(), 1e-5f32);
+        assert_eq!(fixed_rank_relative_tol::<f64>(), 1e-10f64);
+    }
 }