Skip to content

Fix/correlation csr guanlab pbmc multiome - #83

Open
VladimirShitov wants to merge 4 commits into
mainfrom
fix/correlation-csr-guanlab-pbmc-multiome
Open

VladimirShitov wants to merge 4 commits into
mainfrom
fix/correlation-csr-guanlab-pbmc-multiome

Conversation

@VladimirShitov

Copy link
Copy Markdown
Contributor

Describe your changes

This PR fixes the wrong Spearman scores and the three failed runs of 0.2.0, and seeds every method so a rerun reproduces its predictions.

1. correlation: Spearman scores of CSR matrices

dynutils::spearman_rank_sparse() ranks along the compressed dimension, which is a column only for CSC. A CSR matrix (a dgRMatrix in R) was therefore ranked per feature for the per-cell correlation and per cell for the per-gene one. In 0.2.0 this hit:

  • ss_opm, which writes CSR, on every dataset;
  • every method on GEX → ATAC, because the solution file stores ATAC layers as CSR;
  • random_predict (and solution), which write CSR on GEX → ATAC.

The metric now converts both matrices to CSC after reading them. A new component test (test_prediction_layout.py) scores the test prediction as CSR and as CSC, and checks both against Spearman correlations of dense ranks. It fails on main, where the CSR prediction gets 0.137 per-cell Spearman instead of 0.367, and passes here.

Scoring all 176 predictions of the rerun below with both versions of the metric, only mean_spearman_per_cell and mean_spearman_per_gene change:

changed largest difference
ss_opm, every dataset 0.03 – 0.33
every other method, GEX → ATAC only ≤ 7.2e-4
all other scores (Pearson, overall correlations, other datasets) 0

The fixed scores of ss_opm equal those of dense ranks exactly:

dataset per-cell Spearman, 0.2.0 → fixed per-gene Spearman, 0.2.0 → fixed
bmmc_cite GEX → ADT 0.509 → 0.724 0.257 → 0.582
bmmc_cite ADT → GEX 0.136 → 0.408 0.077 → 0.175
bmmc_multiome GEX → ATAC 0.063 → 0.115 0.033 → 0.095
bmmc_multiome ATAC → GEX 0.097 → 0.300 0.078 → 0.133
pbmc_cite GEX → ADT 0.327 → 0.528 0.306 → 0.396
pbmc_cite ADT → GEX 0.353 → 0.611 0.028 → 0.131
pbmc_multiome GEX → ATAC 0.059 → 0.091 0.012 → 0.054
pbmc_multiome ATAC → GEX 0.223 → 0.430 0.041 → 0.128

2. guanlab_dengkw_pm: failed on both pbmc_multiome datasets

From the 0.2.0 run's logs:

  • ATAC → GEX: failed twice with numpy.linalg.LinAlgError: Internal potrf return info = [49313]. The Cholesky of the 32-bit-integer OpenBLAS that scipy's wheels link breaks down once a half of the batches holds more than ~30–46k cells; on pbmc_multiome the halves hold 56k–74k cells. The kernel ridge regression is now solved with numpy.linalg.solve(), through numpy's 64-bit-integer OpenBLAS. Where scipy's Cholesky works, the coefficients of both agree to 6e-13.
  • GEX → ATAC: stopped at hightime's 8 hours (exit 140) during the 7th of its 10 fits.
  • New labels: the component asks for veryhightime and veryhighmem. In the rerun ATAC → GEX peaked at 112 GB and took 2.5 h, and GEX → ATAC took 6.6 h on 30 cores.

3. simple_mlp: failed on pbmc_multiome ATAC → GEX

split() densified whole folds: 143 GB for the ATAC input of pbmc_multiome. The inputs and targets now stay sparse, and only one batch at a time is densified (utils.SparseRows + utils.batch_loader()), drawing the same batches as DataLoader(shuffle=True) did. The unused counts layers are dropped after reading. Training on pbmc_multiome ATAC → GEX now takes 25 min.

4. Seeding: reruns reproduce the predictions

Before this PR, rerunning the benchmark changed many scores by up to 0.02. These steps were unseeded:

  • the truncated SVDs of guanlab_dengkw_pm (its np.random.seed() came after them) and knnr_py;
  • the lmds landmarks of lm and knnr_r;
  • the draws of random_predict;
  • the scvi-tools model and pynndescent search of cellmapper_scvi and cellmapper_linear;
  • the training of simple_mlp, novel, babel and the senkin_tmp networks.

The new src/utils/seed.py seeds python and numpy, plus torch and tensorflow when a method has imported them, and switches them to deterministic GPU algorithms (cudnn.deterministic, a fixed cuBLAS workspace, use_deterministic_algorithms(warn_only=True), tensorflow op determinism). ss_opm, scipenn and scbutterfly were already seeded and now also use deterministic GPU algorithms. The R components call set.seed(). Methods that had no seed argument get --seed (default 0).

Reproducibility is bit-for-bit only on the same GPU model. CPU methods reproduce exactly anywhere. GPU-trained models reproduce exactly only on the same GPU model: A100 and H100 differ by ~1e-3. A rerun on other GPUs, such as de.NBI's T4s, will therefore move the GPU-trained methods by about that much.

