| Expression | Terminology |
|---|---|
| \(\boldsymbol{x}^{\mathsf{T}}\boldsymbol{y}\) | Dot product / Euclidean inner product |
| \(\boldsymbol{x}^{\mathsf{T}}\boldsymbol{A}\boldsymbol{y}\) | Bilinear form (inner product if \(\boldsymbol{A}\) is SPD) |
| \(\boldsymbol{x}^{\mathsf{T}}\boldsymbol{A}\boldsymbol{x}\) | Quadratic form |
| \(\boldsymbol{x}^{\mathsf{T}}\boldsymbol{A}^{-1}\boldsymbol{x}\) | Quadratic form induced by \(\boldsymbol{A}^{-1}\) |
| \((\boldsymbol{x}-\boldsymbol{\mu})^{\mathsf{T}}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\) | Squared Mahalanobis distance |
| \(\sqrt{(\boldsymbol{x}-\boldsymbol{\mu})^{\mathsf{T}}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})}\) | Mahalanobis distance |
Simplex maps and numerically stable softmax / ALR derivatives
2026-08-08
Source:vignettes/softmax-alr-derivatives.qmd
DeCovarT estimates cellular ratios \boldsymbol{p}\in\Delta^{J-1} by maximising a Gaussian-convolution log-likelihood over an unconstrained parameter \boldsymbol{\rho}\in\mathbb{R}^{J-1}. This vignette collects the simplex coordinate maps, the matrix-algebra vocabulary used in the analytic gradients, and the closed-form softmax / additive log-ratio derivatives needed for efficient optimisation. For every map we give the explicit tensor formula, a base-R implementation (checked against numDeriv (Gilbert and Varadhan 2019)), and a batched PyTorch counterpart. Naming follows the DeCovarT API: additive_logistic() for \boldsymbol{\rho}\mapsto\boldsymbol{p} and additive_log_ratio() for \boldsymbol{p}\mapsto\boldsymbol{\rho}.
Simplex coordinate maps
The forward map additive_logistic() sends unconstrained coordinates \boldsymbol{\rho} to the open simplex,
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}. \tag{1}
Eq. 1 is a softmax in which the last category J is pinned as a reference (\rho_J\equiv 0); in Aitchison’s framework it is the additive logistic transform, i.e. the inverse additive log-ratio map (\mathrm{alr}^{-1}).
The inverse map additive_log_ratio() sends a composition back to \mathbb{R}^{J-1},
\rho_j=\log\!\left(\frac{p_j}{p_J}\right),\qquad j=1,\ldots,J-1. \tag{2}
Eq. 2 is the additive log-ratio (\mathrm{alr}) transform with reference part p_J, equivalently the multinomial-logit link with reference category J.
Table 1 summarises the two helpers and their equivalents across statistical, machine-learning, and compositional-data-analysis literatures.
| DeCovarT function | Direction | Transform | Equivalent formulations |
|---|---|---|---|
additive_logistic() |
\boldsymbol{\rho}\mapsto\boldsymbol{p} (to \Delta^{J-1}) | inverse additive log-ratio (\mathrm{alr}^{-1}) | softmax with reference category J; compositions::alrInv(); skbio.stats.composition.alr_inv()
|
additive_log_ratio() |
\boldsymbol{p}\mapsto\boldsymbol{\rho} (to \mathbb{R}^{J-1}) | additive log-ratio (\mathrm{alr}) | multinomial-logit link (reference J); compositions::alr(); skbio.stats.composition.alr()
|
Matrix-induced scalar products
The DeCovarT log-likelihood and its derivatives evaluate several matrix-induced scalar products involving a symmetric positive-definite covariance or precision matrix \boldsymbol{A} (non-degenerate Gaussian laws). Table 2 records the standard terminology; a quadratic form is the bilinear form evaluated on the same vector twice, q(\boldsymbol{x})=b(\boldsymbol{x},\boldsymbol{x}).
Package helpers .bilinear_form() and .squared_mahalanobis_distance() implement the middle rows of Table 2. The second solves \boldsymbol{\Sigma}\boldsymbol{z}=\boldsymbol{\delta} rather than forming \boldsymbol{\Sigma}^{-1} explicitly. Calling \boldsymbol{x}^{\mathsf{T}}\boldsymbol{A}\boldsymbol{y} a “dot product” is misleading unless \boldsymbol{A}=\boldsymbol{I}; see also the dot product and inner product space articles for the Euclidean special case \langle\boldsymbol{x},\boldsymbol{y}\rangle=\boldsymbol{x}^{\mathsf{T}}\boldsymbol{y}.
Log-sum-exp and softmax
The softmax of a score vector \boldsymbol{z}\in\mathbb{R}^{K} is evaluated through the log-sum-exp (LSE) to avoid overflow,
s_i=\operatorname{softmax}(\boldsymbol{z})_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. \tag{3}
import torch
import torch.nn.functional as F
def softmax(z: torch.Tensor, dim: int = -1) -> torch.Tensor:
return F.softmax(z, dim=dim)
def log_sum_exp(z: torch.Tensor, dim: int = -1) -> torch.Tensor:
return torch.logsumexp(z, dim=dim)Softmax Jacobian and Hessian
The softmax Jacobian is the difference between a diagonal and a rank-one term,
\frac{\partial s_i}{\partial z_j}=s_i(\delta_{ij}-s_j), \qquad \mathbf{J}_{\mathrm{softmax}}=\operatorname{diag}(\boldsymbol{s}) -\boldsymbol{s}\boldsymbol{s}^{\top}. \tag{4}
This matrix also equals the Hessian of the log-sum-exp, \nabla^2\operatorname{LSE}(\boldsymbol{z})=\operatorname{diag}(\boldsymbol{s}) -\boldsymbol{s}\boldsymbol{s}^{\top}. Because softmax is vector-valued, its second derivative is a third-order tensor H_{ijk}=\partial^2 s_i/(\partial z_j\partial z_k),
H_{ijk}=s_i\bigl[(\delta_{ij}-s_j)(\delta_{ik}-s_k)-s_j(\delta_{jk}-s_k)\bigr], \tag{5}
which, slice by slice, is the rank-one correction
\mathbf{H}^{(i)}=s_i\Bigl[(\boldsymbol{e}_i-\boldsymbol{s}) (\boldsymbol{e}_i-\boldsymbol{s})^{\top} -\bigl(\operatorname{diag}(\boldsymbol{s})-\boldsymbol{s}\boldsymbol{s}^{\top}\bigr) \Bigr], \tag{6}
with \boldsymbol{e}_i the i-th canonical basis vector. Listing 3 materialises Eq. 4 and Eq. 6 in O(K^2) per slice.
softmax_jacobian <- function(z) {
s <- softmax(z)
diag(s) - tcrossprod(s)
}
softmax_hessian <- function(z) {
s <- softmax(z)
k <- length(s)
jac <- diag(s) - tcrossprod(s)
hessian <- array(0, dim = c(k, k, k))
for (i in seq_len(k)) {
di <- (seq_len(k) == i) - s # e_i - s
hessian[i, , ] <- s[i] * (tcrossprod(di) - jac)
}
hessian
}def softmax_jacobian(z: torch.Tensor) -> torch.Tensor:
"""z: [..., K] -> [..., K, K]."""
s = softmax(z, dim=-1)
return torch.diag_embed(s) - s.unsqueeze(-1) * s.unsqueeze(-2)
def softmax_hessian(z: torch.Tensor) -> torch.Tensor:
"""z: [K] -> H[i, j, k] of shape [K, K, K]."""
s = softmax(z, dim=0)
eye = torch.eye(s.numel(), dtype=s.dtype, device=s.device)
diff = eye - s # (e_i - s) stacked over i
outer = diff.unsqueeze(-1) * diff.unsqueeze(-2)
jac = torch.diag(s) - torch.outer(s, s)
return s[:, None, None] * (outer - jac.unsqueeze(0))Additive log-ratio coordinates
Let \boldsymbol{p}=(p_1,\ldots,p_J) be a composition with p_j>0 and \sum_j p_j=1, and take the last part J as reference. The additive log-ratio (alr) transform sends the composition to \mathbb{R}^{J-1},
\rho_i=\log\!\Bigl(\frac{p_i}{p_J}\Bigr),\qquad i=1,\ldots,J-1, \tag{7}
and its inverse, the additive logistic transform (\mathrm{alr}^{-1}), pins a reference logit at zero, \tilde{\boldsymbol{\rho}}=(\rho_1,\ldots,\rho_{J-1},0), and applies a softmax,
\boldsymbol{p}=\operatorname{softmax}(\tilde{\boldsymbol{\rho}}), \qquad p_i=\frac{e^{\rho_i}}{1+\sum_{j<J}e^{\rho_j}}\ (i<J), \qquad p_J=\frac{1}{1+\sum_{j<J}e^{\rho_j}}. \tag{8}
def additive_logistic(rho: torch.Tensor) -> torch.Tensor:
"""rho: [..., J-1] -> p: [..., J] (reference category = last)."""
reference = torch.zeros_like(rho[..., :1])
return torch.softmax(torch.cat([rho, reference], dim=-1), dim=-1)
def additive_log_ratio(p: torch.Tensor) -> torch.Tensor:
"""p: [..., J] -> rho: [..., J-1]."""
log_p = torch.log(p)
return log_p[..., :-1] - log_p[..., -1:]Derivatives of the additive logistic map (\boldsymbol{\rho}\mapsto\boldsymbol{p})
The additive logistic map is the softmax restricted to the first J-1 inputs, so its Jacobian \mathbf{J}_{\boldsymbol{\psi}}\in\mathcal{M}_{J\times(J-1)} is the softmax Jacobian with the reference column dropped,
\frac{\partial p_i}{\partial\rho_a}=p_i(\delta_{ia}-p_a), \qquad \mathbf{J}_{\boldsymbol{\psi}} =\bigl[\operatorname{diag}(\boldsymbol{p}) -\boldsymbol{p}\boldsymbol{p}^{\top}\bigr]_{:,\,1:(J-1)}, \qquad a=1,\ldots,J-1. \tag{9}
Its Hessian is a tensor \mathbf{H}_{\boldsymbol{\psi}}\in\mathcal{M}_{J\times(J-1)\times(J-1)}, the softmax Hessian Eq. 6 restricted to the first J-1 input indices,
\frac{\partial^2 p_i}{\partial\rho_a\partial\rho_b} =p_i\bigl[(\delta_{ia}-p_a)(\delta_{ib}-p_b)-p_a(\delta_{ab}-p_b)\bigr], \qquad a,b=1,\ldots,J-1. \tag{10}
additive_logistic_jacobian <- function(rho) {
p <- additive_logistic(rho)
j <- length(p)
(diag(p) - tcrossprod(p))[, -j, drop = FALSE]
}
additive_logistic_hessian <- function(rho) {
p <- additive_logistic(rho)
j <- length(p)
jac <- diag(p) - tcrossprod(p)
hessian <- array(0, dim = c(j, j - 1, j - 1))
for (i in seq_len(j)) {
di <- (seq_len(j) == i) - p
slice <- p[i] * (tcrossprod(di) - jac)
hessian[i, , ] <- slice[-j, -j]
}
hessian
}def additive_logistic_jacobian(rho: torch.Tensor) -> torch.Tensor:
"""rho: [J-1] -> Jacobian [J, J-1]."""
p = additive_logistic(rho)
return (torch.diag(p) - torch.outer(p, p))[:, :-1]
def additive_logistic_hessian(rho: torch.Tensor) -> torch.Tensor:
"""rho: [J-1] -> Hessian [J, J-1, J-1]."""
reference = torch.zeros_like(rho[:1])
augmented = torch.cat([rho, reference])
return softmax_hessian(augmented)[:, :-1, :-1]Derivatives of the additive log-ratio map (\boldsymbol{p}\mapsto\boldsymbol{\rho})
Treating \boldsymbol{p} as ambient coordinates in \mathbb{R}^{J} (admissible tangent directions satisfy \sum_i \mathrm{d}p_i=0), the alr Jacobian \mathbf{J}_{\mathrm{alr}}\in\mathcal{M}_{(J-1)\times J} is sparse,
\frac{\partial\rho_i}{\partial p_j}=\frac{\delta_{ij}}{p_i}-\frac{\delta_{jJ}}{p_J}, \qquad \mathbf{J}_{\mathrm{alr}} =\Bigl[\operatorname{diag}\!\bigl(p_1^{-1},\ldots,p_{J-1}^{-1}\bigr) \ \big|\ -p_J^{-1}\boldsymbol{1}_{J-1}\Bigr], \tag{11}
and its Hessian \mathbf{H}_{\mathrm{alr}}\in\mathcal{M}_{(J-1)\times J\times J} is diagonal in the last two modes,
\frac{\partial^2\rho_i}{\partial p_j\partial p_k} =-\frac{\delta_{ij}\delta_{ik}}{p_i^{2}} +\frac{\delta_{jJ}\delta_{kJ}}{p_J^{2}}. \tag{12}
additive_log_ratio_jacobian <- function(p) {
j <- length(p)
idx <- seq_len(j - 1)
jac <- matrix(0, j - 1, j)
jac[cbind(idx, idx)] <- 1 / p[idx]
jac[, j] <- -1 / p[j]
jac
}
additive_log_ratio_hessian <- function(p) {
j <- length(p)
hessian <- array(0, dim = c(j - 1, j, j))
for (i in seq_len(j - 1)) {
hessian[i, i, i] <- -1 / p[i]^2
hessian[i, j, j] <- 1 / p[j]^2
}
hessian
}def additive_log_ratio_jacobian(p: torch.Tensor) -> torch.Tensor:
"""p: [J] -> Jacobian [J-1, J]."""
j = p.numel()
jac = torch.zeros(j - 1, j, dtype=p.dtype, device=p.device)
idx = torch.arange(j - 1, device=p.device)
jac[idx, idx] = 1.0 / p[:-1]
jac[:, -1] = -1.0 / p[-1]
return jac
def additive_log_ratio_hessian(p: torch.Tensor) -> torch.Tensor:
"""p: [J] -> Hessian [J-1, J, J]."""
j = p.numel()
hessian = torch.zeros(j - 1, j, j, dtype=p.dtype, device=p.device)
idx = torch.arange(j - 1, device=p.device)
hessian[idx, idx, idx] = -1.0 / p[:-1].square()
hessian[:, -1, -1] = 1.0 / p[-1].square()
return hessianThe compositions package and closed-form checks
The compositions package (van den Boogaart et al. 2025) provides a general toolbox for Aitchison geometry on the simplex. Beyond the additive log-ratio pair alr() / alrInv() used here, it implements the centred log-ratio clr() / clrInv() and the isometric log-ratio ilr() / ilrInv(), together with composition classes (acomp(), rcomp()) and the perturbation and powering operations of the Aitchison simplex. DeCovarT relies only on the additive log-ratio pair because it yields the sparsest (J-1)-dimensional parametrisation with an interpretable reference category, but any of the alternative bases would define an equally valid unconstrained coordinate system. DeCovarT’s additive_logistic() and additive_log_ratio() match alrInv() / alr() with the last part fixed as reference (Listing 11).
p <- c(0.2, 0.3, 0.5)
rho <- additive_log_ratio(p)
if (requireNamespace("compositions", quietly = TRUE)) {
rho_pkg <- as.numeric(compositions::alr(compositions::acomp(p)))
cat("alr agreement:", isTRUE(all.equal(as.numeric(rho), rho_pkg)), "\n")
p_pkg <- as.numeric(compositions::alrInv(rho_pkg))
cat("alrInv round-trip:", isTRUE(all.equal(p, p_pkg)), "\n")
}The analytic Jacobians and Hessians are validated against Richardson extrapolation in R (Listing 12).
rho <- c(0.4, -1.1)
p <- c(0.2, 0.3, 0.5)
# additive logistic: Jacobian and per-output Hessian
jac_ok <- all.equal(
numDeriv::jacobian(additive_logistic, rho),
additive_logistic_jacobian(rho)
)
hess_ok <- all.equal(
numDeriv::hessian(function(r) additive_logistic(r)[1], rho),
additive_logistic_hessian(rho)[1, , ]
)
# additive log-ratio: ambient Jacobian
alr_jac_ok <- all.equal(
numDeriv::jacobian(additive_log_ratio, p),
additive_log_ratio_jacobian(p)
)
c(logistic_jacobian = isTRUE(jac_ok),
logistic_hessian = isTRUE(hess_ok),
alr_jacobian = isTRUE(alr_jac_ok))Numerical speed-ups and solver safeguards
The analytic maps above are only half of a usable optimiser: each Newton-type iteration also evaluates the Gaussian-convolution log-likelihood and its derivatives with respect to \boldsymbol{p} (or \boldsymbol{\rho}). Several practical bottlenecks and numerical inconsistencies showed up when running the hybrid scenario in Deconvolution use cases; the changes below live in R/03_03_DeCovarT_estimate_ratios_frequentist.R.
Raise the Newton–Raphson evaluation budget
deconvolute_ratios_Newton_Raphson() wraps stats::nlminb(). An earlier control list set eval.max = 1, which caps the total number of objective evaluations for the whole run (not per iteration). The solver therefore evaluated the log-likelihood once at the equi-balanced start and stopped, returning the untouched initial guess on every sample. Removing that entry (keeping iter.max, rel.tol, ) restores genuine Newton steps; see the before/after contrast (Newton_Raphson versus Newton_Raphson_fixed) in Table solver diagnostics.
Cache a Cholesky factorisation of \boldsymbol{\Sigma}(\boldsymbol{p})
Within one iteration, optim() / nlminb() / marqLevAlg() treat the log-likelihood, gradient and Hessian as independent callbacks, yet they all hit the same trial \boldsymbol{p}. Assembling \boldsymbol{\Sigma}(\boldsymbol{p})=\sum_j p_j^2\boldsymbol{\Sigma}_j and factorising it separately in each callback (and again inside the constrained Hessian chain rule) paid for an O(G^3) factorisation up to four times per iteration. The internal helper .sigma_p_factorisation() caches a single Cholesky factor keyed on exact equality of (p, Sigma) and returns
\boldsymbol{\Sigma}(\boldsymbol{p}) = \mathbf{R}^{\mathsf{T}}\mathbf{R}, \qquad \log\det\boldsymbol{\Sigma}(\boldsymbol{p}) = 2\sum_{g=1}^{G}\log R_{gg}, \qquad \boldsymbol{\Sigma}(\boldsymbol{p})^{-1} = \mathbf{R}^{-1}(\mathbf{R}^{\mathsf{T}})^{-1} \tag{13}
via chol() / chol2inv() (Eq. 13). The unconstrained log-likelihood then uses the cached log-determinant and inverse,
\ell_{\boldsymbol{y}\,|\,\boldsymbol{\zeta}}(\boldsymbol{p}) = -\log\det\boldsymbol{\Sigma}(\boldsymbol{p}) -\tfrac{1}{2} \boldsymbol{r}^{\mathsf{T}} \boldsymbol{\Sigma}(\boldsymbol{p})^{-1} \boldsymbol{r}, \qquad \boldsymbol{r} = \boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p}, \tag{14}
and both gradient_loglik_unconstrained() and hessian_loglik_unconstrained() reuse the same \boldsymbol{\Sigma}(\boldsymbol{p})^{-1} (Eq. 14). Analytic gradients and Hessians still match numDeriv to \sim 10^{-8}–10^{-9} after the refactor.
Guard the box-constrained L-BFGS-B path
deconvolute_ratios_L_BFGS_B() optimises directly in \boldsymbol{p} with box constraints [0,1]^J. Those boxes alone do not enforce \mathbf{1}^{\mathsf{T}}\boldsymbol{p}=1, so line searches can drive \boldsymbol{\Sigma}(\boldsymbol{p}) singular. Both the objective and the analytic gradient are wrapped: near \sum_j p_j\approx 0 (or on a failed chol()), the log-likelihood returns a finite penalty and the gradient a zero vector rather than aborting optim():
safe_loglik <- function(p, y, mean_signature_matrix, Sigma) {
if (sum(p) < 1e-8) {
return(-1e12)
}
tryCatch(
loglik_multivariate(p, y, mean_signature_matrix, Sigma),
error = function(e) -1e12
)
}
safe_gradient <- function(p, y, mean_signature_matrix, Sigma) {
if (sum(p) < 1e-8) {
return(rep(0, length(p)))
}
tryCatch(
gradient_loglik_unconstrained(p, y, mean_signature_matrix, Sigma),
error = function(e) rep(0, length(p))
)
}Supply the analytic gradient to L-BFGS-B
The same solver previously relied on finite-difference gradients. Passing gr = safe_gradient (the unconstrained analytic score) removes that cost and keeps the safeguarded behaviour above.
marqLevAlg Hessian sign under minimize = FALSE
Caution 1:
marqLevAlg(minimize = FALSE)does not fliphess
marqLevAlg(Philipps et al. 2023) implements the Marquardt–Levenberg algorithm with the relative-distance-to-minimum (RDM) stopping rule of Commenges et al. (Commenges et al. 2006):C_k \approx m^{-1} \boldsymbol{U}(\boldsymbol{\theta}_k)^{\mathsf{T}} \mathbf{G}(\boldsymbol{\theta}_k)^{-1} \boldsymbol{U}(\boldsymbol{\theta}_k). \tag{15}
With
minimize = FALSE, the package only negates the user-suppliedfnandgr. The analytichessis passed through unchanged, while the internal Cholesky / inversion routines (dsinv,dchole) assume a positive-definite matrix at the optimum (minimisation convention). When maximising a log-likelihood,hessian_loglik_constrained()is negative-definite at the MLE, sodsinvfails (ier = -1) on every iteration,rdmsticks at the unevaluated sentinelepsd + 1(e.g.1.0001), and the algorithm always exhaustsmaxiter(istop = 2) even when the iterate is already at the maximum.Workaround used in DeCovarT: pass
hess = function(...) -hessian_loglik_constrained(...)and keepminimize = FALSE. Reported upstream as VivianePhilipps/marqLevAlgParallel#3.
Table 3 summarises the empirical confirmation on the hybrid J=3 scenario (m=J-1=2 free ALR coordinates). Point estimates barely moved (the wrongly scaled Newton step was often rescued by the internal line search), but iterations dropped roughly 20–40\times and istop / rdm became honest relative to Eq. 15.
| Hessian sign | istop | iterations | rdm |
|---|---|---|---|
| Wrong (log-likelihood Hessian as-is) | 2 (hit maxiter) | 200 / 200 | 1.0001 (sentinel) |
| Negated (`-hess`) | 1 (converged) | 4–11 | ~1e-13 to 1e-21 |
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.
Gilbert, Paul, and Ravi Varadhan. 2019. numDeriv: Accurate Numerical Derivatives. http://optimizer.r-forge.r-project.org/.
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.
van den Boogaart, K. Gerald, Raimon Tolosana-Delgado, and Matevz Bren. 2025. Compositions: Compositional Data Analysis. http://www.stat.boogaart.de/compositions/.