From 2f51739405a9ff422c269b318e9a6196df784a24 Mon Sep 17 00:00:00 2001 From: fishidaho Date: Thu, 17 Sep 2026 15:46:32 -0700 Subject: [PATCH 1/2] Stop a column with no variance amplifying its own rounding noise `test_recipe_matches_reference[scanpy-VCSRArray]` has been failing on main: a 6x1 all-ones matrix came back as 1.05e-08 where the dense reference said 1.0. Neither number is right. Both are rounding noise multiplied by the reciprocal of rounding noise, and the two implementations simply amplified different noise. Depth normalization sends every row of a one-column matrix to the same value, so `g` is constant and the column's true variance is 0. `variance` is computed one-pass as `E[x^2] - E[x]^2`, which here subtracts two numbers of size `mean ** 2` that cancel almost entirely, leaving noise of order `n_rows * eps * mean ** 2` rather than exactly 0. `std > 0.0` cannot tell that from signal, so it took the dividing branch and scaled the column by `1 / 1.37e-07`. It is an ill-conditioned input rather than a wrong formula on either side, so widening the test's tolerance would only have hidden it: the value being compared is arbitrary, and at 5 rows the same input already gave 0. Judge "no variance" against the noise floor instead, which scales both with the column's magnitude and with how many terms were summed to reach it (`sqrt(n_rows * eps) * |mean|`). A constant column now centers to 0, which is what centering it means and what scanpy does with a zero-variance gene. An all-zero column has `mean == 0`, so its threshold stays 0 and its behaviour is unchanged. The two dense references carried the same `std > 0` test and so produced the same amplified noise; they now share the convention, still computing it independently through `numpy`. Measured: constant columns agree with the reference across 216 combinations of row count, magnitude and recipe, all centering to 0. Columns with real variance are untouched -- 20x5 through 1000x40 agree to 1.7e-13 against values up to 6.2. The threshold is set by the one-pass formula, not by the data. A two-pass or Welford variance would leave noise of order `eps` instead of `n_rows * eps` and admit a tighter one, at the cost of a second pass over `nnz`; that is a larger change than this fix needs. Co-Authored-By: Claude Opus 5 (1M context) --- src/vsparse/_norm_common.py | 33 +++++++++++++++++++++++- tests/test_property_normalization.py | 7 ++++- tests/test_vcs_norm_recipes.py | 38 +++++++++++++++++++++++++++- 3 files changed, 75 insertions(+), 3 deletions(-) diff --git a/src/vsparse/_norm_common.py b/src/vsparse/_norm_common.py index 4d5a78f..8b07475 100644 --- a/src/vsparse/_norm_common.py +++ b/src/vsparse/_norm_common.py @@ -608,6 +608,35 @@ def _compute_row_scale(arr: Any, recipe: Recipe) -> np.ndarray: return row_scale +def _zero_variance_tol(mean: np.ndarray, n_rows: int) -> np.ndarray: + """Below what ``std`` a column counts as having no variance at all. + + ``variance`` above is the one-pass ``E[x^2] - E[x]^2``. For a column whose + values are all (nearly) equal, that subtracts two numbers of size + ``mean ** 2`` which cancel almost completely, and what survives is rounding + noise rather than signal. Summing ``n_rows`` terms accumulates a relative + error of order ``n_rows * eps``, so the noise left in ``variance`` is of + order ``n_rows * eps * mean ** 2`` and the noise in ``std`` of order + ``sqrt(n_rows * eps) * |mean|``. + + A plain ``std > 0`` test cannot see that, and dividing by such a ``std`` + amplifies noise into the output instead of scaling anything: a constant + column comes back as an arbitrary O(1) value rather than the 0 that + centering it should give. So judge "no variance" against that noise floor, + which scales both with the column's magnitude and with how many terms were + summed to reach it. + + Columns that are genuinely all zero have ``mean == 0``, hence a threshold + of 0, and still take the ``std > 0`` branch. + + The floor is set by the one-pass formula, not by the data: a stable + two-pass or Welford variance would leave noise of order ``eps`` rather than + ``n_rows * eps`` and would admit a tighter threshold, at the cost of a + second pass over ``nnz``. + """ + return np.sqrt(n_rows * np.finfo(np.float64).eps) * np.abs(mean) + + class NormalizedViewBase: """Shared implementation for the normalized VCSC/VCSR views. @@ -708,7 +737,9 @@ def __init__( col_mean = mean if self.recipe.center else np.zeros(n_cols, dtype=np.float64) if self.recipe.post_scale: with np.errstate(divide="ignore", invalid="ignore"): - col_post_scale = np.where(std > 0.0, 1.0 / std, 1.0) + col_post_scale = np.where( + std > _zero_variance_tol(mean, n_rows), 1.0 / std, 1.0 + ) else: col_post_scale = np.ones(n_cols, dtype=np.float64) else: diff --git a/tests/test_property_normalization.py b/tests/test_property_normalization.py index a7db114..f6d44dd 100644 --- a/tests/test_property_normalization.py +++ b/tests/test_property_normalization.py @@ -83,7 +83,12 @@ def _recipe_reference(dense: np.ndarray, recipe: str) -> np.ndarray: if recipe in ("scanpy", "pearson"): std = g.std(axis=0) with np.errstate(divide="ignore", invalid="ignore"): - s = np.where(std > 0, 1.0 / std, 1.0) + # A column with no real variance leaves only rounding noise in + # `std`; dividing by it amplifies noise instead of scaling. Judge + # "no variance" relative to the column's own magnitude, as + # `_zero_variance_tol` does in the library. + tol = np.sqrt(dense.shape[0] * np.finfo(np.float64).eps) * np.abs(g.mean(axis=0)) + s = np.where(std > tol, 1.0 / std, 1.0) else: s = np.ones(dense.shape[1]) diff --git a/tests/test_vcs_norm_recipes.py b/tests/test_vcs_norm_recipes.py index e9a4a3e..b366faa 100644 --- a/tests/test_vcs_norm_recipes.py +++ b/tests/test_vcs_norm_recipes.py @@ -66,7 +66,12 @@ def _reference(dense: np.ndarray, recipe: str) -> np.ndarray: if recipe in ("scanpy", "pearson"): std = g.std(axis=0) with np.errstate(divide="ignore", invalid="ignore"): - s = np.where(std > 0, 1.0 / std, 1.0) + # A column with no real variance leaves only rounding noise in + # `std`; dividing by it amplifies noise instead of scaling. Judge + # "no variance" relative to the column's own magnitude, as + # `_zero_variance_tol` does in the library. + tol = np.sqrt(dense.shape[0] * np.finfo(np.float64).eps) * np.abs(g.mean(axis=0)) + s = np.where(std > tol, 1.0 / std, 1.0) else: s = np.ones(dense.shape[1]) @@ -356,3 +361,34 @@ def test_anndata_cache_does_not_pin_a_dropped_view(): assert ref() is None # Still reusable -- from obs/varm/uns if not from the retained statistics. assert adata.normalized("scanpy", recalculate=False) is not None + + +@pytest.mark.parametrize("n_rows", [2, 5, 6, 17, 64, 501]) +@pytest.mark.parametrize("value", [1.0, 7.0, 9999.0]) +@pytest.mark.parametrize("recipe", ["scanpy", "pearson"]) +def test_a_column_with_no_variance_centers_to_zero(vcls, n_rows, value, recipe): + """A constant column has no variance to scale by, so centering must leave 0. + + ``variance`` is computed one-pass as ``E[x^2] - E[x]^2``, which for such a + column cancels to rounding noise rather than to exactly 0. Testing + ``std > 0`` would therefore take the dividing branch and amplify that noise + into an arbitrary O(1) value -- 1.05e-08 for the 6x1 case below, against + the 0 the same input gives at 5 rows. The scale of the noise moves with + the column's magnitude and with how many terms were summed, so the + threshold has to move with both. + """ + dense = np.full((n_rows, 1), value) + v = vcls.from_scipy(_scipy_for(vcls, dense)) + out = v.normalized(recipe).toarray() + np.testing.assert_allclose(out, np.zeros_like(dense), atol=1e-12) + + +def test_a_constant_column_does_not_suppress_its_neighbours(vcls): + """Zeroing a no-variance column must not touch the columns beside it.""" + rng = np.random.default_rng(0) + dense = rng.integers(1, 50, size=(40, 5)).astype(float) + dense[:, 2] = 4.0 + v = vcls.from_scipy(_scipy_for(vcls, dense)) + out = v.normalized("pearson").toarray() + varying = np.delete(out, 2, axis=1) + assert np.abs(varying).max() > 0.5 From 53afb3e79fb05a9e05e539701b9ca8e95682690c Mon Sep 17 00:00:00 2001 From: fishidaho Date: Thu, 17 Sep 2026 17:03:17 -0700 Subject: [PATCH 2/2] Compute the column variance two-pass, and shrink the constant bound to fit The previous commit stopped a column with no variance amplifying its own rounding noise, but did it by widening the threshold rather than by fixing what produced the noise. That left the threshold sized by the arithmetic instead of by the data, and scaling as `sqrt(n_rows)`: at a million cells it called any column with a coefficient of variation under 1.7e-05 constant, against a smallest real signal here of about 1e-04 -- a margin of 6x, and 2x at ten million. Fine today, eroding exactly as the cohorts grow. Compute the variance about the mean instead, as the corrected two-pass algorithm (Chan, Golub & LeVeque), which is what scikit-learn's sparse `_csr_mean_variance_axis0` does -- the code scanpy calls for this. A column with no spread now has every deviation exactly 0, so its variance is 0 rather than the residue of subtracting two numbers of size `mean ** 2`. The implicit zeros never reach a kernel and never need to: each sits at `g(0) == 0`, so its deviation is exactly `-mean` and its contribution is closed-form. `_finish_variance` folds them in for both formats, which is also what keeps the two layouts bit-comparable -- verified to 8e-15. The second pass costs what the layout makes it cost. VCSC already knows a column's mean once its own nonzeros are summed, so the deviation loop fuses into the same `prange` body; VCSR cannot, since a column's mean is not final until every row has been scattered, so it gets a genuine second scatter. Measured on 4M nonzeros, per view and then cached: recipe VCSC VCSR scanpy 7.0 -> 10.5 ms 5.3 -> 8.9 ms pearson 5.3 -> 6.6 ms 3.7 -> 6.1 ms parafac2 8.1 -> 11.5 ms 6.7 -> 10.9 ms cp10k_log1p unchanged unchanged raw unchanged unchanged Only the recipes that need the variance pay, and the cost is re-evaluating `g` rather than re-reading memory. Extrapolated to a 2.3B-nonzero cohort that is roughly two seconds, once, against a BiCV run measured in minutes. With the error down at `eps` rather than `sqrt(eps)`, the threshold becomes scikit-learn's `_is_constant_feature` bound, and stops growing with the data: n_rows was now 6 3.65e-08 1.33e-15 10,000 1.49e-06 2.22e-12 1,300,000 1.70e-05 2.89e-10 10,000,000 4.71e-05 2.22e-09 Detecting constant columns is still needed, and is not what this replaces: a column with no spread has no unit-variance scaling in exact arithmetic either. The two are separate jobs, and only one of them was ever numerical. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01Rvu3cf8ZL7F5EPX22eo6Je --- src/vsparse/_norm_common.py | 161 +++++++++++++++++---------- tests/test_property_normalization.py | 13 ++- tests/test_vcs_norm_recipes.py | 53 ++++++--- 3 files changed, 146 insertions(+), 81 deletions(-) diff --git a/src/vsparse/_norm_common.py b/src/vsparse/_norm_common.py index 8b07475..97b2d1a 100644 --- a/src/vsparse/_norm_common.py +++ b/src/vsparse/_norm_common.py @@ -234,20 +234,20 @@ def _g_np(x: np.ndarray, g_code: int) -> np.ndarray: @numba.njit(cache=True, parallel=True) def _column_stats_major_is_col( - major_ptr, values, value_ptr, indices, row_scale, need_b, need_gstats, g_code + major_ptr, values, value_ptr, indices, row_scale, need_b, need_gstats, g_code, n_rows ): - """Per-column ``gsum`` (raw material for ``b``) and sum/sum-of-squares of ``g(scaled)``. + """Per-column ``gsum`` (raw material for ``b``) and the ``g(scaled)`` variance inputs. - Fused into one pass per column (rather than two separate dispatches): - unlike the VCSR scatter passes below, a VCSC column's ``gsum`` depends - only on that column's own nonzeros, so it's already final by the time - the second (``g``-transform) loop over the same nonzeros needs it -- - no need to wait for every other column to finish first. + Sums cover the stored entries only; :func:`_finish_variance` folds in the + implicit zeros. A VCSC column's statistics depend only on its own + nonzeros, so all three loops fuse into one pass per column. """ n_major = major_ptr.shape[0] - 1 gsum = np.ones(n_major, dtype=np.float64) col_sum = np.zeros(n_major, dtype=np.float64) - col_sumsq = np.zeros(n_major, dtype=np.float64) + col_m2 = np.zeros(n_major, dtype=np.float64) + col_corr = np.zeros(n_major, dtype=np.float64) + col_nnz = np.zeros(n_major, dtype=np.float64) for j in numba.prange(n_major): # ty: ignore[not-iterable] gs = 1.0 if need_b: @@ -259,17 +259,30 @@ def _column_stats_major_is_col( gsum[j] = gs if need_gstats and gs > 0.0: s0 = 0.0 - s1 = 0.0 + count = 0.0 for u in range(major_ptr[j], major_ptr[j + 1]): v = values[u] for k in range(value_ptr[u], value_ptr[u + 1]): - scaled = v / row_scale[indices[k]] / gs - gy = _g(scaled, g_code) - s0 += gy - s1 += gy * gy + s0 += _g(v / row_scale[indices[k]] / gs, g_code) + count += 1.0 col_sum[j] = s0 - col_sumsq[j] = s1 - return gsum, col_sum, col_sumsq + col_nnz[j] = count + # Deviations about the mean, not `sumsq - mean ** 2`: for a + # column with no spread every deviation is 0 exactly, where the + # latter subtracts two numbers of size `mean ** 2` and keeps only + # their rounding error. + mean = s0 / n_rows if n_rows > 0 else 0.0 + m2 = 0.0 + corr = 0.0 + for u in range(major_ptr[j], major_ptr[j + 1]): + v = values[u] + for k in range(value_ptr[u], value_ptr[u + 1]): + d = _g(v / row_scale[indices[k]] / gs, g_code) - mean + m2 += d * d + corr += d + col_m2[j] = m2 + col_corr[j] = corr + return gsum, col_sum, col_m2, col_corr, col_nnz # -- statistics: major=rows -- scatter-add passes ---------------------------- @@ -297,15 +310,16 @@ def _scaled_col_sums_vcs(major_ptr, values, value_ptr, indices, row_scale, n_col def _gstats_col_sums_vcs( major_ptr, values, value_ptr, indices, row_scale, gene_scale, g_code, n_cols, nthreads ): + """Per-column sum of ``g(scaled)`` over stored entries, and how many there were.""" n_major = major_ptr.shape[0] - 1 chunk = (n_major + nthreads - 1) // nthreads partial_sum = np.zeros((nthreads, n_cols), dtype=np.float64) - partial_sumsq = np.zeros((nthreads, n_cols), dtype=np.float64) + partial_nnz = np.zeros((nthreads, n_cols), dtype=np.float64) for t in numba.prange(nthreads): # ty: ignore[not-iterable] start = t * chunk end = min(n_major, start + chunk) local_sum = partial_sum[t] - local_sumsq = partial_sumsq[t] + local_nnz = partial_nnz[t] for i in range(start, end): rs = row_scale[i] for u in range(major_ptr[i], major_ptr[i + 1]): @@ -314,11 +328,57 @@ def _gstats_col_sums_vcs( c = indices[k] gs = gene_scale[c] if gs > 0.0: - scaled = v / rs / gs - gy = _g(scaled, g_code) - local_sum[c] += gy - local_sumsq[c] += gy * gy - return partial_sum.sum(axis=0), partial_sumsq.sum(axis=0) + local_sum[c] += _g(v / rs / gs, g_code) + local_nnz[c] += 1.0 + return partial_sum.sum(axis=0), partial_nnz.sum(axis=0) + + +@numba.njit(cache=True, parallel=True) +def _gstats_col_deviations_vcs( + major_ptr, values, value_ptr, indices, row_scale, gene_scale, g_code, col_mean, n_cols, nthreads +): + """Per-column ``sum(y - mean)`` and ``sum((y - mean) ** 2)`` over stored entries. + + A second walk of the values: on a row-major layout a column's mean is not + final until every row has been scattered into it. + """ + n_major = major_ptr.shape[0] - 1 + chunk = (n_major + nthreads - 1) // nthreads + partial_m2 = np.zeros((nthreads, n_cols), dtype=np.float64) + partial_corr = np.zeros((nthreads, n_cols), dtype=np.float64) + for t in numba.prange(nthreads): # ty: ignore[not-iterable] + start = t * chunk + end = min(n_major, start + chunk) + local_m2 = partial_m2[t] + local_corr = partial_corr[t] + for i in range(start, end): + rs = row_scale[i] + for u in range(major_ptr[i], major_ptr[i + 1]): + v = values[u] + for k in range(value_ptr[u], value_ptr[u + 1]): + c = indices[k] + gs = gene_scale[c] + if gs > 0.0: + d = _g(v / rs / gs, g_code) - col_mean[c] + local_m2[c] += d * d + local_corr[c] += d + return partial_m2.sum(axis=0), partial_corr.sum(axis=0) + + +def _finish_variance(col_sum, m2_stored, corr_stored, col_nnz, n_rows): + """Per-column mean and variance, given either layout's stored-only sums.""" + if n_rows <= 0: + zeros = np.zeros_like(col_sum) + return zeros, zeros + mean = col_sum / n_rows + # Implicit zeros sit at `g(0) == 0`, so each deviates by exactly `-mean` + # and their contribution is closed-form rather than iterated. + n_zero = n_rows - col_nnz + m2 = m2_stored + n_zero * mean**2 + corr = corr_stored - n_zero * mean + # `corr ** 2 / n` is the corrected two-pass term (Chan, Golub & LeVeque): + # the deviations are about the computed mean, not the exact one. + return mean, np.clip((m2 - corr**2 / n_rows) / n_rows, 0.0, None) # -- full materialization ----------------------------------------------------- @@ -608,33 +668,15 @@ def _compute_row_scale(arr: Any, recipe: Recipe) -> np.ndarray: return row_scale -def _zero_variance_tol(mean: np.ndarray, n_rows: int) -> np.ndarray: - """Below what ``std`` a column counts as having no variance at all. - - ``variance`` above is the one-pass ``E[x^2] - E[x]^2``. For a column whose - values are all (nearly) equal, that subtracts two numbers of size - ``mean ** 2`` which cancel almost completely, and what survives is rounding - noise rather than signal. Summing ``n_rows`` terms accumulates a relative - error of order ``n_rows * eps``, so the noise left in ``variance`` is of - order ``n_rows * eps * mean ** 2`` and the noise in ``std`` of order - ``sqrt(n_rows * eps) * |mean|``. - - A plain ``std > 0`` test cannot see that, and dividing by such a ``std`` - amplifies noise into the output instead of scaling anything: a constant - column comes back as an arbitrary O(1) value rather than the 0 that - centering it should give. So judge "no variance" against that noise floor, - which scales both with the column's magnitude and with how many terms were - summed to reach it. - - Columns that are genuinely all zero have ``mean == 0``, hence a threshold - of 0, and still take the ``std > 0`` branch. - - The floor is set by the one-pass formula, not by the data: a stable - two-pass or Welford variance would leave noise of order ``eps`` rather than - ``n_rows * eps`` and would admit a tighter threshold, at the cost of a - second pass over ``nnz``. - """ - return np.sqrt(n_rows * np.finfo(np.float64).eps) * np.abs(mean) +def _is_constant_column(variance: np.ndarray, mean: np.ndarray, n_rows: int) -> np.ndarray: + """Whether each column's variance is indistinguishable from zero.""" + # A constant column has no unit-variance scaling to be put on, and + # dividing by whatever the arithmetic left behind amplifies rounding noise + # into an arbitrary O(1) value. Bound from scikit-learn's + # `_is_constant_feature` (Chan, Golub & LeVeque), sized to the error + # `_finish_variance` can leave. + eps = np.finfo(np.float64).eps + return variance <= n_rows * eps * variance + (n_rows * mean * eps) ** 2 class NormalizedViewBase: @@ -687,7 +729,7 @@ def __init__( if self._format == "csc": # One fused pass per column for both -- see _column_stats_major_is_col. if need_b or need_gstats: - gene_scale, col_sum, col_sumsq = _column_stats_major_is_col( + gene_scale, col_sum, col_m2, col_corr, col_nnz = _column_stats_major_is_col( arr.major_ptr, arr.values, arr.value_ptr, @@ -696,10 +738,11 @@ def __init__( need_b, need_gstats, self.recipe.g_code, + n_rows, ) else: gene_scale = np.ones(n_cols, dtype=np.float64) - col_sum = col_sumsq = np.zeros(n_cols, dtype=np.float64) + col_sum = col_m2 = col_corr = col_nnz = np.zeros(n_cols, dtype=np.float64) else: # VCSR can't fuse these: gene_scale[c] isn't final until every row # has been scattered into it, so the g-transform pass has to wait @@ -713,7 +756,7 @@ def __init__( gene_scale = np.ones(n_cols, dtype=np.float64) if need_gstats: nthreads = numba.get_num_threads() - col_sum, col_sumsq = _gstats_col_sums_vcs( + kernel_args = ( arr.major_ptr, arr.values, arr.value_ptr, @@ -721,24 +764,26 @@ def __init__( row_scale, gene_scale, self.recipe.g_code, + ) + col_sum, col_nnz = _gstats_col_sums_vcs(*kernel_args, n_cols, nthreads) + col_m2, col_corr = _gstats_col_deviations_vcs( + *kernel_args, + col_sum / n_rows if n_rows > 0 else np.zeros(n_cols, dtype=np.float64), n_cols, nthreads, ) else: - col_sum = col_sumsq = np.zeros(n_cols, dtype=np.float64) + col_sum = col_m2 = col_corr = col_nnz = np.zeros(n_cols, dtype=np.float64) self.gene_scale = gene_scale if need_gstats: - mean = col_sum / n_rows if n_rows > 0 else np.zeros(n_cols, dtype=np.float64) - variance = np.clip( - col_sumsq / n_rows - mean**2 if n_rows > 0 else np.zeros(n_cols), 0.0, None - ) + mean, variance = _finish_variance(col_sum, col_m2, col_corr, col_nnz, n_rows) std = np.sqrt(variance) col_mean = mean if self.recipe.center else np.zeros(n_cols, dtype=np.float64) if self.recipe.post_scale: with np.errstate(divide="ignore", invalid="ignore"): col_post_scale = np.where( - std > _zero_variance_tol(mean, n_rows), 1.0 / std, 1.0 + _is_constant_column(variance, mean, n_rows), 1.0, 1.0 / std ) else: col_post_scale = np.ones(n_cols, dtype=np.float64) diff --git a/tests/test_property_normalization.py b/tests/test_property_normalization.py index f6d44dd..2cd78ab 100644 --- a/tests/test_property_normalization.py +++ b/tests/test_property_normalization.py @@ -83,12 +83,13 @@ def _recipe_reference(dense: np.ndarray, recipe: str) -> np.ndarray: if recipe in ("scanpy", "pearson"): std = g.std(axis=0) with np.errstate(divide="ignore", invalid="ignore"): - # A column with no real variance leaves only rounding noise in - # `std`; dividing by it amplifies noise instead of scaling. Judge - # "no variance" relative to the column's own magnitude, as - # `_zero_variance_tol` does in the library. - tol = np.sqrt(dense.shape[0] * np.finfo(np.float64).eps) * np.abs(g.mean(axis=0)) - s = np.where(std > tol, 1.0 / std, 1.0) + # Constant columns have no unit-variance scaling; same bound as + # the library's `_is_constant_column`, over numpy's own std. + eps = np.finfo(np.float64).eps + n = dense.shape[0] + var = std**2 + constant = var <= n * eps * var + (n * g.mean(axis=0) * eps) ** 2 + s = np.where(constant, 1.0, 1.0 / std) else: s = np.ones(dense.shape[1]) diff --git a/tests/test_vcs_norm_recipes.py b/tests/test_vcs_norm_recipes.py index b366faa..51c4ea1 100644 --- a/tests/test_vcs_norm_recipes.py +++ b/tests/test_vcs_norm_recipes.py @@ -15,7 +15,7 @@ import scipy.sparse as sp from vsparse import RECIPES, Recipe, VCSCAnnData, VCSCArray, VCSRArray -from vsparse._norm_common import NORM_CACHE_MAXSIZE +from vsparse._norm_common import NORM_CACHE_MAXSIZE, _is_constant_column @pytest.fixture(params=[VCSCArray, VCSRArray]) @@ -66,12 +66,13 @@ def _reference(dense: np.ndarray, recipe: str) -> np.ndarray: if recipe in ("scanpy", "pearson"): std = g.std(axis=0) with np.errstate(divide="ignore", invalid="ignore"): - # A column with no real variance leaves only rounding noise in - # `std`; dividing by it amplifies noise instead of scaling. Judge - # "no variance" relative to the column's own magnitude, as - # `_zero_variance_tol` does in the library. - tol = np.sqrt(dense.shape[0] * np.finfo(np.float64).eps) * np.abs(g.mean(axis=0)) - s = np.where(std > tol, 1.0 / std, 1.0) + # Constant columns have no unit-variance scaling; same bound as + # the library's `_is_constant_column`, over numpy's own std. + eps = np.finfo(np.float64).eps + n = dense.shape[0] + var = std**2 + constant = var <= n * eps * var + (n * g.mean(axis=0) * eps) ** 2 + s = np.where(constant, 1.0, 1.0 / std) else: s = np.ones(dense.shape[1]) @@ -367,16 +368,7 @@ def test_anndata_cache_does_not_pin_a_dropped_view(): @pytest.mark.parametrize("value", [1.0, 7.0, 9999.0]) @pytest.mark.parametrize("recipe", ["scanpy", "pearson"]) def test_a_column_with_no_variance_centers_to_zero(vcls, n_rows, value, recipe): - """A constant column has no variance to scale by, so centering must leave 0. - - ``variance`` is computed one-pass as ``E[x^2] - E[x]^2``, which for such a - column cancels to rounding noise rather than to exactly 0. Testing - ``std > 0`` would therefore take the dividing branch and amplify that noise - into an arbitrary O(1) value -- 1.05e-08 for the 6x1 case below, against - the 0 the same input gives at 5 rows. The scale of the noise moves with - the column's magnitude and with how many terms were summed, so the - threshold has to move with both. - """ + """Depth normalization flattens a one-column matrix, so centering leaves 0.""" dense = np.full((n_rows, 1), value) v = vcls.from_scipy(_scipy_for(vcls, dense)) out = v.normalized(recipe).toarray() @@ -392,3 +384,30 @@ def test_a_constant_column_does_not_suppress_its_neighbours(vcls): out = v.normalized("pearson").toarray() varying = np.delete(out, 2, axis=1) assert np.abs(varying).max() > 0.5 + + +@pytest.mark.parametrize("n_rows", [6, 10_000, 1_300_000]) +def test_small_but_real_variance_is_not_called_constant(n_rows): + """A coefficient of variation of 1e-06 is signal, not noise, at any scale here.""" + mean = np.array([9.21]) + std = 1e-6 * mean[0] + assert not _is_constant_column(np.array([std**2]), mean, n_rows)[0] + + +@pytest.mark.parametrize("n_rows", [6, 10_000, 1_300_000]) +def test_constant_detection_still_catches_a_flat_column(n_rows): + """Variance down at the arithmetic's own noise floor is treated as zero.""" + mean = np.array([9.21]) + noise = n_rows * np.finfo(np.float64).eps * mean[0] + assert _is_constant_column(np.array([(noise * 0.1) ** 2]), mean, n_rows)[0] + + +def test_both_layouts_compute_the_same_variance(dense): + """Both layouts must produce the same statistics.""" + if dense.sum() == 0: + pytest.skip("all-zero matrix: median row total is 0") + for recipe in ("scanpy", "pearson", "parafac2"): + r = VCSRArray.from_scipy(sp.csr_array(dense)).normalized(recipe) + c = VCSCArray.from_scipy(sp.csc_array(dense)).normalized(recipe) + np.testing.assert_allclose(r.col_post_scale, c.col_post_scale, rtol=1e-12) + np.testing.assert_allclose(r.col_mean, c.col_mean, rtol=1e-12, atol=1e-15)