Checklist before requesting a review

  • I have performed a self-review of my code

  • Check the correct box. Does this PR contain:

    • Breaking changes
    • New functionality
    • Major changes
    • Minor changes
    • Bug fixes
  • Proposed changes are described in the CHANGELOG.md

  • CI Tests succeed and look good!

VladimirShitov and others added 4 commits October 25, 2026 09:50
…rman scores

dynutils::spearman_rank_sparse() ranks along the compressed dimension, which is a
column only for CSC. ss_opm writes its predictions as CSR, so its per-cell and
per-gene Spearman correlations were computed on ranks along the wrong axis on every
dataset. Convert the solution and the prediction to CSC after reading them, and add
a component test that scores a prediction as CSR and as CSC against dense ranks.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…bit-integer LAPACK

scipy's 32-bit-integer OpenBLAS Cholesky reports the kernel matrix as not positive
definite once a half of the batches holds more than ~30-46k cells, so the method
failed on pbmc_multiome ATAC -> GEX in 0.2.0. Also ask for veryhighmem (118 GB peak)
and veryhightime (GEX -> ATAC hit hightime's 8 hours).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
split() densified whole training and validation folds (143 GB for the ATAC input of
pbmc_multiome), so ATAC -> GEX failed on it in every run since run_2026-08-14.
utils.SparseRows serves one batch at a time as dense tensors through a BatchSampler,
drawing the same batches as DataLoader(shuffle=True) did; the unused counts layers
are dropped after reading.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Rerunning the benchmark did not reproduce its predictions: guanlab_dengkw_pm's and
knnr_py's truncated SVDs, the lmds landmarks of lm and knnr_r, random_predict, the
scvi-tools and pynndescent steps of cellmapper, and the training of simple_mlp, novel,
babel and senkin_tmp's networks were unseeded. src/utils/seed.py seeds python, numpy,
torch and tensorflow and makes them use deterministic GPU algorithms (also for the
already seeded ss_opm, scipenn and scbutterfly); the R components call set.seed().
Methods without one get a --seed argument (default 0).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Comment on lines +51 to +54
- name: "--seed"
type: integer
default: 0
description: Seed for the weight initialisation, the validation split and the order of the training batches.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if you want to set the seed in all methods, add it to src/api/comp_method.yaml. that way we can also override it in the workflow more easily, for when we do wantr different seeds.

Comment thread CHANGELOG.md
Comment on lines +5 to +13
* Seed every stochastic method and control, so that rerunning the benchmark reproduces its predictions. Unseeded in 0.2.0 were the truncated SVDs of `guanlab_dengkw_pm` (its `np.random.seed()` came after them) and `knnr_py`, the `lmds` landmarks of `lm` and `knnr_r`, the draws of `random_predict`, the scvi-tools model and pynndescent search of `cellmapper_scvi` and `cellmapper_linear`, and the training of `simple_mlp`, `novel`, `babel` and the `senkin_tmp` networks. The new `src/utils/seed.py` seeds python, numpy, torch and tensorflow and makes them use deterministic GPU algorithms, which `ss_opm`, `scipenn` and `scbutterfly`, already seeded, now use too; the R components call `set.seed()`. Methods without one get a `--seed` argument (default 0).

## BUG FIXES

* `correlation`: Convert the solution and the prediction to CSC before correlating them. `dynutils::spearman_rank_sparse()` ranks along the compressed dimension, which is a column only for CSC, so a CSR matrix (read as a `dgRMatrix`) was ranked per feature for the per-cell Spearman correlation and per cell for the per-gene one. In 0.2.0 this hit `mean_spearman_per_cell` and `mean_spearman_per_gene` of `ss_opm`, which writes CSR, on every dataset, and of every method on GEX -> ATAC, whose solution holds the ATAC layer as CSR (`random_predict` also writes CSR there). All other scores are unchanged. A new component test scores a prediction as CSR and as CSC and checks both against dense ranks.

* `guanlab_dengkw_pm`: Solve the kernel ridge regression with `numpy.linalg.solve()` instead of `scipy.linalg.cho_factor()`/`cho_solve()`. The Cholesky of the 32-bit-integer OpenBLAS that scipy's wheels link reports the kernel matrix as not positive definite once a half of the batches holds more than ~30-46k cells (`LinAlgError: Internal potrf return info = [49313]` on `pbmc_multiome` ATAC -> GEX in 0.2.0, where the halves hold 56k-74k cells); numpy's wheels link a 64-bit-integer OpenBLAS whose LU solve handles these systems. Also ask for `veryhighmem` and `veryhightime` rather than `highmem` and `hightime`: ATAC -> GEX peaks at 112-118 GB, and GEX -> ATAC, which 0.2.0 stopped at `hightime`'s 8 hours, takes 6.5-8 hours on 30 cores.

* `simple_mlp_train`, `simple_mlp_predict`: Keep the inputs and targets sparse and densify one batch at a time (`utils.SparseRows` + `utils.batch_loader()`), and drop the unused `counts` layers after reading. `split()` densified whole folds, 143 GB for the ATAC input of `pbmc_multiome`, so ATAC -> GEX failed on it in every run since run_2026-08-14, 0.2.0 included. Batches are drawn as before.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

make this less wordy.

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

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants