fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG - #4138

Open
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch
Open

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG#4138
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch

Conversation

@Arshammik

Copy link
Copy Markdown

Hi,

In _highly_variable_pearson_residuals, the clip parameter is reassigned inside the per-batch loop when it comes in as None: if clip is None: clip = np.sqrt(n). After the first iteration, clip is no longer None, so every subsequent batch silently reuses the first batch's threshold instead of computing its own sqrt(n_batch).

This matters because Lause, Berens & Kobak (2021) define the clip as sqrt(N) over the cells entering the residual computation, and with batch_key each batch is its own residual computation. For a 5000/200 split, the small batch gets sqrt(5000) ≈ 70.7 instead of sqrt(200) ≈ 14.1, so residuals that should be clipped in the small batch are not, its per-gene variances are inflated, and HVG ranking shifts. There is also a quieter reproducibility consequence: because np.unique sorts the labels, the threshold that leaks into all batches is whichever one is alphabetically first, so renaming "A"↔"B" on the same data changes the HVG output.

The fix is a one-line scoping change — a loop-local clip_batch that is recomputed each iteration from the current batch's n. The user-facing clip argument is untouched, so any caller that passed an explicit value still gets exactly that value applied to every batch.

Two tests in tests/test_highly_variable_genes.py cover this. The first monkeypatches _calculate_res_sparse to capture the clip value passed on each batch and asserts two distinct values appear for a 5000/200 split; the second swaps batch labels A↔B on identical data and asserts residual_variances is unchanged. Both fail on main and pass on this branch, and the rest of the pearson-residual suite is green.

rapids-singlecell ports this routine verbatim and inherits the same bug; I have a parallel fix queued there that will reference this PR number once assigned.

Thanks,
Arsham

@codecov

codecovBot commented May 23, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.98%. Comparing base (a656a33) to head (8a6bf5f).
✅ All tests successful. No failed tests found.

Additional details and impacted files
@@ Coverage Diff @@## main #4138 +/- ##
=======================================
Coverage 81.98% 81.98% =======================================
Files 134 134 Lines 13235 13236 +1 =======================================
+ Hits 10851 10852 +1 
Misses 2384 2384 
FlagCoverage Δ
hatch-test.low-vers79.07% <100.00%> (+<0.01%)⬆️
hatch-test.pre81.86% <100.00%> (+<0.01%)⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing linesCoverage Δ
...c/scanpy/experimental/pp/_highly_variable_genes.py94.28% <100.00%> (+0.05%)⬆️

@Zethson

Copy link
Copy Markdown
Member

Thanks! Could you please ensure that the pre-commit checks pass?

@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from 47b118c to c588f69CompareMay 25, 2026 14:18
@ArshammikArshammik changed the title Fix per-batch clip threshold in pearson_residuals HVGfix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVGMay 25, 2026
…esiduals
`_highly_variable_pearson_residuals` reassigned the function-scoped `clip`
parameter from `None` to `sqrt(n_batch_0)` on the first batch iteration,
which made the `if clip is None:` check False on every subsequent batch.
The first batch's clip was therefore reused for every other batch, even
though Lause-Berens-Kobak (2021) specifies a per-batch threshold of
`sqrt(n_cells_in_batch)`.
For an unbalanced split — e.g. two batches of 5000 and 200 cells — the
small batch silently inherited `sqrt(5000) ≈ 70.7` instead of using its
own `sqrt(200) ≈ 14.1`. Residuals in the small batch that should be
clipped were not, inflating per-gene residual variance and shifting HVG
ranking. The bug also made the result depend on alphabetical batch-label
ordering (`np.unique` picks the first label), which is a reproducibility
hazard.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter is
untouched.
Two new tests in `tests/test_highly_variable_genes.py` cover the
regression: one instruments `_calculate_res_sparse` to capture the
per-batch clip value, the other verifies that swapping batch labels A<->B
on identical data produces identical `residual_variances`.
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from c588f69 to 402b167CompareMay 25, 2026 14:43
@Arshammik

Copy link
Copy Markdown
Author

@Zethson confirming that the pre-commit checks have passed.

@Zethson

Copy link
Copy Markdown
Member

@jlause do you have any thoughts on this matter?

