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(). Sec. 3 compares the softmax, additive, centred and isometric log-ratio charts and explains why the solvers use ILR.

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.

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)$$

Applying Eq. 2 and Eq. 3 with ∂Σ/∂p_j = 2p_jΣ_j gives the ambient score.

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

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

The cellular ratios satisfy p_j ≥ 0 and ∑_jp_j = 1, the compositional or unit-simplex constraint. An optimiser can handle it in three ways: eliminate it with a reparametrisation, enforce it softly with penalties or barriers, or enforce it exactly with projections and Karush–Kuhn–Tucker (KKT) conditions. DeCovarT uses the first route: a C² diffeomorphism from free Euclidean coordinates onto the open simplex, so every Newton step stays feasible by construction. Four such charts are common: the softmax, the additive log-ratio (ALR), the centred log-ratio (CLR) and the isometric log-ratio (ILR). This section defines the notation they share (Sec. 3.1), presents each chart (Sec. 3.2–Sec. 3.5), explains why DeCovarT uses ILR (Sec. 3.6) and compares the four side by side (Sec. 3.7). The derivatives of the ILR chart (Sec. 3.8, Sec. 3.9), their link to the Lagrangian (Sec. 3.10) and perspectives on Bayesian charts and exact zeros (Sec. 3.11) follow.

The compositional background is chapter 4 of Modelling and Analysis of Compositional Data (Pawlowsky-Glahn et al. 2015), cited below by page.

Notation for compositions and charts

The closure operator. 𝒞 rescales any positive vector to unit sum. It discards the total and keeps the relative sizes: 𝒞(cx) = 𝒞(x) for every c > 0. A composition is therefore an equivalence class of proportional positive vectors, and 𝒞 picks its representative on the simplex (Pawlowsky-Glahn et al. 2015, 33–35). Two consequences matter here.

• Only scale-free functions of p are meaningful, and the log-ratios log (p_j/p_ℓ) are the simplest. A constant added to log p does not change the composition, so every log-ratio chart must cancel the direction 1_J.
• Closure creates dependence. Library-size normalisation of counts is a closure, and it induces negative covariances between parts that need not exist in the absolute abundances.

Closure also gives the simplex a Euclidean vector-space structure, the Aitchison geometry: perturbation x ⊕ y = 𝒞(x ∘ y) plays the role of addition, powering α ⊙ x = 𝒞(x^α) the role of scalar multiplication, and the inner product is ⟨x, y⟩_a = clr (x)^Tclr (y) (Pawlowsky-Glahn et al. 2015, Property 4.2, p. 35). The log-ratio charts below are linear maps for this structure; they differ in whether they also preserve its inner product.

Softmax and log-sum-exp

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

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

Its Jacobian is J_softmax = S(s) = diag (s) − ss^T, which is also ∇²LSE (z). The second derivative is the third-order tensor

H_jℓm = s_j[(δ_jℓ − s_ℓ)(δ_jm − s_m) − s_ℓ(δ_ℓm − s_m)].   (11)

The softmax with J free scores is over-parametrised: softmax (z + c1_J) = softmax (z), because closure removes the common factor e^c. Accordingly S(s)1_J = 0, the Jacobian has rank J − 1, and the Hessian of ℓ(softmax (z)) has a zero eigenvalue along 1_J: a Newton system in z is singular. The three log-ratio charts are three ways to remove this redundant direction. Each writes p = softmax (Bθ) with θ ∈ ℝ^J − 1 and a matrix B ∈ ℝ^{J × (J − 1)} whose columns, together with 1_J, span ℝ^J. ALR pins one score to zero, CLR constrains the scores to sum to zero, and ILR spans 1_J^⟂ with an orthonormal basis.

Additive log-ratio chart

Pin a reference cell type r, by default the last one (r = J). The ALR coordinates and their inverse, the additive logistic transform, are (Aitchison 1982)

$$ \rho^{\mathrm{alr}}_{j} = \log\!\left(\frac{p_{j}}{p_{J}}\right), \qquad p_{j} = \frac{e^{\rho^{\mathrm{alr}}_{j}}}{1+\sum_{\ell<J}e^{\rho^{\mathrm{alr}}_{\ell}}}\ (j<J), \qquad p_{J} = \frac{1}{1+\sum_{\ell<J}e^{\rho^{\mathrm{alr}}_{\ell}}}. \qquad(12)$$

