# Nonparametric Sparse Tensor Factorization with Hierarchical Gamma Processes

Conor Tillinghast
Department of Mathematics
University of Utah
Salt Lake City, UT 84112
<EMAIL_ADDRESS>
&Zheng Wang
School of Computing
University of Utah
Salt Lake City, UT 84112
<EMAIL_ADDRESS>
&Shandian Zhe
School of Computing
University of Utah
Salt Lake City, UT 84112
<EMAIL_ADDRESS>

###### Abstract

We propose a nonparametric factorization approach for sparsely observed
tensors. The sparsity does not mean zero-valued entries are massive or
dominated. Rather, it implies the observed entries are very few, and even
fewer with the growth of the tensor; this is ubiquitous in practice. Compared
with the existent works, our model not only leverages the structural
information underlying the observed entry indices, but also provides extra
interpretability and flexibility — it can simultaneously estimate a set of
location factors about the intrinsic properties of the tensor nodes, and
another set of sociability factors reflecting their extrovert activity in
interacting with others; users are free to choose a trade-off between the two
types of factors. Specifically, we use hierarchical Gamma processes and
Poisson random measures to construct a tensor-valued process, which can freely
sample the two types of factors to generate tensors and always guarantees an
asymptotic sparsity. We then normalize the tensor process to obtain
hierarchical Dirichlet processes to sample each observed entry index, and use
a Gaussian process to sample the entry value as a nonlinear function of the
factors, so as to capture both the sparse structure properties and complex
node relationships. For efficient inference, we use Dirichlet process
properties over finite sample partitions, density transformations, and random
features to develop an stochastic variational estimation algorithm. We
demonstrate the advantage of our method in several benchmark datasets.

##  1 Introduction

Tensor factorization is a popular tool to analyze multiway data, such as
purchase records among (users, products, sellers) in online shopping. While
numerous algorithms have been proposed, e.g., (Tucker,, 1966; Harshman,, 1970;
Chu and Ghahramani,, 2009; Kang et al.,, 2012), most of these methods are
essentially modeling dense tensors. For instance, many of them (Kolda and
Bader,, 2009; Kang et al.,, 2012; Choi and Vishwanathan,, 2014; Xu et al.,,
2012) use tensor and matrix algebras (Kolda,, 2006) to update the latent
factors and hence demand all the tensor entries be observed (although the
majority are often assumed to be zero-valued). While other methods (Rai et
al.,, 2014; Zhe et al., 2016b, ; Du et al.,, 2018), especially with a Bayesian
formulation, can perform element-wise factorization to handle a small subset
of observed entries, from modeling perspective, they are equivalent to
sampling a full tensor first and then marginalizing out the unobserved entry
values. One might consider adding a Bernoulli likelihood (or a thresholding
function) to model the presence of each entry. However, this still leads to
dense tensors. In theory, these models are random function prior models (Lloyd
et al.,, 2012), and the generated mask tensors are exchangeable arrays.
According to the Aldous-Hoover theorem (Aldous,, 1981; Hoover,, 1979), these
tensors are “either trivially empty or dense”. That is, the present entries
(i.e., with mask one) increases linearly with the size of the tensor.

However, in practice, tensors are often sparsely observed. The ratio of
observed entries is very small, e.g., 0.01%. This ratio can even decrease with
more nodes added. For example, in online shopping, with the growth of users,
products and sellers, the number of actual purchases — i.e., the present
entries — while growing, can take an even smaller percentage of all possible
purchases (entries), i.e., all (user, product, seller) combinations, because
the latter grows much faster. Hence, the (asymptotic) sparse nature in these
real-world applications implies that the existent dense models are
misspecified. Note that the sparsity here does not refer to that the zero-
valued entries are massive or dominated (Choi and Vishwanathan,, 2014; Hu et
al.,, 2015)

To address the model misspecification, Tillinghast and Zhe, (2021) recently
used Dirichlet processes (DPs) (normalized Gamma processes) and GEM
distributions to build two asymptotically sparse tensor models, which have
shown significant advantages over the popular dense models in entry value and
index (link) prediction. However, their model structures severely restrict the
type of factors to be learned — either all the factors are sociability factors
to account for the generation of the interactions, or only one factor is the
sociability factor and the remaining are location factors about the intrinsic
properties of the tensor nodes. The former misses the chance to capture the
intrinsic node properties while the latter can be insufficient to estimate the
sociabilities. Accordingly, both models can limit the flexibly of the factor
estimation and lower the interpretability.

To address these limitations, we propose a new nonparametric factorization
model, which not only guarantees an asymptotic sparsity, but also provides
more flexibility and interpretability. Our model can estimate an arbitrary
number of sociability and location factors so as to fully capture the two
types of tensor node properties. Specifically, we use hierarchical Gamma
processes (H$\Gamma$P) to construct a summation of product measures, witch
which we construct a tensor-variate process. We show that the proportion of
the sampled entries is asymptotically tending to zero with the growth of
tensor size. More important, the $\Gamma$Ps in the top level enable the
sharing of locations across the low-level $\Gamma$Ps. As a result, for each
node, our model can sample not only an arbitrary number of location factors
representing the node’s intrinsic properties, but also an arbitrary number of
sociability factors accounting for its extrovert interaction activity in
different communities. Hence, the results can provide more interpretability
and users are flexible to select a trade-off. Next, for convenient modeling
and inference, we normalize our tensor process to obtain hierarchical
Dirichlet processes so as to sample the observed entry indices; we then use a
Gaussian process to sample the entry value as a nonlinear function of the
latent factors of the associated nodes so as to embed both the sparse
structural knowledge and complex relationships between the tensor nodes. For
efficient inference, we use the Dirichlet process property over finite, non-
overlapping sample space partitions to aggregate the sociabilities of non-
active nodes, probability density transformations, and random Fourier features
to develop a scalable stochastic variational estimation algorithm.

For evaluation, we first ran simulations to confirm that our method can indeed
generate increasingly sparser tensors. In three real-world applications, we
examined the performance in predicting missing entry indices and entry values.
In both tasks, our approach nearly always obtains the best performance, as
compared with the state-of-the-part dense models and (Tillinghast and Zhe,,
2021). The factors estimated by our method reveal clear, interesting
clustering patterns while those estimated from the popular dense models
cannot.

##  2 Background

Notation and Settings. We denote a full $K$-mode tensor by
${\mathcal{Y}}\in\mathbb{R}^{D_{1}\times\ldots\times D_{K}}$, where mode $k$
includes $D_{k}$ nodes. Each entry is indexed by ${\bf
i}=(i_{1},\ldots,i_{K})$ and the entry value by $y_{{\bf i}}$. Practical
tensors are often very sparse. That means most entries of ${\mathcal{Y}}$ are
nonexistent. Given a small set of present entries, $\mathcal{D}=\\{({\bf
i}_{1},y_{{\bf i}_{1}}),\ldots,({\bf i}_{N},y_{{\bf i}_{N}})\\}$ where
$N\ll\prod_{k=1}^{D}D_{k}$, we aim to estimate a set of factors ${\bf
u}_{k}^{j}$ to represent each node $j$ in mode $k$. We denote by
${\mathcal{U}}=\\{{\bf U}^{1},\ldots,{\bf U}^{K}\\}$ all the factors, where
each ${\bf U}^{k}=[{\bf u}^{k}_{1},\ldots,{\bf u}^{k}_{D_{k}}]^{\top}$.

Dense Tensor Models. The classical Tucker (Tucker,, 1966) and
CANDECOMP/PARAFAC (CP) decomposition (Harshman,, 1970) require ${\mathcal{Y}}$
to be fully observed. The Tucker decomposition assumes
${\mathcal{Y}}={\mathcal{W}}\times_{1}{\bf
U}^{1}\times_{2}\ldots\times_{K}{\bf U}^{K}$ where
${\mathcal{W}}\in\mathbb{R}^{r_{1}\times\ldots\times r_{K}}$ is a parametric
core tenor and $\times_{k}$ is the tensor-matrix product at mode $k$ (Kolda,,
2006). CP decomposition is a special case of Tucker decomposition where
${\mathcal{W}}$ is restricted to be diagonal. While many methods can perform
element-wise factorization over a subset of observations, e.g., (Zhe et al.,
2016b, ; Rai et al.,, 2015; Schein et al.,, 2015; Schein et al., 2016b, ; Zhe
and Du,, 2018), from the Bayesian perspective, these methods can be viewed as
first sampling all the entry values of ${\mathcal{Y}}$ given the factors
${\mathcal{U}}$, and then marginalize out the unobserved ones. In other words,
they never consider the sparsity of the observed entry indices (links).
Therefore, these methods are in essence still modeling a dense tensor. In
theory, this has been justified by the random function prior models (Lloyd et
al.,, 2012), where a Bernoulli distribution is used to sample the presence (or
existence) of each entry given the factors,
$p({\mathcal{Z}}|{\mathcal{U}})=\prod_{{\bf i}}p(z_{{\bf
i}}|{\mathcal{U}})=\prod_{\bf i}\text{Bern}(z_{{\bf i}}|f({\bf x}_{{\bf
i}}))$, where ${\mathcal{Z}}$ is a mask tensor, $z_{\bf i}=1$ means the entry
${\bf i}$ is generated, and $z_{\bf i}=0$ means ${\bf i}$ is nonexistent, $f$
is a function of latent factors, e.g., element-wise CP/Tucker form or a
Gaussian process (GP). The mask tensor ${\mathcal{Z}}$ is known to be
exchangeable (Lloyd et al.,, 2012). According to the Aldous-Hoover theorem
(Aldous,, 1981; Hoover,, 1979), the asymptotic111by asymptotic, it means when
we sample a sequence of ${\mathcal{Z}}$’s, with increasing sizes (till
infinity). proportion of the present entries (i.e., the entries with mask 1)
is $\eta=\int\text{Bern}(z_{{\bf i}}|f({\bf x}_{{\bf i}})){\rm
d}p({\mathcal{U}})$. This implies the sampled tensors are either empty
($\eta=0$) or dense ($\eta>0$) almost surely, i.e., the number of sampled
entries grows linearly with the tensor size
($\Theta(\eta\prod_{k=1}^{K}D_{K})$).

