Derivatives of the DeCovarT generative model under simplex transforms

DeCovarT maximises a Gaussian-convolution log-likelihood for cellular ratios p ∈ ΔJ − 1. This vignette derives the analytic first- and second-order derivatives used by the frequentist solvers, first with respect to unconstrained ambient proportions p, then after the isometric log-ratio (ILR) map that enforces the simplex. It closes with the expected Fisher information and the ILR delta-method covariance that feed vcov.decovart_fit() and confint.decovart_fit(). Package entry points are gradient_loglik_unconstrained(), hessian_loglik_unconstrained(), gradient_loglik_constrained(), hessian_loglik_constrained(), expected_fisher_unconstrained() and vcov_ilr_delta(). The additive log-ratio chart is retained as an appendix (Sec. 8).

Note 1: Notation

  • μ ∈ ℝG × J: mean signature.
  • Σj ∈ SPD(G): cell-type covariance; Σ(p) = ∑jpj2Σj and Ω(p) = Σ(p)−1.
  • p: proportions on the open simplex (ambient score equations treat them as free coordinates in (0, 1)J).
  • ρ ∈ ℝJ − 1: unconstrained log-ratio coordinates (ILR / Helmert chart V); p = softmax (Vρ). The manuscript DAG uses the same letter for the ALR chart; here ρ denotes the ILR coordinates used by the solvers.

Matrix calculus reminders

The score and Hessian of the convolution log-likelihood follow from a handful of matrix-derivative identities. The statements below mirror the matrix-calculus remarks in the manuscript Methods appendix.

Remark 1 (First-order matrix calculus). Let A = A(p) and B = B(p) be invertible matrix-valued functions of a scalar p, and let U, V be constant matrices of compatible size. Then

$$ \begin{aligned} \frac{\partial\det(\boldsymbol{A})}{\partial p} &= \det(\boldsymbol{A})\, \mathrm{tr}\!\left( \boldsymbol{A}^{-1} \frac{\partial\boldsymbol{A}}{\partial p} \right), \\ \frac{\partial(\boldsymbol{U}\boldsymbol{A}\boldsymbol{V})}{\partial p} &= \boldsymbol{U} \frac{\partial\boldsymbol{A}}{\partial p} \boldsymbol{V}, \\ \frac{\partial\boldsymbol{A}^{-1}}{\partial p} &= - \boldsymbol{A}^{-1} \frac{\partial\boldsymbol{A}}{\partial p} \boldsymbol{A}^{-1}. \end{aligned} \qquad(1)$$

The chain rule on log det  yields

$$ \frac{\partial\log\det(\boldsymbol{A})}{\partial p} = \mathrm{tr}\!\left( \boldsymbol{A}^{-1} \frac{\partial\boldsymbol{A}}{\partial p} \right), \qquad \frac{\partial\log\det(\boldsymbol{A}^{-1})}{\partial p} = - \mathrm{tr}\!\left( \boldsymbol{A}^{-1} \frac{\partial\boldsymbol{A}}{\partial p} \right). \qquad(2)$$

For a residual quadratic form with symmetric precision Ω and mean map linear in p,

$$ \frac{\partial}{\partial p} \bigl( (\boldsymbol{y}-\boldsymbol{x}p)^{\top} \boldsymbol{\Omega} (\boldsymbol{y}-\boldsymbol{x}p) \bigr) = -2\, \boldsymbol{x}^{\top} \boldsymbol{\Omega} (\boldsymbol{y}-\boldsymbol{x}p). \qquad(3)$$

Remark 2 (Second-order matrix calculus). Differentiating the inverse once more gives

$$ \frac{\partial^{2}\boldsymbol{A}^{-1}}{\partial p_{i}\partial p_{j}} = \boldsymbol{A}^{-1} \Biggl( \frac{\partial\boldsymbol{A}}{\partial p_{i}} \boldsymbol{A}^{-1} \frac{\partial\boldsymbol{A}}{\partial p_{j}} - \frac{\partial^{2}\boldsymbol{A}}{\partial p_{i}\partial p_{j}} + \frac{\partial\boldsymbol{A}}{\partial p_{j}} \boldsymbol{A}^{-1} \frac{\partial\boldsymbol{A}}{\partial p_{i}} \Biggr) \boldsymbol{A}^{-1}, \qquad(4)$$

and the trace is linear, ∂tr(A)/∂pi = tr(∂A/∂pi). Combining Eq. 1 with that linearity yields

$$ \frac{\partial^{2}\log\det(\boldsymbol{A}^{-1})}{\partial p_{i}\partial p_{j}} = - \mathrm{tr}\!\left[ \boldsymbol{A}^{-1} \frac{\partial^{2}\boldsymbol{A}}{\partial p_{i}\partial p_{j}} \right] + \mathrm{tr}\!\left[ \left( \boldsymbol{A}^{-1} \frac{\partial\boldsymbol{A}}{\partial p_{i}} \right) \left( \boldsymbol{A}^{-1} \frac{\partial\boldsymbol{A}}{\partial p_{j}} \right) \right]. \qquad(5)$$

Tip 2: Inner products in the score

Several terms are Ω-inner products of signature columns, μ⋅j⊤Ωμ⋅k. Package helper .inner_product() evaluates that bilinear form; see also the squared Mahalanobis distance helper in the same family.

Unconstrained score and Hessian

Conditional on the reference moments ζ = (μ, {Σj}), the bulk follows y ∣ ζ, p ∼ 𝒩G(μp, Σ(p)). Up to an additive constant independent of p,

$$ \ell_{\boldsymbol{y}\mid\boldsymbol{\zeta}}(\boldsymbol{p}) = \tfrac{1}{2}\log\det\boldsymbol{\Omega}(\boldsymbol{p}) - \tfrac{1}{2} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega}(\boldsymbol{p}) (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p}). \qquad(6)$$

Important 3: The determinant term carries a factor $\tfrac{1}{2}$

Both terms of Eq. 6 are halved, because log det Ω = −log det Σ enters the Gaussian log-density as $-\tfrac{1}{2}\log\det\boldsymbol{\Sigma}(\boldsymbol{p})$. DeCovarT releases before 2.3.0 implemented −log det Σ(p), which doubled the determinant contribution. The objective was then not the stated Gaussian model, and it disagreed with the expected Fisher information of Sec. 4, which had always been the one implied by the correct density. With the factor restored, loglik_multivariate() matches mvtnorm::dmvnorm(..., log = TRUE) up to the constant $\tfrac{G}{2}\log(2\pi)$, and 𝔼[−H] = I(p) holds exactly (see Note 7). Point estimates move only slightly, because the two terms trade off differently, but likelihood values, AIC, likelihood-ratio statistics, and any curvature-based standard error do change.

Warning 4: What “unconstrained” means here

Sec. 2 treats each pj as an ambient coordinate in (0, 1) and ignores ∑jpj = 1. The resulting gradient and Hessian are the building blocks of the ILR chain rule in Sec. 3. They are not the coordinates used by Marquardt–Levenberg or Newton–Raphson in the package (those optimise in ρ-space). The box-constrained L-BFGS-B path is the exception: it works directly in p with [0, 1]J boxes and does not enforce the simplex exactly.

Applying Eq. 2 and Eq. 3 with ∂Σ/∂pj = 2pjΣj gives the ambient score. Colour coding matches the manuscript: purple determinant / precision terms, blue mean residual, orange covariance quadratic.

Theorem 1 (Unconstrained gradient) $$ \begin{aligned} \frac{\partial\ell}{\partial p_{j}} &= \underbrace{-p_{j}\,\mathrm{tr}\bigl( \boldsymbol{\Omega}\boldsymbol{\Sigma}_{j} \bigr)}_{\text{purple}} + \underbrace{ (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega} \boldsymbol{\mu}_{\cdot j} }_{\text{blue}} \\ &\quad+ \underbrace{ p_{j}\, (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega} \boldsymbol{\Sigma}_{j} \boldsymbol{\Omega} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p}) }_{\text{orange}}. \end{aligned} \qquad(7)$$

Differentiating Eq. 7 once more yields the ambient Hessian H ∈ ℳJ × J. Diagonal and off-diagonal blocks keep the same colour pairing.

Theorem 2 (Unconstrained Hessian) For i ∈ {1, …, J},

$$ \begin{aligned} \mathbf{H}_{ii} &= \frac{\partial^{2}\ell}{\partial p_{i}^{2}} = \underbrace{ -\,\mathrm{tr}(\boldsymbol{\Omega}\boldsymbol{\Sigma}_{i}) + 2p_{i}^{2}\,\mathrm{tr}\bigl( (\boldsymbol{\Omega}\boldsymbol{\Sigma}_{i})^{2} \bigr) }_{\text{purple}} \\ &\quad \underbrace{ -2p_{i} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega} \boldsymbol{\Sigma}_{i} \boldsymbol{\Omega} \boldsymbol{\mu}_{\cdot i} - \boldsymbol{\mu}_{\cdot i}^{\top} \boldsymbol{\Omega} \boldsymbol{\mu}_{\cdot i} }_{\text{blue}} \\ &\quad \underbrace{ -2p_{i} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega} \boldsymbol{\Sigma}_{i} \boldsymbol{\Omega} \boldsymbol{\mu}_{\cdot i} - (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega} \bigl( 4p_{i}^{2}\boldsymbol{\Sigma}_{i}\boldsymbol{\Omega}\boldsymbol{\Sigma}_{i} - \boldsymbol{\Sigma}_{i} \bigr) \boldsymbol{\Omega} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p}) }_{\text{orange}}, \end{aligned} \qquad(8)$$

and for i ≠ j,

$$ \begin{aligned} \mathbf{H}_{ij} &= \frac{\partial^{2}\ell}{\partial p_{i}\partial p_{j}} = \underbrace{ 2p_{i}p_{j}\, \mathrm{tr}\bigl( \boldsymbol{\Omega}\boldsymbol{\Sigma}_{j} \boldsymbol{\Omega}\boldsymbol{\Sigma}_{i} \bigr) }_{\text{purple}} \\ &\quad \underbrace{ -2p_{i} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega} \boldsymbol{\Sigma}_{i} \boldsymbol{\Omega} \boldsymbol{\mu}_{\cdot j} - \boldsymbol{\mu}_{\cdot i}^{\top} \boldsymbol{\Omega} \boldsymbol{\mu}_{\cdot j} }_{\text{blue}} \\ &\quad \underbrace{ -2p_{j} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega} \boldsymbol{\Sigma}_{j} \boldsymbol{\Omega} \boldsymbol{\mu}_{\cdot i} - 4p_{i}p_{j} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top} \boldsymbol{\Omega} \boldsymbol{\Sigma}_{i} \boldsymbol{\Omega} \boldsymbol{\Sigma}_{j} \boldsymbol{\Omega} (\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p}) }_{\text{orange}}. \end{aligned} \qquad(9)$$

The stationary equation ∇ℓ(p) = 0 has no closed root, so DeCovarT uses iterative solvers. Analytic gradients and Hessians are checked against Richardson extrapolation in numDeriv (Gilbert and Varadhan 2019).

Constrained optimisation on the simplex

Constrained optimisation may eliminate constraints by reparametrisation, enforce them softly with penalties or barriers, or enforce them exactly with projections / KKT methods. DeCovarT uses a C2 diffeomorphism from free Euclidean coordinates ρ ∈ ℝJ − 1 onto the open simplex, so every Newton step stays feasible by construction.

Log-sum-exp and softmax

The softmax of scores z ∈ ℝK is evaluated through the log-sum-exp (LSE) to avoid overflow,

$$ s_{i} = \frac{e^{z_{i}}}{\sum_{k}e^{z_{k}}} = e^{z_{i}-\operatorname{LSE}(\boldsymbol{z})}, \qquad \operatorname{LSE}(\boldsymbol{z}) = m+\log\sum_{k}e^{z_{k}-m}, \quad m=\max_{k}z_{k}. \qquad(10)$$

Its Jacobian is Jsoftmax = diag (s) − ss⊤, which is also ∇2LSE (z). The second derivative is the third-order tensor

Hijk = si[(δij − sj)(δik − sk) − sj(δjk − sk)].   (11)

Isometric log-ratio map

Choose a Helmert sub-matrix V ∈ ℝJ × (J − 1) with

VTV = IJ − 1,   VT1 = 0.   (12)

ILR coordinates and their inverse are (Pawlowsky-Glahn and Buccianti 2011)

ρ = VTlog p = VTclr (p),   p = softmax (Vρ) = 𝒞{exp (Vρ)},   (13)

where 𝒞 denotes closure to unit sum. No cell type is pinned as a reference (contrast Sec. 8). Package helpers helmert_basis(), isometric_logistic() and isometric_log_ratio() implement Eq. 12–Eq. 13 with a log-sum-exp softmax, matching the compositions::ilr() / ilrInv() convention on a Helmert basis (van den Boogaart et al. 2025; Egozcue et al. 2003).

Proposition 1 (ILR is a C2 diffeomorphism) ψ : ℝJ − 1 → (0, 1)J with ∑jpj = 1 is bijective and twice continuously differentiable, so every critical point in ρ-space maps to a unique open simplex composition and conversely. Any other valid ILR basis is V⋆ = VQ with QTQ = I, which rotates coordinates as ρ⋆ = QTρ and leaves eigenvalues, condition numbers and traces of properly transformed quadratic forms unchanged.

Note 5: Why ILR is the default chart

Changing an ALR denominator changes the coordinates, not the composition. If derivatives, Fisher information and the delta method are transformed correctly, $\hat{\boldsymbol{p}}$, $\ell(\hat{\boldsymbol{p}})$ and $\operatorname{Var}(\hat{\boldsymbol{p}})$ must agree after mapping back to the simplex. A material change in $\hat{\boldsymbol{p}}$ when the reference cell is permuted is therefore numerical or optimisation sensitivity, not a statistical reference dependence.

ILR still improves the chart because distinct orthonormal bases differ by an orthogonal rotation, whereas changing an ALR denominator is a generally non-orthogonal reparametrisation. Euclidean gradient norms, Hessian condition numbers and Newton step geometry are then natural up to rotation. A fixed Helmert basis introduces an ordering convention, not a denominator. The ALR formulae are kept in Sec. 8 for comparison and for the structural-zero discussion (Sec. 8.1).

Efficient derivatives of ψ

Write S(p) = diag (p) − ppT. The Jacobian is

$$ \mathbf{J}_{\boldsymbol{\psi}}(\boldsymbol{\rho}) = \frac{\partial\boldsymbol{p}}{\partial\boldsymbol{\rho}^{\mathsf{T}}} = \mathbf{S}(\boldsymbol{p})\mathbf{V}. \qquad(14)$$

jacobian_isometric_logistic() returns this J × (J − 1) matrix. Let $\bar{\boldsymbol{v}}=\mathbf{V}^{\mathsf{T}}\boldsymbol{p}$, $\boldsymbol{q}_{i}=\mathbf{V}_{i\cdot}-\bar{\boldsymbol{v}}$ and CV = VTS(p)V. Slice i of the map Hessian is

Hψi = pi(qiqiT − CV).   (15)

hessian_isometric_logistic() stores the tensor of shape (J − 1) × (J − 1) × J.

Change of variables: ambient p to ILR ρ

Let ℓρ(ρ) = ℓy ∣ ζ(ψ(ρ)). The chain rule maps the unconstrained derivatives of Sec. 2 into free Euclidean coordinates.

Theorem 3 (Constrained score and Hessian) ∇ρℓ = JψT∇pℓ = VTS(p)∇pℓ,   (16)

$$ \mathbf{H}_{\boldsymbol{\rho}} = \mathbf{J}_{\boldsymbol{\psi}}^{\mathsf{T}} \mathbf{H}_{\boldsymbol{p}} \mathbf{J}_{\boldsymbol{\psi}} + \sum_{i=1}^{J} \frac{\partial\ell}{\partial p_{i}} \mathbf{H}_{\psi_{i}}. \qquad(17)$$

At an interior KKT point every active score equals the same multiplier, ∂ℓ/∂pi = λ. Because ∑ipi(ρ) = 1 for every ρ, ∑iHψi = 0, so the second summand of Eq. 17 vanishes and Hρ = JψTHpJψ. Away from stationarity that contraction must be kept.

Important 6: Implementation map

Object Package function
Ambient score / Hessian gradient_loglik_unconstrained(), hessian_loglik_unconstrained()
Map Jacobian / Hessian jacobian_isometric_logistic(), hessian_isometric_logistic()
ILR score / Hessian gradient_loglik_constrained(), hessian_loglik_constrained()

Eq. 17 is the sum of a covariant pull-back of H and a contraction of the ambient score against the map Hessian. Omitting the second summand is a common bug when porting first-order code to Newton methods.

CLR coordinates and Bayesian formulations sit on the same footing (Pawlowsky-Glahn and Buccianti 2011): a Dirichlet or logistic-normal prior on p pairs naturally with the Gaussian bulk likelihood (the experimental MAP CTS path in scripts/experimental/03_04_DeCovarT_estimate_CTS_MAP_Bayesian.R, derived in the perspectives vignette, is one starting point). Regardless of coordinate system, the workflow is the same: optimise in an unconstrained parameterisation, then map back to the simplex for interpretation.

Expected Fisher information and Wald inference

Point estimation alone does not report uncertainty. Under a regular large-sample regime the MLE is asymptotically normal with covariance given by the inverse Fisher information. DeCovarT builds Wald standard errors from the expected Fisher information of the multivariate-normal convolution, then maps that information through the ILR chart.

Unconstrained expected Fisher information

For y ∼ 𝒩G(μp, Σ(p)) the expected Fisher information has entries

I(p)jk = μ⋅j⊤Ω(p)μ⋅k + 2pjpk tr(Ω(p)ΣjΩ(p)Σk).   (18)

The first summand is the mean contribution (a Ω-inner product of signature columns); the second is the covariance contribution of the quadratic map p ↦ Σ(p). See the multivariate-normal formula on Wikipedia: Fisher information. expected_fisher_unconstrained() implements Eq. 18.

Note 7: Observed versus expected information

The analytic Hessian of Sec. 2 is the observed information (up to sign). Wald intervals in decovart_fit use the expected information Eq. 18, which does not depend on the residual y − μp and is positive definite under the usual full-rank assumptions on (μ, {Σj}).

The two are consistent by construction. Substituting 𝔼[r] = 0 and 𝔼[rr⊤] = Σ(p) into Eq. 8 and Eq. 9 cancels the determinant traces against the residual traces and returns Eq. 18 exactly, i.e. 𝔼[−H] = I(p). That identity only holds with the factor $\tfrac{1}{2}$ of Note 3 in place.

Constrained case: ILR delta method

Fisher information transforms as a covariant quadratic form under the ILR chart. With Jψ = S(p)V,

Iρ = JψTI(p)Jψ = VTS(p)I(p)S(p)V.   (19)

The mean/covariance split of Eq. 18 pulls back separately, Iρ = Iρ, mean + Iρ, cov. Under regularity, $\hat{\boldsymbol{\rho}} \overset{a}{\sim} \mathcal{N}(\boldsymbol{\rho}_{0},I_{\boldsymbol{\rho}}^{-1})$ (Vaart 2000). The first-order delta method (Oehlert 1992) then returns the simplex covariance used by vcov.decovart_fit(),

$$ \mathrm{Var}(\hat{\boldsymbol{p}}) \approx \mathbf{J}_{\boldsymbol{\psi}} I_{\boldsymbol{\rho}}^{-1} \mathbf{J}_{\boldsymbol{\psi}}^{\mathsf{T}}. \qquad(20)$$

Diagonal square roots are asymptotic standard errors; Wald intervals at level 1 − α are p̂j ± z1 − α/2 SEj with zq = Φ−1(q), as in confint.decovart_fit(). vcov_ilr_delta() implements Eq. 19–Eq. 20. Orthogonal rotations of V leave Eq. 20 unchanged. The construction is undefined on the simplex boundary (the log-ratio chart blows up); the helper then returns NA with a warning.

Warning 8: Wald is the cheapest, not the safest, interval

Eq. 20 linearises the ILR chart at $\hat{\boldsymbol{p}}$, so it is neither invariant to a non-orthogonal reparametrisation nor confined to [0, 1], and it degenerates entirely when a proportion reaches a simplex face. The companion vignette MLE properties and asymptotic inference derives the alternatives: profile likelihood-ratio intervals, the chi-bar-square calibration of Chernoff and of Self and Liang for nulls on the boundary, a parametric bootstrap, and a reference-sample bootstrap of donors or cells (or a Dirichlet composition sweep). It also shows why these limits need replicate samples that share one composition. Shuffling genes or cell-type names of an averaged signature is not a bootstrap: the maximised likelihood is equivariant under relabelling.

