Repository navigation
Glm solver rework - #55
Merged
Merged
Conversation
Three defects in NBGLMBatchFitter, all found while benchmarking wide designs. min_mu was applied as a clamp on the IRLS weights, on the variance, and on several division guards, as well as the floor on the fitted mean that DESeq2 defines. Clamping the weights solves a different estimating equation: on genes below one count per cell the returned coefficients left a score of 50-100 where the maximum likelihood estimate has a score of zero, and 40-100 units of excess deviance. They now leave a score of ~1e-6. min_mu floors the fitted mean and nothing else; its default is 0 on both fitters, and de.py passes 0.5 explicitly so the DESeq2-compatible path is unchanged. The same conflation floored NBGLMFitter's standard errors at sqrt(min_mu) -- forcing se >= 0.707 at the old default -- and the Cook's distance leverage denominator. The per-gene Hessian was formed with a three-operand numpy.einsum, which stops routing through BLAS once its intermediate exceeds an internal budget and falls back to a nested loop. Formed with a chunked gemm, fit_batch on 1,000 cells and 500 genes went from 16.05s to 0.15s at design width 41 and from 35.38s to 0.23s at width 61, agreeing to 7e-12; below width 21 there is no material difference. The IRLS loop accepted Newton steps that increased a gene's deviance, iterated every gene until the slowest converged, re-estimated the dispersion inside each iteration, and solved the normal equations in the caller's column basis. It now halves such steps, freezes converged genes, estimates the dispersion around the loop with a consistency refit, and solves in a unit-RMS basis. Convergence requires both the relative deviance change and the largest coefficient change to fall below tol: deviance alone is quadratic in the coefficients near the optimum, so it only pins them to about sqrt(tol). Adds a real Poisson family and fixed_dispersion on fit_batch, and StructuredGLMBatchFitter for [covariates | one-hot groups] designs, whose disjoint group supports make the per-gene Hessian arrowhead-structured. Each Newton step goes through the Schur complement of the diagonal block: exact against a dense per-gene solve to 1.4e-15, and 14.2x faster at 31 covariates and 200 groups. It also converges where the dense path does not (744 of 800 genes against 31 of 800 on that design), because an empty group's coefficient is unbounded and the group ridge and clip make it finite and defined. Shared numerics live in the new crispyx._irls so the two fitters cannot drift. Correctness on sparse genes is checked against the score equation rather than another implementation: cross-solver agreement is a weak instrument there, since a disagreement says the two differ without saying which is right. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…stop reporting non-estimable effects build_design_matrix one-hot encodes categorical covariates, so nb_glm_test with a batch, donor or lane covariate produces a [covariates | one-hot groups] design -- the shape StructuredGLMBatchFitter exists for. fit_nb_glm_batch_auto takes that path when the group columns outnumber the remaining covariates and the dense path otherwise, since the advantage grows with that ratio and was measured at only 1.5-2.5x for 20 groups beside 12 covariates. Results come back in the column order of the design either way. The intercept and the perturbation column are held out of the group block, because clip_group bounds group coefficients, which suits a nuisance level effect and not the effect under test. End-to-end log-fold-changes agree with the dense path to 1e-3 and the statistic correlation is above 0.999. A (gene, perturbation) pair with no counts in one arm no longer receives an effect estimate. A log-link GLM estimates the effect as a difference of log means, so such a pair has no finite effect: the likelihood has no interior maximum and the coefficient runs to the boundary. What was reported came from wherever the fitter stopped -- on one gene the effect was -18.9 at min_mu=0 and -5.1 at min_mu=0.5, with the Wald statistic moving from 0.08 to 32.6 on identical counts, that is from no evidence to wildly significant purely as an artefact of the mean floor. (The collapse at min_mu=0 is the Hauck-Donner effect: as the coefficient diverges its standard error grows faster still.) Those pairs are now reported as untested, as genes with no counts anywhere already were, and are excluded from the multiple-testing correction rather than entering it with artefactual p-values. They are deliberately not reported as an effect of zero, which would describe a completely silenced gene as unchanged; the observation stays visible in pts and pts_rest, and apeglm shrinkage still gives a bounded estimate where one is wanted. min_cells_ctrl and min_cells_pert are separate because the informative direction depends on the screen: a knockdown screen expects hits abundant in control and depleted in the perturbed arm, an activation screen the reverse. Both default to 1, which is symmetric and also the exact boundary between an effect that exists and one that does not. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…0.1.5 Estimates change for lowly-expressed genes -- toward the reference implementations, not away -- so results from 0.1.4 on covariate-adjusted or low-count genes will not reproduce exactly. nb_glm_test keeps its own min_mu default of 0.5, so the DESeq2-compatible path is unchanged. Renames the API section from "Negative binomial GLM", which no longer describes it: the module now fits Poisson as well, and carries a second fitter for designs with a one-hot block. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…rsion Four defects in the 0.1.5 solver, all of which let a fit be reported as something it is not. Step halving aliased its own interpolation endpoint. The accepted half-step was written back into the array that the next retry interpolated against, so the distance shrank geometrically-squared rather than by halves: the third retry landed at 2^-6 of the Newton step rather than 2^-3. A gene needing repeated damping therefore stopped moving, and was then reported as converged *because* nothing had moved, at a point that is not a stationary point of the deviance -- which de.py gates its output on. Each retry now interpolates towards the full Newton point, and a shortened step is accepted only if it improves on the deviance the iteration started from, not merely on the failed full step. A gene for which nothing in the halving range improves is left where it was rather than moved uphill, and is judged by the steps it refused instead of by how far it moved. Cells on the fitted-mean floor are exempt from that judgement: their mean does not move with the coefficients so they carry no gradient, yet the normal equations still count them, and the line search then describes the floor rather than the fit. A gene at the eta divergence clips gets no such exemption. The Numba path fitted at a hard-coded dispersion of 0.1 and reported a separately estimated one beside it. On an intercept-plus-perturbation design the coefficients are the two group means whatever the dispersion is, but the standard errors are not: on an overdispersed fixture they came out up to 3.9x too small, so every Wald statistic built on them was up to 3.9x too large. The kernel is now run twice, as the NumPy path already was. The NB deviance was formed as the difference of (y + r) log(y + r) and (y + r) log(mu + r), each of order r n log r. At the alpha clip floor r is 1e8, so a near-Poisson gene's deviance was the small difference of two very large sums: measured against a 60-digit reference, the absolute error was 0.109 at 3,000 cells and 182 at 100,000, far above the 1e-6 relative tolerance the convergence test uses. Such a gene could never converge and its damping fired on noise. Evaluated with log1p the same quantities are accurate to 5e-10 and 7e-7. The ridge penalty was added in the preconditioned basis, which penalises the caller's coefficient j by ridge * scale_j**2. It is now divided by the column scale. No effect at the 1e-6 default. Finally, StructuredGLMBatchFitter._irls set eta beside the coefficients rather than from them, so a design without an exact column of ones started from a point its coefficients did not describe and compared damped steps against a deviance belonging to a different model. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…rded work StructuredGLMBatchFitter.fit_batch had no gene batching, unlike NBGLMBatchFitter: it held about a dozen (n_samples, n_genes) work arrays live per Newton step, which at 3,000 cells and 8,563 genes is around 2 GB, against the 100 MB of work arrays the dense fitter targets. It now takes the same gene_batch_size argument, defaulting to "auto". Measured at 1,500 cells and 1,500 genes, peak allocation falls from 239 MB to 116 MB at the same runtime; the fit is unchanged beyond the summation order BLAS picks for a different matrix shape. The negative-binomial path also ran a full Poisson fit and an intermediate NB fit only for their fitted means, and computed -- then discarded -- the per-gene standard errors of both, each of which inverts a d_X matrix per gene. Neither intermediate fit now forms them. The Poisson fitter is built once rather than per batch. schur_solve re-formed the two (n_genes, n_features, n_groups) arrays that schur_complement had just built, holding three where one would do, and duplicating a gemm, on every Newton step. Both now come from one helper. Both IRLS loops copied the counts of the active gene set twice per iteration -- once for the working response, once for the deviance -- including on the first iteration, where the active set is every gene and the copies are of arrays already in hand. The full-batch case now uses the originals and the deviance subset reuses the slice the loop already made. Also removes NBGLMBatchFitter._xtx, which nothing has read since the preconditioned solve replaced it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A (gene, perturbation) pair with no counts in one arm is reported as untested, on the grounds -- stated in the documentation -- that the observation itself survives in pts and pts_rest, where a complete knockdown expressed in 95% of control cells and none of the perturbed ones is plainly visible. It did not: the uncached path built both behind the validity mask, so the evidence was zeroed away on precisely the pairs it was there for, and a silenced gene became indistinguishable from one that was never measured. The cached path did not do this, so the two disagreed. Fractions of expressing cells describe the data rather than a fit, so both paths now report them for every gene. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ty rationale Folds the fixes into the 0.1.5 entry rather than opening a new version, since 0.1.5 has not shipped from this branch, and drops its release date line. The argument for reporting a pair with an empty arm as untested -- the measurement, the Hauck-Donner note, and which asymmetry suits CRISPRi against CRISPRa -- was written out four times, in _nonestimable_glm_mask, in its own docstring below that, in nb_glm_test's parameter list and in the changelog, which is four copies to keep in step the next time the rule changes. It now lives with the function that implements it, and the other three point at it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
``min_mu`` floors the fitted mean to steady the IRLS iteration. It is not part of the model whose uncertainty is being reported, and DESeq2 keeps it out of the Wald standard error: irls_solver deliberately "returns an UNthresholded mu (as in the R code)" and wald_test rebuilds W = mu / (1 + mu * disp) from exp(X @ LFC) * size_factors, with no floor anywhere. crispyx built those weights from the floored mean instead, which overstates how much every floored cell knows. Measured against PyDESeq2's own wald_test on identical coefficients, the standard errors on floor-binding genes were a median of 23.7% and at worst 34.9% too small; genes the floor never touches already agreed to 0.000%. Rebuilt from the unfloored mean they agree to 0.00%. Isolated on a real covariate-adjusted design (Adamson, 9,496 genes) with the coefficients and dispersions held fixed, the floored weights give standard errors a median of 5% and a 10th percentile of 14% too small. The weights are still clipped, at the generic ETA_MIN/ETA_MAX divergence bounds rather than at the mean floor, so a diverging gene stays finite. Two things deliberately keep the floored weights. The weights inside the IRLS iteration are the DESeq2 iteration itself; changing those changes the coefficients. So is the leverage behind Cook's distance, which DESeq2 computes with the floored W inside irls_solver. The floor belongs to the iteration, and that is the whole of where it now appears. The Numba kernel no longer computes standard errors at all -- it built them from its floored weights, and the caller now forms them from the unfloored mean like every other path, so ~25 lines of duplicated 2-feature SE algebra come out of the kernel with it. Not every reported standard error was affected: nb_glm_test on a two-group comparison without covariates discards the fitter's and recomputes with _statistics._compute_se_batched, which never floored. The routes that do use it are the covariate-adjusted one (fit_nb_glm_batch_auto), the structured solver, and the public fitter API. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… binds A cell at the fitted-mean floor has d mu / d beta = 0: it contributes nothing to the gradient of the deviance being evaluated while still contributing its full weight to X'WX and its full residual to X'Wz. The Newton step therefore solves one model while the deviance measures another, and can point uphill at a perfectly good fit. No step length repairs that, so the deviance-guarded line search added in 0.1.5 was the wrong instrument for those genes, and 854aaee's exemption -- passing them through the convergence test unchecked -- was the wrong repair: it reported a gene as converged without having tested it. This is not an edge case. On a screen-shaped fixture the floor binds on 73 of 120 genes, with a median of 291 of 300 cells floored on those. DESeq2 has no line search at all: its IRLS is a fixed-point iteration stopped on |dev - dev_old| / (|dev| + 0.1) < tol, and crispyx's own Numba kernel has always been one too. Only the two NumPy loops judged a non-descent iteration by a descent test. A gene resting on the floor now takes the full Newton step and is judged by that ratio, paired with the same coefficient-change test the other branch uses -- on its own the ratio fires on the first iteration when a warm start happens to leave the deviance flat while the coefficients are still moving. The split is between steps and tests, not between genes. A gene with some floored cells that is descending normally keeps the strict test; only once the full step has failed to descend does the floor decide whether that is a stall (revert, report not converged) or the floor holding the gene where it is (keep the step, use the ratio). Judging by the current iterate alone stopped well-expressed genes early, at a score of 1220 where they otherwise reach 1e-6. Against PyDESeq2 on five fixtures, coefficients on floor-binding genes go from 2.8e-02 away to 6.4e-08..3.3e-07; every gene converges where one to four used to fall short; the worst gene takes 23 iterations against a 200 cap it used to hit; and min_mu=0.5 is no longer slower than min_mu=0. On real data (Adamson, 9,496 genes, covariate-adjusted) this returns 103 and 89 genes that were being dropped, at 3.4 s against 5.7 s. They are not marginal: their median pts_rest is 0.69 and 0.77 against 0.30 for the genes already reported, and 44 of them reach padj < 0.05. Genes reported before are untouched -- median |dLFC| 2e-05. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Both fold into the 0.1.5 entry, which has not shipped from this branch. The convergence change is user-visible in two ways worth stating plainly: low-count genes that were reported as NaN now carry estimates, which enlarges the BH denominator, and p-values on floor-binding genes grow because their standard errors were too small. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The entry quoted the isolated figures -- standard errors a median of 24% too small -- next to a claim that p-values grow, which reads as though hit lists move by that much. They do not: nb_glm_test re-estimates the dispersion after the fit and absorbs most of it. Measured on the Adamson run, p-values for genes reported both before and after move by a median factor of 1.0000, with 7 losing significance and 19 gaining it. Also says plainly that a non-converged gene previously had no result at all, not merely no log-fold-change, so the convergence fix changes what is tested rather than only what is estimated. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
nb_glm_test already wrote NaN to every column of a (gene, perturbation) pair it excluded. t_test and wilcoxon_test did not, and disagreed with each other and with themselves: depending on which of the three Wilcoxon paths ran, a gene dropped by min_cells_expressed came back as p = 1.0 beside a fold change of exactly 0.0 -- "tested, no change", asserted about a gene nobody tested -- or as a NaN p-value beside that same 0.0, while t_test returned a NaN p-value beside a finite fold change and effect size. Worse, the 1.0 leaked across rows. The kernels leave their own defaults in the columns they skip and the whole block was written out, so a gene excluded for one perturbation but tested for another picked up that other perturbation's p = 1.0: on the Adamson subset with min_cells_expressed=20, 48 genes in the AARS row carried a p-value only because AMIGO3 happened to be in the same run. A row is no longer a function of the run's composition. All four paths now mask every column derived from the comparison by the validity mask, and pts / pts_rest -- descriptions of the data, not of the comparison -- stay populated, as they already did in t_test and nb_glm_test. Zeroing them removed the only evidence that an excluded gene was a 0%-vs-95% knockdown. A NaN p-value does not enter the Benjamini-Hochberg denominator, so adjusted p-values on the genes that were tested get slightly smaller. That is the correct denominator. One exception, deliberate: in the batch-stratified path a perturbation sharing no batch with the control keeps its pooled fold change beside a NaN p-value, because only the rank test is stratified and the pooled columns follow the pooled filter. Also fixes t_test folding untested cells into the last perturbation. With an explicit perturbations= subset, cells outside it code to -1 and -1 indexes the last row, so testing one perturbation out of three compared control against that perturbation plus both untested ones -- its sums, means, variances, expressing counts and cell count. wilcoxon_test and nb_glm_test select cells by label and were unaffected. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
logfoldchanges reproduces scanpy's formula, so a gene with no expressing cells in one arm has the 1e-9 constant as that arm's term: the magnitude reported says how small the constant is, not how large the effect is, with a ceiling of log2(1 / 1e-9) = 29.9. Splitting 500 real control cells at random into two arms of 250 -- where every difference is noise -- puts roughly 119 of 11,630 genes above |log2FC| = 2, the largest at 24.9. The p-values are sound: 4.1-4.9% rejection against a nominal 5% across sparsity strata, and the smallest padj over five null runs is 0.73. So this is a caveat about ranking on the effect-size column, not about the test, and the FAQ says so rather than leaving a user to discover it from a hit list. pts and pts_rest separate a real complete knockdown (pts = 0 beside a high pts_rest) from a gene that was never detectable (both near zero), and the FAQ shows the mask that follows from it: on a real Adamson perturbation it takes the largest reported |log2FC| from 23.6 to 3.0 and brings the targeted gene itself into the top five, at the cost of setting aside 2,612 of 11,630 genes as not measurable in both arms. That mask excludes complete knockdowns by construction, so the text says to look for them in pts rather than at the top of the ranking. comparison.rst now states that parity with scanpy includes this convention, rather than leaving r > 0.9999 to imply the column is safe to sort on. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The stored outputs came from 0.1.3. A stale src/crispyx.egg-info left by an older editable install still advertised 0.1.4 and shadowed the installed metadata -- the notebooks put src/ first on sys.path, so importlib.metadata resolved the version there rather than in site-packages -- so it is removed and the version banner now reads 0.1.5. Results are unchanged. The tutorials run with the default filters, where no gene is excluded, so the reporting fix is invisible in them; the rest of the diff is timestamps, execution counts, figures and two progress lines that earlier versions did not print. Both notebooks execute without error and the HTML build is clean apart from one pre-existing docstring warning in glm.py. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
An excluded (gene, perturbation) pair is dropped before any fit, so it has no MLE to shrink. Two places handed it one anyway. The per-gene shrinkage returns a fallback standard error whenever it does not fit the gene, and that fallback was the literal 1.0 whenever the MLE standard error was not finite -- which is exactly the untested case. The NaN the rest of the result reports for such a gene was overwritten with a number, and lfc_base="log2" scaled it, so an untested gene surfaced with se = 1.44 beside a NaN p-value. The fallback is now the MLE standard error as it stands, NaN included. The cached worker then built mle_coef from the whole coefficient block. That block is zero-initialised, so a gene outside the validity mask arrived with a finite [0, 0] start; on the global-dispersion fast path its dispersion is finite too, so nothing short-circuited and the gene was fully re-fitted under the prior. On the separation fixture the silenced gene came back with logfc = -14.07 beside a NaN p-value. The uncached worker has always masked mle_coef, so the two workers disagreed on identical input. Both now mask it. The docstring and the changelog pointed users at lfc_shrinkage_type="apeglm" for a bounded estimate on non-estimable pairs. That was never true on the path that masks -- the pair is excluded before the fit, so shrinkage has nothing to bound -- and with the fix above it is not true anywhere. Both now say so. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Both batch fitters judge a gene resting on the mean floor by DESeq2's deviance ratio paired with the coefficient-change test, since the line search says nothing about a fit the floor has flattened. The ratio was inlined in the np.select that follows the in-place writes, and the deviance it compares against aliases the array those writes land in whenever the whole batch is active -- which is every iteration until the first gene converges or stalls. So it compared the new deviance with itself: identically zero, always below tol, and the paired test silently degenerated to the coefficient-change term alone. Both fitters now form deviance_ratio beside relative_change, before the writes, and the comment at each site says why it has to live there. The term is live again: tracing 400 low-count genes at min_mu=0.5 and tol=1e-8, 502 of 606 floor-resting gene-iterations now carry a ratio at or above tol, up to 1.9e-3. No verdict changes on the fixtures -- across 8 seeds, three iteration budgets and two tolerances, not one convergence flag, iteration count or coefficient differs, because the coefficient term dominates in every case measured. This is a latent defect fixed by inspection, not an observed misfit: what it prevents is a gene whose Newton step raises the deviance while its coefficients move less than tol being called converged at a non-stationary point, and nb_glm_test gates its p-values on that flag. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
OneHotGroups is where the arrowhead structure is asserted: disjoint one-hot columns make G'WG diagonal, which is what schur_solve exploits and what makes the structured solve exact rather than approximate. The dense branch raised on a cell belonging to two groups. The sparse branch assigned labels in ascending column order, so the same input silently kept the higher-indexed group and the exactness guarantee quietly stopped holding. It also read the index array directly, so an explicitly stored zero counted as membership. Both are now handled: stored zeros are dropped before the supports are compared, a shared cell raises the same ValueError the dense branch raises, and the labels are assigned in one vectorised write instead of a loop over columns. Reachable through fit_glm_onehot and StructuredGLMBatchFitter with a scipy.sparse groups matrix; the internal callers pass dense or already-validated input and were unaffected. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
![claude[bot]](https://avatars.githubusercontent.com/in/1236702?s=60&v=4){kind=small}
There was a problem hiding this comment.
If your organization's extra usage balance is empty, an organization admin can add extra usage credits at claude.ai/admin-settings/usage. If its monthly spend limit was reached, an admin can raise it on the same page. If neither applies, contact Anthropic support.
Once extra usage is available, comment @claude review on this pull request, or push a new commit, to trigger a review.
The tutorials resolved their paths to absolute ones, so every path crispyx echoed back -- "Saving -> ...", the backed AnnData repr, the disk-space report -- committed the author's home directory into outputs that ship to readers verbatim (nbsphinx_execute = "never"). Address data and outputs relative to the notebook directory instead, and capture the aggressive-filtering warning rather than letting Python's default banner print an absolute path to qc.py. Both notebooks re-executed. DiskEstimate now labels itself with the path it was given instead of a resolved copy; the free-space walk-up still resolves internally, so only the message changes. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.