Sparse Tensor Model. Recently, Tillinghast and Zhe, (2021) proposed
nonparametric models that guarantee to generate asymptotically sparse tensors,
i.e., the proportion of sampled entries tend to zero with the growth of the
tensor size ($o(\prod_{k=1}^{K}D_{K})$). The idea is to use Gamma processes
($\Gamma$Ps) to construct a tensor-variate random process. This is equivalent
to sampling a Dirichlet process (Ferguson,, 1973) (normalized $\Gamma$P) for
each mode $k$, and combine the weights of each node to sample the present
entries,

$\displaystyle\text{DP}^{k}=\sum\nolimits_{j}w_{j}^{k}\delta_{\bm{\theta}^{k}_{j}},\;\;\;p({\bf
i}_{n})=\prod\nolimits_{k}w_{i_{nk}}^{k},$ (1)

where $1\leq k\leq K$, $1\leq n\leq N$, $w_{j}^{k}$ and $\bm{\theta}^{k}_{j}$
are the DP weights and locations for node $j$ in mode $k$, which naturally
represent the sociability and location factors. The sociability factor
$w_{j}^{k}$ represents the activity of the node in interacting with others to
generate the entry index, i.e., link. The location factors naturally
represents the intrinsic properties of the node. Then $w_{j}^{k}$ and
$\bm{\theta}^{k}_{j}$ are jointly used to the sample the value of the observed
entries $y_{{\bf i}_{n}}$. To incorporate multiple sociability factors,
Tillinghast and Zhe, (2021) proposed a second model that draws multiple
weights for each node via the GEM distribution (Griffiths,, 1980; Engen,,
1975; McCloskey,, 1965). Note that GEM essentially samples DP weights (no
locations). These weights are integrated to sample both the observed entry
indices and entry values.

##  3 Model

Despite being successful, the popular dense tensor models can be misspecified
for the sparse data that are ubiquitous in practice. Although the work of
(Tillinghast and Zhe,, 2021) have overcome this issue, their model structures
severely restrict the type of the factors to be estimated. From (1), we can
only estimate one sociability factor for each node to account for their
complex interactions (i.e., links), which can be quite insufficient. However,
to incorporate multiple sociability factors, we have to use GEM distributions
and drop the location factors (the second model). Accordingly, we cannot
capture any intrinsic properties of the nodes. Note that we cannot simply draw
multiple DPs, because the indices of the nodes in different DPs are not
necessarily the same and we are unable to align their location factors.

To address these limitations, we propose a new sparse tensor factorization
model, which not only guarantees the asymptotic sparsity in tensor generation,
but also is free to estimate an arbitrary number of sociability and location
factors so as to fully capture the extrovert and intrinsic properties of the
tensor nodes. Accordingly, our model can provide more interpretability and
additional flexibly to select the trade-off between the two types of the
factors (e.g., cross-validation).

###  3.1 H$\mathbf{\Gamma}$P Based Sparse Tensor-Variate Process

Specifically, we use hierarchical Gamma processes (H$\Gamma$Ps) to construct a
sparse tensor-variate process. At the top level, we sample a Gamma process
($\Gamma$P) (Hougaard,, 1986) in each mode $k$,

$\displaystyle L^{\alpha}_{k}$
$\displaystyle\sim\Gamma\text{P}(\lambda_{\alpha})\;\;(1\leq k\leq K),$ (2)

where $\lambda_{\alpha}$ is a Lebesgue base measure restricted to
$[0,\alpha]^{R_{1}}=\underbrace{[0,\alpha]\times\ldots\times[0,\alpha]}_{R_{1}\text{times}}(\alpha>0)$.
Accordingly, $L^{\alpha}_{k}$ is a discrete measure,

$\displaystyle L^{\alpha}_{k}$
$\displaystyle=\sum\nolimits_{j=1}^{\infty}w^{\alpha}_{kj}\cdot\delta_{\bm{\theta}^{k}_{j}}$
(3)

where $\delta_{[\cdot]}$ is the Dirac Delta measure, $\\{w^{\alpha}_{kj}>0\\}$
and $\\{\bm{\theta}^{k}_{j}\in[0,\alpha]^{R_{1}}\\}$ are the weights and
locations, which correspond to an infinite number of tensor nodes in mode $k$.
We refer to each $\bm{\theta}^{k}_{j}$ as $R_{1}$ location factors to
represent the intrinsic properties of node $j$ in mode $k$. In the second
level, we use $L_{k}^{\alpha}$ as the base measure to sample $R_{2}$ Gamma
processes in each mode $k$,

$\displaystyle W^{\alpha}_{k,r}\mid
L^{\alpha}_{k}\sim\Gamma\text{P}(L^{\alpha}_{k}),\;\;W^{\alpha}_{k,r}=\sum\nolimits_{j=1}^{\infty}v^{\alpha}_{krj}\cdot\delta_{\bm{\theta}^{k}_{j}},$
(4)

where $1\leq r\leq R_{2}$. Note that due to the common base measure
$L^{\alpha}_{k}$, the locations are shared across all
$\\{W^{\alpha}_{k,r}\\}_{r=1}^{R_{2}}$. Hence, for each node $j$ in mode $k$,
we not only generate $R_{1}$ location factors $\bm{\theta}_{j}^{k}$, but also
$R_{2}$ weights, $\\{v^{\alpha}_{krj}\\}_{r=1}^{R_{2}}$. We refer to the
normalized weights as sociability factors, since we will use them to sample
the observed entry indices (links). More detailed discussions will be given in
Sec. 3.2. Next, we use the H$\Gamma$Ps to construct a product measure sum as
the rate (or mean) measure to sample a Poisson random measure (PRM) (Kingman,,
1992), which represents the sampled tensor entries,

$\displaystyle A$
$\displaystyle=\sum_{r=1}^{R_{2}}W^{\alpha}_{1,r}\times\ldots\times
W^{\alpha}_{K,r},$ $\displaystyle T$
$\displaystyle\mid\\{W^{\alpha}_{k,r}\\}_{1\leq k\leq K,1\leq r\leq
R_{2}}\sim\text{PRM}(A).$ (5)

Accordingly, $T$ has the following form, $T=\sum_{{\bf
i}\in{\mathcal{S}}}c_{{\bf
i}}\cdot\delta_{(\bm{\theta}^{1}_{i_{1}},\ldots,\bm{\theta}^{K}_{i_{K}})}$,
where each point represents an entry index, ${\mathcal{S}}$ is the set of all
the sampled points in the PRM, i.e., entry indices, $c_{{\bf i}}>0$ is the
count of the point ${\bf i}$, and
$(\bm{\theta}^{1}_{i_{1}},\ldots,\bm{\theta}^{K}_{i_{K}})$ is the location of
that point.

In a nut shell, the H$\Gamma$Ps sample an infinite number of nodes in each
mode, then the PRM picks the nodes from each mode to generate the entries
(i.e., links).

######  Lemma 3.1.

For any fixed $\alpha>0$, the number of sampled entries
$N^{\alpha}=|{\mathcal{S}}|$ via (2)(4)(5) is finite almost surely. When
$\alpha\rightarrow\infty$, $N^{\alpha}\rightarrow\infty$ a.s.

We leave the proof in the Appendix. To examine the sparsity, we are interested
in the active nodes that show up in the sampled entry indices ${\mathcal{S}}$.
For example, if an entry $(1,3,5)$ is sampled, then node 1, 3, 5 in mode 1, 2,
3 are active nodes. Denote by $D_{k}^{\alpha}$ the number of (distinct) active
nodes in mode $k$. If we use these active nodes to construct a full tensor,
which we refer to as the active tensor, the size will be
$\prod_{k=1}^{K}D_{k}^{\alpha}$. The asymptotic sparsity means the proportion
of the present entries in the active tensor is small, and the increase of
present entries is much slower than the growth of the active tensor size. This
is guaranteed by

######  Lemma 3.2.

$N^{\alpha}=o(\prod_{k=1}^{K}D^{\alpha}_{k})$ almost surely as
$\alpha\rightarrow\infty$, i.e.,
$\underset{\alpha\rightarrow\infty}{\lim}\frac{N^{\alpha}}{\prod_{k=1}^{K}D^{\alpha}_{k}}=0\;\;\;a.s.$

The proof framework is similar to that in (Tillinghast and Zhe,, 2021), and we
leave the details in the Appendix. We refer to the model defined by (2)(4)(5)
as our Sparse Tensor-variate Process (STP).

###  3.2 Projection on Finite Data

Directly using the STP for modeling and inference is inconvenient, because the
infinite discrete measures of $\Gamma$Ps and PRMs are computationally
challenging. On the other hand, in practice, we always observe a finite number
of entries, $\mathcal{D}=\\{({\bf i}_{1},y_{1}),\ldots,({\bf
i}_{N},y_{N})\\}$. Hence, we can use the standard PRM construction (Kingman,,
1992) to equivalently model the sampling procedure. That is, we normalize the
rate measure $A=\sum_{r=1}^{R_{2}}W^{\alpha}_{1r}\times\ldots\times
W^{\alpha}_{Kr}$ in (5) to obtain a probability measure, and use it to sample
$N$ entry indices (points) independently. To normalize $A$, we need to first
normalize each $\Gamma$P $W^{\alpha}_{kr}(1\leq k\leq K,1\leq r\leq R_{2})$,
which gives a Dirichlet process (DP) (Ferguson,, 1973) with the base measure
as the normalized base measure of $W^{\alpha}_{kr}$. Since the base measure of
$W^{\alpha}_{kr}$ is another $\Gamma$P, $L^{\alpha}_{k}$ (see (2)), its
normalization is a DP again. Hence, we obtain a set of hierarchical Dirichlet
processes (HDPs) (Teh et al.,, 2006),

$\displaystyle G_{k}$
$\displaystyle\sim\text{DP}\big{(}\widetilde{\alpha},\text{Uniform}([0,\alpha]^{R_{1}})\big{)},$
$\displaystyle H_{r}^{k}$ $\displaystyle\sim\text{DP}(\gamma_{r}^{k},G_{k})$
(6)

where $1\leq k\leq K$, $1\leq r\leq R_{2}$,
$\widetilde{\alpha}=\alpha^{R_{1}}$, and
$\gamma_{r}^{k}\sim\text{Gamma}(\cdot|1,\widetilde{\alpha})$. Then the
normalized $A$ is
${\widehat{A}}=\frac{1}{R_{2}}\sum_{r=1}^{R_{2}}H^{1}_{r}\times\ldots\times
H^{K}_{r}$. The sampling of the $N$ observed entry indices is as follows. We
first sample each $G_{k}$ and $H_{r}^{k}$ from (6), which gives