Equivalently p = softmax (ρ₁^alr, …, ρ_J − 1^alr, 0) and ρ^alr = D_Jlog p. The “1” in the denominator is the reference part; the word additive refers to the sum of exponentials in that denominator, not to the Aitchison perturbation (Pawlowsky-Glahn et al. 2015, 47). Package helpers additive_logistic() and additive_log_ratio() implement Eq. 12 and 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, ∂p_j/∂ρ_a^alr = p_j(δ_ja − p_a) for a = 1, …, J − 1 (jacobian_additive_logistic(); hessian_additive_logistic()).

What the reference changes, and what it does not. All ALR charts describe the same family of compositions; only the coordinates differ. Sec. 3.6 shows that the estimate $\hat{\boldsymbol{p}}$, the maximised log-likelihood and the delta-method covariance of $\hat{\boldsymbol{p}}$ do not depend on the reference. The numerical geometry does: ALR coordinates refer to an oblique, non-orthogonal basis of the simplex (Pawlowsky-Glahn et al. 2015, 47), so Euclidean lengths, angles and Hessian condition numbers in ρ^alr change with the reference.

Centred log-ratio coefficients

The CLR transform divides each part by the geometric mean (Aitchison 1982):

$$ \boldsymbol{\xi} = \operatorname{clr}(\boldsymbol{p}) = \log\frac{\boldsymbol{p}}{g(\boldsymbol{p})} = \mathbf{H}_{J}\log\boldsymbol{p}, \qquad \boldsymbol{p} = \mathcal{C}(\exp\boldsymbol{\xi}) = \operatorname{softmax}(\boldsymbol{\xi}). \qquad(13)$$

CLR treats every cell type symmetrically and preserves the Aitchison geometry: distances and inner products between compositions equal Euclidean ones between their CLR vectors (Pawlowsky-Glahn et al. 2015, Property 4.2, p. 35). Pairwise log-ratios are differences, log (p_j/p_ℓ) = ξ_j − ξ_ℓ. The price is a constraint on the coefficients: 1_J^Tξ = 0. CLR gives J coefficients for J − 1 degrees of freedom, so they are not coordinates in a basis, a random CLR vector has a singular covariance matrix, and the coefficients change when a subcomposition is analysed (Pawlowsky-Glahn et al. 2015, 34–35). Optimising over a free ξ ∈ ℝ^J is the over-parametrised softmax of Sec. 3.2 again.

CLR is nevertheless the hub between the other charts. Since D_r1_J = 0 and V^T1_J = 0, both ALR and ILR are linear images of CLR: ρ^alr = D_rξ and ρ = V^Tξ.

Isometric log-ratio chart

ILR expresses ξ in an orthonormal basis of the hyperplane 1_J^⟂ (Egozcue et al. 2003). Choose V ∈ ℝ^{J × (J − 1)} with

V^TV = I_J − 1,   V^T1_J = 0,   VV^T = H_J.   (14)

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

ρ = V^Tlog p = V^Tclr (p),   p = softmax (Vρ) = 𝒞{exp (Vρ)}.   (15)

In the book’s notation the contrast matrix is Ψ = V^T, and Eq. 14 is its Eq. (4.5); Eq. 15 is Eqs. (4.7)–(4.8) (Pawlowsky-Glahn et al. 2015, 36–37). ILR keeps the strengths of both previous charts: like ALR it gives J − 1 unconstrained coordinates, and like CLR it preserves the Aitchison geometry, now with d_a(x, y) = ‖ilr (x) − ilr (y)‖ (Property 4.5, p. 37). No cell type is a reference. Package helpers helmert_basis(), isometric_logistic() and isometric_log_ratio() implement Eq. 14–Eq. 15 with a log-sum-exp softmax, matching the compositions::ilr() / ilrInv() convention on a Helmert basis (van den Boogaart et al. 2025).

The Helmert basis: definition, role and use

Definition. The ILR chart needs one orthonormal basis of 1_J^⟂; there are infinitely many. The Helmert basis is the one with a closed form. Column k = 1, …, J − 1 of helmert_basis() is

$$ \boldsymbol{v}_{k} = \frac{1}{\sqrt{k(k+1)}} \bigl(\underbrace{1,\ldots,1}_{k},\,-k,\,\underbrace{0,\ldots,0}_{J-k-1}\bigr)^{\mathsf{T}}. \qquad(16)$$