Numerical speed-ups and solver safeguards

The analytic maps above are only half of a usable optimiser. Practical bottlenecks that showed up on the hybrid scenario in Covariance-driven hybrid scenario live in R/03_03_DeCovarT_estimate_ratios_frequentist.R.

Cache a Cholesky factorisation of Σ(p)

Within one iteration, objective / gradient / Hessian callbacks hit the same trial p. The helper .sigma_p_factorisation() caches a single Cholesky factor R with Σ(p) = RTR and returns log det Σ(p) = 2∑glog Rgg together with the precision Ω(p) = R−1R−T (via chol2inv) without repeating an O(G3) factorisation.

loglik_multivariate() then evaluates the Mahalanobis term exactly as mvtnorm::dmvnorm (Genz et al. 2026): it solves RTz = y − μp by backsolve(..., transpose = TRUE) and takes ‖z‖2, rather than forming Ω(p) and a dense quadratic form. The two routes agree to machine precision; the backsolve path is the one used for the objective, while the cached inverse is retained for the analytic score and Hessian, which need Ω(p) in several trace and inner-product terms. A QR factorisation of the covariance itself would recover the same log-determinant and quadratic form at a larger O(G3) constant and is not used: Cholesky is the natural factorisation of a symmetric positive-definite matrix.

The implemented log-density omits the additive $-\tfrac{G}{2}\log(2\pi)$ of the full Gaussian, which does not depend on p and therefore cannot change the MLE, the score, or the Hessian. Tests check that subtracting that constant recovers dmvnorm(..., log = TRUE) exactly.

Guard the box-constrained L-BFGS-B path

deconvolute_ratios_L_BFGS_B() optimises in p with boxes [0, 1]J that do not enforce 1⊤p = 1. Near ∑jpj ≈ 0 (or on a failed chol()), safeguarded wrappers return a finite penalty and a zero gradient rather than aborting optim().

marqLevAlg Hessian sign under minimize = FALSE

marqLevAlg (Philipps et al. 2023) implements Marquardt–Levenberg with the relative-distance-to-minimum (RDM) rule of Commenges et al. (Commenges et al. 2006). With minimize = FALSE, the package negates only fn and gr; the analytic hess is passed through unchanged while internal Cholesky routines assume a positive-definite matrix at the optimum. When maximising a log-likelihood, pass hess = function(...) -hessian_loglik_constrained(...). Reported upstream as VivianePhilipps/marqLevAlgParallel#3.

Model objects and S3 accessors

fit_decovart() returns a decovart_fit object with standard S3 methods for coefficients, fitted values, residuals, and (under regularity) a delta-method covariance. The chunk below mirrors the appendix equivariance setup:

library(DeCovarT)
genes <- paste0("g", 1:3)
cts <- paste0("ct", 1:2)
mu <- matrix(
  c(20, 40, 25, 35, 30, 22),
  nrow = 3,
  dimnames = list(genes, cts)
)
Sigma <- array(0, dim = c(3, 3, 2), dimnames = list(genes, genes, cts))
Sigma[, , 1] <- diag(c(1.0, 1.2, 0.8))
Sigma[, , 2] <- diag(c(1.5, 0.9, 1.1))
p_true <- c(0.65, 0.35)
y <- matrix(
  drop(mu %*% p_true + rnorm(3, sd = 0.05)),
  ncol = 1L,
  dimnames = list(genes, "s1")
)

fit <- fit_decovart(
  mu,
  y,
  Sigma = Sigma,
  method = "Marquardt-Levenberg",
  itmax = 80
)
coef(fit)
nobs(fit)
c(
  n_genes = attr(nobs(fit), "n_genes"),
  n_celltypes = attr(nobs(fit), "n_celltypes")
)
fitted(fit)
residuals(fit)
vcov(fit)

formula() and predict() are not implemented: DeCovarT does not forecast bulk expression.

Numerical consistency of the likelihood and its derivatives

Analytic formulae can be differentiated correctly and still implement the wrong objective. The checks below separate those two questions. They are the same identities that tests/testthat/test-03_03_DeCovarT.R and tests/testthat/test-03_06_decovart_inference.R enforce on every package build; the chunks here are a readable, small-(G) illustration.

Density against mvtnorm

The first check is against an independent implementation of the multivariate Gaussian log-density, not against our own inverse. With the $\tfrac{1}{2}\log\det$ factor restored, DeCovarT and mvtnorm::dmvnorm (Genz et al. 2026) differ by exactly $\tfrac{G}{2}\log(2\pi)$ (DeCovarT omits the negative constant):