$\displaystyle
G_{k}=\sum\nolimits_{j=1}^{\infty}\beta^{k}_{j}\cdot\delta_{\bm{\theta}^{k}_{j}},\;\;\;H_{r}^{k}=\sum\nolimits_{j=1}^{\infty}\omega^{k}_{rj}\cdot\delta_{\bm{\theta}^{k}_{j}},$
(7)

where each location $\bm{\theta}^{k}_{j}$ is independently sampled from
$\text{Uniform}([0,\alpha]^{R_{1}})$. Then we obtain the sample of the
normalized rate measure as

$\displaystyle{\widehat{A}}=\sum\nolimits_{{\bf
i}=(1,\ldots,1)}^{(\infty,\ldots,\infty)}w_{{\bf
i}}\cdot\delta_{(\bm{\theta}^{1}_{i_{1}},\ldots,\bm{\theta}^{K}_{i_{K}})},$
(8)

where $w_{{\bf
i}}=\frac{1}{R_{2}}\sum_{r=1}^{R_{2}}\prod_{k=1}^{K}\omega_{ri_{k}}^{k}$ which
is essentially a probability measure over all possible tensor entries, i.e.,
$\sum_{{\bf i}}w_{{\bf i}}=1$. We then sample each observed entry ${\bf
i}_{n}\sim{\widehat{A}}$, and hence the probability is

$\displaystyle p({\mathcal{S}})=\prod\nolimits_{n=1}^{N}p({\bf
i}_{n})=\prod\nolimits_{n=1}^{N}w_{{\bf i}_{n}}.$ (9)

Now, it can be seen clearly that in (6), for each node $j$ in mode $k$, we
sample the location $\bm{\theta}_{j}^{k}$ and a set of weights
$\bm{\nu}^{k}_{j}=[\omega^{k}_{1j},\ldots,\omega^{k}_{R_{2}j}]$ from $R_{2}$
DPs. From (8), we can see that these weights represent the activity of the
node interacting with other nodes. Each weight naturally corresponds to the
activity in one community/group. Hence, we refer to $\bm{\nu}^{k}_{j}$ as the
sociability factors of the node in $R_{2}$ overlapping communities. The
location factors $\bm{\theta}_{j}^{k}$ naturally represent the node’s
intrinsic properties. Thanks to the HDP (H$\Gamma$P) structure, the locations
can be shared in an arbitrary number of low-level DPs ($\Gamma$Ps). That
means, we are free to choose the number of sociability factors and location
factors. Hence, it provides not only more interpretability but also an extra
flexibility to select their trade-off.

Given the sampled entries ${\mathcal{S}}$, we then sample the entry values
${\bf y}=[y_{{\bf i}_{1}},\ldots,y_{{\bf i}_{N}}]$. In this work, we mainly
consider continuous values. It is straightforward to extend our model for
other value types. We sample each value from a Gaussian noise model,

$\displaystyle p(y_{{\bf i}_{n}}|{\mathcal{S}})=\mathcal{N}(y_{{\bf
i}_{n}}|f({\bf x}_{{\bf i}_{n}}),\sigma^{2}),$ (10)

where ${\bf x}_{\bf i}=[{\bf u}^{1}_{i_{n1}};\ldots;{\bf u}^{K}_{i_{nK}}]$ are
the factors of all the nodes associated with entry ${\bf i}_{n}$, each ${\bf
u}^{k}_{i_{nk}}=[\bm{\theta}^{k}_{i_{nk}};\bm{\nu}^{k}_{i_{nk}}]$, i.e.,
including $R_{1}$ location and $R_{2}$ sociability factors, and $f(\cdot)$ is
a latent factorization function. To flexibly estimate $f(\cdot)$ so as to
capture the complex relationships between the tensor nodes (in terms of their
factor representations), we assign a Gaussian process (Rasmussen and
Williams,, 2006) prior over $f(\cdot)$. Accordingly, the function values ${\bf
f}=[f({\bf x}_{{\bf i}_{1}}),\ldots,f({\bf x}_{{\bf i}_{N}})]$ follow a joint
Gaussian distribution, $p({\bf f})=\mathcal{N}({\bf f}|{\bf 0},{\bf K})$,
where ${\bf K}$ is an $N\times N$ kernel (covariance) matrix, each element
$[{\bf K}]_{st}=\kappa({\bf x}_{{\bf i}_{s}},{\bf x}_{{\bf i}_{t}})$, and
$\kappa(\cdot,\cdot)$ is a kernel function.

From (9) and (10), we can see via coupling the HDPs and GP, both the
structural knowledge in sparse entry indices and complex relationships of the
tensor nodes can be captured and encoded into the sociability and location
factors.

##  4 Algorithm

The estimation of our model is still challenging. First, we have to deal with
infinite discrete measures (see (7)) from HDPs. Second, the exact GP prior
includes an $N\times N$ covariance matrix. When $N$ is large, the computation
is very costly or even infeasible. To address these issues, we use the DP
property over sample partitions (Ferguson,, 1973), probability density
transformations, and random Fourier features (Rahimi et al.,, 2007; Lázaro-
Gredilla et al.,, 2010) to develop an efficient, scalable variational learning
algorithm.

Specifically, since we only need to estimate the latent factors for finite
active nodes, i.e., the nodes appearing in the observed entries (the infinite
inactive nodes do not have any observed data), we can marginalize (7) into
finite measures. Suppose we have $D_{k}$ active nodes in each mode $k$.
Without loss of generality, we can index these nodes by $1,\ldots,D_{k}$, and
the remaining inactive nodes by $D_{k+1},\ldots,\infty$. In the first level,
for each $G_{k}$, we consider
$\bm{\beta}^{k}=[\beta^{k}_{1},\ldots,\beta^{k}_{D_{k}},\beta^{k}_{>D_{k}}]$
where $\beta^{k}_{>D_{k}}=\sum_{j=D_{k}+1}^{\infty}\beta^{k}_{j}$, i.e., the
DP weights for every active node and the aggregated weight of all the inactive
nodes. In the second level, we partition the sample space of each $H^{k}_{r}$
into $D_{k}+1$ subsets,
$\\{\\{\bm{\theta}_{1}^{k}\\},\ldots,\\{\bm{\theta}_{D_{k}}^{k}\\},\\{\text{Others}\\}\\}$.
Since $H^{k}_{r}\sim\text{DP}(\gamma_{r}^{k},G_{k})$, its measures of the
$K+1$ subsets will jointly follow a Dirichlet distribution. The values of
these measures are
$\bm{\omega}^{k}_{r}=\\{\omega_{r1}^{k},\ldots,\omega_{rD_{k}}^{k},\omega_{r>D_{k}}^{k}\\}$,
where $\omega^{k}_{r>D_{k}}=\sum_{j=D_{k}+1}^{\infty}\omega^{k}_{rj}$, namely,
the sociability of each active node and the aggregated sociability of all the
inactive nodes. The Dirichlet distribution is parameterized by
$[\gamma_{r}^{k}G_{k}(\\{\bm{\theta}_{1}^{k}\\}),\ldots,\gamma_{r}^{k}G_{k}(\\{\bm{\theta}_{D_{k}}^{k}\\}),\gamma_{r}^{k}G_{k}(\\{\text{Others}\\})]=\gamma_{r}^{k}\bm{\beta}^{k}$.
Therefore, we have

$\displaystyle
p(\bm{\omega}^{k}_{r}|\gamma^{k}_{r},\bm{\beta}^{k})=\text{Dir}(\bm{\omega}^{k}_{r}|\gamma^{k}_{r}\bm{\beta}^{k})$
$\displaystyle=\frac{\Gamma(\gamma_{r}^{k})}{\prod_{j}\Gamma(\gamma_{r}^{k}\beta^{k}_{j})}\prod\nolimits_{j}\left(\omega^{k}_{rj}\right)^{\gamma_{r}^{k}\beta^{k}_{j}-1},$
(11)

where $\Gamma(\cdot)$ is a Gamma function. Given the sociabilities of the
active nodes, we have already been able to sample the observed entry indices
from (9). The remaining problem is how to obtain the prior of $\bm{\beta}^{k}$
for each $G_{k}$. From the stick-breaking construction (Ishwaran and James,,
2001), we have
$\xi^{k}_{j}\sim\text{Beta}(1,\widetilde{\alpha}),\beta^{k}_{j}=\xi^{k}_{j}\prod\nolimits_{t=1}^{j-1}(1-\xi^{k}_{t})\;(1\leq
j\leq\infty)$, from which we can obtain a reverse mapping,

$\displaystyle(\xi^{k}_{1},\ldots,\xi^{k}_{D_{k}})=(\frac{\beta^{k}_{1}}{\Lambda^{k}_{1}},\ldots,\frac{\beta^{k}_{D_{k}}}{\Lambda^{k}_{D_{k}}}),$
(12)

where
$\Lambda^{k}_{j}=1-\sum_{t=1}^{j-1}\beta^{k}_{t}=\prod_{t=1}^{j-1}(1-\xi^{k}_{t})$.
Therefore, we can use probability density transformation to obtain the prior
distribution of $\bm{\beta}^{k}$. Since each $\Lambda_{j}^{k}$ is only
determined by $\\{\beta_{1}^{k},\ldots,\beta_{j-1}^{k}\\}$, the Jacobian ${\bf
J}$ is a lower-triangular matrix and $|{\bf
J}|=\prod_{j=1}^{D_{k}}\frac{1}{\Lambda^{k}_{j}}$. Accordingly, we can derive

$\displaystyle
p(\bm{\beta}^{k})=\prod\nolimits_{j=1}^{D_{k}}\text{Beta}(\frac{\beta^{k}_{j}}{\Lambda^{k}_{j}}|1,\widetilde{\alpha})\frac{1}{\Lambda^{k}_{j}}.$
(13)

Note that since $\beta^{k}_{>D_{k}}=1-\sum_{j=1}^{D_{k}}\beta^{k}_{j}$, it
does not need to be explicitly computed in the prior. To conveniently estimate
each $\bm{\beta}^{k}$ and $\bm{\omega}^{k}_{r}$ and to avoid incorporating
additional constraints (nonnegative, summation is one), we parameterize
$\bm{\beta}^{k}=\text{softmax}({\widetilde{\bm{\beta}}}^{k})$ and
$\bm{\omega}^{k}_{r}=\text{softmax}({\widetilde{\bm{\omega}}}^{k}_{r})$ and
estimate the free parameters
$\\{{\widetilde{\bm{\beta}}}^{k},{\widetilde{\bm{\omega}}}^{k}_{r}\\}$
instead.