Each column sums to zero and has unit norm, and two columns are orthogonal because the later one is constant on the support of the earlier one. The basis is the Gram–Schmidt orthonormalisation of the generating system of the simplex (Egozcue et al. 2003), and its contrast matrix is a Helmert matrix without the constant row (Pawlowsky-Glahn et al. 2015, Example 4.6, p. 40). The book lists the same basis in reverse order, with the first coordinate contrasting the first J − 1 parts with the last one.

Role: every coordinate is a balance. Inserting Eq. 16 into Eq. 15 gives

$$ \rho_{k} = \boldsymbol{v}_{k}^{\mathsf{T}}\log\boldsymbol{p} = \sqrt{\frac{k}{k+1}}\, \log\frac{g(p_{1},\ldots,p_{k})}{p_{k+1}}, \qquad(17)$$

the normalised log-ratio between the geometric mean of the first k cell types and type k + 1. This is the balance of a sequential binary partition: level k splits {1, …, k + 1} into {1, …, k} (sign +) and {k + 1} (sign −) (Pawlowsky-Glahn et al. 2015, Eq. (4.9), pp. 38–39). For J = 4 the sign matrix is

The factor $\sqrt{k/(k+1)}$ makes balances comparable across levels. Balances are geometrically orthogonal, so each carries information the others do not; geometric orthogonality does not imply that their estimates are uncorrelated (Pawlowsky-Glahn et al. 2015, 41). The Helmert basis thus does three jobs at once:

• it cancels the closure direction, because every column sums to zero: adding a constant to log p leaves ρ unchanged;
• it removes the redundant softmax direction while keeping Euclidean geometry, because the columns are orthonormal;
• it is explicit and sparse, so Vρ costs O(J) operations per column and needs no QR or Gram–Schmidt step at run time.

Use in practice. DeCovarT builds V once from J and proceeds as follows.

1. Map a starting composition to coordinates with isometric_log_ratio(p), which computes V^Tlog p.
2. Iterate on ρ ∈ ℝ^J − 1 with the score and Hessian of Theorem 3; no step can leave the simplex.
3. Map back with isometric_logistic(rho), a log-sum-exp softmax of Vρ.
4. Push the covariance of $\hat{\boldsymbol{\rho}}$ to $\hat{\boldsymbol{p}}$ with the Jacobian S(p)V (Sec. 4.2).

The order of the cell types fixes which balances are reported, not the estimate. Both helpers accept another basis through their V argument: for instance, a partition that first contrasts immune with stromal types gives coordinates with a biological reading. Any such basis is VQ for an orthogonal Q (Proposition 2).

Figure 1 summarises the construction. Gram–Schmidt applied to the successive differences e_k − e_k + 1, which span 1_J^⟂, returns exactly the columns of Eq. 16 for every J (panel D); the resulting chart is isometric, so equal Aitchison lengths become equal Euclidean lengths (panels B–C).

Why DeCovarT uses the ILR chart

The question is which chart to optimise in, given that all of them describe the same model. Two kinds of invariance answer it: statistical invariance, which every chart has, and numerical invariance, which only ILR has under a change of convention.

Statistical invariance holds for every chart. Let θ^′ = Aθ be any invertible linear change of coordinates, for instance between ILR and ALR (ρ^alr = A_rρ with A_r = D_rV) or between two ALR references. The likelihood is a function of p alone, so the maximiser is equivariant and $\hat{\boldsymbol{p}}$ and $\ell(\hat{\boldsymbol{p}})$ are the same in every chart. The delta method is invariant too: the Jacobian becomes J_ψA⁻¹ and the information A^−TℐA⁻¹, so

(J_ψA⁻¹)(A^−TℐA⁻¹)⁻¹(J_ψA⁻¹)^T = J_ψℐ⁻¹J_ψ^T.   (18)

In this sense ALR is reference-invariant: whichever cell type is the denominator, a correct implementation returns the same $\hat{\boldsymbol{p}}$ and the same Wald covariance for $\hat{\boldsymbol{p}}$. A material change in $\hat{\boldsymbol{p}}$ when the reference is permuted signals a numerical or convergence problem, not a statistical dependence on the reference. Pure Newton steps share this invariance, because the Newton step is invariant under affine changes of coordinates (Boyd et al. 2004, sec. 9.5).

