diff --git a/pyproject.toml b/pyproject.toml
index f7089d4c..ea329e0b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -13,7 +13,7 @@ dependencies = [
     "pandas>=3.0.0",
     "pyarrow>=19.0",
     "tensorly>=0.9.0",
-    "parafac2>=1.6.0",
+    "parafac2>=1.7.0",
     "anndata>=0.13",
     "pacmap>=0.9",
     "tqdm>=4.66.1",
diff --git a/scrise/rank_selection.py b/scrise/rank_selection.py
index a53b5a2a..2621d8c4 100644
--- a/scrise/rank_selection.py
+++ b/scrise/rank_selection.py
@@ -18,6 +18,8 @@
 import numpy as np
 import pandas as pd
 import scipy.sparse as sps
+from parafac2.compress import compress_dataset
+from parafac2.parafac2 import parafac2_nd
 from parafac2.utils import calc_W, condition_slices, project_data, rmatmul
 from tqdm import tqdm
 
@@ -195,9 +197,89 @@ def _test_block_moments(
     return sums[gene_idx], squares[gene_idx]
 
 
+def _fit_at_ranks(
+    X_in: anndata.AnnData,
+    ranks: Sequence[int],
+    tolerance: float,
+    max_iter: int,
+    compress: int | tuple[int, int | None] | str | bool | None,
+    compression_kwarg: dict[str, Any] | None,
+    parafac2_kwarg: dict[str, Any] | None,
+    rng: np.random.Generator,
+) -> list[tuple[tuple, float]]:
+    """Fit PARAFAC2 against ``X_in`` at every rank in ``ranks``.
+
+    When ``compression_kwarg`` is given, ``X_in`` is compressed *once*,
+    explicitly, sized for the *largest* rank in ``ranks`` (CANDELINC
+    compression to ``L`` dimensions is valid for fitting any rank ``<= L``,
+    not just the rank it was sized for -- see
+    :func:`parafac2.compress.compress_dataset`/:class:`~parafac2.compress.CompressedData`),
+    and every rank's fit reuses that same compressed representation. Compression
+    is the expensive, ``O(nnz)`` raw-data pass; the fit itself
+    (:func:`~parafac2.parafac2.parafac2_nd` against an already-compressed
+    ``CompressedData``) only touches the small dense compressed cores, so
+    this turns what used to be one full compression *per rank* into one
+    compression for the whole set.
+
+    Without ``compression_kwarg``, each rank's fit instead goes through
+    :func:`~scrise._pf2_utils.run_parafac2` exactly as before -- unchanged,
+    since ``parafac2_nd``'s own internal per-call compression shortcut
+    doesn't hand back a reusable ``CompressedData`` object to hoist out of
+    the loop.
+
+    Returns a list of ``((weights, (A, B, C), projections), r2x)`` tuples,
+    one per rank, in the same order as ``ranks``.
+    """
+    if not compression_kwarg:
+        return [
+            run_parafac2(
+                X_in,
+                rank=rank,
+                random_state=int(rng.integers(np.iinfo(np.int32).max)),
+                tol=tolerance,
+                n_iter_max=max_iter,
+                compress=compress,
+                compression_kwarg=None,
+                parafac2_kwarg=parafac2_kwarg,
+            )
+            for rank in ranks
+        ]
+
+    if not compress:
+        raise ValueError("compression_kwarg requires compress to be set.")
+
+    parafac2_kwarg = dict(parafac2_kwarg) if parafac2_kwarg else {}
+    normalize_slices = parafac2_kwarg.pop("normalize_slices", False)
+    backend = parafac2_kwarg.pop("backend", None)
+
+    compressed = compress_dataset(
+        X_in,
+        L=compress,
+        rank=max(ranks),
+        random_state=int(rng.integers(np.iinfo(np.int32).max)),
+        normalize_slices=normalize_slices,
+        backend=backend,
+        **compression_kwarg,
+    )
+    return [
+        parafac2_nd(
+            compressed,
+            rank=rank,
+            random_state=int(rng.integers(np.iinfo(np.int32).max)),
+            tol=tolerance,
+            n_iter_max=max_iter,
+            normalize_slices=normalize_slices,
+            backend=backend,
+            compress=None,  # already compressed above
+            **parafac2_kwarg,
+        )
+        for rank in ranks
+    ]
+
+
 def _bicv_trial(
     X: anndata.AnnData,
-    rank: int,
+    ranks: Sequence[int],
     held_out_cell_frac: float,
     held_out_gene_frac: float,
     seed: int,
@@ -206,12 +288,14 @@ def _bicv_trial(
     compress: int | tuple[int, int | None] | str | bool | None = "auto",
     compression_kwarg: dict[str, Any] | None = None,
     parafac2_kwarg: dict[str, Any] | None = None,
-) -> dict[str, float]:
-    """Run a single bi-cross-validation trial and return its scores and shape.
+) -> list[dict[str, float]]:
+    """Run a single bi-cross-validation trial, evaluated at every rank in ``ranks``.
 
-    Splits cells (stratified by condition) and genes into train/test blocks,
-    fits PARAFAC2 on the train-cell x train-gene block, then predicts the
-    held-out test-cell x test-gene block:
+    Splits cells (stratified by condition) and genes into train/test blocks
+    *once*, fits PARAFAC2 on the train-cell x train-gene block at each rank
+    (compressing that train block once, not once per rank -- see
+    :func:`_fit_at_ranks`), then predicts the held-out test-cell x test-gene
+    block from each rank's fit:
 
     - Gene loadings for the held-out genes are estimated by regressing the
       train cells' expression of those genes onto the fitted (train-cell)
@@ -222,12 +306,16 @@ def _bicv_trial(
 
     R2X is then computed by reconstructing the held-out block from these
     estimates and comparing against the (mean-centered) observed values.
+    That reference block (and its moments) depends only on the split, not on
+    the rank, so it too is computed once and reused across every rank below.
 
     Returns
     -------
-    dict[str, float]
-        ``BiCV R2X`` (the held-out block), ``Train Block R2X`` (in-sample, on
-        the block the model was fit to), the four block sizes, and ``Seed``.
+    list[dict[str, float]]
+        One dict per rank (in the same order as ``ranks``), each with
+        ``Rank``, ``BiCV R2X`` (the held-out block), ``Train Block R2X``
+        (in-sample, on the block the model was fit to), the four block
+        sizes, and ``Seed``.
     """
     rng = np.random.default_rng(seed)
     cond_idx = X.obs["condition_unique_idxs"].to_numpy().astype(int)
@@ -239,17 +327,16 @@ def _bicv_trial(
     train_gene_mask, test_gene_mask = _split_genes(X.n_vars, held_out_gene_frac, rng)
 
     X_train = X[train_cell_mask][:, train_gene_mask].copy()
-    (weights, (A, B, C), P_train), train_block_r2x = run_parafac2(
+    fits = _fit_at_ranks(
         X_train,
-        rank=rank,
-        random_state=int(rng.integers(np.iinfo(np.int32).max)),
-        tol=tolerance,
-        n_iter_max=max_iter,
-        compress=compress,
-        compression_kwarg=compression_kwarg,
-        parafac2_kwarg=parafac2_kwarg,
+        ranks,
+        tolerance,
+        max_iter,
+        compress,
+        compression_kwarg,
+        parafac2_kwarg,
+        rng,
     )
-    A = A * weights
 
     # Everything below reaches the raw data through products against the raw
     # matrix rather than materialising a block of it. Cells are restricted by
@@ -267,38 +354,17 @@ def _bicv_trial(
     n_genes = X.n_vars
     test_gene_idx = np.flatnonzero(test_gene_mask)
     means_test_genes = means[test_gene_mask]
-
-    # Estimate gene loadings for the held-out genes from the train cells.
-    #   Z^T (X[train, test] - 1 mu^T) = (Z^T X[train])[:, test] - (Z^T 1) mu^T
     cond_train = cond_idx[train_cell_mask]
-    Z = _cell_loadings(P_train, B, A, cond_train, n_cond)
-    ZtY = np.asarray(
-        rmatmul(np.ascontiguousarray(Z.T), _restrict_rows(X_mat, train_cell_mask)),
-        dtype=np.float64,
-    )[:, test_gene_idx]
-    ZtY -= np.outer(Z.sum(axis=0), means_test_genes)
-    # `lstsq` on the rank x rank system keeps the minimum-norm
-    # behaviour when the fit is rank deficient.
-    C_test = np.linalg.lstsq(Z.T @ Z, ZtY, rcond=None)[0].T
-
-    # Estimate projections for the held-out cells from the train genes.
-    C_full = np.zeros((n_genes, C.shape[1]))
-    C_full[train_gene_mask] = C
     cond_test = cond_idx[test_cell_mask]
-    X_test = _restrict_rows(X_mat, test_cell_mask)
-    W_test = calc_W(X_test, means, C_full)
     cond_slices_test = condition_slices(cond_test, n_cond)
-    P_test, _ = project_data(W_test, [A, B, C], cond_slices_test)
-
-    # Score the held-out block. `A` carries the training slice's energy, so the
-    # held-out loadings need rescaling for the held-out slice's size.
-    L = _cell_loadings(P_test, B, A, cond_test, n_cond)
-    L *= _holdout_scale(cond_train, cond_test, n_cond)
-    LtY = np.asarray(rmatmul(np.ascontiguousarray(L.T), X_test), dtype=np.float64)[
-        :, test_gene_idx
-    ]
-    LtY -= np.outer(L.sum(axis=0), means_test_genes)
-
+    scale = _holdout_scale(cond_train, cond_test, n_cond)
+
+    # These depend only on the (fixed, once-per-trial) split, not on rank, so
+    # -- unlike compression/fitting above -- they were already shared across
+    # ranks even before this change; now made explicit and computed once here
+    # rather than once per rank.
+    X_train_lazy = _restrict_rows(X_mat, train_cell_mask)
+    X_test_lazy = _restrict_rows(X_mat, test_cell_mask)
     col_sums, col_squares = _test_block_moments(X_mat, test_cell_mask, test_gene_idx)
     n_test_cells = int(test_cell_mask.sum())
     ss_tot = float(
@@ -308,19 +374,57 @@ def _bicv_trial(
             + n_test_cells * means_test_genes**2
         )
     )
-    cross = float(np.sum(C_test.T * LtY))
-    ss_fit = float(np.sum((L.T @ L) * (C_test.T @ C_test)))
-    ss_res = ss_tot - 2.0 * cross + ss_fit
 
-    return {
-        "BiCV R2X": 1.0 - ss_res / ss_tot,
-        "Train Block R2X": float(train_block_r2x),
-        "NTrainGenes": int(train_gene_mask.sum()),
-        "NTestGenes": int(test_gene_mask.sum()),
-        "NTrainCells": int(train_cell_mask.sum()),
-        "NTestCells": int(test_cell_mask.sum()),
-        "Seed": int(seed),
-    }
+    results = []
+    for rank, ((weights, (A, B, C), P_train), train_block_r2x) in zip(
+        ranks, fits, strict=True
+    ):
+        A = A * weights
+
+        # Estimate gene loadings for the held-out genes from the train cells.
+        #   Z^T (X[train, test] - 1 mu^T) = (Z^T X[train])[:, test] - (Z^T 1) mu^T
+        Z = _cell_loadings(P_train, B, A, cond_train, n_cond)
+        ZtY = np.asarray(
+            rmatmul(np.ascontiguousarray(Z.T), X_train_lazy), dtype=np.float64
+        )[:, test_gene_idx]
+        ZtY -= np.outer(Z.sum(axis=0), means_test_genes)
+        # `lstsq` on the rank x rank system keeps the minimum-norm
+        # behaviour when the fit is rank deficient.
+        C_test = np.linalg.lstsq(Z.T @ Z, ZtY, rcond=None)[0].T
+
+        # Estimate projections for the held-out cells from the train genes.
+        C_full = np.zeros((n_genes, C.shape[1]))
+        C_full[train_gene_mask] = C
+        W_test = calc_W(X_test_lazy, means, C_full)
+        P_test, _ = project_data(W_test, [A, B, C], cond_slices_test)
+
+        # Score the held-out block. `A` carries the training slice's energy,
+        # so the held-out loadings need rescaling for the held-out slice's
+        # size.
+        L = _cell_loadings(P_test, B, A, cond_test, n_cond)
+        L = L * scale
+        LtY = np.asarray(
+            rmatmul(np.ascontiguousarray(L.T), X_test_lazy), dtype=np.float64
+        )[:, test_gene_idx]
+        LtY -= np.outer(L.sum(axis=0), means_test_genes)
+
+        cross = float(np.sum(C_test.T * LtY))
+        ss_fit = float(np.sum((L.T @ L) * (C_test.T @ C_test)))
+        ss_res = ss_tot - 2.0 * cross + ss_fit
+
+        results.append(
+            {
+                "Rank": rank,
+                "BiCV R2X": 1.0 - ss_res / ss_tot,
+                "Train Block R2X": float(train_block_r2x),
+                "NTrainGenes": int(train_gene_mask.sum()),
+                "NTestGenes": int(test_gene_mask.sum()),
+                "NTrainCells": int(train_cell_mask.sum()),
+                "NTestCells": int(test_cell_mask.sum()),
+                "Seed": int(seed),
+            }
+        )
+    return results
 
 
 def _resolve_bicv_inputs(
@@ -394,11 +498,22 @@ def bicv(
 
     For each candidate rank, computes both the ordinary in-sample fit R2X
     (using the full dataset, as in :func:`scrise.factorization.rise_pca_r2x`)
-    and the BiCV R2X (repeated ``n_repeats`` times with independent random
-    cell/gene splits). The fit R2X increases monotonically with rank; the
-    BiCV R2X penalizes overfitting and typically peaks near the rank that
-    best generalizes to held-out data. Plot both with
-    :func:`scrise.plotting.plot_bicv_r2x` to select a rank.
+    and the BiCV R2X (``n_repeats`` independent random cell/gene splits,
+    each evaluated at every rank -- see :func:`_bicv_trial`). The fit R2X
+    increases monotonically with rank; the BiCV R2X penalizes overfitting and
+    typically peaks near the rank that best generalizes to held-out data.
+    Plot both with :func:`scrise.plotting.plot_bicv_r2x` to select a rank.
+
+    Each of the ``n_repeats`` splits (and the full dataset, for the in-sample
+    fit) is compressed once -- sized for the *largest* rank requested -- and
+    every rank's PARAFAC2 fit reuses that same compression, rather than
+    recompressing per rank (see :func:`_fit_at_ranks`). Compression is by far
+    the most expensive, ``O(nnz)`` step against the raw data; fitting a
+    smaller rank from an already-compressed representation only touches the
+    small dense compressed cores. This applies whenever ``compression_kwarg``
+    is given (as it must be to control e.g. ``n_power_iter``); without it,
+    each rank still goes through its own call into ``parafac2_nd``'s internal
+    compression shortcut, unchanged.
 
     Parameters
     ----------
@@ -409,8 +524,9 @@ def bicv(
     ranks : sequence of int
         Candidate rank values to evaluate (e.g., [5, 10, 15, 20, 25, 30]).
     n_repeats : int, optional (default: 3)
-        Number of independent random cell/gene splits per rank. Higher
-        values give a less noisy BiCV estimate but take longer.
+        Number of independent random cell/gene splits, each evaluated at
+        every rank in ``ranks``. Higher values give a less noisy BiCV
+        estimate but take longer.
     held_out_cell_frac : float, optional (default: 0.5)
         Fraction of cells held out per condition in each BiCV trial.
     held_out_gene_frac : float, optional (default: 0.5)
@@ -470,35 +586,38 @@ def bicv(
 
     rng = np.random.default_rng(random_state)
     rows = []
-    for rank in tqdm(ranks, desc="BiCV rank selection"):
-        _, fit_r2x = run_parafac2(
-            X,
-            rank=rank,
-            random_state=int(rng.integers(np.iinfo(np.int32).max)),
-            tol=tolerance,
-            n_iter_max=max_iter,
-            compress=compress,
-            compression_kwarg=compression_kwarg,
-            parafac2_kwarg=parafac2_kwarg,
-        )
+
+    # The full (unsplit) dataset is the same for every rank, so -- like the
+    # per-trial train blocks below -- it's compressed once (sized for the
+    # largest rank) and reused, rather than once per rank.
+    fit_results = _fit_at_ranks(
+        X, ranks, tolerance, max_iter, compress, compression_kwarg, parafac2_kwarg, rng
+    )
+    for rank, (_, fit_r2x) in zip(ranks, fit_results, strict=True):
         # The full-data fit has no split, so the per-trial columns are absent
         # here rather than zero.
         rows.append({"Rank": rank, "Repeat": 0, "Metric": "Fit R2X", "R2X": fit_r2x})
 
-        for repeat in range(n_repeats):
-            trial_seed = int(rng.integers(np.iinfo(np.int32).max))
-            trial = _bicv_trial(
-                X,
-                rank,
-                held_out_cell_frac,
-                held_out_gene_frac,
-                trial_seed,
-                tolerance,
-                max_iter,
-                compress,
-                compression_kwarg,
-                parafac2_kwarg,
-            )
+    # One held-out split per repeat, evaluated at every rank (see
+    # `_bicv_trial`), rather than one independent split per (rank, repeat)
+    # pair -- the split, and the compression of its train block, no longer
+    # need to be redone for each rank.
+    for repeat in tqdm(range(n_repeats), desc="BiCV repeats"):
+        trial_seed = int(rng.integers(np.iinfo(np.int32).max))
+        trial_rows = _bicv_trial(
+            X,
+            ranks,
+            held_out_cell_frac,
+            held_out_gene_frac,
+            trial_seed,
+            tolerance,
+            max_iter,
+            compress,
+            compression_kwarg,
+            parafac2_kwarg,
+        )
+        for trial in trial_rows:
+            rank = trial.pop("Rank")
             rows.append(
                 {
                     "Rank": rank,
diff --git a/scrise/tests/test_rank_selection.py b/scrise/tests/test_rank_selection.py
index 16da0e25..6e1a0474 100644
--- a/scrise/tests/test_rank_selection.py
+++ b/scrise/tests/test_rank_selection.py
@@ -157,6 +157,71 @@ def test_bicv_parafac2_kwarg_and_compression_kwarg():
         )
 
 
+def test_bicv_compresses_once_per_trial_not_once_per_rank(monkeypatch):
+    """The whole point of the redesign: one compression per held-out split
+    (and one for the in-sample fit), reused across every rank -- not one
+    compression per (rank, repeat) pair as before.
+    """
+    import scrise.rank_selection as rs
+
+    X = _make_test_data()
+    ranks = [2, 4, 6]
+    n_repeats = 2
+
+    calls: list[int] = []
+    real_compress_dataset = rs.compress_dataset
+
+    def _counting_compress_dataset(*args, **kwargs):
+        calls.append(kwargs["rank"])
+        return real_compress_dataset(*args, **kwargs)
+
+    monkeypatch.setattr(rs, "compress_dataset", _counting_compress_dataset)
+
+    bicv(
+        X,
+        ranks,
+        n_repeats=n_repeats,
+        random_state=0,
+        max_iter=20,
+        compress="auto",
+        compression_kwarg={"n_power_iter": 1},
+    )
+
+    # One compression for the in-sample fit, plus one per repeat -- not
+    # one per (rank, repeat) pair (which would be len(ranks) * (1 + n_repeats)).
+    assert len(calls) == 1 + n_repeats
+    # Every compression is sized for the *largest* requested rank.
+    assert calls == [max(ranks)] * len(calls)
+
+
+def test_bicv_without_compression_kwarg_still_compresses_per_rank(monkeypatch):
+    """Without compression_kwarg there's no explicit compress_dataset call to
+    hoist -- parafac2_nd's own internal shortcut still runs once per rank,
+    exactly as before this change.
+    """
+    import scrise.rank_selection as rs
+
+    X = _make_test_data()
+    ranks = [2, 4]
+    n_repeats = 1
+
+    calls: list[int] = []
+    real_run_parafac2 = rs.run_parafac2
+
+    def _counting_run_parafac2(*args, **kwargs):
+        calls.append(kwargs["rank"])
+        return real_run_parafac2(*args, **kwargs)
+
+    monkeypatch.setattr(rs, "run_parafac2", _counting_run_parafac2)
+
+    bicv(X, ranks, n_repeats=n_repeats, random_state=0, max_iter=20, compress="auto")
+
+    # One run_parafac2 call per rank for the in-sample fit, plus one per
+    # (rank, repeat) pair for BiCV trials -- unchanged from before.
+    assert len(calls) == len(ranks) + len(ranks) * n_repeats
+    assert sorted(calls) == sorted(ranks * (1 + n_repeats))
+
+
 def test_bicv_adata_alias():
     X = _make_test_data()
     results = bicv(adata=X, ranks=[2], n_repeats=1, random_state=0, max_iter=10)
@@ -325,13 +390,13 @@ def test_reported_seed_reproduces_its_own_trial():
     target = rows.iloc[1]
     replayed = _bicv_trial(
         X,
-        rank=3,
+        ranks=[3],
         held_out_cell_frac=0.5,
         held_out_gene_frac=0.5,
         seed=int(target["Seed"]),
         tolerance=1e-6,
         max_iter=60,
-    )
+    )[0]
     assert replayed["BiCV R2X"] == pytest.approx(target["R2X"], rel=1e-9)
     assert replayed["Train Block R2X"] == pytest.approx(
         target["Train Block R2X"], rel=1e-9
@@ -415,13 +480,13 @@ def test_streamed_scoring_matches_the_dense_formulation(sparse, seed):
 
     got = _bicv_trial(
         X,
-        rank=3,
+        ranks=[3],
         held_out_cell_frac=0.2,
         held_out_gene_frac=0.2,
         seed=seed,
         tolerance=1e-6,
         max_iter=60,
-    )["BiCV R2X"]
+    )[0]["BiCV R2X"]
     want = _dense_reference_trial(X, rank=3, seed=seed)
     assert got == pytest.approx(want, rel=1e-6, abs=1e-9)
 
@@ -436,13 +501,13 @@ def test_streamed_scoring_matches_on_interleaved_conditions():
 
     got = _bicv_trial(
         X,
-        rank=3,
+        ranks=[3],
         held_out_cell_frac=0.2,
         held_out_gene_frac=0.2,
         seed=1,
         tolerance=1e-6,
         max_iter=60,
-    )["BiCV R2X"]
+    )[0]["BiCV R2X"]
     want = _dense_reference_trial(X, rank=3, seed=1)
     assert got == pytest.approx(want, rel=1e-6, abs=1e-9)
 
@@ -705,7 +770,7 @@ def test_exactly_low_rank_data_scores_near_one_at_its_own_rank():
     held-out score still said it was worse than predicting the mean.
     """
     adata = _exact_low_rank(rank_true=4)
-    trial = _bicv_trial(adata, 4, 0.5, 0.5, seed=0, tolerance=1e-8, max_iter=300)
+    trial = _bicv_trial(adata, [4], 0.5, 0.5, seed=0, tolerance=1e-8, max_iter=300)[0]
 
     assert trial["Train Block R2X"] > 0.99
     assert trial["BiCV R2X"] > 0.95
@@ -719,7 +784,7 @@ def test_score_barely_moves_with_the_held_out_fraction(frac):
     same data at the same rank, because the error was sqrt(n_train/n_test).
     """
     adata = _exact_low_rank(rank_true=4)
-    trial = _bicv_trial(adata, 4, frac, frac, seed=0, tolerance=1e-8, max_iter=300)
+    trial = _bicv_trial(adata, [4], frac, frac, seed=0, tolerance=1e-8, max_iter=300)[0]
     assert trial["BiCV R2X"] > 0.9
 
 
@@ -743,7 +808,7 @@ def test_pure_noise_does_not_score_positive():
     adata.obs["condition_unique_idxs"] = cond
     adata.var["means"] = np.zeros(60)
 
-    trial = _bicv_trial(adata, 5, 0.5, 0.5, seed=0, tolerance=1e-8, max_iter=200)
+    trial = _bicv_trial(adata, [5], 0.5, 0.5, seed=0, tolerance=1e-8, max_iter=200)[0]
     assert trial["BiCV R2X"] < 0.02
 
 
diff --git a/uv.lock b/uv.lock
index dd9b6e6f..c3015e90 100644
--- a/uv.lock
+++ b/uv.lock
@@ -271,6 +271,48 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/ec/82/32e3bd191d498e64f6f911ad55d14006a0861e54869d2d32452326399e65/coverage-7.15.2-py3-none-any.whl", hash = "sha256:eb6bcae8d1a9d305351ecb108232441d11c5cfe9de840a04388ba5d2db8d735c", size = 213375, upload-time = "2026-07-15T18:56:17.305Z" },
 ]
 
+[[package]]
+name = "cuda-bindings"
+version = "13.4.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "cuda-pathfinder" },
+]
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/38/cf/165b4d449f94956c2a60930cf5dfeb27132ead60a7e7f2c37819df1cba07/cuda_bindings-13.4.1-cp313-cp313-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:b601c0cbf0dffb648f68e56a60b320738a20210293f33896a1964a6438cc65f1", size = 6313772, upload-time = "2026-09-10T01:17:03.969Z" },
+    { url = "https://files.pythonhosted.org/packages/d9/f9/cf021d1560541caa1f35f3e7e311d2678dbacb4fa6a4573b63470fe1ae00/cuda_bindings-13.4.1-cp313-cp313-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:0e2c357698588b06ebd65811ee2013b6650dd0a10d16d899924aabad0d606d76", size = 6924300, upload-time = "2026-09-10T01:17:06.23Z" },
+    { url = "https://files.pythonhosted.org/packages/66/93/d41e182cde0f1dd3c6c6aac018cb550ffedf55ba2a67c427320e91048592/cuda_bindings-13.4.1-cp313-cp313-win_amd64.whl", hash = "sha256:11bce22fe70a19eb6e727ae7126940245160e0223b3f1cb3e209a313449f75ea", size = 6289385, upload-time = "2026-09-10T01:17:08.372Z" },
+    { url = "https://files.pythonhosted.org/packages/d0/cd/e2100b3d5e2b2e023a603a6463df4582c73d057daed9e7cf6802cd3b8455/cuda_bindings-13.4.1-cp313-cp313-win_arm64.whl", hash = "sha256:73d86523ee2f749cf312fdada7be83160c84ddf4845232769f044695ec1c0e68", size = 5293378, upload-time = "2026-09-10T01:17:10.356Z" },
+    { url = "https://files.pythonhosted.org/packages/58/17/74346b49114779920929ec0ea1361f0a0357262d6020cfdd017f670da638/cuda_bindings-13.4.1-cp314-cp314-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:d7df2dddb81feb15787e8c4a13b7aa3f4c23eabeec7746037e94dbbed065fc6c", size = 6406821, upload-time = "2026-09-10T01:17:12.344Z" },
+    { url = "https://files.pythonhosted.org/packages/f6/77/2f9a38be7399a34e3703b6ce1be60bb2c56611d324750f8d544ff90b0473/cuda_bindings-13.4.1-cp314-cp314-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:62df23df11074e9833bf348bbcf0b8eec2fcbded4f305c6fbaa3ed067433e97d", size = 6972525, upload-time = "2026-09-10T01:17:15.098Z" },
+    { url = "https://files.pythonhosted.org/packages/fb/88/28195fa9f0c7aebf9978b5508e463289091a574240c4c50f9035b1d96764/cuda_bindings-13.4.1-cp314-cp314-win_amd64.whl", hash = "sha256:fb66fa7ea3f257a0451d54edc4b5a4a052eb29c32ba15b104f465ecb24907ffb", size = 6434166, upload-time = "2026-09-10T01:17:17.097Z" },
+    { url = "https://files.pythonhosted.org/packages/2d/f5/41388671cbe81477aa30fbccbf195ca6f207eeaca5b1268282f1aba827cd/cuda_bindings-13.4.1-cp314-cp314-win_arm64.whl", hash = "sha256:506457cc214f3dac300044431a36b3f608f7a8c495c5bb40b9008d162ff1d5d3", size = 5450466, upload-time = "2026-09-10T01:17:19.099Z" },
+    { url = "https://files.pythonhosted.org/packages/7d/de/41197aebf91f6c5f82b35e06e3b4bbe08edddb335c4d6e53c1f1fe542e0f/cuda_bindings-13.4.1-cp314-cp314t-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:0947a6c491a622b076e9afcbc0444eefd9cbe001630558bbcba8d0a90dab7fe7", size = 6243282, upload-time = "2026-09-10T01:17:21.158Z" },
+    { url = "https://files.pythonhosted.org/packages/be/d6/1b697092f53cfd4d721fbbe13e7d66dc72a9d508156af9be4c86607e4373/cuda_bindings-13.4.1-cp314-cp314t-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:51e730e80d997b6037566033a79aae3995e1981654054536f512e5c091679c70", size = 6803245, upload-time = "2026-09-10T01:17:23.309Z" },
+    { url = "https://files.pythonhosted.org/packages/f4/19/c4e0cb61555f553bc01a533f558b2ba775c88b0df5ac3dfba70125948712/cuda_bindings-13.4.1-cp314-cp314t-win_amd64.whl", hash = "sha256:2a95d8e231d219c20beed1c3c97f58603375a3f2d1883d076a6584190606c9b0", size = 7331806, upload-time = "2026-09-10T01:17:25.641Z" },
+    { url = "https://files.pythonhosted.org/packages/15/dd/80fb0ea15ed15ac1b2b27d2c1d0eef4e725f37effbc4946ddd0795cfc8e3/cuda_bindings-13.4.1-cp314-cp314t-win_arm64.whl", hash = "sha256:76c10f8ad5cf55676518b2cb3d1e3d867826389b29d8418b402d23ba4df7d655", size = 5421917, upload-time = "2026-09-10T01:17:27.82Z" },
+]
+
+[[package]]
+name = "cuda-core"
+version = "1.2.0"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "cuda-pathfinder" },
+    { name = "numpy" },
+]
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/c5/88/145fdb860aa0cafbe197af627567f2f5347d33e61651b38a7361487eb64c/cuda_core-1.2.0-cp313-cp313-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:2a77d40f613e4d151f4466818d1baab271237122c7240450fe66fe8a071ff86f", size = 6154425, upload-time = "2026-09-03T00:12:36.418Z" },
+    { url = "https://files.pythonhosted.org/packages/aa/bf/ccfe842dc59003a721a37103c1964bb797b4b3362d07f82b6bbc71926aec/cuda_core-1.2.0-cp313-cp313-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:0dddba476d990c62a17b58b2c198a3b2768d575987577febda8f3e88140df192", size = 6562184, upload-time = "2026-09-03T00:12:38.441Z" },
+    { url = "https://files.pythonhosted.org/packages/9f/fd/1ebc50d239f61211f8732a7b6e8a35d363e2828b5eeca1930dc3a89925f0/cuda_core-1.2.0-cp313-cp313-win_amd64.whl", hash = "sha256:f19593699a49a52d08300d44c1c84db8629dc01702b36336401c0ba87fe5067f", size = 5923555, upload-time = "2026-09-03T00:12:40.511Z" },
+    { url = "https://files.pythonhosted.org/packages/e0/7c/fd915c7f8cdb954744c1fda75fefa75a39e73508b6aae8ec7369800c8491/cuda_core-1.2.0-cp314-cp314-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:4206c19f3c787c6dd4fdce3625e842686dc73ba665d46d16d6c57659a94db0ac", size = 6233228, upload-time = "2026-09-03T00:12:42.122Z" },
+    { url = "https://files.pythonhosted.org/packages/76/69/2e255b5db595188ba12911fcd9efa45ccf71b02cd96809bc950923dc76d4/cuda_core-1.2.0-cp314-cp314-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:e138be12af795c69c1cb61e38562ea70e5a38616b2753e2a265231f1941394ef", size = 6595561, upload-time = "2026-09-03T00:12:44.001Z" },
+    { url = "https://files.pythonhosted.org/packages/16/c5/7d8a1d90071f5aac2a3bf11f1829535a60519a673809b5c3bc3cdb30df54/cuda_core-1.2.0-cp314-cp314-win_amd64.whl", hash = "sha256:419e884cbe26d85936c8ffea27865a6f2589677f34a81eecb0fe1ec3715c06a0", size = 6051744, upload-time = "2026-09-03T00:12:45.754Z" },
+    { url = "https://files.pythonhosted.org/packages/e7/f3/22b67d9b035fb22d055a39ca21affc4246008da7e72aa024a8fd771b9990/cuda_core-1.2.0-cp314-cp314t-manylinux_2_24_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:c4228d27af1993f8516b51db6c74846ebd1f772bb34132ee721da61369b7febd", size = 6408313, upload-time = "2026-09-03T00:12:48.128Z" },
+    { url = "https://files.pythonhosted.org/packages/a1/82/5edfcbff353cb4dedd13415e5793ca710b302a26f9ee9144d27e9a680bc1/cuda_core-1.2.0-cp314-cp314t-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:ccfd9905e43dc5ef0f88e26a3e3f1509d84c0306710dec351c5ab1db0da5752e", size = 6683167, upload-time = "2026-09-03T00:12:50.178Z" },
+    { url = "https://files.pythonhosted.org/packages/38/eb/30258a7fe759902f1603b960aaedd843666179e686d6ad47d04ba9914a1c/cuda_core-1.2.0-cp314-cp314t-win_amd64.whl", hash = "sha256:29ec78a8a45ad029e8d31b0ee67fe658964efc97a0176eb9d7ad354677116923", size = 6972170, upload-time = "2026-09-03T00:12:52.247Z" },
+]
+
 [[package]]
 name = "cuda-pathfinder"
 version = "1.6.0"
@@ -1345,6 +1387,29 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/dd/7e/2eecb277d8a98184d881fb98a738363fd4f14577a4d2d7f8264266e82623/nvidia_nvjitlink_cu12-12.9.86-py3-none-win_amd64.whl", hash = "sha256:cc6fcec260ca843c10e34c936921a1c426b351753587fdd638e8cff7b16bb9db", size = 35584936, upload-time = "2025-06-05T20:16:08.525Z" },
 ]
 
+[[package]]
+name = "nvmath-python"
+version = "1.0.0"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "cuda-bindings" },
+    { name = "cuda-core" },
+    { name = "cuda-pathfinder" },
+    { name = "numpy" },
+    { name = "pywin32", marker = "sys_platform == 'win32'" },
+]
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/3b/43/1bf0b6d6e97ad3ea5040a026c5aff03f1c60fceb977bb943c2a6bbc25294/nvmath_python-1.0.0-cp313-cp313-manylinux_2_28_aarch64.whl", hash = "sha256:f550e28e87c253ab8747c559c1c4faeb15b79f7c2684bb23e8d4c42df4f5852d", size = 5056945, upload-time = "2026-07-09T20:16:00.903Z" },
+    { url = "https://files.pythonhosted.org/packages/83/f2/7ac589040aeaa92bf72506b480f5a6b9065cb160e10513db07420f0a0eff/nvmath_python-1.0.0-cp313-cp313-manylinux_2_28_x86_64.whl", hash = "sha256:f1e1bc567e56e7ff72515ac34c5af328478e7b8226c6f92b6b2f324925bb38a8", size = 5379834, upload-time = "2026-07-09T20:16:24.375Z" },
+    { url = "https://files.pythonhosted.org/packages/53/52/1cd2fb1fdc2ce3c29341a192ca47209593c7a4e8e1f79f2512c2980ed03f/nvmath_python-1.0.0-cp313-cp313-win_amd64.whl", hash = "sha256:d3a4085a9f98faa223fbb13d6244f8db4d61d48f3f5785e84f6c115b6ee7e8a4", size = 4087396, upload-time = "2026-07-09T20:16:47.375Z" },
+    { url = "https://files.pythonhosted.org/packages/7a/a2/512f5d0a9e2bdfe0b52c7afe7dbdd3c140fa18ab41c1c863f59eee8808df/nvmath_python-1.0.0-cp314-cp314-manylinux_2_28_aarch64.whl", hash = "sha256:7f78fb221df17fbb8e1303d536f5b945122c10febd74b37dea9cb8e44cc586a3", size = 5094362, upload-time = "2026-07-09T20:17:09.604Z" },
+    { url = "https://files.pythonhosted.org/packages/4c/74/a07afd16866ee9f3afcf60e7d5dd1ad3ddc5fb9d56998fce0629eb2671a9/nvmath_python-1.0.0-cp314-cp314-manylinux_2_28_x86_64.whl", hash = "sha256:b6ebb9ae1a2d16dc0b7431d0e6c3c72aa7b11f3e3e83d4dc4f046b086fff8fb8", size = 5395750, upload-time = "2026-07-09T20:17:33.616Z" },
+    { url = "https://files.pythonhosted.org/packages/43/da/d94c5904b0c7878be75994ba541b7937c9ecef50dc633fa10d149885cf4f/nvmath_python-1.0.0-cp314-cp314-win_amd64.whl", hash = "sha256:89807bd2e8a8ef79bf73868638c83fe5f94082f7d8d7b87e08cc94d072959ec4", size = 4204235, upload-time = "2026-07-09T20:17:57.785Z" },
+    { url = "https://files.pythonhosted.org/packages/85/20/f4627d20fef8a2df6dabe8758717d254535529104e6a17c80078276b3ba3/nvmath_python-1.0.0-cp314-cp314t-manylinux_2_28_aarch64.whl", hash = "sha256:e572183b26b850e9e7fd11abbab2b3cc74ab5e13afb04025049e293fee84b4ad", size = 5009618, upload-time = "2026-07-09T20:18:19.401Z" },
+    { url = "https://files.pythonhosted.org/packages/a2/bc/cc88e62b71c37e1b569680507e9eca5f8ea27a8a31a58db3aa6f8250b9ed/nvmath_python-1.0.0-cp314-cp314t-manylinux_2_28_x86_64.whl", hash = "sha256:96fe75363269113d933c9e158da5ce2d585adfab1b7bfb8fdda1a16a202cba52", size = 5292397, upload-time = "2026-07-09T20:18:39.616Z" },
+    { url = "https://files.pythonhosted.org/packages/35/65/958fb26ab46f11b4e8503792f9dc86ef3a18f705206f71181793ca6b3be1/nvmath_python-1.0.0-cp314-cp314t-win_amd64.whl", hash = "sha256:b6ab77f61fca32a77c50075b78f366c7dc05733f46f744068644d63a9530b611", size = 4818972, upload-time = "2026-07-09T20:19:01.313Z" },
+]
+
 [[package]]
 name = "packaging"
 version = "26.2"
@@ -1417,7 +1482,7 @@ wheels = [
 
 [[package]]
 name = "parafac2"
-version = "1.6.0"
+version = "1.7.0"
 source = { registry = "https://pypi.org/simple" }
 dependencies = [
     { name = "anndata" },
@@ -1426,15 +1491,16 @@ dependencies = [
     { name = "tensorly" },
     { name = "tqdm" },
 ]
-sdist = { url = "https://files.pythonhosted.org/packages/5b/11/f69fc560bf17fbdcf269fdb63d9106724ded27dee395e74afb40e061e93f/parafac2-1.6.0.tar.gz", hash = "sha256:b7e7817a28fc59ba1f1d80469c65d74d4a6044bf9f383364eb07cef1f4d3516e", size = 28064, upload-time = "2026-09-02T01:02:14.823Z" }
+sdist = { url = "https://files.pythonhosted.org/packages/40/55/35d771bcd672150a84993ea89764da418d7685df1925c17728e9fc35edfe/parafac2-1.7.0.tar.gz", hash = "sha256:074907a94c64ec0fce4c22cdf73396bb511f7cabd04af353f831302f14e4b1ee", size = 43927, upload-time = "2026-09-14T17:07:34.24Z" }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/06/58/654144ff025927b36895e515b7b01083959afa7ada96909acf49d364a08f/parafac2-1.6.0-py3-none-any.whl", hash = "sha256:dd5b1575bc7c1a779f794ec1c0d2e9c5dc0c32feb8ad6e46e60b0c592b57223a", size = 33469, upload-time = "2026-09-02T01:02:13.391Z" },
+    { url = "https://files.pythonhosted.org/packages/72/8e/634f3f244c9cb7f6b277f1240d09710e4f5fe3c62843242a21239d3c0d74/parafac2-1.7.0-py3-none-any.whl", hash = "sha256:c49ecc4eb4d6e971bce65cbb8ae3f9979ae0450a49b492026bce64633661dd74", size = 53172, upload-time = "2026-09-14T17:07:32.935Z" },
 ]
 
 [package.optional-dependencies]
 gpu = [
     { name = "cupy-cuda12x", extra = ["ctk"], marker = "sys_platform != 'darwin'" },
     { name = "mlx", marker = "platform_machine == 'arm64' and sys_platform == 'darwin'" },
+    { name = "nvmath-python", marker = "sys_platform != 'darwin'" },
 ]
 
 [[package]]
@@ -1828,6 +1894,22 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/0b/d7/1959b9648791274998a9c3526f6d0ec8fd2233e4d4acce81bbae76b44b2a/python_dotenv-1.2.2-py3-none-any.whl", hash = "sha256:1d8214789a24de455a8b8bd8ae6fe3c6b69a5e3d64aa8a8e5d68e694bbcb285a", size = 22101, upload-time = "2026-03-01T16:00:25.09Z" },
 ]
 
+[[package]]
+name = "pywin32"
+version = "312"
+source = { registry = "https://pypi.org/simple" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/2d/41/12fbfd7f36ed2146d8bc9de96c2741296bf0d490b98508496cff322e274c/pywin32-312-cp313-cp313-win32.whl", hash = "sha256:7a27df850933d16a8eabfbaeb73d52b273e2da667f80d70b01a89d1f6828d02c", size = 6370184, upload-time = "2026-06-04T07:49:36.253Z" },
+    { url = "https://files.pythonhosted.org/packages/ba/db/36a78e3403099d31d9746d13fdcde5accc43c1155f375a34d15983a479a7/pywin32-312-cp313-cp313-win_amd64.whl", hash = "sha256:c53e878d15a1c44788082bfe712a905433473aa38f86375b7cf8b45e3acbaaf9", size = 6914298, upload-time = "2026-06-04T07:49:38.876Z" },
+    { url = "https://files.pythonhosted.org/packages/84/37/c1697194092b76de9ed47ca124323f02c57ffc8a45c06f88a3d5acaf01eb/pywin32-312-cp313-cp313-win_arm64.whl", hash = "sha256:59aba5d5940842075343a5ddc6b11f1cdf0d1567fe745290359dfbcc7c2eb831", size = 6727640, upload-time = "2026-06-04T07:49:41.083Z" },
+    { url = "https://files.pythonhosted.org/packages/fc/2b/1f3cded5822fd49c02f40544cbb5f58c7cfd6b1694869fd476cb6170ee97/pywin32-312-cp314-cp314-win32.whl", hash = "sha256:a77a90fbb6881238d2ca9c6fd797b25817f3768fe78d214a90137ff055a75f5b", size = 6468928, upload-time = "2026-06-04T07:49:43.188Z" },
+    { url = "https://files.pythonhosted.org/packages/21/82/3bf86d2e2808902013132e1ce905a7da0da53790f3836c64bf44d55e24f3/pywin32-312-cp314-cp314-win_amd64.whl", hash = "sha256:a4dd3a848290ef724347b19f301045831d8e802fa4464f491b98b1e0a081432e", size = 7024157, upload-time = "2026-06-04T07:49:45.34Z" },
+    { url = "https://files.pythonhosted.org/packages/a4/0e/73f6d6800b4f27655abd9e9f6aaeaefcddb2b946e4674efa2bab184a7f7b/pywin32-312-cp314-cp314-win_arm64.whl", hash = "sha256:9fce94568364e0155e6dfb781ac5d95903be8baf28670632beab1b523f300daa", size = 6839598, upload-time = "2026-06-04T07:49:47.613Z" },
+    { url = "https://files.pythonhosted.org/packages/eb/61/caa39686032d2ebdd04ff0ab5cbe163126c0066d98e00c9018646e42393b/pywin32-312-cp315-cp315-win32.whl", hash = "sha256:5c1fbe4a937a73ae9297384a3da38518cbc694c68ad8a809b2e19acd350f03ed", size = 6471159, upload-time = "2026-06-04T07:49:50.035Z" },
+    { url = "https://files.pythonhosted.org/packages/0f/cd/7e1de64a4a6f69c04214169657ccab0d93a670ea50e35eb8f489d7378249/pywin32-312-cp315-cp315-win_amd64.whl", hash = "sha256:c2f03a0f73f804a13c2735b99392b0cd426bb4f2c4d0178e5ac966a0f21618d5", size = 7025293, upload-time = "2026-06-04T07:49:54.857Z" },
+    { url = "https://files.pythonhosted.org/packages/23/ed/4532e9388e65fa16b46776ef47ad631a64eda1631884488af707666350ed/pywin32-312-cp315-cp315-win_arm64.whl", hash = "sha256:a8597d28f267b39074aef51fa593530082b39cbe5a074226096857b1fed2dfb9", size = 6840337, upload-time = "2026-06-04T07:49:57.531Z" },
+]
+
 [[package]]
 name = "pyyaml"
 version = "6.0.3"
@@ -2168,7 +2250,7 @@ requires-dist = [
     { name = "numpy", specifier = ">=2.2" },
     { name = "pacmap", specifier = ">=0.9" },
     { name = "pandas", specifier = ">=3.0.0" },
-    { name = "parafac2", specifier = ">=1.6.0" },
+    { name = "parafac2", specifier = ">=1.7.0" },
     { name = "parafac2", extras = ["gpu"], marker = "extra == 'gpu'" },
     { name = "pyarrow", specifier = ">=19.0" },
     { name = "scikit-learn", specifier = ">=1.6" },