Next, to scale to a large number of observed entries, we use random Fourier
features to construct a sparse GP approximation (Rahimi et al.,, 2007; Lázaro-
Gredilla et al.,, 2010). Specifically, according to Bochner’s theorem (Rudin,,
1962), any stationary kernel $\kappa({\bf b}_{1},{\bf b}_{2})=\kappa({\bf
b}_{1}-{\bf b}_{2})$ can be viewed as an expectation, $\kappa({\bf b}_{1}-{\bf
b}_{2})=\kappa({\bf 0})\mathbb{E}_{p({\bf z})}[e^{i{\bf z}^{\top}{\bf
a}_{1}}\big{(}e^{i{\bf z}^{\top}{\bf a}_{2}}\big{)}^{\dagger}]$, where
$\dagger$ is the complex conjugate, $p({\bf z})=\mathcal{N}({\bf z}|0,\tau{\bf
I})$ is computed from the Fourier transform of $\kappa(\cdot)$. Hence, we can
draw $M$ frequencies ${\bf Z}=\\{{\bf z}_{1},\ldots,{\bf z}_{M}\\}$ from
$p({\bf z})$ to construct an Monte-Carlo approximation, $\kappa({\bf
b}_{1},{\bf b}_{2})\approx\frac{\kappa({\bf 0})}{M}\sum_{m=1}^{M}e^{i{\bf
s}_{m}^{\top}{\bf b}_{1}}\big{(}e^{is_{m}^{\top}{\bf
b}_{2}}\big{)}^{\dagger}=\frac{\kappa({\bf 0})}{M}\bm{\phi}({\bf
b}_{1})^{\top}\bm{\phi}({\bf b}_{2})$, where

$\bm{\phi}({\bf b})=[\cos({\bf z}_{1}^{\top}{\bf b}),\sin({\bf
z}_{1}^{\top}{\bf b}),\ldots,\cos({\bf z}_{M}^{\top}{\bf b}),\sin({\bf
z}_{M}^{\top}{\bf b})].$

Here, $\bm{\phi}(\cdot)$ is a $2M$ dimensional nonlinear feature mapping,
which is called random Fourier features. We can then use a Bayesian linear
regression model over the random Fourier features to approximate the GP.
Specifically, we use the RBF kernel, $\kappa({\bf b}_{1},{\bf
b}_{2})=\exp(-\frac{1}{2}\tau\|{\bf b}_{1}-{\bf b}_{2}\|^{2})$. The
corresponding frequency distribution $p({\bf z})=\mathcal{N}({\bf z}|{\bf
0},\tau{\bf I})$. We sample $M$ frequencies, ${\bf Z}=\\{{\bf
z}_{1},\ldots,{\bf z}_{M}\\}$, from $p({\bf
Z})=\prod_{m=1}^{M}\mathcal{N}\big{(}{\bf z}_{m}|{\bf 0},\tau{\bf I}\big{)}$,
and a weight vector ${\bf g}$ from $p({\bf g})=\mathcal{N}({\bf g}|{\bf
0},\frac{1}{M}{\bf I})$. We then represent the function $f(\cdot)$ in (10) by

$\displaystyle f({\bf x}_{{\bf i}_{n}})=\bm{\phi}({\bf x}_{{\bf
i}_{n}})^{\top}{\bf g},$ (14)

where the input ${\bf x}_{{\bf i}_{n}}$ includes the factors of all the nodes
associated with ${\bf i}_{n}$. Since the sociabilities are often very small
and close to zero, we use their softmax parameters
$\\{{\widetilde{\omega}}^{k}_{i_{nk}}\\}$ in ${\bf x}_{{\bf i}_{n}}$. When we
marginalize out the weight vector ${\bf g}$, the joint Gaussian distribution
is recovered, and the kernel takes the form of the Monte-Carlo approximation.
We will keep ${\bf g}$ to prevent computing the full covariance matrix, and to
scale to a large number of observed entries.

Combining (13)(11)(9)(14)(10), we obtain the joint probability of our model,

$\displaystyle
p(\text{Joint})=\prod\nolimits_{k=1}^{K}p(\bm{\beta}^{k})\left(\prod\nolimits_{j=1}^{D_{k}}\text{Uniform}(\bm{\theta}^{k}_{j}|[0,\alpha]^{R_{1}})\right)$
$\displaystyle\cdot\left(\prod\nolimits_{r=1}^{R_{2}}\text{Gamma}(\gamma^{k}_{r}|1,\widetilde{\alpha})\text{Dir}(\bm{\omega}^{k}_{r}|\gamma^{k}_{r}\bm{\beta}^{k})\right)\mathcal{N}({\bf
g}|{\bf 0},\frac{1}{M}{\bf I})$
$\displaystyle\cdot\prod\nolimits_{m=1}^{M}\mathcal{N}({\bf z}_{m}|{\bf
0},\tau{\bf I})\prod\nolimits_{n=1}^{N}w_{{\bf i}_{n}}\cdot\mathcal{N}(y_{{\bf
i}_{n}}|\bm{\phi}({\bf x}_{{\bf i}_{n}})^{\top}{\bf g},\sigma^{2}).$

We use variational inference (Wainwright and Jordan,, 2008) for model
estimation. We introduce a variational posterior $q({\bf g})=\mathcal{N}({\bf
g}|\bm{\mu},{\bf L}{\bf L}^{\top})$, and construct a variational evidence
lower bound (ELBO), ${\mathcal{L}}=\mathbb{E}_{q}\left[{\log
p(\text{Joint})}-\log{q({\bf g})}\right]$. We maximize the ELBO to estimate
$q({\bf g})$ and all the other parameters, including
$\\{{\widetilde{\bm{\beta}}}^{k}\\}$, $\\{\bm{\theta}_{j}^{k}\\}$,
$\\{\gamma_{r}^{k},{\widetilde{\bm{\omega}}}^{k}_{r}\\}$, ${\bf Z}$, etc. Due
to the additive structure over the observed entries in ${\mathcal{L}}$, we can
use mini-batch stochastic optimization for efficient and scalable estimation.

Complexity. The time complexity of our inference is
$O(NM^{2}+R\sum_{k=1}^{K}D_{k})$ where $R=R_{1}+R_{2}$, and $N$ is the number
of observed entries. Since $M$ is fixed and $M\ll N$, the time complexity is
linear in $N$. The space complexity is $O(M^{2}+M+\sum_{k}D_{k}R)$, which is
to store the latent factors, the frequencies, and the variational posterior.

##  5 Related Work

Numerous tensor factorization methods have been developed, e.g., (Shashua and
Hazan,, 2005; Chu and Ghahramani,, 2009; Sutskever et al.,, 2009; Acar et
al.,, 2011; Hoff,, 2011; Kang et al.,, 2012; Yang and Dunson,, 2013; Rai et
al.,, 2014; Choi and Vishwanathan,, 2014; Hu et al.,, 2015; Zhao et al.,,
2015; Rai et al.,, 2015; Du et al.,, 2018). While most methods adopt the
multilinear forms of the classical Tucker (Tucker,, 1966) and CP
decompositions (Harshman,, 1970), a few GP based factorization models (Xu et
al.,, 2012; Zhe et al.,, 2015; Zhe et al., 2016a, ; Zhe et al., 2016b, ; Pan
et al., 2020b, ) were recently proposed to capture the nonlinear relationships
between the tensor nodes.

Most methods model dense tensors. From the Bayesian viewpoint, they all belong
to random function prior models (Lloyd et al.,, 2012) on exchangeable arrays.
Recently, Tillinghast and Zhe, (2021) used Gamma processes ($\Gamma$Ps), a
commonly used completely random measure (CRM) (Kingman,, 1967, 1992; Lijoi et
al.,, 2010), to construct a Poisson random measure that represents
asymptotically sparse tensors. Their work can be viewed as an extension of the
pioneer works Caron and Fox, (2014); Williamson, (2016); Caron and Fox, (2017)
in sparse random graph modeling. However, their single-level $\Gamma$Ps and
accordingly single-level Dirichlet processes (DPs) for finite projection
impose a severe restriction on the type of the learned factors. For each
tensor node, either all the factors are exclusively sociability factors or
only one factor can be the sociability factor. The former misses the location
factors that represent the intrinsic properties of the node, and the latter is
often inadequate to capture the extrovert interaction activity across
different communities. To overcome this limitation, our model uses
hierarchical Gamma processes (H$\Gamma$Ps) to construct the sparse tensor
model. Due to the sharing of the locations of $\Gamma$Ps in the second level,
our model does not have the node alignment issue and is able to estimate an
arbitrary set of sociability and location factors to fully capture the
extrovert activities and inward attributes respectively, hence giving more
interpretability and an additional flexibility to select the number of two
types of factors. We couple the normalized H$\Gamma$Ps and GPs to sample both
the sparse entry indices (hyperlinks) and entry values, so as to absorb both
the sparse structural knowledge and complex relationships into the factor
estimates. We point out that the normalized H$\Gamma$P is HDP (Teh et al.,,
2006), an important extension of DP (Ferguson,, 1973). HDPs are popular in
nonparametric mixture modeling in text mining (Hong and Davison,, 2010; Wang
et al.,, 2011; McFarland et al.,, 2013), which is often referred to as topic
models. For inference, Tillinghast and Zhe, (2021) directly used the stick-
breaking construction of DPs, which, however, are not available for our model,
because it is intractable for HDP inference. Instead, we used DP measures over
a finite partition of the sample space to integrate the measures of non-active
nodes, and density transformation to obtain the density of finite point
masses. The strategy was also used in HDP mixture models, e.g., (Liang et
al.,, 2007; Bryant and Sudderth,, 2012). Schein et al., 2016a  used Gamma and
Poisson distributions to develop a Bayesian Poisson Tucker decomposition
model. They discussed that if one uses a Gamma process instead of the Gamma
distribution to sample the latent core tensor, the core tensor will be
asymptotically sparse. However, it is unclear if the tensor itself will also
be sparse.

Another recent line of research (Zhe and Du,, 2018; Pan et al., 2020a, ; Wang
et al.,, 2020) uses GPs to construct different point processes for temporal
events modeling and embedding the participants of these events. While these
problems are formulated as decomposing a special type of tensors, i.e., event-
tensors, where each entry is a sequence of temporal events (rather than
numerical values), these models focus on building expressive non-Poisson
processes to capture complex temporal dependencies between the events, e.g.,
the Hawkes processes with local or global triggering effects (Zhe and Du,,
2018; Pan et al., 2020a, ), the non-Poisson, non-Hawkes process (Wang et al.,,
2020) that captures the long-term, short-term, triggering and inhibition
effects.

