pp.regress_out returns wrong residuals when a covariate barely varies

[selmanozleyen]

Hi, here I reproduced some numeric errors. But my main question is do you prefer this to be kept for reproducibility with older versions or fixed? I am asking because I am working on numbafication of this function, and I ran into this.

The vectorized path in numpy_regress_out solves the normal equations with an explicit inverse:

inv_gram_matrix = np.linalg.inv(regressor.T @ regressor)
coeff = inv_gram_matrix @ (regressor.T @ data)
data = get_resid(data, regressor, coeff)

Forming XᵗX squares the condition number of the design. A nearly-constant covariate is nearly collinear with the intercept column, so cond(X) grows and float64 runs out of digits — while the det(regressors.T @ regressors) != 0 guard still reports the design as full rank and routes it here.

regress_out only ever returns residuals, never coefficients. The residual is well-conditioned even when the coefficients are not, so this precision loss is introduced by the algorithm rather than being inherent to the problem — pinv on the same designs is accurate to ~1e-9.

Result

On main (3dadf246), scanpy 1.14.0.dev19+g3dadf246c, numpy 2.4.6:

pct_mito column	cond(X)	det != 0?	max err	% of resid SD	== inv?	err via pinv
healthy spread (0–15%)	2.5e+01	True	2.77e-14	0.00%	True	2.33e-15
tight QC (2.0% ± 1e-4)	1.8e+05	True	1.11e-06	0.00%	True	3.01e-12
near-constant (2.0 ± 1e-7)	1.7e+08	True	3.85e-01	133.42%	True	2.95e-09
identical to 1e-9	7.2e+09	True	3.28e-01	113.84%	True	2.04e-07

• % of resid SD — the error as a percentage of the returned residuals' own standard deviation. In the bottom two rows the output is off by more than the entire spread of the values it is reporting. Healthy inputs are unaffected (0.00%), which is why no test catches this.
• == inv? is True throughout, confirming scanpy's output is exactly the normal-equations result — these are the current implementation's own numbers, not a modelled approximation.
• err via pinv is what the same designs give through np.linalg.pinv: ~1e-9 and ~1e-7. The information is present in the data; the current method discards it.
• det != 0 is True in every row, including the worst.

Reference residuals come from modified Gram–Schmidt in np.longdouble (~18 significant digits vs float64's ~16), which never forms XᵗX and so never squares the condition number.

Reproducer

Click to expand

"""`sc.pp.regress_out` returns wrong residuals when a covariate barely varies."""

from importlib.metadata import version

import anndata as ad
import numpy as np
import pandas as pd

import scanpy as sc


def reference_residuals(y, X):
    """OLS residuals via modified Gram-Schmidt in longdouble (~18 digits vs float64's ~16).

    Orthonormalize the design's columns, then project: no matrix is inverted, so the
    condition number is never squared. Columns that are (nearly) combinations of earlier
    ones are dropped at a tolerance relative to their own norm -- the same decision
    `pinv` makes via its singular-value cutoff.
    """
    basis = []
    for j in range(X.shape[1]):
        col = X[:, j].astype(np.longdouble)
        v = col.copy()
        for q in basis:
            v -= (q @ v) * q
        norm = np.sqrt(v @ v)
        if norm > 1e-13 * np.sqrt(col @ col):
            basis.append(v / norm)
    r = y.astype(np.longdouble).copy()
    for q in basis:
        r -= (q @ r) * q
    return np.asarray(r, dtype=np.float64)


rng = np.random.default_rng(0)
n_obs, n_vars = 2000, 50

x = rng.random((n_obs, n_vars))
n_counts = rng.lognormal(8, 0.4, n_obs)
n_counts /= n_counts.max()

covariates = {
    "healthy spread (0-15%)": rng.random(n_obs) * 0.15,
    "tight QC (2.0% +/- 1e-4)": 2.0 + 1e-4 * rng.random(n_obs),
    "near-constant (2.0 +/- 1e-7)": 2.0 + 1e-7 * rng.random(n_obs),
    "identical to 1e-9": 1.0 + 1e-9 * rng.random(n_obs),
}

print(f"scanpy {version('scanpy')}, numpy {np.__version__}\n")
head = (
    f"{'pct_mito column':<30}{'cond(X)':>9}{'det!=0':>8}"
    f"{'err':>10}{'% resid SD':>12}{'==inv?':>8}{'err pinv':>10}"
)
print(head)
print("-" * len(head))

for label, pct_mito in covariates.items():
    adata = ad.AnnData(
        x.copy(),
        obs=pd.DataFrame(
            {"n_counts": n_counts, "pct_mito": pct_mito},
            index=[f"cell_{i}" for i in range(n_obs)],
        ),
    )
    sc.pp.regress_out(adata, keys=["n_counts", "pct_mito"])
    got = np.asarray(adata.X)

    # the design regress_out builds internally: intercept, then one column per key
    design = np.c_[np.ones(n_obs), n_counts, pct_mito]
    expected = np.column_stack([
        reference_residuals(x[:, j], design) for j in range(n_vars)
    ])

    # what the current normal-equations path computes, and what `pinv` would give
    via_inv = x - design @ (np.linalg.inv(design.T @ design) @ (design.T @ x))
    via_pinv = x - design @ (np.linalg.pinv(design) @ x)

    err = np.abs(got - expected).max()
    print(
        f"{label:<30}{np.linalg.cond(design):>9.1e}"
        f"{str(np.linalg.det(design.T @ design) != 0):>8}"
        f"{err:>10.2e}{100 * err / expected.std():>11.2f}%"
        # confirms scanpy really took the inv() path, so these are its own numbers
        f"{str(np.allclose(got, via_inv, atol=1e-12)):>8}"
        f"{np.abs(via_pinv - expected).max():>10.2e}"
    )

Output:

scanpy 1.14.0.dev19+g3dadf246c, numpy 2.4.6

pct_mito column                 cond(X)  det!=0       err  % resid SD  ==inv?  err pinv
---------------------------------------------------------------------------------------
healthy spread (0-15%)          2.5e+01    True  2.77e-14       0.00%    True  2.33e-15
tight QC (2.0% +/- 1e-4)        1.8e+05    True  1.11e-06       0.00%    True  3.01e-12
near-constant (2.0 +/- 1e-7)    1.7e+08    True  3.85e-01     133.42%    True  2.95e-09
identical to 1e-9               7.2e+09    True  3.28e-01     113.84%    True  2.04e-07

On the det guard

det is not a conditioning measure. Worse: by cond ≈ 1.7e8 the computed det(XᵗX) has gone negative — mathematically impossible for a positive semi-definite matrix — so the quantity being compared against zero is pure rounding noise, and != 0 still says "full rank, proceed."

Swapping in np.linalg.matrix_rank fixes exact-singularity detection but not this: its default tolerance (max(M,N)·eps ≈ 8.8e-14) still calls these designs full rank. The fix is to remove the need for the guard by routing everything through pinv, which handles full-rank, ill-conditioned and singular designs by the same stable path.

Cost of fixing

Well-conditioned results move by ~1e-12, below every tolerance in the test suite except test_regress_out_reproducible's ordinal case, which is pinned at atol=0.0 — so regress_test_small.npy would need regenerating. Worth noting that atol=0.0 currently pins a computation method rather than a result.