@ZethsonZethson added this to the 1.12.2 milestone May 26, 2026
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch 2 times, most recently from 2754759 to 402b167CompareMay 26, 2026 17:58
Intron7 added a commit to scverse/rapids-singlecell that referenced this pull request May 27, 2026
…arson_residuals` (#674)
* hvg(pearson_residuals): scope per-batch clip threshold
`_highly_variable_pearson_residuals` reassigned the function-scoped
`clip` parameter inside the per-batch loop via `if clip is None: clip =
cp.sqrt(n, ...)`. After the first batch, `clip is None` was False, so
every subsequent batch reused the first batch's threshold instead of
computing its own `sqrt(n_cells_in_batch)` per Lause-Berens-Kobak 2021.
For a 5000/200 split with `clip=None`, the small batch silently inherited
`sqrt(5000) ≈ 70.7` instead of using its own `sqrt(200) ≈ 14.1`. Residuals
in the small batch that should be clipped were not, inflating their
per-gene variance and shifting HVG ranking. The result also depended on
which batch label was alphabetically first (since `np.unique` sorts), so
renaming "A"↔"B" silently changed the HVG output.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter
is untouched.
Two new tests in `tests/test_hvg.py`:
- `test_pearson_residuals_batch_clip_scoped_per_batch` monkeypatches the
underlying `_pr_cuda.csc_hvg_res` / `dense_hvg_res` bindings and
asserts that two distinct clip values reach the kernel for a 5000/200
split.
- `test_pearson_residuals_batch_order_invariant` swaps batch labels A↔B
on identical data and asserts that `residual_variances` is unchanged.
Both fail on the unpatched code and pass on the fix; the full 58-test
pearson HVG subset of `tests/test_hvg.py` still passes.
scanpy ports the same routine and has the same bug — see scverse/scanpy#4138.
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* tests/hvg: slim per review
Drop the whitebox test_pearson_residuals_batch_clip_scoped_per_batch (the
order-invariant test already covers the observable behavior), the
file:line reference docstring on the surviving test, and the inline
comment above the clip_batch assignment.
* docs: add release note for #674
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Severin Dicks <37635888+Intron7@users.noreply.github.com>
@ilan-goldilan-gold modified the milestones: 1.12.2, 1.12.3Jun 29, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.3, 1.12.4Jul 24, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.4, 1.12.5Aug 27, 2026
@flying-sheep

Copy link
Copy Markdown
Member

@giovp you initialized the “experimental” idea. What’s the stabilization plan? What’s the criteria to move these out?

Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants

@Arshammik@Zethson@flying-sheep@ilan-gold
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Add copy buttons to all
 blocks\n(function() {\n function addCopyButtons() {\n document.querySelectorAll('pre code').forEach(function(codeBlock) {\n if (codeBlock.parentElement.hasAttribute('data-copy-added')) return;\n codeBlock.parentElement.setAttribute('data-copy-added', 'true');\n \n var btn = document.createElement('button');\n btn.textContent = 'Copy';\n btn.style.cssText = 'position:absolute;top:4px;right:4px;padding:2px 8px;font-size:11px;background:#4ecdc4;border:none;border-radius:4px;color:#1a1a2e;cursor:pointer;opacity:0.7;transition:opacity 0.2s;';\n btn.onmouseover = function() { this.style.opacity = '1'; };\n btn.onmouseout = function() { this.style.opacity = '0.7'; };\n btn.onclick = function() {\n navigator.clipboard.writeText(codeBlock.textContent).then(function() {\n btn.textContent = 'Copied!';\n setTimeout(function() { btn.textContent = 'Copy'; }, 1500);\n });\n };\n codeBlock.parentElement.style.position = 'relative';\n codeBlock.parentElement.appendChild(btn);\n });\n }\n \n addCopyButtons();\n \n // Re-run on dynamic content\n var observer = new MutationObserver(addCopyButtons);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Add Copy Buttons to Code Blocks");
}
} catch(__e) { console.warn('[Userscript:Add Copy Buttons to Code Blocks]', __e); }
})();
(function(){
try {
var __m = "github.com";
var __re = new RegExp('^' + "github\\.com" + '
Skip to content

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG - #4138

Open
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch
Open

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG#4138
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch

Conversation

@Arshammik

Copy link
Copy Markdown

Hi,

In _highly_variable_pearson_residuals, the clip parameter is reassigned inside the per-batch loop when it comes in as None: if clip is None: clip = np.sqrt(n). After the first iteration, clip is no longer None, so every subsequent batch silently reuses the first batch's threshold instead of computing its own sqrt(n_batch).

This matters because Lause, Berens & Kobak (2021) define the clip as sqrt(N) over the cells entering the residual computation, and with batch_key each batch is its own residual computation. For a 5000/200 split, the small batch gets sqrt(5000) ≈ 70.7 instead of sqrt(200) ≈ 14.1, so residuals that should be clipped in the small batch are not, its per-gene variances are inflated, and HVG ranking shifts. There is also a quieter reproducibility consequence: because np.unique sorts the labels, the threshold that leaks into all batches is whichever one is alphabetically first, so renaming "A"↔"B" on the same data changes the HVG output.

The fix is a one-line scoping change — a loop-local clip_batch that is recomputed each iteration from the current batch's n. The user-facing clip argument is untouched, so any caller that passed an explicit value still gets exactly that value applied to every batch.

Two tests in tests/test_highly_variable_genes.py cover this. The first monkeypatches _calculate_res_sparse to capture the clip value passed on each batch and asserts two distinct values appear for a 5000/200 split; the second swaps batch labels A↔B on identical data and asserts residual_variances is unchanged. Both fail on main and pass on this branch, and the rest of the pearson-residual suite is green.

rapids-singlecell ports this routine verbatim and inherits the same bug; I have a parallel fix queued there that will reference this PR number once assigned.

Thanks,
Arsham

@codecov

codecovBot commented May 23, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.98%. Comparing base (a656a33) to head (8a6bf5f).
✅ All tests successful. No failed tests found.

Additional details and impacted files
@@ Coverage Diff @@## main #4138 +/- ##
=======================================
Coverage 81.98% 81.98% =======================================
Files 134 134 Lines 13235 13236 +1 =======================================
+ Hits 10851 10852 +1 
Misses 2384 2384 
FlagCoverage Δ
hatch-test.low-vers79.07% <100.00%> (+<0.01%)⬆️
hatch-test.pre81.86% <100.00%> (+<0.01%)⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing linesCoverage Δ
...c/scanpy/experimental/pp/_highly_variable_genes.py94.28% <100.00%> (+0.05%)⬆️

@Zethson

Copy link
Copy Markdown
Member

Thanks! Could you please ensure that the pre-commit checks pass?

@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from 47b118c to c588f69CompareMay 25, 2026 14:18
@ArshammikArshammik changed the title Fix per-batch clip threshold in pearson_residuals HVGfix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVGMay 25, 2026
…esiduals
`_highly_variable_pearson_residuals` reassigned the function-scoped `clip`
parameter from `None` to `sqrt(n_batch_0)` on the first batch iteration,
which made the `if clip is None:` check False on every subsequent batch.
The first batch's clip was therefore reused for every other batch, even
though Lause-Berens-Kobak (2021) specifies a per-batch threshold of
`sqrt(n_cells_in_batch)`.
For an unbalanced split — e.g. two batches of 5000 and 200 cells — the
small batch silently inherited `sqrt(5000) ≈ 70.7` instead of using its
own `sqrt(200) ≈ 14.1`. Residuals in the small batch that should be
clipped were not, inflating per-gene residual variance and shifting HVG
ranking. The bug also made the result depend on alphabetical batch-label
ordering (`np.unique` picks the first label), which is a reproducibility
hazard.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter is
untouched.
Two new tests in `tests/test_highly_variable_genes.py` cover the
regression: one instruments `_calculate_res_sparse` to capture the
per-batch clip value, the other verifies that swapping batch labels A<->B
on identical data produces identical `residual_variances`.
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from c588f69 to 402b167CompareMay 25, 2026 14:43
@Arshammik

Copy link
Copy Markdown
Author

@Zethson confirming that the pre-commit checks have passed.

@Zethson

Copy link
Copy Markdown
Member

@jlause do you have any thoughts on this matter?

@ZethsonZethson added this to the 1.12.2 milestone May 26, 2026
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch 2 times, most recently from 2754759 to 402b167CompareMay 26, 2026 17:58
Intron7 added a commit to scverse/rapids-singlecell that referenced this pull request May 27, 2026
…arson_residuals` (#674)
* hvg(pearson_residuals): scope per-batch clip threshold
`_highly_variable_pearson_residuals` reassigned the function-scoped
`clip` parameter inside the per-batch loop via `if clip is None: clip =
cp.sqrt(n, ...)`. After the first batch, `clip is None` was False, so
every subsequent batch reused the first batch's threshold instead of
computing its own `sqrt(n_cells_in_batch)` per Lause-Berens-Kobak 2021.
For a 5000/200 split with `clip=None`, the small batch silently inherited
`sqrt(5000) ≈ 70.7` instead of using its own `sqrt(200) ≈ 14.1`. Residuals
in the small batch that should be clipped were not, inflating their
per-gene variance and shifting HVG ranking. The result also depended on
which batch label was alphabetically first (since `np.unique` sorts), so
renaming "A"↔"B" silently changed the HVG output.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter
is untouched.
Two new tests in `tests/test_hvg.py`:
- `test_pearson_residuals_batch_clip_scoped_per_batch` monkeypatches the
underlying `_pr_cuda.csc_hvg_res` / `dense_hvg_res` bindings and
asserts that two distinct clip values reach the kernel for a 5000/200
split.
- `test_pearson_residuals_batch_order_invariant` swaps batch labels A↔B
on identical data and asserts that `residual_variances` is unchanged.
Both fail on the unpatched code and pass on the fix; the full 58-test
pearson HVG subset of `tests/test_hvg.py` still passes.
scanpy ports the same routine and has the same bug — see scverse/scanpy#4138.
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* tests/hvg: slim per review
Drop the whitebox test_pearson_residuals_batch_clip_scoped_per_batch (the
order-invariant test already covers the observable behavior), the
file:line reference docstring on the surviving test, and the inline
comment above the clip_batch assignment.
* docs: add release note for #674
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Severin Dicks <37635888+Intron7@users.noreply.github.com>
@ilan-goldilan-gold modified the milestones: 1.12.2, 1.12.3Jun 29, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.3, 1.12.4Jul 24, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.4, 1.12.5Aug 27, 2026
@flying-sheep

Copy link
Copy Markdown
Member

@giovp you initialized the “experimental” idea. What’s the stabilization plan? What’s the criteria to move these out?

Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants

@Arshammik@Zethson@flying-sheep@ilan-gold
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Force GitHub README to respect dark mode\n(function() {\n var style = document.createElement('style');\n style.textContent = '\n .markdown-body {\n color-scheme: dark light;\n }\n .markdown-body pre { background: #161b22 !important; }\n .markdown-body code { background: rgba(110, 118, 129, 0.4) !important; }\n .markdown-body table th, .markdown-body table td { border-color: #30363d !important; }\n .markdown-body img { background: #0d1117; }\n .markdown-body blockquote { border-left-color: #8b949e; }\n .markdown-body hr { border-color: #30363d; }\n ';\n document.head.appendChild(style);\n})();", "GitHub Dark Mode README Fix"); } } catch(__e) { console.warn('[Userscript:GitHub Dark Mode README Fix]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + '
Skip to content

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG - #4138

Open
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch
Open

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG#4138
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch

Conversation

@Arshammik

Copy link
Copy Markdown

Hi,

In _highly_variable_pearson_residuals, the clip parameter is reassigned inside the per-batch loop when it comes in as None: if clip is None: clip = np.sqrt(n). After the first iteration, clip is no longer None, so every subsequent batch silently reuses the first batch's threshold instead of computing its own sqrt(n_batch).

This matters because Lause, Berens & Kobak (2021) define the clip as sqrt(N) over the cells entering the residual computation, and with batch_key each batch is its own residual computation. For a 5000/200 split, the small batch gets sqrt(5000) ≈ 70.7 instead of sqrt(200) ≈ 14.1, so residuals that should be clipped in the small batch are not, its per-gene variances are inflated, and HVG ranking shifts. There is also a quieter reproducibility consequence: because np.unique sorts the labels, the threshold that leaks into all batches is whichever one is alphabetically first, so renaming "A"↔"B" on the same data changes the HVG output.

The fix is a one-line scoping change — a loop-local clip_batch that is recomputed each iteration from the current batch's n. The user-facing clip argument is untouched, so any caller that passed an explicit value still gets exactly that value applied to every batch.

Two tests in tests/test_highly_variable_genes.py cover this. The first monkeypatches _calculate_res_sparse to capture the clip value passed on each batch and asserts two distinct values appear for a 5000/200 split; the second swaps batch labels A↔B on identical data and asserts residual_variances is unchanged. Both fail on main and pass on this branch, and the rest of the pearson-residual suite is green.

rapids-singlecell ports this routine verbatim and inherits the same bug; I have a parallel fix queued there that will reference this PR number once assigned.

Thanks,
Arsham

@codecov

codecovBot commented May 23, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.98%. Comparing base (a656a33) to head (8a6bf5f).
✅ All tests successful. No failed tests found.

Additional details and impacted files
@@ Coverage Diff @@## main #4138 +/- ##
=======================================
Coverage 81.98% 81.98% =======================================
Files 134 134 Lines 13235 13236 +1 =======================================
+ Hits 10851 10852 +1 
Misses 2384 2384 
FlagCoverage Δ
hatch-test.low-vers79.07% <100.00%> (+<0.01%)⬆️
hatch-test.pre81.86% <100.00%> (+<0.01%)⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing linesCoverage Δ
...c/scanpy/experimental/pp/_highly_variable_genes.py94.28% <100.00%> (+0.05%)⬆️

@Zethson

Copy link
Copy Markdown
Member

Thanks! Could you please ensure that the pre-commit checks pass?

@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from 47b118c to c588f69CompareMay 25, 2026 14:18
@ArshammikArshammik changed the title Fix per-batch clip threshold in pearson_residuals HVGfix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVGMay 25, 2026
…esiduals
`_highly_variable_pearson_residuals` reassigned the function-scoped `clip`
parameter from `None` to `sqrt(n_batch_0)` on the first batch iteration,
which made the `if clip is None:` check False on every subsequent batch.
The first batch's clip was therefore reused for every other batch, even
though Lause-Berens-Kobak (2021) specifies a per-batch threshold of
`sqrt(n_cells_in_batch)`.
For an unbalanced split — e.g. two batches of 5000 and 200 cells — the
small batch silently inherited `sqrt(5000) ≈ 70.7` instead of using its
own `sqrt(200) ≈ 14.1`. Residuals in the small batch that should be
clipped were not, inflating per-gene residual variance and shifting HVG
ranking. The bug also made the result depend on alphabetical batch-label
ordering (`np.unique` picks the first label), which is a reproducibility
hazard.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter is
untouched.
Two new tests in `tests/test_highly_variable_genes.py` cover the
regression: one instruments `_calculate_res_sparse` to capture the
per-batch clip value, the other verifies that swapping batch labels A<->B
on identical data produces identical `residual_variances`.
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from c588f69 to 402b167CompareMay 25, 2026 14:43
@Arshammik

Copy link
Copy Markdown
Author

@Zethson confirming that the pre-commit checks have passed.

@Zethson

Copy link
Copy Markdown
Member

@jlause do you have any thoughts on this matter?

@ZethsonZethson added this to the 1.12.2 milestone May 26, 2026
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch 2 times, most recently from 2754759 to 402b167CompareMay 26, 2026 17:58
Intron7 added a commit to scverse/rapids-singlecell that referenced this pull request May 27, 2026
…arson_residuals` (#674)
* hvg(pearson_residuals): scope per-batch clip threshold
`_highly_variable_pearson_residuals` reassigned the function-scoped
`clip` parameter inside the per-batch loop via `if clip is None: clip =
cp.sqrt(n, ...)`. After the first batch, `clip is None` was False, so
every subsequent batch reused the first batch's threshold instead of
computing its own `sqrt(n_cells_in_batch)` per Lause-Berens-Kobak 2021.
For a 5000/200 split with `clip=None`, the small batch silently inherited
`sqrt(5000) ≈ 70.7` instead of using its own `sqrt(200) ≈ 14.1`. Residuals
in the small batch that should be clipped were not, inflating their
per-gene variance and shifting HVG ranking. The result also depended on
which batch label was alphabetically first (since `np.unique` sorts), so
renaming "A"↔"B" silently changed the HVG output.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter
is untouched.
Two new tests in `tests/test_hvg.py`:
- `test_pearson_residuals_batch_clip_scoped_per_batch` monkeypatches the
underlying `_pr_cuda.csc_hvg_res` / `dense_hvg_res` bindings and
asserts that two distinct clip values reach the kernel for a 5000/200
split.
- `test_pearson_residuals_batch_order_invariant` swaps batch labels A↔B
on identical data and asserts that `residual_variances` is unchanged.
Both fail on the unpatched code and pass on the fix; the full 58-test
pearson HVG subset of `tests/test_hvg.py` still passes.
scanpy ports the same routine and has the same bug — see scverse/scanpy#4138.
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* tests/hvg: slim per review
Drop the whitebox test_pearson_residuals_batch_clip_scoped_per_batch (the
order-invariant test already covers the observable behavior), the
file:line reference docstring on the surviving test, and the inline
comment above the clip_batch assignment.
* docs: add release note for #674
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Severin Dicks <37635888+Intron7@users.noreply.github.com>
@ilan-goldilan-gold modified the milestones: 1.12.2, 1.12.3Jun 29, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.3, 1.12.4Jul 24, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.4, 1.12.5Aug 27, 2026
@flying-sheep

Copy link
Copy Markdown
Member

@giovp you initialized the “experimental” idea. What’s the stabilization plan? What’s the criteria to move these out?

Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants

@Arshammik@Zethson@flying-sheep@ilan-gold
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Highlight search terms from Google/DuckDuckGo/Bing referrer\n(function() {\n var ref = document.referrer;\n var terms = [];\n \n if (ref.includes('google.com') || ref.includes('duckduckgo.com') || ref.includes('bing.com')) {\n var url = new URL(ref);\n var q = url.searchParams.get('q') || url.searchParams.get('p');\n if (q) {\n terms = q.split(/\\s+/).filter(function(t) { return t.length > 2; });\n }\n }\n \n if (terms.length === 0) return;\n \n var style = document.createElement('style');\n style.textContent = '.userscript-highlight { background: #fbbf24; color: #1a1a2e; padding: 1px 3px; border-radius: 2px; }';\n document.head.appendChild(style);\n \n function highlight(node) {\n if (node.nodeType === 3) { // text node\n var text = node.textContent;\n var found = false;\n terms.forEach(function(term) {\n var regex = new RegExp('(' + term.replace(/[.*+?^${}()|[\\]\\\\]/g, '\\\\') + ')', 'gi');\n if (regex.test(text)) {\n found = true;\n var frag = document.createDocumentFragment();\n var parts = text.split(regex);\n parts.forEach(function(part, i) {\n if (i % 2 === 0) {\n frag.appendChild(document.createTextNode(part));\n } else {\n var span = document.createElement('span');\n span.className = 'userscript-highlight';\n span.textContent = part;\n frag.appendChild(span);\n }\n });\n node.parentNode.replaceChild(frag, node);\n }\n });\n } else if (node.nodeType === 1 && node.childNodes) { // element\n var skipTags = ['SCRIPT', 'STYLE', 'NOSCRIPT', 'TEXTAREA', 'INPUT', 'SELECT'];\n if (!skipTags.includes(node.tagName)) {\n Array.from(node.childNodes).forEach(highlight);\n }\n }\n }\n \n highlight(document.body);\n \n // Re-highlight on dynamic content\n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1 || node.nodeType === 3) highlight(node);\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Highlight Search Terms"); } } catch(__e) { console.warn('[Userscript:Highlight Search Terms]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + '
Skip to content

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG - #4138

Open
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch
Open

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG#4138
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch

Conversation

@Arshammik

Copy link
Copy Markdown

Hi,

In _highly_variable_pearson_residuals, the clip parameter is reassigned inside the per-batch loop when it comes in as None: if clip is None: clip = np.sqrt(n). After the first iteration, clip is no longer None, so every subsequent batch silently reuses the first batch's threshold instead of computing its own sqrt(n_batch).

This matters because Lause, Berens & Kobak (2021) define the clip as sqrt(N) over the cells entering the residual computation, and with batch_key each batch is its own residual computation. For a 5000/200 split, the small batch gets sqrt(5000) ≈ 70.7 instead of sqrt(200) ≈ 14.1, so residuals that should be clipped in the small batch are not, its per-gene variances are inflated, and HVG ranking shifts. There is also a quieter reproducibility consequence: because np.unique sorts the labels, the threshold that leaks into all batches is whichever one is alphabetically first, so renaming "A"↔"B" on the same data changes the HVG output.

The fix is a one-line scoping change — a loop-local clip_batch that is recomputed each iteration from the current batch's n. The user-facing clip argument is untouched, so any caller that passed an explicit value still gets exactly that value applied to every batch.

Two tests in tests/test_highly_variable_genes.py cover this. The first monkeypatches _calculate_res_sparse to capture the clip value passed on each batch and asserts two distinct values appear for a 5000/200 split; the second swaps batch labels A↔B on identical data and asserts residual_variances is unchanged. Both fail on main and pass on this branch, and the rest of the pearson-residual suite is green.

rapids-singlecell ports this routine verbatim and inherits the same bug; I have a parallel fix queued there that will reference this PR number once assigned.

Thanks,
Arsham

@codecov

codecovBot commented May 23, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.98%. Comparing base (a656a33) to head (8a6bf5f).
✅ All tests successful. No failed tests found.

Additional details and impacted files
@@ Coverage Diff @@## main #4138 +/- ##
=======================================
Coverage 81.98% 81.98% =======================================
Files 134 134 Lines 13235 13236 +1 =======================================
+ Hits 10851 10852 +1 
Misses 2384 2384 
FlagCoverage Δ
hatch-test.low-vers79.07% <100.00%> (+<0.01%)⬆️
hatch-test.pre81.86% <100.00%> (+<0.01%)⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing linesCoverage Δ
...c/scanpy/experimental/pp/_highly_variable_genes.py94.28% <100.00%> (+0.05%)⬆️

@Zethson

Copy link
Copy Markdown
Member

Thanks! Could you please ensure that the pre-commit checks pass?

@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from 47b118c to c588f69CompareMay 25, 2026 14:18
@ArshammikArshammik changed the title Fix per-batch clip threshold in pearson_residuals HVGfix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVGMay 25, 2026
…esiduals
`_highly_variable_pearson_residuals` reassigned the function-scoped `clip`
parameter from `None` to `sqrt(n_batch_0)` on the first batch iteration,
which made the `if clip is None:` check False on every subsequent batch.
The first batch's clip was therefore reused for every other batch, even
though Lause-Berens-Kobak (2021) specifies a per-batch threshold of
`sqrt(n_cells_in_batch)`.
For an unbalanced split — e.g. two batches of 5000 and 200 cells — the
small batch silently inherited `sqrt(5000) ≈ 70.7` instead of using its
own `sqrt(200) ≈ 14.1`. Residuals in the small batch that should be
clipped were not, inflating per-gene residual variance and shifting HVG
ranking. The bug also made the result depend on alphabetical batch-label
ordering (`np.unique` picks the first label), which is a reproducibility
hazard.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter is
untouched.
Two new tests in `tests/test_highly_variable_genes.py` cover the
regression: one instruments `_calculate_res_sparse` to capture the
per-batch clip value, the other verifies that swapping batch labels A<->B
on identical data produces identical `residual_variances`.
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from c588f69 to 402b167CompareMay 25, 2026 14:43
@Arshammik

Copy link
Copy Markdown
Author

@Zethson confirming that the pre-commit checks have passed.

@Zethson

Copy link
Copy Markdown
Member

@jlause do you have any thoughts on this matter?

@ZethsonZethson added this to the 1.12.2 milestone May 26, 2026
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch 2 times, most recently from 2754759 to 402b167CompareMay 26, 2026 17:58
Intron7 added a commit to scverse/rapids-singlecell that referenced this pull request May 27, 2026
…arson_residuals` (#674)
* hvg(pearson_residuals): scope per-batch clip threshold
`_highly_variable_pearson_residuals` reassigned the function-scoped
`clip` parameter inside the per-batch loop via `if clip is None: clip =
cp.sqrt(n, ...)`. After the first batch, `clip is None` was False, so
every subsequent batch reused the first batch's threshold instead of
computing its own `sqrt(n_cells_in_batch)` per Lause-Berens-Kobak 2021.
For a 5000/200 split with `clip=None`, the small batch silently inherited
`sqrt(5000) ≈ 70.7` instead of using its own `sqrt(200) ≈ 14.1`. Residuals
in the small batch that should be clipped were not, inflating their
per-gene variance and shifting HVG ranking. The result also depended on
which batch label was alphabetically first (since `np.unique` sorts), so
renaming "A"↔"B" silently changed the HVG output.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter
is untouched.
Two new tests in `tests/test_hvg.py`:
- `test_pearson_residuals_batch_clip_scoped_per_batch` monkeypatches the
underlying `_pr_cuda.csc_hvg_res` / `dense_hvg_res` bindings and
asserts that two distinct clip values reach the kernel for a 5000/200
split.
- `test_pearson_residuals_batch_order_invariant` swaps batch labels A↔B
on identical data and asserts that `residual_variances` is unchanged.
Both fail on the unpatched code and pass on the fix; the full 58-test
pearson HVG subset of `tests/test_hvg.py` still passes.
scanpy ports the same routine and has the same bug — see scverse/scanpy#4138.
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* tests/hvg: slim per review
Drop the whitebox test_pearson_residuals_batch_clip_scoped_per_batch (the
order-invariant test already covers the observable behavior), the
file:line reference docstring on the surviving test, and the inline
comment above the clip_batch assignment.
* docs: add release note for #674
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Severin Dicks <37635888+Intron7@users.noreply.github.com>
@ilan-goldilan-gold modified the milestones: 1.12.2, 1.12.3Jun 29, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.3, 1.12.4Jul 24, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.4, 1.12.5Aug 27, 2026
@flying-sheep

Copy link
Copy Markdown
Member

@giovp you initialized the “experimental” idea. What’s the stabilization plan? What’s the criteria to move these out?

Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants

@Arshammik@Zethson@flying-sheep@ilan-gold
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Strip utm_, fbclid, gclid, etc. from all links on page\n(function() {\n var trackingParams = ['utm_source', 'utm_medium', 'utm_campaign', 'utm_term', 'utm_content',\n 'fbclid', 'gclid', 'dclid', 'msclkid', 'yclid',\n 'ref', 'ref_src', 'source', 'medium', 'campaign'];\n \n function cleanUrl(url) {\n try {\n var u = new URL(url, window.location.origin);\n var changed = false;\n trackingParams.forEach(function(p) {\n if (u.searchParams.has(p)) {\n u.searchParams.delete(p);\n changed = true;\n }\n });\n return changed ? u.toString() : url;\n } catch (e) {\n return url;\n }\n }\n \n function cleanLinks() {\n document.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n \n cleanLinks();\n \n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1) {\n if (node.tagName === 'A') cleanLinks();\n node.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Remove Tracking Parameters from Links"); } } catch(__e) { console.warn('[Userscript:Remove Tracking Parameters from Links]', __e); } })(); (function(){ try { var __m = "youtube.com"; var __re = new RegExp('^' + "youtube\\.com" + '
Skip to content

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG - #4138

Open
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch
Open

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG#4138
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch

Conversation

@Arshammik

Copy link
Copy Markdown

Hi,

In _highly_variable_pearson_residuals, the clip parameter is reassigned inside the per-batch loop when it comes in as None: if clip is None: clip = np.sqrt(n). After the first iteration, clip is no longer None, so every subsequent batch silently reuses the first batch's threshold instead of computing its own sqrt(n_batch).

This matters because Lause, Berens & Kobak (2021) define the clip as sqrt(N) over the cells entering the residual computation, and with batch_key each batch is its own residual computation. For a 5000/200 split, the small batch gets sqrt(5000) ≈ 70.7 instead of sqrt(200) ≈ 14.1, so residuals that should be clipped in the small batch are not, its per-gene variances are inflated, and HVG ranking shifts. There is also a quieter reproducibility consequence: because np.unique sorts the labels, the threshold that leaks into all batches is whichever one is alphabetically first, so renaming "A"↔"B" on the same data changes the HVG output.

The fix is a one-line scoping change — a loop-local clip_batch that is recomputed each iteration from the current batch's n. The user-facing clip argument is untouched, so any caller that passed an explicit value still gets exactly that value applied to every batch.

Two tests in tests/test_highly_variable_genes.py cover this. The first monkeypatches _calculate_res_sparse to capture the clip value passed on each batch and asserts two distinct values appear for a 5000/200 split; the second swaps batch labels A↔B on identical data and asserts residual_variances is unchanged. Both fail on main and pass on this branch, and the rest of the pearson-residual suite is green.

rapids-singlecell ports this routine verbatim and inherits the same bug; I have a parallel fix queued there that will reference this PR number once assigned.

Thanks,
Arsham

@codecov

codecovBot commented May 23, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.98%. Comparing base (a656a33) to head (8a6bf5f).
✅ All tests successful. No failed tests found.

Additional details and impacted files
@@ Coverage Diff @@## main #4138 +/- ##
=======================================
Coverage 81.98% 81.98% =======================================
Files 134 134 Lines 13235 13236 +1 =======================================
+ Hits 10851 10852 +1 
Misses 2384 2384 
FlagCoverage Δ
hatch-test.low-vers79.07% <100.00%> (+<0.01%)⬆️
hatch-test.pre81.86% <100.00%> (+<0.01%)⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing linesCoverage Δ
...c/scanpy/experimental/pp/_highly_variable_genes.py94.28% <100.00%> (+0.05%)⬆️

@Zethson

Copy link
Copy Markdown
Member

Thanks! Could you please ensure that the pre-commit checks pass?

@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from 47b118c to c588f69CompareMay 25, 2026 14:18
@ArshammikArshammik changed the title Fix per-batch clip threshold in pearson_residuals HVGfix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVGMay 25, 2026
…esiduals
`_highly_variable_pearson_residuals` reassigned the function-scoped `clip`
parameter from `None` to `sqrt(n_batch_0)` on the first batch iteration,
which made the `if clip is None:` check False on every subsequent batch.
The first batch's clip was therefore reused for every other batch, even
though Lause-Berens-Kobak (2021) specifies a per-batch threshold of
`sqrt(n_cells_in_batch)`.
For an unbalanced split — e.g. two batches of 5000 and 200 cells — the
small batch silently inherited `sqrt(5000) ≈ 70.7` instead of using its
own `sqrt(200) ≈ 14.1`. Residuals in the small batch that should be
clipped were not, inflating per-gene residual variance and shifting HVG
ranking. The bug also made the result depend on alphabetical batch-label
ordering (`np.unique` picks the first label), which is a reproducibility
hazard.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter is
untouched.
Two new tests in `tests/test_highly_variable_genes.py` cover the
regression: one instruments `_calculate_res_sparse` to capture the
per-batch clip value, the other verifies that swapping batch labels A<->B
on identical data produces identical `residual_variances`.
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from c588f69 to 402b167CompareMay 25, 2026 14:43
@Arshammik

Copy link
Copy Markdown
Author

@Zethson confirming that the pre-commit checks have passed.

@Zethson

Copy link
Copy Markdown
Member

@jlause do you have any thoughts on this matter?

@ZethsonZethson added this to the 1.12.2 milestone May 26, 2026
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch 2 times, most recently from 2754759 to 402b167CompareMay 26, 2026 17:58
Intron7 added a commit to scverse/rapids-singlecell that referenced this pull request May 27, 2026
…arson_residuals` (#674)
* hvg(pearson_residuals): scope per-batch clip threshold
`_highly_variable_pearson_residuals` reassigned the function-scoped
`clip` parameter inside the per-batch loop via `if clip is None: clip =
cp.sqrt(n, ...)`. After the first batch, `clip is None` was False, so
every subsequent batch reused the first batch's threshold instead of
computing its own `sqrt(n_cells_in_batch)` per Lause-Berens-Kobak 2021.
For a 5000/200 split with `clip=None`, the small batch silently inherited
`sqrt(5000) ≈ 70.7` instead of using its own `sqrt(200) ≈ 14.1`. Residuals
in the small batch that should be clipped were not, inflating their
per-gene variance and shifting HVG ranking. The result also depended on
which batch label was alphabetically first (since `np.unique` sorts), so
renaming "A"↔"B" silently changed the HVG output.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter
is untouched.
Two new tests in `tests/test_hvg.py`:
- `test_pearson_residuals_batch_clip_scoped_per_batch` monkeypatches the
underlying `_pr_cuda.csc_hvg_res` / `dense_hvg_res` bindings and
asserts that two distinct clip values reach the kernel for a 5000/200
split.
- `test_pearson_residuals_batch_order_invariant` swaps batch labels A↔B
on identical data and asserts that `residual_variances` is unchanged.
Both fail on the unpatched code and pass on the fix; the full 58-test
pearson HVG subset of `tests/test_hvg.py` still passes.
scanpy ports the same routine and has the same bug — see scverse/scanpy#4138.
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* tests/hvg: slim per review
Drop the whitebox test_pearson_residuals_batch_clip_scoped_per_batch (the
order-invariant test already covers the observable behavior), the
file:line reference docstring on the surviving test, and the inline
comment above the clip_batch assignment.
* docs: add release note for #674
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Severin Dicks <37635888+Intron7@users.noreply.github.com>
@ilan-goldilan-gold modified the milestones: 1.12.2, 1.12.3Jun 29, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.3, 1.12.4Jul 24, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.4, 1.12.5Aug 27, 2026
@flying-sheep

Copy link
Copy Markdown
Member

@giovp you initialized the “experimental” idea. What’s the stabilization plan? What’s the criteria to move these out?

Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants

@Arshammik@Zethson@flying-sheep@ilan-gold
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Auto-enable theater mode on YouTube\n(function() {\n function tryTheater() {\n var btn = document.querySelector('button[aria-label=\"Theater mode\"], ytd-player #player button[title=\"Theater mode\"]');\n if (btn && !btn.classList.contains('activated')) {\n btn.click();\n }\n }\n \n // Try immediately\n tryTheater();\n \n // Try after navigation (SPA)\n var lastUrl = location.href;\n setInterval(function() {\n if (location.href !== lastUrl) {\n lastUrl = location.href;\n setTimeout(tryTheater, 500);\n }\n }, 1000);\n \n // Also try on player load\n var observer = new MutationObserver(tryTheater);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "YouTube Theater Mode Default"); } } catch(__e) { console.warn('[Userscript:YouTube Theater Mode Default]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + '
Skip to content

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG - #4138

Open
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch
Open

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG#4138
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch

Conversation

@Arshammik

Copy link
Copy Markdown

Hi,

In _highly_variable_pearson_residuals, the clip parameter is reassigned inside the per-batch loop when it comes in as None: if clip is None: clip = np.sqrt(n). After the first iteration, clip is no longer None, so every subsequent batch silently reuses the first batch's threshold instead of computing its own sqrt(n_batch).

This matters because Lause, Berens & Kobak (2021) define the clip as sqrt(N) over the cells entering the residual computation, and with batch_key each batch is its own residual computation. For a 5000/200 split, the small batch gets sqrt(5000) ≈ 70.7 instead of sqrt(200) ≈ 14.1, so residuals that should be clipped in the small batch are not, its per-gene variances are inflated, and HVG ranking shifts. There is also a quieter reproducibility consequence: because np.unique sorts the labels, the threshold that leaks into all batches is whichever one is alphabetically first, so renaming "A"↔"B" on the same data changes the HVG output.

The fix is a one-line scoping change — a loop-local clip_batch that is recomputed each iteration from the current batch's n. The user-facing clip argument is untouched, so any caller that passed an explicit value still gets exactly that value applied to every batch.

Two tests in tests/test_highly_variable_genes.py cover this. The first monkeypatches _calculate_res_sparse to capture the clip value passed on each batch and asserts two distinct values appear for a 5000/200 split; the second swaps batch labels A↔B on identical data and asserts residual_variances is unchanged. Both fail on main and pass on this branch, and the rest of the pearson-residual suite is green.

rapids-singlecell ports this routine verbatim and inherits the same bug; I have a parallel fix queued there that will reference this PR number once assigned.

Thanks,
Arsham

@codecov

codecovBot commented May 23, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.98%. Comparing base (a656a33) to head (8a6bf5f).
✅ All tests successful. No failed tests found.

Additional details and impacted files
@@ Coverage Diff @@## main #4138 +/- ##
=======================================
Coverage 81.98% 81.98% =======================================
Files 134 134 Lines 13235 13236 +1 =======================================
+ Hits 10851 10852 +1 
Misses 2384 2384 
FlagCoverage Δ
hatch-test.low-vers79.07% <100.00%> (+<0.01%)⬆️
hatch-test.pre81.86% <100.00%> (+<0.01%)⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing linesCoverage Δ
...c/scanpy/experimental/pp/_highly_variable_genes.py94.28% <100.00%> (+0.05%)⬆️

@Zethson

Copy link
Copy Markdown
Member

Thanks! Could you please ensure that the pre-commit checks pass?

@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from 47b118c to c588f69CompareMay 25, 2026 14:18
@ArshammikArshammik changed the title Fix per-batch clip threshold in pearson_residuals HVGfix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVGMay 25, 2026
…esiduals
`_highly_variable_pearson_residuals` reassigned the function-scoped `clip`
parameter from `None` to `sqrt(n_batch_0)` on the first batch iteration,
which made the `if clip is None:` check False on every subsequent batch.
The first batch's clip was therefore reused for every other batch, even
though Lause-Berens-Kobak (2021) specifies a per-batch threshold of
`sqrt(n_cells_in_batch)`.
For an unbalanced split — e.g. two batches of 5000 and 200 cells — the
small batch silently inherited `sqrt(5000) ≈ 70.7` instead of using its
own `sqrt(200) ≈ 14.1`. Residuals in the small batch that should be
clipped were not, inflating per-gene residual variance and shifting HVG
ranking. The bug also made the result depend on alphabetical batch-label
ordering (`np.unique` picks the first label), which is a reproducibility
hazard.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter is
untouched.
Two new tests in `tests/test_highly_variable_genes.py` cover the
regression: one instruments `_calculate_res_sparse` to capture the
per-batch clip value, the other verifies that swapping batch labels A<->B
on identical data produces identical `residual_variances`.
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from c588f69 to 402b167CompareMay 25, 2026 14:43
@Arshammik

Copy link
Copy Markdown
Author

@Zethson confirming that the pre-commit checks have passed.

@Zethson

Copy link
Copy Markdown
Member

@jlause do you have any thoughts on this matter?

@ZethsonZethson added this to the 1.12.2 milestone May 26, 2026
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch 2 times, most recently from 2754759 to 402b167CompareMay 26, 2026 17:58
Intron7 added a commit to scverse/rapids-singlecell that referenced this pull request May 27, 2026
…arson_residuals` (#674)
* hvg(pearson_residuals): scope per-batch clip threshold
`_highly_variable_pearson_residuals` reassigned the function-scoped
`clip` parameter inside the per-batch loop via `if clip is None: clip =
cp.sqrt(n, ...)`. After the first batch, `clip is None` was False, so
every subsequent batch reused the first batch's threshold instead of
computing its own `sqrt(n_cells_in_batch)` per Lause-Berens-Kobak 2021.
For a 5000/200 split with `clip=None`, the small batch silently inherited
`sqrt(5000) ≈ 70.7` instead of using its own `sqrt(200) ≈ 14.1`. Residuals
in the small batch that should be clipped were not, inflating their
per-gene variance and shifting HVG ranking. The result also depended on
which batch label was alphabetically first (since `np.unique` sorts), so
renaming "A"↔"B" silently changed the HVG output.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter
is untouched.
Two new tests in `tests/test_hvg.py`:
- `test_pearson_residuals_batch_clip_scoped_per_batch` monkeypatches the
underlying `_pr_cuda.csc_hvg_res` / `dense_hvg_res` bindings and
asserts that two distinct clip values reach the kernel for a 5000/200
split.
- `test_pearson_residuals_batch_order_invariant` swaps batch labels A↔B
on identical data and asserts that `residual_variances` is unchanged.
Both fail on the unpatched code and pass on the fix; the full 58-test
pearson HVG subset of `tests/test_hvg.py` still passes.
scanpy ports the same routine and has the same bug — see scverse/scanpy#4138.
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* tests/hvg: slim per review
Drop the whitebox test_pearson_residuals_batch_clip_scoped_per_batch (the
order-invariant test already covers the observable behavior), the
file:line reference docstring on the surviving test, and the inline
comment above the clip_batch assignment.
* docs: add release note for #674
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Severin Dicks <37635888+Intron7@users.noreply.github.com>
@ilan-goldilan-gold modified the milestones: 1.12.2, 1.12.3Jun 29, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.3, 1.12.4Jul 24, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.4, 1.12.5Aug 27, 2026
@flying-sheep

Copy link
Copy Markdown
Member

@giovp you initialized the “experimental” idea. What’s the stabilization plan? What’s the criteria to move these out?

Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants

@Arshammik@Zethson@flying-sheep@ilan-gold
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Remove or un-stick sticky/fixed headers that block content\n(function() {\n function unstick() {\n document.querySelectorAll('header, nav, [role=\"banner\"], .header, .navbar, .sticky, .fixed-top, [style*=\"position: fixed\"], [style*=\"position:sticky\"]').forEach(function(el) {\n if (el.style.position === 'fixed' || el.style.position === 'sticky' || \n getComputedStyle(el).position === 'fixed' || getComputedStyle(el).position === 'sticky') {\n el.style.position = 'static';\n el.style.top = 'auto';\n el.style.zIndex = 'auto';\n }\n });\n }\n \n unstick();\n \n var observer = new MutationObserver(unstick);\n observer.observe(document.body, { childList: true, subtree: true, attributes: true, attributeFilter: ['style', 'class'] });\n})();", "Kill Sticky Headers"); } } catch(__e) { console.warn('[Userscript:Kill Sticky Headers]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + '
Skip to content

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG - #4138

Open
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch
Open

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG#4138
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch

Conversation

@Arshammik

Copy link
Copy Markdown

Hi,

In _highly_variable_pearson_residuals, the clip parameter is reassigned inside the per-batch loop when it comes in as None: if clip is None: clip = np.sqrt(n). After the first iteration, clip is no longer None, so every subsequent batch silently reuses the first batch's threshold instead of computing its own sqrt(n_batch).

This matters because Lause, Berens & Kobak (2021) define the clip as sqrt(N) over the cells entering the residual computation, and with batch_key each batch is its own residual computation. For a 5000/200 split, the small batch gets sqrt(5000) ≈ 70.7 instead of sqrt(200) ≈ 14.1, so residuals that should be clipped in the small batch are not, its per-gene variances are inflated, and HVG ranking shifts. There is also a quieter reproducibility consequence: because np.unique sorts the labels, the threshold that leaks into all batches is whichever one is alphabetically first, so renaming "A"↔"B" on the same data changes the HVG output.

The fix is a one-line scoping change — a loop-local clip_batch that is recomputed each iteration from the current batch's n. The user-facing clip argument is untouched, so any caller that passed an explicit value still gets exactly that value applied to every batch.

Two tests in tests/test_highly_variable_genes.py cover this. The first monkeypatches _calculate_res_sparse to capture the clip value passed on each batch and asserts two distinct values appear for a 5000/200 split; the second swaps batch labels A↔B on identical data and asserts residual_variances is unchanged. Both fail on main and pass on this branch, and the rest of the pearson-residual suite is green.

rapids-singlecell ports this routine verbatim and inherits the same bug; I have a parallel fix queued there that will reference this PR number once assigned.

Thanks,
Arsham

@codecov

codecovBot commented May 23, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.98%. Comparing base (a656a33) to head (8a6bf5f).
✅ All tests successful. No failed tests found.

Additional details and impacted files
@@ Coverage Diff @@## main #4138 +/- ##
=======================================
Coverage 81.98% 81.98% =======================================
Files 134 134 Lines 13235 13236 +1 =======================================
+ Hits 10851 10852 +1 
Misses 2384 2384 
FlagCoverage Δ
hatch-test.low-vers79.07% <100.00%> (+<0.01%)⬆️
hatch-test.pre81.86% <100.00%> (+<0.01%)⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing linesCoverage Δ
...c/scanpy/experimental/pp/_highly_variable_genes.py94.28% <100.00%> (+0.05%)⬆️

@Zethson

Copy link
Copy Markdown
Member

Thanks! Could you please ensure that the pre-commit checks pass?

@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from 47b118c to c588f69CompareMay 25, 2026 14:18
@ArshammikArshammik changed the title Fix per-batch clip threshold in pearson_residuals HVGfix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVGMay 25, 2026
…esiduals
`_highly_variable_pearson_residuals` reassigned the function-scoped `clip`
parameter from `None` to `sqrt(n_batch_0)` on the first batch iteration,
which made the `if clip is None:` check False on every subsequent batch.
The first batch's clip was therefore reused for every other batch, even
though Lause-Berens-Kobak (2021) specifies a per-batch threshold of
`sqrt(n_cells_in_batch)`.
For an unbalanced split — e.g. two batches of 5000 and 200 cells — the
small batch silently inherited `sqrt(5000) ≈ 70.7` instead of using its
own `sqrt(200) ≈ 14.1`. Residuals in the small batch that should be
clipped were not, inflating per-gene residual variance and shifting HVG
ranking. The bug also made the result depend on alphabetical batch-label
ordering (`np.unique` picks the first label), which is a reproducibility
hazard.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter is
untouched.
Two new tests in `tests/test_highly_variable_genes.py` cover the
regression: one instruments `_calculate_res_sparse` to capture the
per-batch clip value, the other verifies that swapping batch labels A<->B
on identical data produces identical `residual_variances`.
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from c588f69 to 402b167CompareMay 25, 2026 14:43
@Arshammik

Copy link
Copy Markdown
Author

@Zethson confirming that the pre-commit checks have passed.

@Zethson

Copy link
Copy Markdown
Member

@jlause do you have any thoughts on this matter?

@ZethsonZethson added this to the 1.12.2 milestone May 26, 2026
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch 2 times, most recently from 2754759 to 402b167CompareMay 26, 2026 17:58
Intron7 added a commit to scverse/rapids-singlecell that referenced this pull request May 27, 2026
…arson_residuals` (#674)
* hvg(pearson_residuals): scope per-batch clip threshold
`_highly_variable_pearson_residuals` reassigned the function-scoped
`clip` parameter inside the per-batch loop via `if clip is None: clip =
cp.sqrt(n, ...)`. After the first batch, `clip is None` was False, so
every subsequent batch reused the first batch's threshold instead of
computing its own `sqrt(n_cells_in_batch)` per Lause-Berens-Kobak 2021.
For a 5000/200 split with `clip=None`, the small batch silently inherited
`sqrt(5000) ≈ 70.7` instead of using its own `sqrt(200) ≈ 14.1`. Residuals
in the small batch that should be clipped were not, inflating their
per-gene variance and shifting HVG ranking. The result also depended on
which batch label was alphabetically first (since `np.unique` sorts), so
renaming "A"↔"B" silently changed the HVG output.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter
is untouched.
Two new tests in `tests/test_hvg.py`:
- `test_pearson_residuals_batch_clip_scoped_per_batch` monkeypatches the
underlying `_pr_cuda.csc_hvg_res` / `dense_hvg_res` bindings and
asserts that two distinct clip values reach the kernel for a 5000/200
split.
- `test_pearson_residuals_batch_order_invariant` swaps batch labels A↔B
on identical data and asserts that `residual_variances` is unchanged.
Both fail on the unpatched code and pass on the fix; the full 58-test
pearson HVG subset of `tests/test_hvg.py` still passes.
scanpy ports the same routine and has the same bug — see scverse/scanpy#4138.
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* tests/hvg: slim per review
Drop the whitebox test_pearson_residuals_batch_clip_scoped_per_batch (the
order-invariant test already covers the observable behavior), the
file:line reference docstring on the surviving test, and the inline
comment above the clip_batch assignment.
* docs: add release note for #674
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Severin Dicks <37635888+Intron7@users.noreply.github.com>
@ilan-goldilan-gold modified the milestones: 1.12.2, 1.12.3Jun 29, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.3, 1.12.4Jul 24, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.4, 1.12.5Aug 27, 2026
@flying-sheep

Copy link
Copy Markdown
Member

@giovp you initialized the “experimental” idea. What’s the stabilization plan? What’s the criteria to move these out?

Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants

@Arshammik@Zethson@flying-sheep@ilan-gold
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Universal Dark Mode - works on any site\n(function() {\n var enabled = true;\n \n function applyDarkMode() {\n if (!enabled) return;\n \n // Create style element if it doesn't exist\n var style = document.getElementById('universal-dark-mode-style');\n if (!style) {\n style = document.createElement('style');\n style.id = 'universal-dark-mode-style';\n document.head.appendChild(style);\n }\n \n // Dark mode CSS - inverts colors but preserves images/video\n style.textContent = '\n /* Invert everything except media */\n html {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #1a1a2e !important;\n }\n \n /* Restore images, videos, iframes, canvas */\n img, video, iframe, canvas, svg, picture, [style*=\"background-image\"] {\n filter: invert(1) hue-rotate(180deg) !important;\n }\n \n /* Preserve specific elements that should not be inverted */\n .no-dark-mode, .no-dark-mode *,\n [data-theme=\"light\"], [data-theme=\"light\"],\n .ace_editor, .ace_editor *,\n .CodeMirror, .CodeMirror *,\n .monaco-editor, .monaco-editor *,\n .markdown-body pre, .markdown-body pre *,\n .highlight, .highlight *,\n pre code, pre code * {\n filter: none !important;\n }\n \n /* Fix common UI elements */\n .modal, .popup, .dropdown-menu, .tooltip, .popover {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #2d2d44 !important;\n border-color: #444 !important;\n }\n \n /* Scrollbars */\n ::-webkit-scrollbar { background: #1a1a2e !important; }\n ::-webkit-scrollbar-thumb { background: #444 !important; }\n ::-webkit-scrollbar-thumb:hover { background: #555 !important; }\n \n /* Selection */\n ::selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ::-moz-selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ';\n }\n \n function removeDarkMode() {\n var style = document.getElementById('universal-dark-mode-style');\n if (style) style.remove();\n }\n \n // Toggle with Alt+Shift+D\n document.addEventListener('keydown', function(e) {\n if (e.altKey && e.shiftKey && e.key === 'D') {\n e.preventDefault();\n enabled = !enabled;\n if (enabled) {\n applyDarkMode();\n console.log('[Universal Dark Mode] Enabled');\n } else {\n removeDarkMode();\n console.log('[Universal Dark Mode] Disabled');\n }\n }\n });\n \n // Apply on load\n applyDarkMode();\n \n // Re-apply on dynamic content\n var observer = new MutationObserver(function(mutations) {\n if (enabled && !document.getElementById('universal-dark-mode-style')) {\n applyDarkMode();\n }\n });\n observer.observe(document.head, { childList: true });\n \n console.log('[Universal Dark Mode] Loaded - Press Alt+Shift+D to toggle');\n})();", "Universal Dark Mode"); } } catch(__e) { console.warn('[Userscript:Universal Dark Mode]', __e); } })(); })();
Skip to content

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG - #4138

Open
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch
Open

fix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVG#4138
Arshammik wants to merge 2 commits into
scverse:mainfrom
Arshammik:fix/hvg-pearson-clip-per-batch

Conversation

@Arshammik

Copy link
Copy Markdown

Hi,

In _highly_variable_pearson_residuals, the clip parameter is reassigned inside the per-batch loop when it comes in as None: if clip is None: clip = np.sqrt(n). After the first iteration, clip is no longer None, so every subsequent batch silently reuses the first batch's threshold instead of computing its own sqrt(n_batch).

This matters because Lause, Berens & Kobak (2021) define the clip as sqrt(N) over the cells entering the residual computation, and with batch_key each batch is its own residual computation. For a 5000/200 split, the small batch gets sqrt(5000) ≈ 70.7 instead of sqrt(200) ≈ 14.1, so residuals that should be clipped in the small batch are not, its per-gene variances are inflated, and HVG ranking shifts. There is also a quieter reproducibility consequence: because np.unique sorts the labels, the threshold that leaks into all batches is whichever one is alphabetically first, so renaming "A"↔"B" on the same data changes the HVG output.

The fix is a one-line scoping change — a loop-local clip_batch that is recomputed each iteration from the current batch's n. The user-facing clip argument is untouched, so any caller that passed an explicit value still gets exactly that value applied to every batch.

Two tests in tests/test_highly_variable_genes.py cover this. The first monkeypatches _calculate_res_sparse to capture the clip value passed on each batch and asserts two distinct values appear for a 5000/200 split; the second swaps batch labels A↔B on identical data and asserts residual_variances is unchanged. Both fail on main and pass on this branch, and the rest of the pearson-residual suite is green.

rapids-singlecell ports this routine verbatim and inherits the same bug; I have a parallel fix queued there that will reference this PR number once assigned.

Thanks,
Arsham

@codecov

codecovBot commented May 23, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.98%. Comparing base (a656a33) to head (8a6bf5f).
✅ All tests successful. No failed tests found.

Additional details and impacted files
@@ Coverage Diff @@## main #4138 +/- ##
=======================================
Coverage 81.98% 81.98% =======================================
Files 134 134 Lines 13235 13236 +1 =======================================
+ Hits 10851 10852 +1 
Misses 2384 2384 
FlagCoverage Δ
hatch-test.low-vers79.07% <100.00%> (+<0.01%)⬆️
hatch-test.pre81.86% <100.00%> (+<0.01%)⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing linesCoverage Δ
...c/scanpy/experimental/pp/_highly_variable_genes.py94.28% <100.00%> (+0.05%)⬆️

@Zethson

Copy link
Copy Markdown
Member

Thanks! Could you please ensure that the pre-commit checks pass?

@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from 47b118c to c588f69CompareMay 25, 2026 14:18
@ArshammikArshammik changed the title Fix per-batch clip threshold in pearson_residuals HVGfix(experimental.pp.hvg): scope per-batch clip threshold in pearson_residuals HVGMay 25, 2026
…esiduals
`_highly_variable_pearson_residuals` reassigned the function-scoped `clip`
parameter from `None` to `sqrt(n_batch_0)` on the first batch iteration,
which made the `if clip is None:` check False on every subsequent batch.
The first batch's clip was therefore reused for every other batch, even
though Lause-Berens-Kobak (2021) specifies a per-batch threshold of
`sqrt(n_cells_in_batch)`.
For an unbalanced split — e.g. two batches of 5000 and 200 cells — the
small batch silently inherited `sqrt(5000) ≈ 70.7` instead of using its
own `sqrt(200) ≈ 14.1`. Residuals in the small batch that should be
clipped were not, inflating per-gene residual variance and shifting HVG
ranking. The bug also made the result depend on alphabetical batch-label
ordering (`np.unique` picks the first label), which is a reproducibility
hazard.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter is
untouched.
Two new tests in `tests/test_highly_variable_genes.py` cover the
regression: one instruments `_calculate_res_sparse` to capture the
per-batch clip value, the other verifies that swapping batch labels A<->B
on identical data produces identical `residual_variances`.
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch from c588f69 to 402b167CompareMay 25, 2026 14:43
@Arshammik

Copy link
Copy Markdown
Author

@Zethson confirming that the pre-commit checks have passed.

@Zethson

Copy link
Copy Markdown
Member

@jlause do you have any thoughts on this matter?

@ZethsonZethson added this to the 1.12.2 milestone May 26, 2026
@Arshammik
Arshammikforce-pushed the fix/hvg-pearson-clip-per-batch branch 2 times, most recently from 2754759 to 402b167CompareMay 26, 2026 17:58
Intron7 added a commit to scverse/rapids-singlecell that referenced this pull request May 27, 2026
…arson_residuals` (#674)
* hvg(pearson_residuals): scope per-batch clip threshold
`_highly_variable_pearson_residuals` reassigned the function-scoped
`clip` parameter inside the per-batch loop via `if clip is None: clip =
cp.sqrt(n, ...)`. After the first batch, `clip is None` was False, so
every subsequent batch reused the first batch's threshold instead of
computing its own `sqrt(n_cells_in_batch)` per Lause-Berens-Kobak 2021.
For a 5000/200 split with `clip=None`, the small batch silently inherited
`sqrt(5000) ≈ 70.7` instead of using its own `sqrt(200) ≈ 14.1`. Residuals
in the small batch that should be clipped were not, inflating their
per-gene variance and shifting HVG ranking. The result also depended on
which batch label was alphabetically first (since `np.unique` sorts), so
renaming "A"↔"B" silently changed the HVG output.
The fix introduces a loop-local `clip_batch` variable so the per-batch
default is recomputed each iteration. The user-facing `clip` parameter
is untouched.
Two new tests in `tests/test_hvg.py`:
- `test_pearson_residuals_batch_clip_scoped_per_batch` monkeypatches the
underlying `_pr_cuda.csc_hvg_res` / `dense_hvg_res` bindings and
asserts that two distinct clip values reach the kernel for a 5000/200
split.
- `test_pearson_residuals_batch_order_invariant` swaps batch labels A↔B
on identical data and asserts that `residual_variances` is unchanged.
Both fail on the unpatched code and pass on the fix; the full 58-test
pearson HVG subset of `tests/test_hvg.py` still passes.
scanpy ports the same routine and has the same bug — see scverse/scanpy#4138.
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* tests/hvg: slim per review
Drop the whitebox test_pearson_residuals_batch_clip_scoped_per_batch (the
order-invariant test already covers the observable behavior), the
file:line reference docstring on the surviving test, and the inline
comment above the clip_batch assignment.
* docs: add release note for #674
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Severin Dicks <37635888+Intron7@users.noreply.github.com>
@ilan-goldilan-gold modified the milestones: 1.12.2, 1.12.3Jun 29, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.3, 1.12.4Jul 24, 2026
@flying-sheepflying-sheep modified the milestones: 1.12.4, 1.12.5Aug 27, 2026
@flying-sheep

Copy link
Copy Markdown
Member

@giovp you initialized the “experimental” idea. What’s the stabilization plan? What’s the criteria to move these out?

Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants

@Arshammik@Zethson@flying-sheep@ilan-gold