##  6 Experiment

###  6.1 Simulation

First, we checked if our H$\Gamma$P based tensor process can indeed produce
sparse tensors. We used our model to generate a series of tensors with
increasingly more present entries, and examined how the proportion of the
present entries varied accordingly. The technique details about sampling is
given in the Appendix. We fixed $R_{1}=1$ and varied the number of sociability
factors ($R_{2}$) from $\\{1,3,5\\}$, and $\alpha$ from $[1,15]$ . Then we
looked at the active nodes in the sampled entries (the nodes present in the
entry indices), with which we can construct an active tensor that only
includes the active nodes in each mode. We calculated the ratio between the
number of sampled entries and the size of the active tensor. The ratio
indicates the sparsity. For each particular $\alpha$, we ran the simulation
for $100$ times and calculated the average ratio, which gives a reliable
estimate of the sparsity. We show how the ratio varied along with the size of
the active tensor in Fig. 1 a. As we can see, the proportion of the sampled
entries drops rapidly with the growth of the active tensor and is tending to
zero (the gap between the curves and the horizontal axis is decreasing). This
is consistent with the (asymptotic) sparse guarantee in Lemma 3.2.

We contrasted our model with two popular dense factorization models, CP-Bayes
(Zhao et al.,, 2015; Du et al.,, 2018) and GPTF (Zhe et al., 2016b, ; Pan et
al., 2020b, ), where we first sampled the latent factors ${\mathcal{U}}$ from
the standard Gaussian prior distribution, and then independently sampled the
presence of each entry ${\bf i}$ from $p(z_{\bf
i}|{\mathcal{U}})=\text{Bern}\big{(}z_{\bf i}|f({\bf x}_{i})\big{)}$, where
${\bf x}_{i}=[{\bf u}_{i_{1}}^{1};\ldots;{\bf u}_{i_{K}}^{K}]$, and $f(\cdot)$
is the factorization function. For CP-Bayes, $f(\cdot)$ is the element-wise
version of the CP form, while for GPTF, $f(\cdot)$ is assigned a GP prior, and
we used a sparse approximation with random Fourier feature (the same as in our
model inference) to sample a large number of entries. We gradually increased
the number of nodes in each mode so as to grow the tensor size. We examined
the proportion of the present entries accordingly (i.e., $z_{\bf i}=1$). At
each step, we ran the sampling procedure for $100$ times, and computed the
average proportion. As shown in Fig. 1b, CP-Bayes and GPTF both generated
dense tensors, where the proportion is almost constant, reflecting that the
present entry number grows linearly with the tensor size. In Fig. 1c-e, we
showcase exemplar tensors generated by each model. We can see that the tensors
entries sampled from CP-Bayes and GPTF (Fig. 1d and e) are much denser than
ours (Fig. 1c) and more uniform. This is consistent with the fact that their
priors are exchangeable and symmetric, and each entry is sampled with the same
probability (see Sec. 2).

|  |  |  |
---|---|---|---|---

(a) Sparsity: Ours

(b) Sparsity: Dense models

(c) Our model

(d) CP-Bayes

(e) GPTF

Figure 1: The ratio between the sampled entries and tensor size (a, b), and
exemplar tensors generated by each model (c, d, e), of size $10\times 10\times
10$. The legend in (a) indicates the number of sociability factors.

|  |  |  |
---|---|---|---|---
|  |  |  |
|  |  |  |

(a) Alog

(b) MovieLens

(c) SG

(d) Alog-mode 2

(e) Alog

(f) MovieLens

(g) SG

(h) MovieLens-mode 2

(i) Alog-Link

(j) MovieLens-Link

(k) SG-Link

(l) SG-mode 3

Figure 2: Prediction accuracy of entry values (a-c, e-g) and entry indices
(i-k), and the structures of the estimated factors by our method (d, h, l).
Note that, the results of some approaches were much worse than the others and
were not included, e.g., the broken curves of CP-WOPT and CP-ALS in (a);
NEST-1 and NEST-2 overlap in (j, k), and NEST-2 overlaps with Ours in (i).

###  6.2 Predictive Performance in Practical Applications