n_genes <- length(y_chk)
ours <- loglik_multivariate(p_chk, y_chk, mu_chk, Sigma_chk)
reference <- if (requireNamespace("mvtnorm", quietly = TRUE)) {
  as.numeric(mvtnorm::dmvnorm(
    y_chk,
    mean = drop(mu_chk %*% p_chk),
    sigma = cov_chk,
    log = TRUE
  ))
} else {
  residual <- y_chk - drop(mu_chk %*% p_chk)
  as.numeric(
    -0.5 * determinant(cov_chk, logarithm = TRUE)$modulus -
      0.5 * drop(crossprod(residual, solve(cov_chk, residual))) -
      0.5 * n_genes * log(2 * pi)
  )
}
c(
  loglik_multivariate = ours,
  dmvnorm = reference,
  difference = ours - 0.5 * n_genes * log(2 * pi) - reference
)

A finite-difference check of the same objective cannot detect a shared factor-of-two error in the log-likelihood and its derivatives. That is why this density comparison is independent of numDeriv.

Richardson extrapolation of the score and Hessian

Once the objective is the intended Gaussian, the analytic score and Hessian are checked against Richardson extrapolation in numDeriv (Gilbert and Varadhan 2019). Richardson’s tableau extrapolates central differences to cancel successive even powers of the step size, which is more accurate than a single two-point stencil at the same O(ε−1) cost in function evaluations.

if (!requireNamespace("numDeriv", quietly = TRUE)) {
  numderiv_report <- c(gradient_rel = NA_real_, hessian_rel = NA_real_)
} else {
  grad_n <- numDeriv::grad(
    loglik_multivariate,
    p_chk,
    method = "Richardson",
    method.args = list(eps = 1e-4, r = 6),
    y = y_chk,
    mean_signature_matrix = mu_chk,
    Sigma = Sigma_chk
  )
  hess_n <- numDeriv::hessian(
    loglik_multivariate,
    p_chk,
    method = "Richardson",
    method.args = list(eps = 1e-4, r = 4),
    y = y_chk,
    mean_signature_matrix = mu_chk,
    Sigma = Sigma_chk
  )
  grad_a <- gradient_loglik_unconstrained(
    p_chk, y_chk, mu_chk, Sigma_chk
  )
  hess_a <- hessian_loglik_unconstrained(
    p_chk, y_chk, mu_chk, Sigma_chk
  )
  numderiv_report <- c(
    gradient_rel = max(abs(grad_a - grad_n)) / max(abs(grad_n)),
    hessian_rel = max(abs(hess_a - hess_n)) / max(abs(hess_n))
  )
}
numderiv_report

Relative discrepancies of order (10^{-8}) are typical on this three-gene toy; that is the same tolerance used in the test suite.

Expected information versus Monte Carlo ([-])

Matching numDeriv shows that the derivatives belong to the implemented objective. Matching the expected Fisher information Eq. 18 to a Monte Carlo average of (-(; )) shows that the objective is the Gaussian whose information matrix we use for Wald intervals. Draw (^{(b)}_G(, ())) at the true () (the parametric bootstrap of bootstrap_decovart(), used here as a diagnostic rather than as an inferential procedure):

set.seed(2)
n_mc <- 40L
info <- expected_fisher_unconstrained(p_chk, mu_chk, Sigma_chk)
neg_hess <- array(NA_real_, c(3L, 3L, n_mc))
for (b in seq_len(n_mc)) {
  y_b <- MASS::mvrnorm(
    n = 1L,
    mu = drop(mu_chk %*% p_chk),
    Sigma = cov_chk
  )
  neg_hess[, , b] <- -hessian_loglik_unconstrained(
    p_chk, y_b, mu_chk, Sigma_chk
  )
}
mc_info <- apply(neg_hess, c(1L, 2L), mean)
c(
  rel_frobenius = sqrt(sum((mc_info - info)^2)) / sqrt(sum(info^2))
)

The same identity is what closed the factor-of-two inconsistency: with the old (-) objective, ([-]) could not equal (I()). A few dozen draws already bring the relative Frobenius gap below (10^{-2}) on this toy; the test suite uses a tighter Monte Carlo.

A related, weaker check compares the empirical covariance of parametric-bootstrap MLEs to the Cramér–Rao map vcov_ilr_delta(). That comparison is noisy at small replication (see MLE properties for why one bulk column is not an asymptotic sample) and is therefore left to the inference vignette and to bootstrap_decovart().

Structure-aware covariance backends

Every solver call factorises Σ(p) = ∑jpj2Σj. The dense Cholesky in .sigma_p_factorisation() is the universal default and the only method that needs no assumptions. When the covariance inherits structure from the underlying gene network, the same log-determinant and solves can be obtained far more cheaply. Exploiting band, sparse and diagonal-plus-low-rank structure in a Newton-type solver is standard practice (Boyd et al. 2004, sec. 9.7); the backends below make that structure explicit for DeCovarT.

new_decovart_covariance() wraps the covariance array with a declared structure, and the operators sigma_logdet(), sigma_solve(), sigma_quadform() and sigma_trace_precision_times() return exactly the quantities the objective, score and Hessian need. They expose operations, never a materialised dense precision: a sparse or Woodbury solve must not be silently undone by a dense chol2inv.

Decision tree

Structure is supplied explicitly, never inferred from floating-point near-zeros, so a numerical near-zero is never mistaken for an exact conditional independence. covariance_structure_from_graph_model() maps a network topology to the recommended backend. A sparse precision Ωj does not imply a sparse covariance Σj; only structure that survives in Σ(p) itself accelerates the likelihood.