Numerical invariance differs. What changes with the chart is the Euclidean geometry the solver sees: gradient norms, damping, stopping rules and the condition number of the Hessian, which controls the floating-point accuracy of every linear solve and of ℐ⁻¹.

Proof sketch. (i) is a similarity by an orthogonal matrix. For (ii), the rows of D_r are e_j − e_r, so D_rD_r^T = I_J − 1 + 1_J − 11_J − 1^T, with eigenvalues 1 (multiplicity J − 2) and J. Because D_r1_J = 0 and V is an orthonormal basis of 1_J^⟂, A_r has the same singular values as D_r. The right singular vector for $\sqrt{J}$ is proportional to V^TD_r^T1_J − 1 = V^T(1_J − Je_r) = −J V^Te_r. The bounds follow from κ(A^−THA⁻¹) ≤ κ(A)²κ(H) applied in both directions. ▫

The constant $\kappa(\mathbf{A}_{r})=\sqrt{J}$ does not depend on r, but the direction it distorts does, so the condition number of the ALR Hessian depends on the reference. Even a perfectly conditioned problem, H_ρ = −I, has κ(H_ρ^alr) = J in ALR coordinates. For an illustrative random negative-definite Hessian with J = 6 (not a DeCovarT fit), κ = 171 in every orthonormal basis, but between 164 and 420 across the six ALR references. A rare cell type, or one whose signature is nearly collinear with others, already has weak curvature along its own direction; choosing it as the ALR reference divides that curvature by J and makes the problem harder.

These differences affect the solver, not the target:

The default solver is Newton–Raphson on ρ (deconvolute_ratios_Newton_Raphson()), a second-order ascent that uses the analytic ILR score and Hessian of Theorem 3. Its pure Newton direction is the same in every chart; its safeguards, a spherical trust region and Euclidean step-length tolerances in stats::nlminb(), are unchanged by a rotation of the ILR basis but not by a change of ALR reference.

Side-by-side comparison of the four charts

All four charts guarantee feasibility, and none reaches a face of the simplex at finite coordinates. For likelihood inference ILR is usually the best suited: it is the only chart that is at once minimal, symmetric in the cell types and isometric, which is why Egozcue et al. introduced it to avoid the drawbacks of both ALR and CLR (Egozcue et al. 2003; Pawlowsky-Glahn et al. 2015, 33). Its cost is interpretation: a balance is less direct than a ratio to one reference, and the Helmert balances depend on the order of cell types unless a partition with a biological meaning is chosen.

Efficient derivatives of ψ

Write S(p) = diag (p) − pp^T. 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(20)$$

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

H_ψj = p_j(h_jh_j^T − C_V).   (21)

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.

Lagrangian, KKT conditions and what the ILR chart removes

Without a chart, the estimation problem is

max_{p ∈ ℝ^J} ℓ(p)  subject to  1_J^Tp = 1,   p ≥ 0,

with Lagrangian

ℒ(p, λ, ν) = ℓ(p) − λ(1_J^Tp − 1) + ν^Tp.   (24)

The constraints are linear, so a constraint qualification holds and every local maximiser satisfies the KKT conditions (Boyd et al. 2004, sec. 5.5.3):

$$ \nabla_{\boldsymbol{p}}\ell(\hat{\boldsymbol{p}}) -\lambda\mathbf{1}_{J} +\boldsymbol{\nu} =\boldsymbol{0}, \qquad \mathbf{1}_{J}^{\mathsf{T}}\hat{\boldsymbol{p}}=1, \qquad \boldsymbol{\nu}\ge\boldsymbol{0}, \qquad \nu_{j}\hat{p}_{j}=0. \qquad(25)$$

• Interior point (p̂_j > 0 for all j): complementary slackness gives ν = 0, hence ∂ℓ/∂p_j = λ for every j. Every cell type has the same marginal gain, so moving mass from one type to another does not change ℓ to first order. Multiplying by $\hat{\boldsymbol{p}}^{\mathsf{T}}$ gives the multiplier, $\lambda=\hat{\boldsymbol{p}}^{\mathsf{T}}\nabla_{\boldsymbol{p}}\ell(\hat{\boldsymbol{p}})$.
• Boundary point (p̂_j = 0): ∂ℓ/∂p_j = λ − ν_j ≤ λ. Adding some of type j would not increase the likelihood.

