Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
82,313 characters · 17 sections · 53 citation commands
Bayesian Poisson-Randomized Gamma Tensor Factorization with Application to International Trade Flows
Dynamic multiway networks record pairwise interactions together with additional indexing dimensions, such as relation type, category, or time kivela2014multilayer. In many empirical settings, these interactions are observed as nonnegative edge weights on a sparse set of node pairs and indices, yielding a natural tensor representation DeDomenicoSoleRibaltaCozzoKivelaMorenoPorterGomezArenas2013. Such weights may be discrete counts, as in social-science and biological networks schein2015bayesian,ChacoffResascoVazquez2018, or continuous, as in economic and financial tensor time series ChenYangZhang2022TensorTS,BillioCasarinIacopiniKaufmann2023BDTR,BabiiGhyselsPan2025TensorPCA. The low-rank structure of such tensors encodes cross-mode dependence, and tensor decomposition has consequently become a standard tool for their analysis KoldaBader2009, gillis2020nonnegative. Existing probabilistic formulations for tensor-valued data have focused mainly on discrete observations, particularly categorical and count-valued data DunsonXing2009,chi2012tensors, schein2016bayesian, whereas financial and economic networks often involve monetary-valued edge weights, yielding tensors with excess zeros, heavy-tailed positive entries, and product- or transaction-type-specific variation in scale and dispersion helpman2008estimating,barigozzi2010multinetwork.
Our motivating application is international trade, where modern databases record dollar-valued flows indexed by exporting country $i$, importing country $j$, product $a$, and year $t$, denoted $y_{ijat}\in \mathbb{R}_{\geq 0}$. Countries export distinctive product baskets to particular sets of partners, and shocks propagate through shared supply chains, transport corridors, and policy regimes, generating clustering and co-movement that spans all four modes Hidalgo2007ProductSpace. The resulting tensor exhibits precisely the challenges described above: for any given product, most country pairs record no trade at all, while active pairs exchange flows ranging from hundreds to billions of dollars. Dispersion also differs markedly across products fally2018commodity. For example, petroleum oil is associated with much larger and more volatile flows than cocoa products.
This paper develops a probabilistic tensor factorization built on the nonnegative CANDECOMP/PARAFAC (CP) decomposition lee1999learning,KoldaBader2009, which represents the four-way trade tensor through the low-rank signal $y_{ijat} \approx \sum_{k=1}^{K} \theta^{(1)}_{ik}\theta^{(2)}_{jk}\theta^{(3)}_{ak}\theta^{(4)}_{tk}$ via exporter, importer, product, and time “factors”, that are all non-negative thus respecting the support of the non-negative data. As in many applications, the non-negativity of the factors yields a “parts-based” representation that is highly interpretable gillis2020nonnegative. (ref) illustrates such an interpretable structure that is surfaced by our proposed model.
Standard non-negative CP models typically assume conditionally Gaussian or Poisson entries shashua2005nonnegative, cemgil2009bayesian,gillis2020nonnegative and govern only the location parameter, like the mean, through the low-rank signal, leaving two limitations. First, they do not accommodate the mixed zero-positive support of semi-continuous data. Second, they offer no mechanism for capturing heterogeneous variability across slices when observations share similar means. We address both through a compound Poisson--Gamma construction, which is conceptually related to two-part “hurdle” models Mullahy1986 and closely connected to Tweedie compound Poisson models SmythJorgensen2002,DunnSmyth2005. Each tensor entry is driven by a latent Poisson count whose mean follows the CP decomposition, so zeros arise naturally from the probability of no events. Conditional on a positive count, the observed magnitude follows a Gamma distribution whose shape depends on the latent count. This construction can be viewed as closely related to a Tweedie model with a restricted power parameterization. Crucially, dispersion is not tied to the mean but is instead allowed to vary freely across modes, allowing the heterogeneous variability observed across product slices and other tensor margins to be captured through dedicated mode-specific parameters. Such separate dispersion modeling is new in probabilistic CP-type tensor models.
To enable posterior inference, we introduce two augmentations of the generative process that express the observed dollar-value as a sum of per-component dollar-values. Since the latent Poisson count driving each component can itself be zero, some of these per-component values are exactly zero under the model. To connect the two augmentations and update the per-component values in closed form, we introduce binary indicator variables that explicitly identify active non-zero components, and update them through a “partially collapsed” vanDykPark2008 sampler that marginalizes over variables with which the indicators have a deterministic relationship. Conditioning on these indicators, we show that all remaining parameter blocks can be updated in closed-form via coordinate ascent variational inference (CAVI) blei2017variational. The resulting approximate inference algorithm is thus a hybrid of variational and Monte Carlo methods. \looseness=-1
We apply the proposed model to a large international trade tensor constructed from a World Bank database wits2025 of size \(151 \times 151 \times 96 \times 27\), recording bilateral trade values among 151 countries as exporters and importers, across 96 high-level product categories, from 1997 to 2023. The model recovers 125 interpretable components, each identifying a coherent combination of exporters, importers, products, and temporal trends. The recovered components reveal economically meaningful structure: concentration and specialization patterns across country and product modes, opposing temporal trends that reflect market-share reallocation across competing exporters, trade-routing shifts around policy changes, and value-chain reconfiguration within industries such as textiles.
Existing statistical approaches to trade data typically operate on fewer dimensions. Gravity models regress bilateral volume on dyad-specific covariates anderson2011gravity, while network and community detection methods summarize pairwise trade relationships without incorporating the product or temporal dimensions wang2023fast, jian2025restricted. Among tensor-based approaches, BillioCasarinIacopiniKaufmann2023BDTR develop a Bayesian dynamic tensor autoregressive model for a two-layer international network of bilateral trade and capital, though at a much smaller scale without granular product-level flows, and kim2020measuring use product-level bilateral trade data to cluster exporter-importer pairs by the composition of goods exchanged. Our framework operates on the full four-way structure of transaction values.
After reviewing basic tensor notation and Bayesian Poisson tensor factorization in (ref), the paper makes four main contributions.
In this section, we briefly review the foundational concepts of non-negative tensor factorization and the Poisson Bayesian tensor factorization approach. For clarity, (ref) develop the methodology and inference under a general $M$-way tensor formulation.
Let $\mathcal{Y}\in \mathbb{R}_{\ge 0}^{I_1 \times I_2 \times \cdots \times I_M}$ be a nonnegative $M$-way tensor with entries $y_{\mathbf{i}}\ge 0$, where $\mathbf{i}\equiv(i_1, \ldots,i_M)$ denotes a multi-index and $i_m\in\{1,\ldots,I_m\}$ for $m=1,\ldots,M$. We write $\circ$ for the vector outer product. For vectors $\mathbf{a}\in\mathbb{R}^{I_1}$ and $\mathbf{b}\in\mathbb{R}^{I_2}$, $\mathbf{a}\circ \mathbf{b}\in\mathbb{R}^{I_1\times I_2}$ yields a matrix with entries $(\mathbf{a}\circ\mathbf{b})_{ij}=a_i b_j$. This is a rank-one matrix in the usual matrix rank sense. A rank-$K$ CANDECOMP/PARAFAC (CP) representation writes
where $\boldsymbol{\Theta}^{(m)}=[\boldsymbol{\theta}^{(m)}_{:, 1},\ldots,\boldsymbol{\theta}^{(m)}_{:,k}] \in\mathbb{R}_{\ge 0}^{I_m\times K}$ is the mode-$m$ factor matrix. Thus an $M$-way tensor admits exactly $M$ such factor matrices, one for each mode. Equivalently, entry $\mathbf{i}$ is
Nonnegative tensor decomposition estimates the $M$ mode-specific factor matrices $\{\Theta^{(m)}\}_m$ so that the low-rank reconstruction $\widehat{\mathcal{Y}}$ in (ref) approximates $\mathcal{Y}$ well. Two broad approaches are common. One treats $\widehat{\mathcal{Y}}$ as a deterministic approximation and solves $\min_{\boldsymbol{\Theta}^{(1)},\ldots,\boldsymbol{\Theta}^{(M)}\ge 0} \ \ \mathcal{L}\!\left(\mathcal{Y},\widehat{\mathcal{Y}}\right)$ where $\mathcal{L}$ is a loss or divergence (e.g., squared Frobenius loss or generalized KL divergence). The other specifies a probabilistic model for $\mathcal{Y}$ with mean or intensity tensor $\widehat{\mathcal{Y}}$ and learns the factors via maximum likelihood or Bayesian inference. Either way, the CP structure yields a parsimonious yet interpretable representation: an unrestricted tensor over $I_1\times\cdots\times I_M$ cells involves on the order of $\prod_{m=1}^M I_m$ parameters, whereas a rank-$K$ CP parameterization uses only $K\sum_{m=1}^M I_m$ mode-specific loadings.
CP factors are identifiable up to componentwise rescaling and permutation. Specifically, for any component $k$ and constants $c_1,\ldots,c_M>0$ with $\prod_{m=1}^M c_m=1$, replacing $\boldsymbol{\theta}^{(m)}_{:,k}$ by $c_m\,\boldsymbol{\theta}^{(m)}_{:,k}$ leaves the implied reconstruction unchanged. This rescaling does not affect the reconstruction (ref). We therefore do not impose an explicit identifiability constraint. Instead, we focus on within-component comparisons when interpreting results in application in (ref).
For count-valued $M$-way data, Bayesian Poisson tensor factorization (BPTF) cemgil2009bayesian,gopalan2015scalable,schein2015bayesian first models entries $y_{\mathbf{i}}\in\mathbb{Z}_{\ge 0}$ as conditionally independent Poisson random variables whose rates form a low-rank tensor:
The Bayesian formulation then places priors on the nonnegative factors, commonly Gamma for conjugate conditional updates, $\theta^{(m)}_{i_m k}\sim \mathrm{Gam}(\textsc{shape}=a_m,\textsc{rate}=b_m)$.
Posterior computation exploits an additive augmentation known as Poisson thinning, also called Poisson allocation. Specifically, we decompose each observed count $y_{\mathbf{i}}$ into $K$ latent component counts $y^{(k)}_{\mathbf{i}}$ such that
By Poisson additivity, this preserves the marginal model in (ref). Moreover, conditioning on the total $y_{\mathbf{i}}$ yields the equivalent multinomial form
where $\tau^{(k)}_{\mathbf{i}} =\nicefrac{\prod_{m=1}^M \theta^{(m)}_{i_m k}}{\sum_{\ell=1}^K \prod_{m=1}^M \theta^{(m)}_{i_m \ell}}$. This augmentation restores conditional conjugacy for the factors. In particular, with $\theta^{(m)}_{i_m k}\sim\text{Gam}(a_m,b_m)$, the complete conditional is
where $\mathbf{i}_{-m}\equiv(i_1,\ldots,i_{m-1},i_{m+1},\ldots,i_M)$ and $\sum_{\mathbf{i}_{-m}}$ denotes summation over all indices other than $i_m$. (ref) form the backbone of scalable Gibbs or variational inference algorithms for Poisson tensor factorization.
We propose the Poisson-randomized Gamma tensor factorization model, which extends the CP decomposition to semi-continuous arrays with many zeros and positive continuous values. The model preserves the interpretability of low-rank CP structure while addressing two limitations of standard Poisson factorization: the inability to accommodate mixed zero-positive support, and the Poisson mean-variance equality that prevents flexible modeling of over-dispersed, heavy-tailed positive magnitudes. In (ref), we introduce the generative model, and (ref) develops the reparameterization that enables efficient posterior inference.
To model positive continuous outcomes with exact zeros, we assume the following
where $\eta_{\mathbf{i}}$ is a latent count variable that is Poisson-distributed as
Here $\delta_0$ denotes a Dirac delta function at zero and the Gamma distribution is parameterized by positive shape $\eta_{\mathbf{i}} > 0 $ and rate $\beta_{\mathbf{i}} > 0$.
Marginalizing over the latent count $\eta_{\mathbf{i}}$ then yields the Poisson-randomized Gamma (PRG) distribution zhou2016augmentable, which we denote by $y_{\mathbf{i}} \sim \mathrm{PRG}(\lambda_{\mathbf{i}},\beta_{\mathbf{i}})$. Thus, the hierarchical construction in (ref) defines a PRG tensor model with low-rank location parameter $\lambda_{\mathbf{i}}$ and dispersion parameter $\beta_{\mathbf{i}}$, and falls within the general framework of (ref). The Poisson component governs the latent event rate, inducing zeros through the probability that $\eta_{\mathbf{i}}=0$ and setting the overall scale of nonzero observations.
In principle, when repeated tensor observations are available, the rate parameter in (ref) may vary across all cells, yielding a fully cell specific rate tensor with entries $\beta_{\mathbf{i}}$. In practice, however, we typically observe only a single tensor array. Under a single observation, lower-dimensional parameterizations are required to facilitate stable estimation via statistical pooling. This includes, for example, constant rate $\beta$, mode-specific structures such as $\beta_{i_{m^\star}}$, pairwise structures such as $\beta_{i_j,i_a}$, and other lower-dimensional parameterizations. For example, in our international trade application ((ref)), restricting $\beta_{\mathbf{i}}\equiv \beta_{i_{m^\star}}$ allows variance dispersion to vary uniquely across products. To maintain a unified presentation that covers any of these pooling structures, we retain the general notation $\beta_{\mathbf{i}}$ throughout the majority of the paper.
An observed zero $y_{\mathbf{i}}=0$ occurs only when $\eta_{\mathbf{i}}=0$, so during inference, $\eta_{\mathbf{i}}=0$ almost surely for all zero entries. More generally, if an entity along any mode (e.g., a particular $i_m$) is consistently associated with small outcomes, the likelihood pushes its corresponding factor loading $\Theta^{(m)}$ toward zero, shrinking its contribution to $\eta_{\mathbf{i}}$. Since these factors are constrained to be nonnegative, we place independent conditionally conjugate Gamma priors over them:
where $\mathbb{E}[\theta_{i_m k}^{(m)}]= \nicefrac{\alpha}{\alpha \nu^{(m)}}= \nicefrac{1}{\nu^{(m)}}$ and $\boldsymbol{\nu} = (\nu^{(1)}, \dots, \nu^{(M)})^T$ collects these mode-specific rates. In all of our experiments, we fix $\alpha=0.1$, which encourages shrinkage among the latent factors.
Having specified the full hierarchical model, we first establish an entrywise identifiability property of the PRG sampling distribution.
Let $\boldsymbol{\Lambda}, \mathbf{B} \in \mathbb{R}_+^{I_1\times\cdots\times I_M}$ denote the intensity and rate tensors with entries $\lambda_{\mathbf{i}}$ and $\beta_{\mathbf{i}}$. As an entrywise distribution level result, (ref) uniquely determines each parameter pair $(\lambda_{\mathbf{i}}, \beta_{\mathbf{i}})$ and ensures that no simultaneous scaling can produce the same marginal data distribution. Consequently, due to conditional independence, equality of the joint sampling distribution implies entrywise equality of the full parameter tensors, meaning $\boldsymbol{\Lambda} = \tilde{\boldsymbol{\Lambda}}$ and $\mathbf{B} = \tilde{\mathbf{B}}$. To lift this result from $\boldsymbol{\Lambda}$ to its latent factors, we impose the following standard CP uniqueness condition.
Combining the entrywise PRG identifiability with the CP uniqueness condition yields global structural identifiability of the full rate tensor and the latent CP factors:
Our goal in the subsequent (ref) is to perform Bayesian inference for the latent factors $\{\Theta^{(m)}\}_{m}$ under the model in (ref). Direct posterior computation is intractable because the likelihood couples the factors through a sum--product structure, so the full conditionals are not available in closed form. We therefore use an equivalent auxiliary-variable augmentation that preserves the original marginal model while simplifying posterior computation. Our construction combines Poisson–multinomial allocation with a new thinning for the Gamma layer in (ref) that links positive magnitudes to the latent structure.
Scheme 1: Poisson thinning. The first augmentation is Poisson thinning, which provides a latent-source allocation of the Poisson CP mean across the \(K\) rank-one components. This is the tensor analogue of allocative Poisson factorization: each latent count is decomposed into component-specific sub-counts that record how much of the total intensity is attributed to each CP component. In our setting, the same Poisson--multinomial identity used in Bayesian Poisson tensor factorization reviewed in (ref) applies to the latent Poisson layer, except that the total count is now the latent variable \(\eta_{\mathbf{i}}\) rather than an observed count---i.e., \[ \eta_{\mathbf{i}}\sim \mathrm{Pois}(\mu_{\mathbf{i}}), \qquad \mu_{\mathbf{i}} = \sum_{k=1}^K \prod_{m=1}^M \theta^{(m)}_{i_mk}, \] where the per-component sub-counts are then defined as
which preserve the original marginal model by Poisson additivity. Equivalently, conditional on the total \(\eta_{\mathbf{i}}\), the sub-counts follow a multinomial allocation,
This Poisson--multinomial thinning representation is valuable both statistically and computationally: it yields an interpretable component-wise allocation of latent events, restores the conditional structure needed for efficient factor updates, and implies that whenever \(\eta_{\mathbf{i}}=0\), all sub-counts \(\eta^{(k)}_{\mathbf{i}}\) are automatically zero almost surely, so no component-wise allocation is needed for that entry. Hence, this thinning step scales as \(O(\|\mathcal{Y}\|_0 K)\), where \(\|\mathcal{Y}\|_0\) denotes the number of nonzero tensor entries.
Scheme 2: Gamma thinning and binary indicators. In our Poisson-randomized Gamma specification in (ref), the latent count $\eta_{\mathbf{i}}$ plays a dual role: it determines whether $y_{\mathbf{i}}$ is zero and, conditional on being positive, acts as the Gamma shape parameter. Under Scheme 1, we thin $\eta_{\mathbf{i}}$ into component-specific sub-counts $(\eta^{(k)}_{\mathbf{i}})_{k=1}^K$. Exploiting the additivity of the Gamma distribution under a common rate, we introduce a matching decomposition of each observation, \looseness=-1 \[ y_{\mathbf{i}} \;=\; \sum_{k=1}^K y_{\mathbf{i}}^{(k)}, \] where, conditional on $\eta^{(k)}_{\mathbf{i}}$, the latent contributions are independent with
This construction is consistent with the original model: $y_{\mathbf{i}}=0$ occurs if and only if $\eta^{(k)}_{\mathbf{i}}=0$ for all $k$, i.e., if the total count $\eta_{\mathbf{i}}=\sum_{k=1}^K \eta^{(k)}_{\mathbf{i}}$ is zero. Consequently, the probability of a zero observation is preserved under thinning, matching (ref).
This formulation, in turn, makes each nonnegative observation $y_{\mathbf{i}}>0$ now a hybrid sum of components that are partly degenerate at zero and partly Gamma contributions, so the marginal likelihood is not available in a tractable form for direct inference. We introduce auxiliary binary activity indicators $b^{(k)}_{\mathbf{i}}\in\{0,1\}$ that switch component $k$ on or off. Specifically, $b^{(k)}_{\mathbf{i}} =\mb{1} \{\eta^{(k)}_{\mathbf{i}}>0 \}=\mb{1} \{y^{(k)}_{\mathbf{i}}>0 \}$. The active set is then defined as $ \Delta_{\mathbf{i}}:=\{ k\in \{ 1,\dots, K\}: b_{\mathbf{i}}^{(k)}=1 \}$. We defer detailed motivation and discussion of the indicator and activity set to (ref). Conditioning on $\Delta_{\mathbf{i}}$ and other latent variables, each positive $y_{\mathbf{i}}$ is a sum of independent Gamma random variables with the same rate. When the sequence of latent variables $(y_{\mathbf{i}}^{(k)})_{k\in \Delta_{\mathbf{i}}}$ for a given multi-index $\mathbf{i}$ consists of independent Gamma random variables sharing a common rate, their normalized vector follows a Dirichlet distribution:
This yields a second thinning step, which refer to throughout as Gamma--Dirichlet thinning.
Although one could in principle work only with the first thinning scheme, it leaves an intractable product-sum term that obstructs closed-form updates. The second thinning scheme and the binary indicators decouple this term and are key to making inference practical, as will become clear in (ref).
We now present the complete hierarchical model, collecting the observation model, latent structure, and priors into a single specification as follows:
This summary consolidates the reparameterizations and augmentation introduced above and will serve as the basis for posterior inference. Throughout the rest of the paper, we work with this model representation.
Given an observed tensor $\mathcal{Y}$, posterior inference targets the joint distribution of all unknown quantities under the hierarchical model in (ref). We partition these unknowns into the binary activity indicators $\mathcal{A}\stackrel{\mathrm{def}}{=} \{ b^{(k)}_{\mathbf{i}} \}_{\mathbf{i},\,k}$ and the remaining latent variables $\mathcal{Z}\stackrel{\mathrm{def}}{=} \{y_{\mathbf{i}}^{(k)}\}_{\mathbf{i}, k} \;\cup\; \{\eta_{\mathbf{i}}^{(k)}\}_{\mathbf{i}, k} \;\cup\; \{\Theta^{(m)}\}_{m}$. Our goal is therefore to approximate the posterior $p(\mathcal{A},\mathcal{Z} \mid \mathcal{Y},\mathcal{H})$, where $\mathcal{H}\stackrel{\mathrm{def}}{=}\{\alpha\} \;\cup\; \{\boldsymbol{\nu}\} \;\cup\; \{\beta_{\mathbf{i}}\}_{\mathbf{i}}$ denotes the collection of model hyperparameters. In what follows, we treat $\mathcal{H}$ separately from the posterior updates: $\alpha$ is fixed a priori at $0.1$, while $\boldsymbol{\nu}$ and $\{ \beta_{\mathbf{i}}\}$ are updated by an empirical Bayes step. We therefore suppress its dependence in the notation.
In principle, one could perform full MCMC using the conditionals, but this becomes computationally prohibitive for high-dimensional tensors. We therefore adopt a hybrid strategy that alternates between sampling the binary indicators $\mathcal{A}$ and updating the remaining block $\mathcal{Z}$ by coordinate-ascent variational inference (CAVI) conditional on $\mathcal{A}$. We first describe the variational approximation and the general CAVI framework, then derive the model-specific updates, and finally present the sampling step for $\mathcal{A}$.
Conditional on the current binary indicators $\mathcal{A}$, we approximate the conditional posterior $p(\mathcal{Z}\mid \mathcal{Y},\mathcal{A})$ by a tractable variational distribution $q(\mathcal{Z})$. Variational inference replaces the exact posterior by a simpler family of distributions and chooses the member that is closest in Kullback--Leibler divergence \(\mathrm{KL}\!\left(q(\mathcal{Z})\,\|\,p(\mathcal{Z}\mid \mathcal{Y},\mathcal{A})\right)\), or equivalently maximizes the evidence lower bound (ELBO) Jordan1999VI, $\mathrm{ELBO}(q)=\mathbb{E}_{q}\left[\log p(\mathcal{Y},\mathcal{Z}\mid \mathcal{A})\right]-\mathbb{E}_{q}\left[\log q(\mathcal{Z})\right]$. The derivation of the ELBO for our model is given in Supplementary (ref).
In CAVI, one assumes that $q$ factorizes over blocks of $\mathcal{Z}$, and then updates each block holding the others fixed. Under a mean-field approximation, the optimal factor for a block $\mathcal{Z}_j$ satisfies \[ q^\star(\mathcal{Z}_j)\ \propto\ \exp\!\left\{ \mathbb{E}_{q(\mathcal{Z}_{-j})}\left[\log p(\mathcal{Z}_j\mid \mathcal{Z}_{-j},\mathcal{Y},\mathcal{A})\right] \right\}, \] where $\mathcal{Z}_{-j}$ denotes all remaining blocks. Thus, each variational update has the same functional form as the corresponding complete conditional, with its parameters replaced by expectations under the current variational distribution blei2017variational.
We next provide the complete conditionals and corresponding variational updates for the latent variables $\mathcal{Z}$, conditional on the sampled indicators $\mathcal{A}$. We also give the associated empirical Bayes updates for the hyperparameters.
(i) Update for Dirichlet mixing proportions $\pi^{(k)}_{\mathbf{i}}$ and Gamma sub-mass $y^{(k)}_{\mathbf{i}}$. For an entry with $y_{\mathbf{i}}>0$, define proportions over the active set $\pi^{(k)}_{\mathbf{i}}=\nicefrac{y^{(k)}_{\mathbf{i}}}{y_{\mathbf{i}}}$ for $k\in\Delta_{\mathbf{i}}$, so that $\sum_{k\in\Delta_{\mathbf{i}}}\pi^{(k)}_{\mathbf{i}}=1$ and $y^{(k)}_{\mathbf{i}}=\pi^{(k)}_{\mathbf{i}}\,y_{\mathbf{i}}$. The complete conditional is Dirichlet, as in the Gamma thinning step (ref), with $(\pi^{(k)}_{\mathbf{i}})_{k\in\Delta_{\mathbf{i}}} \mid \cdot\ \sim\ \mathrm{Dir}\big((\eta^{(k)}_{\mathbf{i}})_{k\in\Delta_{\mathbf{i}}}\big)$. However, because \((\eta^{(k)}_{\mathbf{i}})_{k\in\Delta_{\mathbf{i}}}\) is itself latent, the variational distribution for \((\pi^{(k)}_{\mathbf{i}})_{k\in\Delta_{\mathbf{i}}}\) is written with free Dirichlet parameters \((\omega^{(k)}_{\mathbf{i}})_{k\in\Delta_{\mathbf{i}}}\) as $q\big((\pi^{(k)}_{\mathbf{i}})_{k\in\Delta_{\mathbf{i}}}\big)=\mathrm{Dir}\big( (\pi^{(k)}_{\mathbf{i}})_{k\in\Delta_{\mathbf{i}}};(\omega^{(k)}_{\mathbf{i}})_{k\in\Delta_{\mathbf{i}}}\big)$. The coordinate update is then obtained by matching these variational parameters to the expected conditional parameters, yielding $\omega^{(k)}_{\mathbf{i}}=\mathbb{E}_{q}\left[\eta^{(k)}_{\mathbf{i}}\right], \ k\in\Delta_{\mathbf{i}}$: \[ \mathbb{E}_{q}\left[y^{(k)}_{\mathbf{i}}\right] = y_{\mathbf{i}}\,\mathbb{E}_{q}\left[\pi^{(k)}_{\mathbf{i}}\right] = y_{\mathbf{i}}\, \frac{\omega^{(k)}_{\mathbf{i}}}{\sum_{\ell\in\Delta_{\mathbf{i}}}\omega^{(\ell)}_{\mathbf{i}}}, \] and the geometric expectation used in subsequent derivations is \[ \mathbb{G}_{q}\left[y^{(k)}_{\mathbf{i}}\right] =\exp\!\big(\mathbb{E}_{q}\left[\log y^{(k)}_{\mathbf{i}}\right]\big) = y_{\mathbf{i}}\, \exp\!\left\{\psi(\omega^{(k)}_{\mathbf{i}}) -\psi\!\Big(\sum_{\ell\in\Delta_{\mathbf{i}}}\omega^{(\ell)}_{\mathbf{i}}\Big)\right\}. \]
(ii) Update for Poisson sub-counts $\eta^{(k)}$. For each component $k$, the conditional model $ y^{(k)}_{\mathbf{i}}\mid \eta^{(k)}_{\mathbf{i}} \sim \mathrm{Gam}(\eta^{(k)}_{\mathbf{i}},\beta_{\mathbf{i}}), $ so that $y^{(k)}_{\mathbf{i}}>0$ implies $\eta^{(k)}_{\mathbf{i}}\ge 1$ almost surely. Since the prior is $\eta^{(k)}_{\mathbf{i}}\sim \mathrm{Pois}\big(\lambda^{(k)}_{\mathbf{i}}\big)$ with $\lambda^{(k)}_{\mathbf{i}}=\prod_{m=1}^M \theta^{(m)}_{i_m k}$ denoted as the $k-$th component CP rate, conditioning only on the event $y^{(k)}_{\mathbf{i}}>0$ yields the zero-truncated Poisson law, \[ \eta^{(k)}_{\mathbf{i}}\mid \big(y^{(k)}_{\mathbf{i}}>0\big) \ \sim\ \mathrm{Pois}_+\!\left(\lambda^{(k)}_{\mathbf{i}}\right). \]
More importantly for inference, conditioning on positive $y^{(k)}_{\mathbf{i}}$ and $\lambda^{(k)}_{\mathbf{i}}$ leads to a closed-form posterior for the integer shape parameter. Up to constants not depending on $\eta^{(k)}_{\mathbf{i}}=n\ge 1$, \[ \mathbb{P}\!\left(\eta^{(k)}_{\mathbf{i}}=n \mid y^{(k)}_{\mathbf{i}},\lambda^{(k)}_{\mathbf{i}},\beta_{\mathbf{i}}\right) \ \propto\ \frac{\big(\lambda^{(k)}_{\mathbf{i}}\big)^n}{n!}\cdot \frac{\beta_{\mathbf{i}}^n \big(y^{(k)}_{\mathbf{i}}\big)^{n-1}}{\Gamma(n)}. \] The normalizing constant is a modified Bessel function of the first kind, so the complete conditional is a zero-truncated Bessel distribution yuan2000bessel:
Following zhou2016augmentable, if $N\sim \mathrm{Bessel}_{-1}(\alpha)$, then for $n\in\{1,2,\ldots\}$, \[ \mathbb{P}(N=n) = \frac{(\nicefrac{\alpha}{2})^{2n-1}}{I_{-1}(\alpha)\,n!\,\Gamma(n)}, \qquad I_{-1}(\alpha) = \sum_{n=1}^\infty \frac{(\nicefrac{\alpha}{2})^{2n-1}}{n!\,\Gamma(n)}. \] Since the order is fixed at $-1$ here, the conditional law is indexed only by the scalar argument $\alpha$. In this fixed-order form, the Bessel distribution is a one-parameter discrete exponential family, and more broadly the Bessel family is a power-series class of underdispersed count distributions. Accordingly, in the variational approximation we write $q\big(\eta^{(k)}_{\mathbf{i}}\big)=\mathrm{Bessel}_{-1}\!\big(\eta^{(k)}_{\mathbf{i}};\,\zeta^{(k)}_{\mathbf{i}}\big)$, where $\zeta^{(k)}_{\mathbf{i}}$ denotes the corresponding Bessel variational parameter.
For $k\in\Delta_{\mathbf{i}}$, the complete conditional in (ref) is zero-truncated Bessel. Let $\lambda^{(k)}_{\mathbf{i}}=\prod_{m=1}^M \theta^{(m)}_{i_mk}$, and denote its geometric expectation by \[ \mathbb{G}_{q}\left[\lambda^{(k)}_{\mathbf{i}}\right] =\exp\!\big(\mathbb{E}_{q}\left[\log \lambda^{(k)}_{\mathbf{i}}]\right]\big) = \prod_{m=1}^M \exp\!\big(\mathbb{E}_{q}\left[\log \theta^{(m)}_{i_m k}\right]\big) = \prod_{m=1}^M \mathbb{G}_{q}\left[\theta^{(m)}_{i_m k}\right]. \] The corresponding CAVI update therefore remains in the same family: \[ q\big(\eta^{(k)}_{\mathbf{i}}\big) = \mathrm{Bessel}_{-1}\!\big(\eta^{(k)}_{\mathbf{i}};\,\zeta^{(k)}_{\mathbf{i}}\big), \qquad \zeta^{(k)}_{\mathbf{i}} = 2\sqrt{\beta_{\mathbf{i}}\,\mathbb{G}_{q}\left[y^{(k)}_{\mathbf{i}}\right]\,\mathbb{G}_{q}\left[\lambda^{(k)}_{\mathbf{i}}\right]}. \] Under the variational distribution $q(\eta^{(k)}_{\mathbf{i}})=\mathrm{Bessel}_{-1}(\eta^{(k)}_{\mathbf{i}};\,\zeta^{(k)}_{\mathbf{i}})$, the expectation of $\eta^{(k)}_{\mathbf{i}}$ is \[ \mathbb{E}_{q}\left[\eta^{(k)}_{\mathbf{i}}\right] = \zeta^{(k)}_{\mathbf{i}} \cdot \frac{I_{0}(\zeta^{(k)}_{\mathbf{i}})}{I_{-1}(\zeta^{(k)}_{\mathbf{i}})} . \] From this update, the necessity of the second thinning scheme and the binary indicators becomes clear in (ref). Under the first thinning scheme alone, the Bessel update depends on a sum-product term inside the square root, $\eta^{(k)}_{\mathbf{i}} \mid \cdot \sim\ \mathrm{Bessel}_{-1}\!\left( 2\sqrt{ \beta_{\mathbf{i}}\, y^{(k)}_{\mathbf{i}} \sum_{k=1}^K \prod_{m=1}^M \theta^{(m)}_{i_m k}} \right)$, so inference requires approximating $\mathbb{G}_{q}\!\left[\sum_{k=1}^K \prod_{m=1}^M \theta^{(m)}_{i_m k}\right]$. In our experiments not reported here, delta-method approximations, Monte Carlo estimates, and arithmetic-mean substitutes were all either inaccurate or computationally expensive. The second thinning scheme decouples this bottleneck and restores closed-form inference.
(iii) Update for factor matrices $\Theta$. Under the Poisson allocation $\eta^{(k)}_{\mathbf{i}}\sim\text{Pois}(\lambda^{(k)}_{\mathbf{i}})$ and the Gamma prior $\theta^{(m)}_{i_m k}\sim \text{Gam}(\alpha,\alpha\nu^{(m)})$, the complete conditional of $\theta^{(m)}_{i_m k}$ is Gamma by Poisson--Gamma conjugacy
Indeed, collecting only the terms that depend on $\theta^{(m)}_{i_m k}$, \[ \log p(\theta^{(m)}_{i_m k}\mid \cdot ) = \Big(\alpha-1+\sum_{\mathbf{i}_{-m}}\eta^{(k)}_{\mathbf{i}}\Big)\log \theta^{(m)}_{i_m k} -\Big(\alpha\nu^{(m)}+\sum_{\mathbf{i}_{-m}}\prod_{\ell\neq m}\theta^{(\ell)}_{i_\ell k}\Big)\theta^{(m)}_{i_m k} +\text{const}. \] Taking expectations with respect to the variational distribution, the CAVI update for $\theta^{(m)}_{i_m k}$ remains in the Gamma family $q\big(\theta^{(m)}_{i_m k}\big) = \mathrm{Gam}\big(\theta^{(m)}_{i_m k};\gamma^{(m)}_{i_m k},\delta^{(m)}_{i_m k}\big)$, with parameters updated as \[ \gamma^{(m)}_{i_m k} = \alpha+\sum_{\mathbf{i}_{-m}}\mathbb{E}_{q}\left[\eta^{(k)}_{\mathbf{i}}\right], \qquad \delta^{(m)}_{i_m k} = \alpha\nu^{(m)} + \sum_{\mathbf{i}_{-m}}\prod_{\ell\neq m}\mathbb{E}_{q}\left[\theta^{(\ell)}_{i_\ell k}\right], \] where the product of expectations follows from the mean-field factorization across modes. With $\psi(\cdot)$ denotes the digamma function, the arithmetic and geometric expectations are \[ \mathbb{E}_{q}\left[\theta^{(m)}_{i_m k}\right] = \frac{\gamma^{(m)}_{i_m k}}{\delta^{(m)}_{i_m k}}, \qquad \mathbb{G}_{q}\left[\theta^{(m)}_{i_m k}\right] = \exp\!\Big(\mathbb{E}_{q}\left[\log \theta^{(m)}_{i_m k}\right]\Big) = \frac{\exp\!\big(\psi(\gamma^{(m)}_{i_m k})\big)}{\delta^{(m)}_{i_m k}}. \]
(iv) Empirical Bayes updates for $\beta_{i_{m^\star}}$ and $\boldsymbol{\nu}$. The Gamma rate parameters $\beta_{\mathbf{i}}$ and $\nu^{(m)}$ are not assigned priors, so we update them by empirical Bayes, i.e., by maximizing the ELBO with respect to these parameters while holding the current variational distribution fixed schein2015bayesian. This yields closed-form updates of moment type:
For the Gamma rate parameter, the entrywise update takes the form $\beta_{\mathbf{i}} =\nicefrac{\sum_{k=1}^K \mathbb{E}_{q}\left[\eta^{(k)}_{\mathbf{i}}\right]} {y_{\mathbf{i}}}$. More generally, when the rate parameter is constrained to be shared across a collection of tensor entries, the corresponding update is obtained by pooling the numerator and denominator over all entries in that collection. In particular, if the rate parameter is assumed to depend only on a designated mode $m^\star$, then the update becomes
where $\sum_{i_{-m^\star}}$ denotes summation over all indices except $i_{m^\star}$.
The binary indicator $b^{(k)}_{\mathbf{i}}$ is updated alternately with the remaining latent variables and hyperparameters. Throughout this subsection, latent quantities such as $\eta^{(k)}_{\mathbf{i}}$ denote their current variational expectations from (ref).
This augmented variable addresses a computational bottleneck at the zero--positive boundary. Under the Poisson-randomized Gamma construction ((ref)), $y_{\mathbf{i}}^{(k)}=0 \iff \eta^{(k)}_{\mathbf{i}}=0$, coupling their supports and rendering the full conditionals degenerate on the boundary. Coordinate-ascent or Gibbs updates that treat these variables separately inherit near-reducibility at this interface: in MCMC, components exhibit poor mixing across the boundary; in mean-field variational inference, the independent factorization assumption conflicts with the coupled support, promoting excessively sparse solutions.
When $y_{\mathbf{i}}=0$, the update is deterministic: $b^{(k)}_{\mathbf{i}}=0$ for all $k$. When $y_{\mathbf{i}}>0$, we sample $b^{(1)}_{\mathbf{i}},\ldots,b^{(K)}_{\mathbf{i}}$ sequentially, each conditioned on $\{b^{(\ell)}_{\mathbf{i}}\}_{\ell\neq k}$, thereby avoiding marginalization over $2^{K-1}$ activation patterns. Crucially, each $b^{(k)}_{\mathbf{i}}$ is drawn from a collapsed conditional that marginalizes out $(y^{(k)}_{\mathbf{i}},\eta^{(k)}_{\mathbf{i}})$; conditioning on either variable directly would render the update degenerate, reintroducing the boundary stickiness the augmentation is designed to eliminate.
Let $y_{\mathbf{i}}^{(-k)}=\{y^{(\ell)}_{\mathbf{i}}\}_{\ell\neq k}$ and $b_{\mathbf{i}}^{(-k)}=\{b^{(\ell)}_{\mathbf{i}}\}_{\ell\neq k}$ denote the collections excluding component $k$, and define $s_{\mathbf{i}}^{(-k)}\;=\;\sum_{\ell\neq k} b^{(\ell)}_{\mathbf{i}}\,\eta^{(\ell)}_{\mathbf{i}}$. Then, for $y_{\mathbf{i}}>0$,
where the two weights correspond to the two competing explanations of $y_{\mathbf{i}}$:
Here $f_{y\mid b^{(k)}=1}$ is the marginal density of $y_{\mathbf{i}}$ when component $k$ is active, obtained by convolving the PRG term of the active component with the Gamma contribution of the others:
where $\mathrm{PRG}_{+}$ denotes the zero-truncated Poisson--Gamma mixture density induced by $\eta^{(k)}\sim \mathrm{Pois}_{+}(\lambda^{(k)}_{\mathbf{i}})$ and $y^{(k)}\mid \eta^{(k)}\sim \mathrm{Gam}(\eta^{(k)},\beta_{\mathbf{i}})$. Its closed form is given in Supplementary (ref).
Together, for $y_{\mathbf{i}}>0$ the collapsed Bernoulli update (ref) admits the closed form
The Bessel function in (ref) is implemented in libraries like scipy.special.ive SciPyIveDocs but can be expensive to compute for large $s_{\mathbf{i}}^{(-k)}$. We thus adopt the hybrid sampling rule in (ref).
The approximate branch of (ref) is justified by (ref) below, which shows that as $s_{\mathbf{i}}^{(-k)}\to\infty$ the exact conditional converges to: \[ b^{(k)}_{\mathbf{i}}=\mathbf{1}(\kappa^{(k)}_{\mathbf{i}} >0),\quad \kappa^{(k)}_{\mathbf{i}} \sim \operatorname{Pois}\!\left( \nicefrac{\beta_{\mathbf{i}} y_{\mathbf{i}} \lambda^{(k)}_{\mathbf{i}}}{s_{\mathbf{i}}^{(-k)} } \right), \] corresponding to the Bernoulli--Poisson link of zhou2015infinite. Although the result is asymptotic, the approximation is accurate for moderate $s_{\mathbf{i}}^{(-k)}$, motivating the threshold $20$ above.
Explicit expressions for the two conditional Bernoulli success probabilities, together with the proof, are given in Supplementary (ref). (ref), together with the CAVI updates for $(y^{(k)},\eta^{(k)},\Theta)$, constitutes the full hybrid inference procedure summarized in (ref).
In this section, we evaluate the proposed model against recent baselines in a series of controlled synthetic experiments designed to assess predictive performance under the data-generating process described in (ref). The same experimental setup, baselines, and reconstruction method are reused in (ref) for out-of-sample prediction on the real trade data.
Data. We generate a four-way tensor $\mathcal{Y} \in \mathbb{R}_{\ge 0}^{I \times J \times A \times T}$ with rank $R=20$ and mode dimensions $(I,J,A,T) = (50,40,30,10)$. We first construct four factor matrices $\Theta^{(1)} {\in} \mathbb{R}_{\ge 0}^{I \times R}$, $\Theta^{(2)} {\in} \mathbb{R}_{\ge 0}^{J \times R}$, $\Theta^{(3)} {\in} \mathbb{R}_{\ge 0}^{A \times R}$, and $\Theta^{(4)} \in \mathbb{R}_{\ge 0}^{T \times R}$ by first drawing each entry from $\mathrm{Gam}(1,1)$ and then imposing sparsity by setting each entry to zero with probability $0.3$. We then define $\lambda_{ijat}{=}\sum_{r=1}^R \theta_{ir}^{(1)}\theta_{jr}^{(2)}\theta_{ar}^{(3)}\theta_{tr}^{(4)}$, sample $\eta_{ijat}\sim\mathrm{Pois}(\lambda_{ijat})$, and set $y_{ijat}{=}0$ if $\eta_{ijat}{=}0$ and $y_{ijat}{\sim}\mathrm{Gam}(\eta_{ijat},\beta_a)$ otherwise, where $(\beta_a)_{a=1}^A$ is set to an equally spaced grid from $0.01$ to $0.1$. For each tensor, we then generate 30 random 80-20 train-test splits. \looseness=-1
Models. We compare against two baselines for nonnegative tensors with zeros and positive continuous values: non-negative tensor factorization minimizing the least squared loss (NTF-LS) and minimizing the $\beta$-divergence with $\beta\!=\!1.5$ (NTF-$\beta$), with the latter corresponding to a Tweedie loss. To fit NTF-LS we use the implementation in Tensorly JMLRtensorly and for NTF-$\beta$ we use NNEinFact hood2026near. We also consider two variants of our method: shared, which uses a common Gamma rate $\beta_{\mathbf{i}}=\beta$ across all entries, and local, which allows the rate to vary along one designated mode $\beta_{\mathbf{i}}=\beta_a$. Each method is fit to each train-test split of each data tensor with $K\in\{15,20,25\}$.
Evaluation. We evaluate performance via mean absolute error over all held-out entries (MAE) and over nonzero held-out entries (MAE-NZ), mean absolute relative error over nonzero entries (MArel-NZ), Hamming loss on zeros (HAM-Z), and the area under the ROC curve (AUC) for zero versus nonzero entries. For NTF-$\beta$, we use multiple learning rates and report the best in terms of test-set MAE. For NTF-LS and NTF-$\beta$, held-out entries are reconstructed from point-estimated factors, while for BPGPTF they are reconstructed with the geometric mean under the variational posterior, $\widehat{y}_{\mathbf{i}}=\nicefrac{\sum_{k=1}^K \prod_{m=1}^M \mathbb{G}_{}\left[\theta^{(m)}_{i_m k}\right]}{\beta_{\mathbf{i}}}$.
(ref) reports the mean and standard error over 30 random train--test splits, where the proposed model with the local Gamma rate is labeled local. The local model achieves the best predictive performance across all metrics and rank settings, with particularly large gains in the Hamming loss on zeros, indicating improved recovery of sparsity. For all methods, predictive performance generally improves as $K$ increases.
We apply the proposed model to international merchandise trade data. In (ref) we describe the dataset and preprocessing. We then evaluate out-of-sample predictive performance under the strong-generalization regime in (ref). In (ref), we illustrate exploratory analysis through a case study of the fitted components.
We analyze annual merchandise-trade flows from the World Integrated Trade Solution (WITS) wits2025, a World Bank platform providing standardized access to trade data reported by national customs authorities. We use export values at the HS-1996 2-digit (HS2) level, yielding all 96 product categories (e.g., 01 Live animals, 84 Machinery, 85 Electrical machinery) for 196 countries from 1997 to 2023, measured in thousand USD. We arrange the data as a four-mode tensor $\mathcal{Y}\in\mathbb{R}_{\ge0}^{I\times J\times A\times T}$, where $i$ indexes exporters, $j$ indexes importers, $a$ indexes HS2 products, and $t$ indexes years. Each entry $y_{ijat}$ records the value of product $a$ exported from country $i$ to country $j$ in year $t$, so $y_{ijat}=0$ may reflect either a true zero flow or an unreported observation. More details are in Supplementary (ref).
Missingness is pervasive in international trade data, arising from reporter-year non-submission, confidentiality suppressions, reporting thresholds, and HS revisions UNComtradeMethodology,WITSUserGuide. We preprocess the data by flagging reporter-year lapses (exporters with no positive flows in a given year), imputing short internal gaps, and retaining a binary observation mask so that structurally missing entries do not contribute to estimation; details are in Supplementary (ref). The resulting tensor covers 151 countries, 27 years, 96 products, and 59,100,192 entries, of which 21.3% are nonzero. To allow product-specific dispersion, we set the Gamma rate to vary across HS2 categories: $y_{ijat}\mid \eta_{ijat} \overset{ind.}{\sim} \mathrm{Gam}(\eta_{ijat},\beta_a)$.
We adopt the strong-generalization regime of marlin2004collaborative, using a leave-one-year-out design. For each year $t_0 \in \{1,\dots,T\}$, we fit the model to the tensor $\mathcal{Y}_{\neg t_0} \in \mathbb{R}_{\geq 0}^{I \times J \times A \times (T-1)}$ defined as the concatenation $\mathcal{Y}_{\neg t_0} \overset{\mathrm{def}}{=} \textrm{concat}(\{\mathcal{Y}_{:,:,:,t}\}_{t \neq t_0})$ of the set of all 3-mode slices for each year $t \neq t_0$. We then split $\mathcal{Y}_{t_0} = \mathcal{Y}_{t_0}^{\textrm{obs}} \cup \mathcal{Y}_{t_0}^{\textrm{eval}}$ into a 75% observed set and a 25% evaluation set, refitting only the time-mode factor $\Theta^{(4)}$ on the observed set while fixing $\Theta^{(1)}$, $\Theta^{(2)}$, and $\Theta^{(3)}$. We consider three masking schemes: random holdout, complement-of-dense-block, and dense-block, following the evaluation design of schein2015bayesian. The dense block is defined by the $N_0$ highest-volume countries in each held-out year; see Supplementary (ref) for details and an illustration.
(ref) summarizes performance across $K\in\{5,10,50,100,125\}$ for the three masking schemes; full numerical results at $K=125$ are in Supplementary (ref). NTF-LS is competitive at small ranks but saturates as $K$ increases. Our proposed model with the local Gamma rate benefits substantially from additional rank, consistently outperforming all baselines at moderate to large $K$ across masking regimes.
We now illustrate how the proposed model can be used for exploratory analysis. Since $K=125$ gives the best out-of-sample predictive performance, we refit the model to the full unmasked dataset. The estimated product-specific Gamma rates $\beta_a$ are strongly negatively associated with both the mean and variance of nonzero entries in each product subtensor (Spearman $\rho=-0.991$ and $-0.947$), consistent with heavier industries exhibiting larger and more volatile flows; see Supplementary (ref) for details.
For each component $k$, the nonnegative factors $\boldsymbol{\theta}_{:,k}^{(1)}$, $\boldsymbol{\theta}_{:,k}^{(2)}$, $\boldsymbol{\theta}_{:,k}^{(3)}$, $\boldsymbol{\theta}_{:,k}^{(4)}$ describe its exporter, importer, product, and temporal structure. Earlier (ref) illustrates a component dominated by Russian exports of mineral, metal, and industrial products, with temporal deviations around the 2009 crisis, 2014 Crimea annexation, 2020 pandemic, and 2022 Russia--Ukraine war. A small number of components reflect reporting conventions or definitional breaks rather than substantive trade structure; see Supplementary (ref) for an example.
Concentration and entropy. To compare components systematically, we normalize each factor vector into shares ($\tilde{\theta}_{ik} = \nicefrac{\theta^{(1)}_{ik}}{\sum_{i'} \theta^{(1)}_{i'k}}$, analogously for importers and products) and measure each by the normalized Shannon entropy $H_k = -\nicefrac{1}{\ln I} \cdot \sum_{i} \tilde{\theta}_{ik}\ln\tilde{\theta}_{ik} \in [0,1]$, where $I$ is the number of entities in the mode and $H_k=1$ indicates a uniform distribution while $H_k\to 0$ indicates full concentration on a single entity. (ref) reveals a sharp cross-mode trade-off: components with more uniform product coverage ($H^P_k$ high) systematically have more concentrated partner profiles ($H_{\min,k} = \min\{H^E_k, H^M_k\}$ low). Components anchored to a single product class tend to involve many trading partners, reflecting either the geographic concentration of natural resources (e.g.\ cocoa, coffee) or a small number of technologically advanced exporters (e.g.\ aircraft, pharmaceuticals), while components spanning many products are driven by a few large, diversified economies.
\newlength{\tradepanelheight} {3.9cm}
Two-way trade balance. For each ordered direction $i\to j$, the decomposition implies a product distribution $q_{ij}(a) = \sum_{k=1}^{K} \pi_{ij}(k) \tilde\theta^{(3)}_{ak}$ with $\pi_{ij}(k) \propto \tilde \theta^{(1)}_{ik} \tilde\theta^{(2)}_{jk}$, i.e.\ the share of $i$'s exports to $j$ in each product $a$ as attributed by the latent channels. We measure compositional symmetry by the Grubel--Lloyd index $\mathrm{GL}_{ij} = 1 - \tfrac{1}{2}\sum_a |q_{ij}(a) - q_{ji}(a)| \in [0, 1]$, with $\mathrm{GL}_{ij}\to 1$ for similar product baskets in both directions (intra-industry trade) and $\mathrm{GL}_{ij}\to 0$ for disjoint baskets (inter-industry specialization). The classical index grubel_lloyd1975 is computed on raw bilateral flows, while ours operates on the rank-$K$ latent reconstruction $q_{ij}$, suppressing flow noise and isolating structural composition. Among pairs whose smaller direction of trade lies in the top 5% of directional volumes ((ref)), the symmetric end is dominated by pairs within EU, especially Visegr\'ad countries, the canonical intra-industry pattern among similar economies helpman_krugman1985, while the asymmetric end pairs commodity exporters (Chile, Brazil, Peru, Argentina) with industrial economies (China, U.S., Japan, Korea, Russia), reproducing the Heckscher–Ohlin pattern, in which factor-endowment differences drive trade between resource-abundant and capital-abundant economies feenstra2016advanced. (ref) shows $q_{ij}$ and $q_{ji}$ at the HS2 product level for four pairs. For Luxembourg and Germany, both directions trade similar product mixes and GL is high, and for Chile and China, each direction concentrates in different products and GL is low. These two extremes correspond to the intra- and inter-industry regimes that kim2020measuring obtain by clustering compositional profiles formed directly from observed dollar trade. More details are in Supplementary (ref).
Opposite time-factor trends. Most time factors exhibit gradual growth consistent with expanding global trade, with deviations around major shocks. A subset of components, however, displays opposing time trends despite sharing similar importer sets and product families. We interpret these through three lenses.
(i) Market share reallocation across latent export components. (ref) shows the U.S.-centered component for manufactures (HS84/85/90) into Japan, Korea, and Canada declining in relative prominence, while China- and Germany-centered components into overlapping markets rise. This reflects a loss of relative weight in the non-negative decomposition rather than an absolute decline in U.S. exports; the German time factor also rises more slowly than China's, consistent with U.S. manufacturing losing competitiveness to Chinese imports pierce_schott_2016 and China's rising export market share in medium- and high-tech goods, where Germany has lost ground since 2017 meinen_nagengast_2025.
(ii) Trade-routing and policy changes. A China$\to$Hong Kong textiles component (Supplementary (ref)) drops sharply around 2005, consistent with reduced re-exporting, the routing of Chinese goods through Hong Kong to third markets hktid_mainland_entrepot.
(iii) Supply-chain restructuring. (ref) shows declining exports of textile intermediates (HS52/54/55/56/58/59/60) from Japan, Korea, and Thailand to China alongside rising Chinese exports of overlapping intermediates to Vietnam, Cambodia, and Indonesia. This is consistent with China's growing role as a regional supplier of fibers and fabrics to garment manufacturers across the textiles and clothing value chain wto_gvc_textiles_clothing_2024.
We proposed Bayesian Poisson-randomized Gamma tensor factorization for semi-continuous multiway data which combines low-rank nonnegative CP structure with a compound Poisson--Gamma likelihood and mode-specific dispersion. The model extends Bayesian Poisson tensor factorization for count data schein2015bayesian to semi-continuous observations, complementing existing extensions to binary data zhou2015infinite, compound Poisson exponential-dispersion families BasbugEngelhardt2016HCPF, and bounded responses on [0,1] albert2024doubly. \looseness=-1
A central computational contribution is the treatment of the binary activity indicators governing the zero-positive boundary. Their exact conditional involves a modified Bessel function of the first kind, which is costly to evaluate at scale. We introduce a partially collapsed sampler that avoids marginalizing over exponentially many activation configurations, and derive an asymptotic Bernoulli--Poisson approximation to the Bessel-family conditional that is accurate in the regime most common in practice. This reduces sampling to a single exponential evaluation per entry, enabling GPU-accelerated inference for tensors with tens of millions of entries. Applied to 60 million international trade flows, the model outperforms competing tensor baselines in out-of-sample prediction and recovers 125 interpretable components spanning exporters, importers, products, and time.
The proposed framework opens several future directions. The rank is treated as fixed, where a natural extension is to infer it from the data in a nonparametric Bayesian manner, as in zhao2015bayesian and zhou2015infinite. The temporal factors are modeled without temporal structure, treated as independent and exchangeable across years, and incorporating smoothness priors or regularization across adjacent time points would capture gradual trends between periods. Dyadic covariates such as bilateral distance, tariff schedules, and preferential trade agreements could enter through the latent Poisson rate tensor, connecting the model to structural gravity frameworks and tensor regression with structured coefficient priors guhaniyogi2017bayesian. Beyond monetary-valued data, the framework applies broadly to any zero-inflated, positive-continuous multiway array, including drug-response tensors in pharmacology, species-abundance data in ecology, and healthcare utilization records in epidemiology.
{\selectfont }