%%{init: {"theme": "sandstone"}}%%
flowchart TD
  A["Sigma(p) = sum_j p_j^2 Sigma_j"] --> B{"Declared covariance structure?"}
  B -->|"none (default)"| D["Dense Cholesky\nO(G^3)"]
  B -->|"disconnected modules\n(stochastic block model)"| E["Block Cholesky\nO(sum_b G_b^3)"]
  B -->|"ordered local dependence\n(band / AR)"| F["Banded Cholesky\nO(G b^2)"]
  B -->|"fixed sparse pattern\n(Erdos-Renyi, scale-free, small-world)"| G["Sparse Cholesky\nsymbolic ordering computed once"]
  B -->|"shared programs / bridge hubs\n(pathways, housekeeping genes)"| H["Woodbury / Schur\nO(G r^2)"]
Figure 1: Choosing a covariance backend; dense Cholesky is the default and structured backends are exact refinements.

Each refinement is exact, not approximate:

  • Block (disconnected gene modules, e.g. a stochastic block model with no between-block edges (Holland et al. 1983)). Σ(p) is block-diagonal, so the log-likelihood separates over blocks and one Cholesky per block replaces one dense Cholesky.
  • Band / AR. Only genes within bandwidth b covary; a banded Cholesky costs O(Gb2) (Golub and Van Loan 2013).
  • Sparse (Erdős–Rényi, scale-free (Barabási and Albert 1999), small-world (Watts and Strogatz 1998)). The fill-reducing ordering is a symbolic factorisation computed once and reused for every numeric refactor as p changes; only the numerical Cholesky repeats (Chen et al. 2008).
  • Diagonal + low rank. Genes interact through r ≪ G shared regulatory programs, Σj = Dj + UCjUT — the factor-analysis form Σ = WWT + Ψ. The Woodbury identity and the matrix-determinant lemma (Hager 1989) move the O(G3) cost to the low rank r. A few bridge / housekeeping hub genes give the arrow-matrix special case handled by the same Schur-complement algebra.
# Declare structure explicitly, then use the operator interface.
cov <- new_decovart_covariance(Sigma, structure = "diag_lowrank",
  diagonal = D, loadings = U, core = C)
ld <- sigma_logdet(cov, p)            # log det Sigma(p)
z  <- sigma_solve(cov, p, residual)   # Sigma(p)^{-1} r, no dense inverse
q  <- sigma_quadform(cov, p, residual)

Computational complexity

For a dense symmetric positive-definite G × G matrix the Cholesky factorisation costs about G3/3 flops before the cheaper triangular solves (Golub and Van Loan 2013). The structured backends reduce that leading term:

Table 1: Leading factorisation cost per backend. G genes, J cell types, B blocks of size Gb, bandwidth b, rank r ≪ G.
Structure When it applies Factorisation cost Backend argument
Dense universal default O(G3) (default)
Block disconnected modules (SBM) O(∑bGb3);  ≈ B2 saving for B equal blocks structure = "block"
Band / AR bandwidth b O(Gb2) structure = "band"
Sparse fixed sparse pattern pattern-dependent; symbolic ordering reused across p structure = "sparse"
Diagonal + low rank r shared programs O(Gr2) solve, O(r3) core structure = "diag_lowrank"

These are worst-case flop counts and only rough predictors of wall-clock time on cache-based hardware with tuned BLAS; the block and low-rank savings ( ≈ B2 and O(G2) → O(Gr2) respectively) are the most robust. The "band" and "sparse" backends use the imported Matrix package (supernodal sparse Cholesky in CHOLMOD, (Chen et al. 2008)). scripts/supp_S2_covariance_inversion.R verifies each backend against the dense path to machine precision. Approximate sparse Cholesky of otherwise dense kernel covariances is a separate research line (Schäfer et al. 2020); the "sparse" backend here is exact for a declared sparse pattern, not that approximation.

Appendix: additive log-ratio chart

Pin the last category as reference (ρJ ≡ 0). The additive logistic transform ψ (inverse ALR) and its inverse are (Aitchison 1982)

$$ p_{j} = \frac{e^{\rho_{j}}}{\sum_{k<J}e^{\rho_{k}}+1}\ (j<J), \qquad p_{J} = \frac{1}{\sum_{k<J}e^{\rho_{k}}+1}, \qquad \rho_{j} = \log\!\left(\frac{p_{j}}{p_{J}}\right). \qquad(21)$$

Equivalently, p = softmax (ρ1, …, ρJ − 1, 0). Package helpers additive_logistic() and additive_log_ratio() implement Eq. 21; they match compositions::alrInv() / alr() with the last part as reference (van den Boogaart et al. 2025). The Jacobian is the softmax Jacobian with the reference column dropped, ∂pi/∂ρa = pi(δia − pa) for a = 1, …, J − 1 (jacobian_additive_logistic(); hessian_additive_logistic()).

Note 9: Why ALR is not the default

ALR is a valid C2 diffeomorphism onto the open simplex. The last cell type is nevertheless a privileged denominator: permuting which type occupies position J changes the Euclidean metric on coordinate space. ILR / Helmert (Sec. 3.2) removes that privilege up to orthogonal rotation, which is why the solvers and vcov.decovart_fit() use ILR. Keep ALR when comparing against a multinomial-logit reference-category parameterisation.

