--- title: "Derivatives of the DeCovarT generative model under simplex transforms" author: "Bastien Chassagnol" date: "`r Sys.Date()`" vignette: > %\VignetteIndexEntry{Derivatives of the DeCovarT generative model under simplex transforms} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} knitr: opts_chunk: collapse: true comment: "#>" eval: false bibliography: - ../inst/packages.bib - ../inst/REFERENCES.bib link-citations: true crossref: eq-prefix: "Eq." tbl-prefix: "Table" lst-prefix: "Listing" sec-prefix: "Sec." format: html: embed-resources: true toc: true --- ```{r} #| label: setup #| include: false #| eval: true # CRAN vignettes use Quarto Markdown tables (no tinytable) so vignette # rebuilds do not depend on Suggests packages that change S4 APIs. # Chunks that call DeCovarT stay at the YAML eval: false default: the # Quarto CLI on Windows cannot see R CMD build's temporary library # (https://github.com/quarto-dev/quarto-r/issues/217). Live checks are # in tests/testthat/. NULL ``` DeCovarT maximises a Gaussian-convolution log-likelihood for cellular ratios $\boldsymbol{p}\in\Delta^{J-1}$. This vignette derives the analytic first- and second-order derivatives used by the frequentist solvers, first with respect to unconstrained ambient proportions $\boldsymbol{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()`](../reference/fit_decovart.html) and [`confint.decovart_fit()`](../reference/fit_decovart.html). Package entry points are [`gradient_loglik_unconstrained()`](../reference/gradient_loglik_unconstrained.html), [`hessian_loglik_unconstrained()`](../reference/hessian_loglik_unconstrained.html), [`gradient_loglik_constrained()`](../reference/gradient_loglik_constrained.html), [`hessian_loglik_constrained()`](../reference/hessian_loglik_constrained.html), [`expected_fisher_unconstrained()`](../reference/expected_fisher_unconstrained.html) and [`vcov_ilr_delta()`](../reference/vcov_ilr_delta.html). The additive log-ratio chart is retained as an appendix (@sec-alr). ::: {#nte-notation .callout-note} ## Notation - $\boldsymbol{\mu}\in\mathbb{R}^{G\times J}$: mean signature. - $\boldsymbol{\Sigma}_j\in\mathrm{SPD}(G)$: cell-type covariance; $\boldsymbol{\Sigma}(\boldsymbol{p})=\sum_j p_j^{2}\boldsymbol{\Sigma}_j$ and $\boldsymbol{\Omega}(\boldsymbol{p})=\boldsymbol{\Sigma}(\boldsymbol{p})^{-1}$. - $\boldsymbol{p}$: proportions on the open simplex (ambient score equations treat them as free coordinates in $(0,1)^{J}$). - $\boldsymbol{\rho}\in\mathbb{R}^{J-1}$: unconstrained log-ratio coordinates (ILR / Helmert chart $\mathbf{V}$); $\boldsymbol{p}=\operatorname{softmax}(\mathbf{V}\boldsymbol{\rho})$. The manuscript DAG uses the same letter for the ALR chart; here $\boldsymbol{\rho}$ denotes the ILR coordinates used by the solvers. ::: # Matrix calculus reminders {#sec-matrix-calculus} 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. ::: {#rem-first-order-matrix} ## First-order matrix calculus Let $\boldsymbol{A}=\boldsymbol{A}(p)$ and $\boldsymbol{B}=\boldsymbol{B}(p)$ be invertible matrix-valued functions of a scalar $p$, and let $\boldsymbol{U}$, $\boldsymbol{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} $$ {#eq-matrix-first-order} 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). $$ {#eq-logdet} For a residual quadratic form with symmetric precision $\boldsymbol{\Omega}$ 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). $$ {#eq-quad-form-deriv} ::: ::: {#rem-second-order-matrix} ## 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}, $$ {#eq-matrix-second-order} and the trace is linear, $\partial\mathrm{tr}(\boldsymbol{A})/\partial p_{i} =\mathrm{tr}(\partial\boldsymbol{A}/\partial p_{i})$. Combining @eq-matrix-first-order 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]. $$ {#eq-logdet-hess} ::: ::: {#nte-inner-products .callout-tip} ## Inner products in the score Several terms are $\boldsymbol{\Omega}$-inner products of signature columns, $\boldsymbol{\mu}_{\cdot j}^{\top}\boldsymbol{\Omega}\boldsymbol{\mu}_{\cdot k}$. Package helper [`.inner_product()`](../reference/dot-inner_product.html) evaluates that bilinear form; see also the squared Mahalanobis distance helper in the same family. ::: # Unconstrained score and Hessian {#sec-unconstrained} Conditional on the reference moments $\boldsymbol{\zeta}=(\boldsymbol{\mu},\{\boldsymbol{\Sigma}_{j}\})$, the bulk follows $\boldsymbol{y}\mid\boldsymbol{\zeta},\boldsymbol{p} \sim\mathcal{N}_{G}\bigl(\boldsymbol{\mu}\boldsymbol{p}, \boldsymbol{\Sigma}(\boldsymbol{p})\bigr)$. Up to an additive constant independent of $\boldsymbol{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}). $$ {#eq-loglik} ::: {#nte-halfdet .callout-important} ## The determinant term carries a factor $\tfrac{1}{2}$ Both terms of @eq-loglik are halved, because $\log\det\boldsymbol{\Omega}=-\log\det\boldsymbol{\Sigma}$ enters the Gaussian log-density as $-\tfrac{1}{2}\log\det\boldsymbol{\Sigma}(\boldsymbol{p})$. DeCovarT releases before 2.3.0 implemented $-\log\det\boldsymbol{\Sigma}(\boldsymbol{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-fisher-wald, 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 $\mathbb{E}[-\mathbf{H}]=I(\boldsymbol{p})$ holds exactly (see @nte-obs-vs-exp). 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. ::: ::: {#nte-unconstrained-meaning .callout-warning} ## What “unconstrained” means here @sec-unconstrained treats each $p_{j}$ as an ambient coordinate in $(0,1)$ and ignores $\sum_{j}p_{j}=1$. The resulting gradient and Hessian are the building blocks of the ILR chain rule in @sec-constrained. They are **not** the coordinates used by Marquardt–Levenberg or Newton–Raphson in the package (those optimise in $\boldsymbol{\rho}$-space). The box-constrained L-BFGS-B path is the exception: it works directly in $\boldsymbol{p}$ with $[0,1]^{J}$ boxes and does not enforce the simplex exactly. ::: Applying @eq-logdet and @eq-quad-form-deriv with $\partial\boldsymbol{\Sigma}/\partial p_{j}=2p_{j}\boldsymbol{\Sigma}_{j}$ gives the ambient score. Colour coding matches the manuscript: purple determinant / precision terms, blue mean residual, orange covariance quadratic. ::: {#thm-score-unconstrained} ## 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} $$ {#eq-grad-unconstrained} ::: Differentiating @eq-grad-unconstrained once more yields the ambient Hessian $\mathbf{H}\in\mathcal{M}_{J\times J}$. Diagonal and off-diagonal blocks keep the same colour pairing. ::: {#thm-hess-unconstrained} ## Unconstrained Hessian For $i\in\{1,\ldots,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} $$ {#eq-hess-diag} and for $i\neq 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} $$ {#eq-hess-off} ::: The stationary equation $\nabla\ell(\boldsymbol{p})=\boldsymbol{0}$ has no closed root, so DeCovarT uses iterative solvers. Analytic gradients and Hessians are checked against Richardson extrapolation in `numDeriv` [@R-numDeriv]. # Constrained optimisation on the simplex {#sec-constrained} 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 $C^{2}$ diffeomorphism from free Euclidean coordinates $\boldsymbol{\rho}\in\mathbb{R}^{J-1}$ onto the open simplex, so every Newton step stays feasible by construction. ## Log-sum-exp and softmax {#sec-lse} The softmax of scores $\boldsymbol{z}\in\mathbb{R}^{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}. $$ {#eq-softmax} Its Jacobian is $\mathbf{J}_{\mathrm{softmax}}=\operatorname{diag}(\boldsymbol{s}) -\boldsymbol{s}\boldsymbol{s}^{\top}$, which is also $\nabla^{2}\operatorname{LSE}(\boldsymbol{z})$. The second derivative is the third-order tensor $$ H_{ijk} = s_{i}\bigl[ (\delta_{ij}-s_{j})(\delta_{ik}-s_{k}) - s_{j}(\delta_{jk}-s_{k}) \bigr]. $$ {#eq-softmax-hess} ## Isometric log-ratio map {#sec-ilr} Choose a Helmert sub-matrix $\mathbf{V}\in\mathbb{R}^{J\times(J-1)}$ with $$ \mathbf{V}^{\mathsf{T}}\mathbf{V}=\mathbf{I}_{J-1}, \qquad \mathbf{V}^{\mathsf{T}}\mathbf{1}=\mathbf{0}. $$ {#eq-helmert} ILR coordinates and their inverse are [@pawlowsky-glahnCompositionalDataAnalysis2011] $$ \boldsymbol{\rho} = \mathbf{V}^{\mathsf{T}}\log\boldsymbol{p} = \mathbf{V}^{\mathsf{T}}\operatorname{clr}(\boldsymbol{p}), \qquad \boldsymbol{p} = \operatorname{softmax}(\mathbf{V}\boldsymbol{\rho}) = \mathcal{C}\bigl\{\exp(\mathbf{V}\boldsymbol{\rho})\bigr\}, $$ {#eq-ilr-maps} where $\mathcal{C}$ denotes closure to unit sum. No cell type is pinned as a reference (contrast @sec-alr). Package helpers [`helmert_basis()`](../reference/helmert_basis.html), [`isometric_logistic()`](../reference/isometric_logistic.html) and [`isometric_log_ratio()`](../reference/isometric_log_ratio.html) implement @eq-helmert--@eq-ilr-maps with a log-sum-exp softmax, matching the `compositions::ilr()` / `ilrInv()` convention on a Helmert basis [@R-compositions; @egozcueIsometricLogratioTransformations2003]. ::: {#prp-ilr-diffeo} ## ILR is a $C^{2}$ diffeomorphism $\boldsymbol{\psi}:\mathbb{R}^{J-1}\to(0,1)^{J}$ with $\sum_{j}p_{j}=1$ is bijective and twice continuously differentiable, so every critical point in $\boldsymbol{\rho}$-space maps to a unique open simplex composition and conversely. Any other valid ILR basis is $\mathbf{V}^{\star}=\mathbf{V}Q$ with $Q^{\mathsf{T}}Q=\mathbf{I}$, which rotates coordinates as $\boldsymbol{\rho}^{\star}=Q^{\mathsf{T}} \boldsymbol{\rho}$ and leaves eigenvalues, condition numbers and traces of properly transformed quadratic forms unchanged. ::: ::: {#nte-ilr-vs-alr .callout-note} ## 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-alr for comparison and for the structural-zero discussion (@sec-structural-zeros). ::: ## Efficient derivatives of $\boldsymbol{\psi}$ {#sec-ilr-deriv} Write $\mathbf{S}(\boldsymbol{p})=\operatorname{diag}(\boldsymbol{p}) -\boldsymbol{p}\boldsymbol{p}^{\mathsf{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}. $$ {#eq-ilr-jac} [`jacobian_isometric_logistic()`](../reference/jacobian_isometric_logistic.html) returns this $J\times(J-1)$ matrix. Let $\bar{\boldsymbol{v}}=\mathbf{V}^{\mathsf{T}}\boldsymbol{p}$, $\boldsymbol{q}_{i}=\mathbf{V}_{i\cdot}-\bar{\boldsymbol{v}}$ and $\mathbf{C}_{V}=\mathbf{V}^{\mathsf{T}}\mathbf{S}(\boldsymbol{p}) \mathbf{V}$. Slice $i$ of the map Hessian is $$ \mathbf{H}_{\psi_{i}} = p_{i}\bigl( \boldsymbol{q}_{i}\boldsymbol{q}_{i}^{\mathsf{T}} -\mathbf{C}_{V} \bigr). $$ {#eq-ilr-hess} [`hessian_isometric_logistic()`](../reference/hessian_isometric_logistic.html) stores the tensor of shape $(J-1)\times(J-1)\times J$. ## Change of variables: ambient $\boldsymbol{p}$ to ILR $\boldsymbol{\rho}$ {#sec-chain-rule} Let $\ell_{\boldsymbol{\rho}}(\boldsymbol{\rho}) =\ell_{\boldsymbol{y}\mid\boldsymbol{\zeta}} \bigl(\boldsymbol{\psi}(\boldsymbol{\rho})\bigr)$. The chain rule maps the unconstrained derivatives of @sec-unconstrained into free Euclidean coordinates. ::: {#thm-chain-rule} ## Constrained score and Hessian $$ \nabla_{\boldsymbol{\rho}}\ell = \mathbf{J}_{\boldsymbol{\psi}}^{\mathsf{T}}\nabla_{\boldsymbol{p}}\ell = \mathbf{V}^{\mathsf{T}}\mathbf{S}(\boldsymbol{p})\nabla_{\boldsymbol{p}}\ell, $$ {#eq-chain-first} $$ \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}}. $$ {#eq-chain-second} ::: At an interior KKT point every active score equals the same multiplier, $\partial\ell/\partial p_{i}=\lambda$. Because $\sum_{i}p_{i}(\boldsymbol{\rho})=1$ for every $\boldsymbol{\rho}$, $\sum_{i}\mathbf{H}_{\psi_{i}}=\mathbf{0}$, so the second summand of @eq-chain-second vanishes and $\mathbf{H}_{\boldsymbol{\rho}} =\mathbf{J}_{\boldsymbol{\psi}}^{\mathsf{T}} \mathbf{H}_{\boldsymbol{p}} \mathbf{J}_{\boldsymbol{\psi}}$. Away from stationarity that contraction must be kept. ::: {#nte-chain-impl .callout-important} ## 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-chain-second is the sum of a covariant pull-back of $\mathbf{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-glahnCompositionalDataAnalysis2011]: a Dirichlet or logistic-normal prior on $\boldsymbol{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](theory-decovart-statistical-perspectives.html#sec-map-blup) 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 {#sec-fisher-wald} 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 {#sec-fisher-p} For $\boldsymbol{y}\sim\mathcal{N}_{G}\bigl(\boldsymbol{\mu}\boldsymbol{p}, \boldsymbol{\Sigma}(\boldsymbol{p})\bigr)$ the expected Fisher information has entries $$ I(\boldsymbol{p})_{jk} = \boldsymbol{\mu}_{\cdot j}^{\top} \boldsymbol{\Omega}(\boldsymbol{p}) \boldsymbol{\mu}_{\cdot k} + 2p_{j}p_{k}\, \mathrm{tr}\bigl( \boldsymbol{\Omega}(\boldsymbol{p})\boldsymbol{\Sigma}_{j} \boldsymbol{\Omega}(\boldsymbol{p})\boldsymbol{\Sigma}_{k} \bigr). $$ {#eq-fisher-p} The first summand is the mean contribution (a $\boldsymbol{\Omega}$-inner product of signature columns); the second is the covariance contribution of the quadratic map $\boldsymbol{p}\mapsto\boldsymbol{\Sigma}(\boldsymbol{p})$. See the multivariate-normal formula on [Wikipedia: Fisher information](https://en.wikipedia.org/wiki/Fisher_information#Multivariate_normal_distribution). [`expected_fisher_unconstrained()`](../reference/expected_fisher_unconstrained.html) implements @eq-fisher-p. ::: {#nte-obs-vs-exp .callout-note} ## Observed versus expected information The analytic Hessian of @sec-unconstrained is the *observed* information (up to sign). Wald intervals in `decovart_fit` use the *expected* information @eq-fisher-p, which does not depend on the residual $\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p}$ and is positive definite under the usual full-rank assumptions on $(\boldsymbol{\mu},\{\boldsymbol{\Sigma}_{j}\})$. The two are consistent by construction. Substituting $\mathbb{E}[\boldsymbol{r}]=\boldsymbol{0}$ and $\mathbb{E}[\boldsymbol{r}\boldsymbol{r}^{\top}] =\boldsymbol{\Sigma}(\boldsymbol{p})$ into @eq-hess-diag and @eq-hess-off cancels the determinant traces against the residual traces and returns @eq-fisher-p exactly, i.e. $\mathbb{E}[-\mathbf{H}]=I(\boldsymbol{p})$. That identity only holds with the factor $\tfrac{1}{2}$ of @nte-halfdet in place. ::: ## Constrained case: ILR delta method {#sec-fisher-rho} Fisher information transforms as a covariant quadratic form under the ILR chart. With $\mathbf{J}_{\boldsymbol{\psi}} =\mathbf{S}(\boldsymbol{p})\mathbf{V}$, $$ I_{\boldsymbol{\rho}} = \mathbf{J}_{\boldsymbol{\psi}}^{\mathsf{T}} I(\boldsymbol{p}) \mathbf{J}_{\boldsymbol{\psi}} = \mathbf{V}^{\mathsf{T}} \mathbf{S}(\boldsymbol{p}) I(\boldsymbol{p}) \mathbf{S}(\boldsymbol{p}) \mathbf{V}. $$ {#eq-fisher-rho} The mean/covariance split of @eq-fisher-p pulls back separately, $I_{\boldsymbol{\rho}}=I_{\boldsymbol{\rho},\mathrm{mean}} +I_{\boldsymbol{\rho},\mathrm{cov}}$. Under regularity, $\hat{\boldsymbol{\rho}} \overset{a}{\sim} \mathcal{N}(\boldsymbol{\rho}_{0},I_{\boldsymbol{\rho}}^{-1})$ [@vaartAsymptoticStatistics2000]. The first-order delta method [@oehlertNoteDeltaMethod1992] then returns the simplex covariance used by [`vcov.decovart_fit()`](../reference/fit_decovart.html), $$ \mathrm{Var}(\hat{\boldsymbol{p}}) \approx \mathbf{J}_{\boldsymbol{\psi}} I_{\boldsymbol{\rho}}^{-1} \mathbf{J}_{\boldsymbol{\psi}}^{\mathsf{T}}. $$ {#eq-delta-p} Diagonal square roots are asymptotic standard errors; Wald intervals at level $1-\alpha$ are $\hat{p}_{j}\pm z_{1-\alpha/2}\,\mathrm{SE}_{j}$ with $z_{q}=\Phi^{-1}(q)$, as in [`confint.decovart_fit()`](../reference/fit_decovart.html). [`vcov_ilr_delta()`](../reference/vcov_ilr_delta.html) implements @eq-fisher-rho--@eq-delta-p. Orthogonal rotations of $\mathbf{V}$ leave @eq-delta-p unchanged. The construction is undefined on the simplex boundary (the log-ratio chart blows up); the helper then returns `NA` with a warning. ::: {#nte-wald-limits .callout-warning} ## Wald is the cheapest, not the safest, interval @eq-delta-p 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](theory-DeCovarT-MLE-properties.html) 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. ::: ::: {.callout-warning} ## Numerical speed-ups and solver safeguards {#sec-numerical-speedups} 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](fig03-covariance-driven.html#sec-hybrid-deconvolution) live in `R/03_03_DeCovarT_estimate_ratios_frequentist.R`. --- ### Cache a Cholesky factorisation of $\boldsymbol{\Sigma}(\boldsymbol{p})$ Within one iteration, objective / gradient / Hessian callbacks hit the same trial $\boldsymbol{p}$. The helper `.sigma_p_factorisation()` caches a single Cholesky factor $\boldsymbol{R}$ with $\boldsymbol{\Sigma}(\boldsymbol{p})=\boldsymbol{R}^{\mathsf{T}}\boldsymbol{R}$ and returns $\log\det\boldsymbol{\Sigma}(\boldsymbol{p})=2\sum_g\log R_{gg}$ together with the precision $\boldsymbol{\Omega}(\boldsymbol{p})=\boldsymbol{R}^{-1}\boldsymbol{R}^{-\mathsf{T}}$ (via `chol2inv`) without repeating an $O(G^{3})$ factorisation. [`loglik_multivariate()`](../reference/loglik_multivariate.html) then evaluates the Mahalanobis term exactly as `mvtnorm::dmvnorm` [@R-mvtnorm]: it solves $\boldsymbol{R}^{\mathsf{T}}\boldsymbol{z} =\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p}$ by `backsolve(..., transpose = TRUE)` and takes $\lVert\boldsymbol{z}\rVert^{2}$, rather than forming $\boldsymbol{\Omega}(\boldsymbol{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 $\boldsymbol{\Omega}(\boldsymbol{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(G^{3})$ 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 $\boldsymbol{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 $\boldsymbol{p}$ with boxes $[0,1]^{J}$ that do not enforce $\mathbf{1}^{\top}\boldsymbol{p}=1$. Near $\sum_{j}p_{j}\approx 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`](https://cran.r-project.org/package=marqLevAlg) [@R-marqLevAlg] implements Marquardt–Levenberg with the relative-distance-to-minimum (RDM) rule of Commenges et al. [@commengesNewtonLikeAlgorithmLikelihood2006]. 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](https://github.com/VivianePhilipps/marqLevAlgParallel/issues/3). ::: # Model objects and S3 accessors {#sec-s3} [`fit_decovart()`](../reference/fit_decovart.html) 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: ```{r} #| label: s3-accessors #| warning: false 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 {#sec-numerical-consistency} 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. ```{r} #| label: consistency-setup #| include: false library(DeCovarT) genes <- paste0("g", 1:3) cts <- paste0("ct", 1:3) mu_chk <- matrix( c(20, 40, 15, 40, 20, 25, 25, 30, 35), nrow = 3, dimnames = list(genes, cts) ) Sigma_chk <- array(0, dim = c(3, 3, 3), dimnames = list(genes, genes, cts)) Sigma_chk[, , 1] <- diag(c(1, 1.5, 1.2)) Sigma_chk[, , 2] <- diag(c(1.1, 0.9, 1.3)) Sigma_chk[, , 3] <- diag(c(0.8, 1.2, 1.0)) p_chk <- c(0.5, 0.3, 0.2) set.seed(1) y_chk <- drop(mu_chk %*% p_chk) + c(0.4, -0.2, 0.1) cov_chk <- .compute_global_variance(p_chk, Sigma_chk) ``` ## 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` [@R-mvtnorm] differ by exactly $\tfrac{G}{2}\log(2\pi)$ (DeCovarT omits the negative constant): ```{r} #| label: consistency-dmvnorm 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` [@R-numDeriv]. 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(\varepsilon^{-1})$ cost in function evaluations. ```{r} #| label: consistency-numderiv 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 \(\mathbb{E}[-\mathbf{H}]\) Matching `numDeriv` shows that the derivatives belong to the implemented objective. Matching the *expected* Fisher information [@eq-fisher-p] to a Monte Carlo average of \(-\mathbf{H}(\boldsymbol{p}; \boldsymbol{Y})\) shows that the objective is the Gaussian whose information matrix we use for Wald intervals. Draw \(\boldsymbol{Y}^{(b)}\sim\mathcal{N}_G(\boldsymbol{\mu}\boldsymbol{p}, \boldsymbol{\Sigma}(\boldsymbol{p}))\) at the *true* \(\boldsymbol{p}\) (the parametric bootstrap of [`bootstrap_decovart()`](../reference/bootstrap_decovart.html), used here as a diagnostic rather than as an inferential procedure): ```{r} #| label: consistency-fisher-mc 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 \(-\log\det\) objective, \(\mathbb{E}[-\mathbf{H}]\) could not equal \(I(\boldsymbol{p})\). 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()`](../reference/vcov_ilr_delta.html). That comparison is noisy at small replication (see [MLE properties](theory-DeCovarT-MLE-properties.html#sec-replication) for why one bulk column is not an asymptotic sample) and is therefore left to the inference vignette and to [`bootstrap_decovart()`](../reference/bootstrap_decovart.html). # Structure-aware covariance backends {#sec-structure-backends} Every solver call factorises $\boldsymbol{\Sigma}(\boldsymbol{p})=\sum_j p_j^{2}\boldsymbol{\Sigma}_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 [@boydConvexOptimization2004, §9.7]; the backends below make that structure explicit for DeCovarT. [`new_decovart_covariance()`](../reference/new_decovart_covariance.html) wraps the covariance array with a declared `structure`, and the operators [`sigma_logdet()`](../reference/sigma_logdet.html), [`sigma_solve()`](../reference/sigma_solve.html), [`sigma_quadform()`](../reference/sigma_quadform.html) and [`sigma_trace_precision_times()`](../reference/sigma_trace_precision_times.html) 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 {#sec-backend-decision} 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()`](../reference/covariance_structure_from_graph_model.html) maps a network topology to the recommended backend. A sparse *precision* $\boldsymbol{\Omega}_j$ does not imply a sparse *covariance* $\boldsymbol{\Sigma}_j$; only structure that survives in $\boldsymbol{\Sigma}(\boldsymbol{p})$ itself accelerates the likelihood. ```{mermaid} %%| label: fig-backend-decision %%| fig-cap: "Choosing a covariance backend; dense Cholesky is the default and structured backends are exact refinements." %%{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)"] ``` Each refinement is exact, not approximate: - **Block** (disconnected gene modules, e.g. a stochastic block model with no between-block edges [@hollandStochasticBlockmodelsFirst1983]). $\boldsymbol{\Sigma}(\boldsymbol{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(G b^{2})$ [@golubMatrixComputations2013]. - **Sparse** (Erdős–Rényi, scale-free [@barabasiEmergenceScalingRandom1999], small-world [@wattsCollectiveDynamicsSmallworld1998]). The fill-reducing ordering is a *symbolic factorisation* computed once and reused for every numeric refactor as $\boldsymbol{p}$ changes; only the numerical Cholesky repeats [@chenAlgorithm887CHOLMOD2008]. - **Diagonal + low rank.** Genes interact through $r\ll G$ shared regulatory programs, $\boldsymbol{\Sigma}_j=\mathbf{D}_j+\mathbf{U}\mathbf{C}_j \mathbf{U}^{\mathsf{T}}$ — the factor-analysis form $\boldsymbol{\Sigma}=\mathbf{W}\mathbf{W}^{\mathsf{T}}+\boldsymbol{\Psi}$. The Woodbury identity and the matrix-determinant lemma [@hagerUpdatingInverseMatrix1989] move the $O(G^{3})$ 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. ```{r} #| label: backend-usage #| eval: false #| echo: true # 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 {#sec-backend-complexity} For a dense symmetric positive-definite $G\times G$ matrix the Cholesky factorisation costs about $G^{3}/3$ flops before the cheaper triangular solves [@golubMatrixComputations2013]. The structured backends reduce that leading term: | Structure | When it applies | Factorisation cost | Backend argument | |:----------|:----------------|:-------------------|:-----------------| | Dense | universal default | $O(G^{3})$ | (default) | | Block | disconnected modules (SBM) | $O(\sum_b G_b^{3})$; $\approx B^{2}$ saving for $B$ equal blocks | `structure = "block"` | | Band / AR | bandwidth $b$ | $O(G b^{2})$ | `structure = "band"` | | Sparse | fixed sparse pattern | pattern-dependent; symbolic ordering reused across $\boldsymbol{p}$ | `structure = "sparse"` | | Diagonal + low rank | $r$ shared programs | $O(G r^{2})$ solve, $O(r^{3})$ core | `structure = "diag_lowrank"` | : Leading factorisation cost per backend. $G$ genes, $J$ cell types, $B$ blocks of size $G_b$, bandwidth $b$, rank $r\ll G$. {#tbl-backend-complexity} 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 ($\approx B^{2}$ and $O(G^{2})\!\to\!O(Gr^{2})$ respectively) are the most robust. The `"band"` and `"sparse"` backends use the imported `Matrix` package (supernodal sparse Cholesky in CHOLMOD, [@chenAlgorithm887CHOLMOD2008]). `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äferCompressionInversionApproximate2020]; the `"sparse"` backend here is exact for a declared sparse pattern, not that approximation. # Appendix: additive log-ratio chart {#sec-alr} Pin the last category as reference ($\rho_{J}\equiv 0$). The *additive logistic* transform $\boldsymbol{\psi}$ (inverse ALR) and its inverse are [@aitchisonStatisticalAnalysisCompositional1982] $$ p_{j} = \frac{e^{\rho_{j}}}{\sum_{k