Next, we evaluated the predictive performance in three real-world benchmark
datasets: (1) Alog (Zhe et al., 2016b, ), extracted from a file access log,
depicting the access frequency among users, actions, and resources, of size
$200\times 100\times 200$. The value of each present entry is the logarithm of
the access frequency. The proportion of present entries is $0.3\%$. (2)
MovieLens (https://grouplens.org/datasets/movielens/100k/), a three-mode
(user, movie, time) tensor, of size $1000\times 1700\times 31$. There are 100K
present entries, taking $0.19\%$ of the whole tensor. The entry values are
movie ratings in $[0,10]$. We divided every value by $10$ to normalized them
in $[0,1]$. (3) SG (Li et al.,, 2015; Liu et al.,, 2019), extracted from
Foursquare Singapore data, representing (user, location, point-of-interest)
check-ins, of size $2321\times 5596\times 1600$. The value of each present
entry is the normalized check-in frequency in $[0,1]$. The number of present
entries is $105,764$, taking $0.0005\%$ of whole tensor.

Methods and Settings. We compared with the following state-of-the-art tensor
factorization methods. (1) CP-Bayes (Zhao et al.,, 2015; Du et al.,, 2018), a
Bayesian CP factorization model with Gaussian likelihood. (2) CP-ALS (Bader et
al.,, 2015), CP factorization using alternating least squares to update the
latent factors. (3) CP-WOPT (Acar et al.,, 2011), a fast CP factorization
method using conjugate gradient descent to optimize the latent factors. (4)
P-Tucker (Oh et al.,, 2018), a highly efficient Tucker factorization algorithm
that updates the factor matrices row-wisely in parallel. (5) GPTF (Zhe et al.,
2016b, ; Pan et al., 2020b, ), nonparametric tensor factorization based on
GPs. It is the same as our modeling the entry value as a latent function of
the latent factors (see (10)), except that it samples the latent factors from
a standard Gaussian prior. The same sparse GP approximation, i.e., random
Fourier features, was employed for scalable model estimation. (6) NEST-1 and
(7) NEST-2, which are the two models proposed by Tillinghast and Zhe, (2021).
NEST-1 uses a single DP in each mode to sample the sparse tensors, and can
only estimate one sociability factor for each node; the remaining factors are
location factors (from DP locations). NEST-2 uses multiple GEM distributions
to draw DP weights for sparse tensor modeling, hence all the factors are
sociability factors and there are no location factors. We implemented our
method with PyTorch (Paszke et al.,, 2019). CP-Bayes, GPTF, NEST-1 and NEST-2
were implemented with TensorFlow (Abadi et al.,, 2016). Since all these
methods used stochastic mini-batch optimization, we chose the learning rate
from $\\{10^{-4},2\times 10^{-4},5\times 10^{-4},10^{-3},5\times
10^{-3},10^{-2}\\}$, and set the mini-batch size to $200$ for Alog and
MovieLens, and $512$ for SG. We ran 700 epochs on all the three datasets,
which is enough for convergence. We used the original MATLAB implementation of
CP-ALS and CP-WOPT (http://www.tensortoolbox.org/), and C++ implementation of
P-Tucker (https://github.com/sejoonoh/P-Tucker), and their default settings.
Note that CP-ALS needs to fill nonexisting entries with zero-values to
manipulate the whole tensor.

Missing Entry Value Completion. We first tested the accuracy of predicting the
missing entry values. To this end, we randomly split the existent entries in
each dataset into $80\%$ for training and $20\%$ for test. We then ran each
method, and evaluated the Mean Square Error (MSE) and Mean Absolute Error
(MAE). We varied the number of factors $R$ from $\\{3,5,7,9,11\\}$, and for
each setting, we ran the experiment for five times. Since our approach can
flexibly estimate the two types of factors, i.e., sociability and location
factors, to identify an appropriate trade-off, we ran an extra validation in
the training set to select $R_{1}$ and $R_{2}$ ($R_{1}+R_{2}=R$) (In the
Appendix, we show how the performance of our method varies along with all
possible combinations of $R_{1}$ and $R_{2}$). Note that while NEST also
estimates the two types of factors, their number choice is strictly fixed: for
NEST-1, $R_{1}=1,R_{2}=R-1$, and NEST-2, $R_{1}=R,R_{2}=0$. We report the
average MSE, average MAE and their standard deviations in Fig. 4 a-c and e-g.
As we can see, our method nearly always obtains the best performance, except
when $R=3,5$, our method is slightly worse than NEST-1 and NEST-2. In general,
our method and NEST outperforms the dense models by a large margin, which
confirm the power of the sparse tensor modeling. Our method further improves
upon NEST in most cases, especially when $R$ is relatively bigger ($R=7,9,11$;
see Fig. 4a,b,e,f and g), which demonstrates the advantage of our method that
can estimate an arbitrary number of sociability and location factors, and
allow the selection of their trade-off.

Missing Entry (Link) Prediction. Next, we examined our model in predicting
missing entry indices. To this end, from each dataset, we randomly sampled
$80\%$ of existent entries, and used their indices for training. Then we used
the remaining $20\%$ existent entries plus ten-times nonexistent entries
(randomly generated) for test. We did not incorporate all the nonexistent
entries for test in order to prevent them from dominating the result — they
are too many. We compared with CP-Bayes, GPTF, P-Tucker, NEST-1 and NEST-2.
Note that for CP-Bayes and GPTF, we used a Bernoulli distribution to sample
the presence (or existence) of an entry, for P-Tucker we used a zero-
thresholding function. We computed the area under ROC curve (AUC) of
predictions on the test entries. We repeated the experiment for five times and
report the average AUC and its standard deviation in Fig. 4 i-k. We can see
that our method shows much better prediction accuracy than all the dense
models. On Alog and SG dataset, the prediction accuracy of our method is close
to or slightly better than NEST (i.e., Fig. 4i and g); However, on MovieLens
dataset, our method significantly outperforms NEST by a large margin. See Fig.
4j; note that the results of NEST-1 and NEST-2 overlap. It is worth noting
that the competing dense methods perform much worse on SG than on Alog and
MovieLens. This might because SG is much sparser (0.0005% vs. 0.3% and 0.19%),
and the model misspecification is more severe, leading to much inferior
accuracy.

###  6.3 Pattern Discovery

Finally, we investigated if our method can discover hidden patterns within the
tensor nodes, as compared with the popular dense factorization models. To this
end, we set $R=11$ where $R_{1}=5$ and $R_{2}=6$, and ran our method on all
the three datasets. We then used Principal Component Analysis (PCA) to project
the learned factors in each mode onto a plain. The positions of the points
represent the first and second principal components. To find the patterns, we
ran the k-means algorithm and filled the cluster regions with different
colors. We used the elbow method (Ketchen and Shook,, 1996) to select the
cluster number. In Fig. 4d, h and l, we show the results of the second mode on
Alog and MovieLens and the third mode on SG. Each point corresponds to a
tensor node. As we can see, these results reflect clear, interesting
clustering structures. By contrast, the factors estimated by CP-Bayes and
GPTF, as shown in Fig. 2 of the Appendix, do not exhibit patterns and their
principal components are distributed like a symmetric Gaussian. While the
datasets have been completely anonymized and we are unable to look into the
meaning of the patterns discovered by our method, together this has
demonstrated that our model, by capturing the sparse structural information
within the entry indices, can potentially discover more interesting and
important knowledge and further enhance the interpretability.

##  7 Conclusion

We have presented a novel nonparametric tensor factorization method based on
hierarchical Gamma processes. Our model not only can sample sparse tensors to
match the nature of data sparsity in practice, but also can estimate an
arbitrary number of sociability and location factors so as to discover
interesting and important patterns.

## References

  * Abadi et al., (2016) Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al. (2016).  Tensorflow: A system for large-scale machine learning.  In 12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16), pages 265–283.
  * Acar et al., (2011) Acar, E., Dunlavy, D. M., Kolda, T. G., and Morup, M. (2011).  Scalable tensor factorizations for incomplete data.  Chemometrics and Intelligent Laboratory Systems, 106(1):41–56.
  * Aldous, (1981) Aldous, D. J. (1981).  Representations for partially exchangeable arrays of random variables.  Journal of Multivariate Analysis, 11(4):581–598.
  * Bader et al., (2015) Bader, B. W., Kolda, T. G., et al. (2015).  Matlab tensor toolbox version 2.6.  Available online.
  * Bryant and Sudderth, (2012) Bryant, M. and Sudderth, E. (2012).  Truly nonparametric online variational inference for hierarchical dirichlet processes.  Advances in Neural Information Processing Systems, 25:2699–2707.
  * Caron and Fox, (2014) Caron, F. and Fox, E. B. (2014).  Sparse graphs using exchangeable random measures.  arXiv preprint arXiv:1401.1137.
  * Caron and Fox, (2017) Caron, F. and Fox, E. B. (2017).  Sparse graphs using exchangeable random measures.  Journal of the Royal Statistical Society. Series B, Statistical Methodology, 79(5):1295.
  * Choi and Vishwanathan, (2014) Choi, J. H. and Vishwanathan, S. (2014).  Dfacto: Distributed factorization of tensors.  In Advances in Neural Information Processing Systems, pages 1296–1304.
  * Chu and Ghahramani, (2009) Chu, W. and Ghahramani, Z. (2009).  Probabilistic models for incomplete multi-dimensional arrays.  AISTATS.
  * Du et al., (2018) Du, Y., Zheng, Y., Lee, K.-c., and Zhe, S. (2018).  Probabilistic streaming tensor decomposition.  In 2018 IEEE International Conference on Data Mining (ICDM), pages 99–108. IEEE.
  * Engen, (1975) Engen, S. (1975).  A note on the geometric series as a species frequency model.  Biometrika, 62(3):697–699.
  * Ferguson, (1973) Ferguson, T. S. (1973).  A bayesian analysis of some nonparametric problems.  The annals of statistics, pages 209–230.
  * Griffiths, (1980) Griffiths, R. C. (1980).  Lines of descent in the diffusion approximation of neutral wright-fisher models.  Theoretical population biology, 17(1):37–50.
  * Harshman, (1970) Harshman, R. A. (1970).  Foundations of the PARAFAC procedure: Model and conditions for an”explanatory”multi-mode factor analysis.  UCLA Working Papers in Phonetics, 16:1–84.
  * Hoff, (2011) Hoff, P. (2011).  Hierarchical multilinear models for multiway data.  Computational Statistics & Data Analysis, 55:530–543.
  * Hong and Davison, (2010) Hong, L. and Davison, B. D. (2010).  Empirical study of topic modeling in twitter.  In Proceedings of the first workshop on social media analytics, pages 80–88.
  * Hoover, (1979) Hoover, D. N. (1979).  Relations on probability spaces and arrays of random variables.  Preprint, Institute for Advanced Study, Princeton, NJ, 2.
  * Hougaard, (1986) Hougaard, P. (1986).  Survival models for heterogeneous populations derived from stable distributions.  Biometrika, 73(2):387–396.
  * Hu et al., (2015) Hu, C., Rai, P., and Carin, L. (2015).  Zero-truncated poisson tensor factorization for massive binary tensors.  In UAI.
  * Ishwaran and James, (2001) Ishwaran, H. and James, L. F. (2001).  Gibbs sampling methods for stick-breaking priors.  Journal of the American Statistical Association, 96(453):161–173.
  * Kang et al., (2012) Kang, U., Papalexakis, E., Harpale, A., and Faloutsos, C. (2012).  Gigatensor: scaling tensor analysis up by 100 times-algorithms and discoveries.  In Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 316–324. ACM.
  * Ketchen and Shook, (1996) Ketchen, D. J. and Shook, C. L. (1996).  The application of cluster analysis in strategic management research: an analysis and critique.  Strategic management journal, 17(6):441–458.
  * Kingman, (1967) Kingman, J. (1967).  Completely random measures.  Pacific Journal of Mathematics, 21(1):59–78.
  * Kingman, (1992) Kingman, J. (1992).  Poisson Processes, volume 3.  Clarendon Press.
  * Kolda, (2006) Kolda, T. G. (2006).  Multilinear operators for higher-order decompositions, volume 2.  United States. Department of Energy.
  * Kolda and Bader, (2009) Kolda, T. G. and Bader, B. W. (2009).  Tensor decompositions and applications.  SIAM Review, 51(3):455–500.
  * Lázaro-Gredilla et al., (2010) Lázaro-Gredilla, M., Quiñonero-Candela, J., Rasmussen, C. E., and Figueiras-Vidal, A. R. (2010).  Sparse spectrum gaussian process regression.  The Journal of Machine Learning Research, 11:1865–1881.
  * Li et al., (2015) Li, X., Cong, G., Li, X.-L., Pham, T.-A. N., and Krishnaswamy, S. (2015).  Rank-geofm: A ranking based geographical factorization method for point of interest recommendation.  In Proceedings of the 38th international ACM SIGIR conference on research and development in information retrieval, pages 433–442.
  * Liang et al., (2007) Liang, P., Petrov, S., Jordan, M. I., and Klein, D. (2007).  The infinite pcfg using hierarchical dirichlet processes.  In Proceedings of the 2007 joint conference on empirical methods in natural language processing and computational natural language learning (EMNLP-CoNLL), pages 688–697.
  * Lijoi et al., (2010) Lijoi, A., Prünster, I., et al. (2010).  Models beyond the Dirichlet process.  Bayesian nonparametrics, 28(80):342.
  * Liu et al., (2019) Liu, H., Li, Y., Tsang, M., and Liu, Y. (2019).  CoSTCo: A Neural Tensor Completion Model for Sparse Tensors, page 324–334.  Association for Computing Machinery, New York, NY, USA.
  * Lloyd et al., (2012) Lloyd, J. R., Orbanz, P., Ghahramani, Z., and Roy, D. M. (2012).  Random function priors for exchangeable arrays with applications to graphs and relational data.  In Advances in Neural Information Processing Systems 24, pages 1007–1015.
  * McCloskey, (1965) McCloskey, J. W. (1965).  A model for the distribution of individuals by species in an environment.  Michigan State University. Department of Statistics.
  * McFarland et al., (2013) McFarland, D. A., Ramage, D., Chuang, J., Heer, J., Manning, C. D., and Jurafsky, D. (2013).  Differentiating language usage through topic models.  Poetics, 41(6):607–625.
  * Oh et al., (2018) Oh, S., Park, N., Lee, S., and Kang, U. (2018).  Scalable Tucker factorization for sparse tensors-algorithms and discoveries.  In 2018 IEEE 34th International Conference on Data Engineering (ICDE), pages 1120–1131. IEEE.
  * (36) Pan, Z., Wang, Z., and Zhe, S. (2020a).  Scalable nonparametric factorization for high-order interaction events.  In International Conference on Artificial Intelligence and Statistics, pages 4325–4335. PMLR.
  * (37) Pan, Z., Wang, Z., and Zhe, S. (2020b).  Streaming nonlinear bayesian tensor decomposition.  In Conference on Uncertainty in Artificial Intelligence, pages 490–499. PMLR.
  * Paszke et al., (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019).  Pytorch: An imperative style, high-performance deep learning library.  arXiv preprint arXiv:1912.01703.
  * Rahimi et al., (2007) Rahimi, A., Recht, B., et al. (2007).  Random features for large-scale kernel machines.  In NIPS, volume 3, page 5. Citeseer.
  * Rai et al., (2015) Rai, P., Hu, C., Harding, M., and Carin, L. (2015).  Scalable probabilistic tensor factorization for binary and count data.  In IJCAI.
  * Rai et al., (2014) Rai, P., Wang, Y., Guo, S., Chen, G., Dunson, D., and Carin, L. (2014).  Scalable Bayesian low-rank decomposition of incomplete multiway tensors.  In Proceedings of the 31th International Conference on Machine Learning (ICML).
  * Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. I. (2006).  Gaussian Processes for Machine Learning.  MIT Press.
  * Rudin, (1962) Rudin, W. (1962).  Fourier analysis on groups, volume 121967.  Wiley Online Library.
  * Schein et al., (2015) Schein, A., Paisley, J., Blei, D. M., and Wallach, H. (2015).  Bayesian poisson tensor factorization for inferring multilateral relations from sparse dyadic event counts.  In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1045–1054. ACM.
  * (45) Schein, A., Zhou, M., Blei, D., and Wallach, H. (2016a).  Bayesian poisson tucker decomposition for learning the structure of international relations.  In International Conference on Machine Learning, pages 2810–2819. PMLR.
  * (46) Schein, A., Zhou, M., Blei, D. M., and Wallach, H. (2016b).  Bayesian poisson tucker decomposition for learning the structure of international relations.  In Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48, ICML’16, pages 2810–2819. JMLR.org.
  * Shashua and Hazan, (2005) Shashua, A. and Hazan, T. (2005).  Non-negative tensor factorization with applications to statistics and computer vision.  In Proceedings of the 22th International Conference on Machine Learning (ICML), pages 792–799.
  * Sutskever et al., (2009) Sutskever, I., Tenenbaum, J. B., and Salakhutdinov, R. R. (2009).  Modelling relational data using bayesian clustered tensor factorization.  In Advances in neural information processing systems, pages 1821–1828.
  * Teh et al., (2006) Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M. (2006).  Hierarchical dirichlet processes.  Journal of the american statistical association, 101(476):1566–1581.
  * Tillinghast and Zhe, (2021) Tillinghast, C. and Zhe, S. (2021).  Nonparametric decomposition of sparse tensors.  In International Conference on Machine Learning, pages 10301–10311. PMLR.
  * Tucker, (1966) Tucker, L. (1966).  Some mathematical notes on three-mode factor analysis.  Psychometrika, 31:279–311.
  * Wainwright and Jordan, (2008) Wainwright, M. J. and Jordan, M. I. (2008).  Graphical models, exponential families, and variational inference.  Now Publishers Inc.
  * Wang et al., (2011) Wang, C., Paisley, J., and Blei, D. (2011).  Online variational inference for the hierarchical dirichlet process.  In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, pages 752–760. JMLR Workshop and Conference Proceedings.
  * Wang et al., (2020) Wang, Z., Chu, X., and Zhe, S. (2020).  Self-modulating nonparametric event-tensor factorization.  In International Conference on Machine Learning, pages 9857–9867. PMLR.
  * Williamson, (2016) Williamson, S. A. (2016).  Nonparametric network models for link prediction.  The Journal of Machine Learning Research, 17(1):7102–7121.
  * Xu et al., (2012) Xu, Z., Yan, F., and Qi, Y. A. (2012).  Infinite tucker decomposition: Nonparametric bayesian models for multiway data analysis.  In ICML.
  * Yang and Dunson, (2013) Yang, Y. and Dunson, D. (2013).  Bayesian conditional tensor factorizations for high-dimensional classification.  Journal of the Royal Statistical Society B, revision submitted.
  * Zhao et al., (2015) Zhao, Q., Zhang, L., and Cichocki, A. (2015).  Bayesian cp factorization of incomplete tensors with automatic rank determination.  IEEE transactions on pattern analysis and machine intelligence, 37(9):1751–1763.
  * Zhe and Du, (2018) Zhe, S. and Du, Y. (2018).  Stochastic nonparametric event-tensor decomposition.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 6857–6867.
  * (60) Zhe, S., Qi, Y., Park, Y., Xu, Z., Molloy, I., and Chari, S. (2016a).  Dintucker: Scaling up gaussian process models on large multidimensional arrays.  In Thirtieth AAAI conference on artificial intelligence.
  * Zhe et al., (2015) Zhe, S., Xu, Z., Chu, X., Qi, Y., and Park, Y. (2015).  Scalable nonparametric multiway data analysis.  In Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, pages 1125–1134.
  * (62) Zhe, S., Zhang, K., Wang, P., Lee, K.-c., Xu, Z., Qi, Y., and Ghahramani, Z. (2016b).  Distributed flexible nonlinear tensor factorization.  In Advances in Neural Information Processing Systems, pages 928–936.

## Appendix

##  8 Tensor Sampling Details

In general, to simulate a Poisson random measure (PRM) with the rate measure
$\mu(\cdot)$ on $(\Omega,{\mathcal{B}})$ where $\Omega$ is the universal set
and ${\mathcal{B}}$ is a $\sigma$-algebra that includes $\Omega$, we conduct
two steps.

  * •

Sample the number of points $N\sim\text{Poisson}\big{(}\mu(\Omega)\big{)}$.

  * •

Sample $N$ points (locations) i.i.d with the normalized rate measure
$\frac{\mu}{\mu(\Omega)}$ (which is a probability measure).

We use this standard procedure to sample sparse tensors with our model. Given
a particular $\alpha$, we first sample the total mass

$A([0,\alpha]^{R_{1}})=\frac{1}{R_{2}}\sum_{r=1}^{R_{2}}W^{\alpha}_{1,r}([0,\alpha]^{R_{1}})\times\ldots\times
W^{\alpha}_{K,r}([0,\alpha]^{R_{1}}).$

To sample the total mass, we need to sample each
$W^{\alpha}_{k,r}([0,\alpha]^{R_{1}})$ where $1\leq k\leq K$ and $1\leq r\leq
R_{2}$, which according to the definition of $\Gamma$Ps, follows
$\text{Gamma}(1,L_{k}^{\alpha}([0,\alpha]^{R_{1}}))$. Recursively,
$L_{k}^{\alpha}([0,\alpha]^{R_{1}})$ is sampled from
$\text{Gamma}(1,\lambda_{\alpha}([0,\alpha]^{R_{1}}))$. We then sample the
number of present entries (i.e., points) $N$ from a Poisson distribution with
the total mass as the mean. Next, we use the normalized H$\Gamma$Ps, namely,
HDPs defined in Sec 3.2 of the main paper to sample each entry. To obtain the
probability measure over all possible entries, namely, the normalized rate
measure in (8) of the main paper, we first sample the weights $\beta^{k}_{j}$
in the first level of DPs (see (6) of the main paper). We use the stick-
breaking construction. Usually, when $j>2000$, the weights are below the
machine precision and automatically truncated to zero, and we do not need to
store an infinite number of weights. Following (Teh et al.,, 2006), the stick-
breaking construction of the weights in a second-level DP $H^{k}_{r}$ is

$\displaystyle\nu^{k}_{rj}\sim\text{Beta}\left(\gamma^{k}_{r}\beta^{k}_{j},\gamma^{k}_{r}\left(1-\sum_{l=1}^{j}\beta^{k}_{j}\right)\right),\;\omega^{k}_{rj}=\nu^{k}_{rj}\prod\nolimits_{t=1}^{j-1}(1-\nu^{k}_{rt})\;(1\leq
j\leq\infty).$

Accordingly, $\nu^{k}_{rj}$ will be truncated to zero at which $\beta^{k}_{j}$
is zero, and so is the weight $\omega^{k}_{rj}$. Hence, we obtain a finite set
of nonzero weights, with which we calculate the probability measure to sample
the present entries.

|  |  |
---|---|---|---

(a) Alog

(b) Alog

(c) SG

(d) SG

Figure 3: The performance of our approach with different combinations of
$R_{1}$ and $R_{2}$ (the number of location and sociability factors
respectively). The total number of factors $R$ is fixed to $11$
($R=R_{1}+R_{2}$).

|  |
---|---|---
|  |
|  |

(a) Alog \- Ours

(b) MovieLens \- Ours

(c) SG\- Ours

(d) Alog \- CP-Bayes

(e) MovieLens \- CP-Bayes

(f) SG\- CP-Bayes

(g) Alog \- GPTF

(h) MovieLens -GPTF

(i) SG\- GPTF

Figure 4: Structures of the estimated latent factors for the second mode of
Alog and MovieLens dataset and third mode of SG dataset, by our method, CP-
Bayes and GPTF.

##  9 Combination of $R_{1}$ and $R_{2}$

We investigated the performance of our method with different combinations of
$R_{1}$ and $R_{2}$ (the number of location and sociability factors) on Alog
and SG datasets. To this end, we fixed the total number of factors to $11$ and
varied $R_{2}$, i.e., the number of sociability factors. We used the same data
splits as in Sec. 6.2 of the main paper. The average MSE and MAE for different
combinations of $R_{1}$ and $R_{2}$ are reported in Fig. 3. We can see that
the best setting on Alog is $R_{2}=6,R_{1}=5$, and on SC is $R_{2}=2,R_{1}=9$
(MSE) and $R_{2}=1,R_{1}=10$ (MAE). While the results of all the settings are
better than all the competing dense tensor models, different combinations of
$R_{1}$ and $R_{2}$ result in quite different prediction accuracies. The best
trade-off is application specific. The results also demonstrate the advantage
of our method over NEST models (Tillinghast and Zhe,, 2021), which restrict
the types of factors to be estimated and do not have the flexibility to choose
different $R_{1}$ and $R_{2}$.

##  10 Running Time

The speed of our method is faster than NEST, comparable to GPTF, and slower
than the other competing methods. For example, for $R=5$, on MovieLens and SG
datasets, the average running time (in seconds) of each method is {Ours:
825.1, GPTF: 706.5, NEST-1: 1218.9, NEST-2: 1193.4, CP-Bayes: 323.2, P-Tucker:
7.9, CP-WOPT: 6.1, CP-ALS: 0.8} and {Ours: 1213.5, GPTF: 1096.1, NEST-1:
2112.2, NEST-2: 1910.5, CP-Bayes: 566.4, P-Tucker: 10.8, CP-WOPT: 10.6, CP-
ALS: 0.7}. This is reasonable, because the GP based methods, although using
sparse approximations to avoid computing the full covariance (kernel) matrix,
still need to estimate the frequencies and perform extra nonlinear feature
transformations (i.e., random Fourier features). Hence, they are much more
complex. In addition, our method introduces nonparametric priors over sparse
tensor entries, and the estimation involves discrete, hierarchical probability
measures. Therefore, it requires more computation.

##  11 Proof of Lemma 3.1 and 3.2

###  11.1 Proof of Lemma 3.1

To ease the notation, we consider the case where $R_{1}=R_{2}=1$. It is
straightforward to extend the proof to arbitrary $R_{1}$ and $R_{2}$. Given a
particular $\alpha<\infty$, denote the total number of points in $T$ (see (4)
in the main paper) by $D^{\alpha}$. We observe that by the properties of the
Poisson random measure,
$D^{\alpha}\sim\text{Poisson}(W_{1}^{\alpha}([0,\alpha])\times\dots\times
W_{K}^{\alpha}([0,\alpha]))$. Since these points might overlap, the number of
the distinct points (i.e., entries) $N^{\alpha}\leq D^{\alpha}$. According to
the definition of Gamma process, $W_{k}^{\alpha}([0,\alpha])|L_{k}^{\alpha}$
is a Gamma distributed random variable with mean given by
$L_{k}^{\alpha}([0,\alpha])$, and $L_{k}^{\alpha}([0,\alpha])$ is another
Gamma random variable with mean $\alpha$. Therefore,
$L_{k}^{\alpha}([0,\alpha])<\infty$ almost surely,
$W_{k}^{\alpha}([0,\alpha])|L_{k}^{\alpha}<\infty$ almost surely, and then
$W_{k}^{\alpha}([0,\alpha])<\infty$ almost surely (because
$p(W_{k}^{\alpha}([0,\alpha])<\infty)=\int
p\left(W_{k}^{\alpha}([0,\alpha])<\infty|L_{k}^{\alpha}\right){\rm
d}p(L_{k}^{\alpha})=1$). Hence, their product
$b^{\alpha}=W_{1}^{\alpha}([0,\alpha])\times\dots\times
W_{K}^{\alpha}([0,\alpha])$ is finite almost surely. Since
$D^{\alpha}\sim\text{Poisson}(b^{\alpha})$, we have $D^{\alpha}<\infty$ almost
surely and $N^{\alpha}\leq D^{\alpha}$ implies $N^{\alpha}<\infty$ almost
surely.

Now we consider $W_{k}^{\alpha}((j-1,j])|L_{k}^{\alpha}$ where
$j\in\mathbb{N}^{+}$. This is a Gamma distributed random variable with mean
given by $L_{k}^{\alpha}((j-1,j])$. However $L_{k}^{\alpha}((j-1,j])$ is a
mean-one Gamma distributed random variable. Thus if we integrate out
$L^{\alpha}_{k}$ for $1\leq k\leq K$, $\\{W^{\alpha}_{k}((j-1,j])\\}_{k,j}$
are independent identically distributed random variables because
$L^{\alpha}_{k}$ as a Gamma process itself is independent on disjoint sets.
Note that given $j$, $\\{W^{\alpha}_{k}((j-1,j])\\}$ across $k$ have already
been independent. Furthermore, we have
$\mathbb{E}[W_{k}^{\alpha}((j-1,j])]=1$, and
$\mathbb{E}[W_{1}^{\alpha}((j-1,j])\times\ldots\times
W_{K}^{\alpha}((j-1,j])]=1$.

We consider
$S_{\lfloor\alpha\rfloor}=\sum_{j=1}^{\lfloor\alpha\rfloor}\mathds{1}(X_{j}>0)$
where $\mathds{1}(\cdot)$ is the indicator function, and

$X_{j}\sim\text{Poisson}\left(W_{1}^{\alpha}((j-1,j])\times\ldots\times
W_{K}^{\alpha}((j-1,j])\right).$

Obviously, $S_{\lfloor\alpha\rfloor}\leq N^{\alpha}$. For convenience, let us
define $\zeta_{j}=W_{1}^{\alpha}((j-1,j])\times\ldots\times
W_{K}^{\alpha}((j-1,j])$. We have $\mathbb{E}[X_{j}|\zeta_{j}]=\zeta_{j}$ and
$\mathbb{E}[X_{j}]=\mathbb{E}[\zeta_{j}]=1$ . Since all $\\{X_{j}\\}_{j}$ are
i.i.d, the indicators $\\{\mathds{1}(X_{j}>0)\\}_{j}$ are i.i.d as well.
According to the strong law of large numbers,

$\lim_{\lfloor\alpha\rfloor\rightarrow\infty}\frac{S_{\lfloor\alpha\rfloor}}{\lfloor\alpha\rfloor}=\mathbb{E}[\mathds{1}(X_{j}>0)]=1-e^{-1}>0\,\,a.s.$

Therefore, when $\alpha\rightarrow\infty$,
$\lfloor\alpha\rfloor\rightarrow\infty$ and
$S_{\lfloor\alpha\rfloor}\rightarrow\infty$ a.s. Since $N^{\alpha}\geq
S_{\lfloor\alpha\rfloor}$, we have
$\lim_{\alpha\rightarrow\infty}N^{{\alpha}}=\infty$ a.s.

###  11.2 Proof of Lemma 3.2

Our tensor-variate stochastic process when $R_{2}=1$ is summarized in the
following.

$\displaystyle L^{\alpha}_{k}$
$\displaystyle\sim\Gamma\text{P}(\lambda_{\alpha}),$ $\displaystyle
W^{\alpha}_{k}\mid L^{\alpha}_{k}$
$\displaystyle\sim\Gamma\text{P}(L^{\alpha}_{k}),$ $\displaystyle
T\mid\\{W^{\alpha}_{k}\\}_{k=1}^{\infty}$
$\displaystyle\sim\text{PRM}(W^{\alpha}_{1}\times\dots\times W^{\alpha}_{K}).$

We will follow the similar strategy in (Tillinghast and Zhe,, 2021) to prove
the lemma. That is, given $\alpha$, we define $M^{\alpha}_{k}$ the number of
active nodes in mode $k$ and $N^{\alpha}$ the number of sampled entries. Then
we have
$M^{\alpha}_{k}=\\#\\{\theta^{k}_{i}\in[0,\alpha]|T(A^{\alpha}_{k,\theta_{i}^{k}})>0\\}$,
where
$A^{\alpha}_{k,\theta_{i}^{k}}=[0,\alpha]\times\dots\times\\{\theta^{k}_{i}\\}\times\dots\times[0,\alpha]$.
The goal of the first step is to show
$\lim_{\alpha\rightarrow\infty}\frac{\alpha}{M_{k}^{\alpha}}=0\;\;a.s.$ for
all $k\in\\{1,\dots,K\\}$, and the second step is
$\limsup_{\alpha\rightarrow\infty}N^{\alpha}/\alpha^{K}<\infty\;\;a.s$. Based
on these, it is trivial to show

$\lim_{\alpha\rightarrow\infty}\dfrac{N^{\alpha}}{\prod_{k=1}^{K}M_{k}^{\alpha}}=0\;\;a.s.$

Step 1. Following the same steps as in (Tillinghast and Zhe,, 2021), we can
show that

$\displaystyle\mathbb{E}[M_{k}^{\alpha}|\\{W_{j}^{\infty}\\}_{j\neq
k},L_{k}^{\alpha}]$
$\displaystyle=\mathbb{E}\left[\sum_{\theta_{i}^{k}\in[0,\alpha]}1-\exp\left(-W_{k}^{\infty}(\\{\theta^{k}_{i}\\}\times\prod_{j\neq
k}W^{\infty}_{j}([0,\alpha])\right)\bigg{|}\\{W^{\infty}_{j}\\}_{j\neq
k},L_{k}^{\alpha}\right].$

Note that different from (Tillinghast and Zhe,, 2021), the expectation here is
conditioned on $L_{k}^{\alpha}$, an additional level of $\Gamma$Ps. We then
apply Campbell’s Theorem (Kingman,, 1992) and obtain

$\displaystyle\mathbb{E}[M_{k}^{\alpha}|\\{W_{j}^{\infty}([0,\alpha])\\}_{j\neq
k},L_{k}^{\alpha}]$
$\displaystyle=L_{k}^{\alpha}([0,\alpha])\int_{0}^{\infty}\left(1-\exp\left(-w\times\prod_{i\neq
k}W^{\infty}_{i}([0,\alpha])\right)\right)w^{-1}e^{-w}dw.$

Since $L_{k}^{\alpha}([0,\alpha])$ is a Gamma random variable with mean
$\alpha$, we further take the expectation with respect to the measure induced
by $L_{k}^{\alpha}$. Then, we obtain

$\mathbb{E}[M_{k}^{\alpha}|\\{W_{j}^{\infty}([0,\alpha])\\}_{j\neq
k}]=\alpha\cdot\int_{0}^{\infty}\left(1-\exp\left(-w\times\prod_{i\neq
k}W^{\infty}_{i}([0,\alpha])\right)\right)w^{-1}e^{-w}dw.$

Now, we arrive at the same intermediate result as in (Tillinghast and Zhe,,
2021). It follows that

$\lim_{\alpha\rightarrow\infty}\frac{\alpha}{M_{k}^{\alpha}}=0\;\;a.s.$

The result points out that in each mode, the number of active nodes
$M_{k}^{\alpha}$ grows faster than $\alpha$.

Step 2. We consider the distribution of $W_{k}([j,j+1])\mid L_{k}^{\alpha}$.
This will be Gamma distributed with mean $L_{k}^{\alpha}([j,j+1])$. But
$L_{k}^{\alpha}(j,j+1)$ is Gamma distributed with mean 1 for all $j<\alpha-1$.
Thus by marginalizing out $L_{k}^{\alpha}(j,j+1)$ we see that $W_{k}([j,j+1])$
is identically distributed for all $j<\alpha-1$. By the independence of the
completely random measure (Kingman,, 1967) on disjoint sets, it follows
immediately by the strong law of large numbers

$\lim_{j\rightarrow\infty}\dfrac{W_{k}^{\infty}([0,j])}{j}=\dfrac{\sum_{i=1}^{j}W_{k}^{\infty}((i-1,i])}{j}=E[W_{k}^{\infty}([0,1])]=1\,a.s.$

as $W_{k}^{\infty}((i-1,i])$ are i.i.d random variables. This implies

$\lim_{j\rightarrow\infty}\dfrac{\prod_{k=1}^{K}W_{k}^{\infty}([0,j])}{j^{K}}=1\,a.s.$
(15)

Let us denote $D^{\alpha}$ by the actual number of points sampled in
$T([0,\alpha]^{K})$. Note $N^{\alpha}\leq D^{\alpha}$. Applying Lemma 17 in
(Caron and Fox,, 2014) implies

$Pr\left(\lim_{j\rightarrow\infty}\dfrac{D^{j}}{\prod_{k=1}^{K}W_{k}^{\infty}([0,j])}=1\bigg{|}\\{W^{\infty}_{i}\\}_{i=1}^{K}\right)=1.$

Taking the expectation of both sides of the above expression and combining
with equation (15) implies

$\lim_{j\rightarrow\infty}\dfrac{D^{j}}{j^{K}}=1\,a.s.$

Now, we can use the same bounding technique as in (Tillinghast and Zhe,, 2021)
to show that

$\lim_{\alpha\rightarrow\infty}\dfrac{N^{\alpha}}{\alpha^{K}}\leq\lim_{\alpha\rightarrow\infty}\dfrac{D^{\alpha}}{\alpha^{K}}=1<\infty.$

The result in the second step points out the number of sampled entries
$N^{\alpha}$ grows slower than or as fast as $\alpha^{K}$.

Finally, we extend the case when $R_{2}>1$ by simply applying Poisson
superposition theorem.