Structural zeros and folded simplex models

Neither ALR nor ILR can represent an exact zero at a finite coordinate, because both contain a logarithm. For a true simplex face the package therefore does not send log 0 through the chart. It fixes the inactive set, optimises ILR coordinates on the remaining active face (restricted_mle_decovart()), and checks KKT inequalities for the zero coordinates (boundary_diagnostics(); see also the MLE-properties vignette).

ALR remains in this appendix as the classical reference-category chart (Aitchison 1982). Folded α-models on the simplex (Tsagris and Stewart 2020) and nonparametric regressions that admit exact zeros without a log-ratio chart (Tsagris et al. 2023) are related data-side constructions. They are not used in the DeCovarT convolution, whose mean map is still linear in p. Rounded zeros (detection limits) remain an imputation problem; true structural zeros belong on a face and should be fitted as an active-set / KKT problem, not by substituting a small positive constant into Eq. 13.

References

Aitchison, J. 1982. “The Statistical Analysis of Compositional Data.” Journal of the Royal Statistical Society: Series B (Methodological) 44 (2): 139–60. https://doi.org/10.1111/j.2517-6161.1982.tb01195.x.
Barabási, Albert-László, and Réka Albert. 1999. “Emergence of Scaling in Random Networks.” Science 286. https://doi.org/10.1126/science.286.5439.509.
Boyd, Stephen, Stephen P. Boyd, and Lieven Vandenberghe. 2004. Convex Optimization. 1st ed. Cambridge University Press; Cambridge University Press. https://doi.org/10.1017/cbo9780511804441.
Chen, Yanqing, Timothy A. Davis, William W. Hager, and Sivasankaran Rajamanickam. 2008. “Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate.” ACM Transactions on Mathematical Software 35 (3): 22:1–14. https://doi.org/10.1145/1391989.1391995.
Commenges, Daniel, Helene Jacqmin-Gadda, Cecile Proust, and Jeremie Guedj. 2006. A Newton-Like Algorithm for Likelihood Maximization: The Robust-Variance Scoring Algorithm. arXiv. https://doi.org/10.48550/arxiv.math/0610402.
Egozcue, J. J., V. Pawlowsky-Glahn, G. Mateu-Figueras, and C. Barceló-Vidal. 2003. “Isometric Logratio Transformations for Compositional Data Analysis.” Mathematical Geology 35 (3): 279–300. https://doi.org/10.1023/a:1023818214614.
Genz, Alan, Frank Bretz, Tetsuhisa Miwa, Xuefei Mi, and Torsten Hothorn. 2026. Mvtnorm: Multivariate Normal and t Distributions. https://codeberg.org/thothorn/mvtnorm.
Gilbert, Paul, and Ravi Varadhan. 2019. numDeriv: Accurate Numerical Derivatives. http://optimizer.r-forge.r-project.org/.
Golub, Gene H., and Charles F. Van Loan. 2013. Matrix Computations. JHU Press.
Hager, William W. 1989. “Updating the Inverse of a Matrix.” SIAM Review 31 (2): 221–39. https://doi.org/10.1137/1031049.
Holland, Paul W., Kathryn Blackmond Laskey, and Samuel Leinhardt. 1983. “Stochastic Blockmodels: First Steps.” Social Networks 5. https://doi.org/10.1016/0378-8733(83)90021-7.
Oehlert, Gary W. 1992. “A Note on the Delta Method.” The American Statistician 46 (1): 27–29. https://doi.org/10.1080/00031305.1992.10475842.
Pawlowsky-Glahn, Vera, and Antonella Buccianti, eds. 2011. Compositional Data Analysis: Theory and Applications. Wiley.
Philipps, Viviane, Cecile Proust-Lima, Melanie Prague, Boris Hejblum, Daniel Commenges, and Amadou Diakite. 2023. marqLevAlg: A Parallelized General-Purpose Optimization Based on Marquardt-Levenberg Algorithm.
Schäfer, Florian, T. J. Sullivan, and Houman Owhadi. 2020. Compression, Inversion, and Approximate PCA of Dense Kernel Matrices at Near-Linear Computational Complexity. arXiv. https://doi.org/10.48550/arxiv.1706.02205.
Tsagris, Michail, Abdulaziz Alenazi, and Connie Stewart. 2023. “Flexible Non-Parametric Regression Models for Compositional Response Data with Zeros.” Statistics and Computing 33 (5): 106. https://doi.org/10.1007/s11222-023-10277-5.
Tsagris, Michail, and Connie Stewart. 2020. “A Folded Model for Compositional Data Analysis.” Australian & New Zealand Journal of Statistics 62 (2): 249–77. https://doi.org/10.1111/anzs.12289.
Vaart, A. W. van der. 2000. Asymptotic Statistics. Cambridge University Press.
van den Boogaart, K. Gerald, Raimon Tolosana-Delgado, and Matevz Bren. 2025. Compositions: Compositional Data Analysis. http://www.stat.boogaart.de/compositions/.
Watts, Duncan J., and Steven H. Strogatz. 1998. “Collective Dynamics of ‘Small-World’ Networks.” Nature 393 (6684): 440–42. https://doi.org/10.1038/30918.