The chart removes the multiplier. Because S(p)1_J = p − p(1_J^Tp) = 0, the ILR score of Eq. 22 ignores any multiple of 1_J added to the ambient score:

∇_ρℓ = V^TS(p)(∇_pℓ − c 1_J)  for every c,   S(p)∇_pℓ = diag (p)(∇_pℓ − λ(p)1_J),   (26)

with λ(p) = p^T∇_pℓ. The ILR score is therefore the KKT stationarity residual of Eq. 25, weighted by p, with the multiplier set to its current value. The vector S(p)∇_pℓ lies in 1_J^⟂, on which V^T is an isometry, so

∥∇_ρℓ∥ = ∥S(p)∇_pℓ∥,   ∇_ρℓ = 0 ⇔ ∇_pℓ = λ1_J.

Stationary points in ρ are exactly the interior KKT points, the multiplier is recovered afterwards as $\hat{\lambda}=\hat{\boldsymbol{p}}^{\mathsf{T}}\nabla_{\boldsymbol{p}}\ell(\hat{\boldsymbol{p}})$, and the score norm that boundary_diagnostics() reports does not depend on the basis. Under ALR the same equivalence holds, but the score norm changes with the reference.

Second order: the Hessian of the Lagrangian. Since ∑_jψ_j(ρ) = 1 for every ρ, ∑_jH_ψj = 0, and the curvature term of Eq. 23 also ignores the multiplier:

$$ \sum_{j=1}^{J}\frac{\partial\ell}{\partial p_{j}}\mathbf{H}_{\psi_{j}} = \sum_{j=1}^{J}\Bigl(\frac{\partial\ell}{\partial p_{j}}-\lambda\Bigr) \mathbf{H}_{\psi_{j}}. $$

At an interior KKT point every bracket vanishes, so

H_ρ = J_ψ^T∇_pp²ℒ J_ψ,   ∇_pp²ℒ = H_p,   (27)

because linear constraints add no curvature to the Lagrangian. The columns of J_ψ = S(p)V span 1_J^⟂, the directions d with 1_J^Td = 0 that keep the sum equal to one. Eq. 27 is thus the Hessian of the Lagrangian reduced to the tangent space of the constraint, and H_ρ ≺ 0 is the classical second-order sufficient condition d^TH_pd < 0 for all d ≠ 0 with 1_J^Td = 0. Away from stationarity the curvature term must be kept: it is the price of a curved chart.

What the chart cannot do. For finite ρ every p_j is strictly positive, so the inequality multipliers ν never become active. A maximiser on a face of the simplex lies at infinity in ρ: the iterates drift, ‖∇_ρℓ‖ decays only through the factor diag (p) in Eq. 26, and the Hessian becomes singular. The boundary is handled separately (Sec. 3.11). The benchmark reports the projected-gradient residual $\lVert\Pi_{\Delta}(\hat{\boldsymbol{p}}+\nabla_{\boldsymbol{p}}\ell)-\hat{\boldsymbol{p}}\rVert$, which is zero exactly at a KKT point, interior or not.

Perspectives: Bayesian charts and the simplex boundary

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)_jℓ = μ_⋅j^⊤Ω(p)μ_⋅ℓ + 2p_jp_ℓ tr(Ω(p)Σ_jΩ(p)Σ_ℓ).   (28)

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. 28.

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_ψ = V^TS(p)I(p)S(p)V.   (29)

The mean/covariance split of Eq. 28 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(30)$$

Diagonal square roots are asymptotic standard errors; Wald intervals at level 1 − α are p̂_j ± z_{1 − α/2} SE_j with z_q = Φ⁻¹(q), as in confint.decovart_fit(). vcov_ilr_delta() implements Eq. 29–Eq. 30. Orthogonal rotations of V leave Eq. 30 unchanged. The construction is undefined on the simplex boundary (the log-ratio chart blows up); the helper then returns NA with a warning.

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,
  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(ε⁻¹) 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. 28 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) = ∑_jp_j²Σ_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.

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(Gb²) (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 = D_j + UC_jU^T — the factor-analysis form Σ = WW^T + Ψ. The Woodbury identity and the matrix-determinant lemma (Hager 1989) move the O(G³) 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 G³/3 flops before the cheaper triangular solves (Golub and Van Loan 2013). The structured backends reduce that leading term:

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 ( ≈ B² and O(G²) → O(Gr²) 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.

References