Stochastic Gradients for Large-Scale Tensor Decomposition

TL;DR

This paper introduces a stochastic gradient method for large-scale tensor decomposition, enabling efficient analysis of big data with various loss functions, and demonstrates its scalability and effectiveness.

Contribution

It proposes a novel stochastic gradient approach for generalized tensor decomposition that handles large-scale and sparse tensors efficiently.

Findings

Method scales to large tensors

Effective for sparse and dense data

Outperforms traditional approaches in experiments

Abstract

Tensor decomposition is a well-known tool for multiway data analysis. This work proposes using stochastic gradients for efficient generalized canonical polyadic (GCP) tensor decomposition of large-scale tensors. GCP tensor decomposition is a recently proposed version of tensor decomposition that allows for a variety of loss functions such as Bernoulli loss for binary data or Huber loss for robust estimation. The stochastic gradient is formed from randomly sampled elements of the tensor and is efficient because it can be computed using the sparse matricized-tensor-times-Khatri-Rao product (MTTKRP) tensor kernel. For dense tensors, we simply use uniform sampling. For sparse tensors, we propose two types of stratified sampling that give precedence to sampling nonzeros. Numerical results demonstrate the advantages of the proposed approach and its scalability to large-scale problems.

Tables

Table 1. Table 1: Empirical bias and variance of stochastic gradients at an initial guess and a final model for GCP with Bernoulli loss using the artificial sparse tensor and sampling methods of Fig.   5 . For each method and model, we use N = 1000 𝑁 1000 N=1000 stochastic gradient realizations.

sampling	Initial guess: ${\Vert 𝐠\Vert }_2=$ 2.00e+07		Final model: ${\Vert 𝐠\Vert }_2=$ 3.10e+06
method	emp. bias	emp. var.	emp. bias	emp. var.
uniform	1.08e+06	1.26e+15	3.90e+07	1.52e+18
stratified	7.48e+05	5.64e+14	3.17e+06	9.84e+15
semi-stratified	7.68e+05	5.70e+14	3.14e+06	9.93e+15

Equations

$$\mathcal{X}\approx \mathcal{M}\text{where}\mathcal{M}=\sum _{j=1}^r{\mathbf{a}}_1(:,j)\circ {\mathbf{a}}_2(:,j)\circ \cdots \circ {\mathbf{a}}_d(:,j).$$

$$\text{minimize }F(\mathcal{X},\mathcal{M})\equiv \sum _{i_1=1}^{n_1}\cdots \sum _{i_d=1}^{n_d}f(x_{i_1\dots i_d},m_{i_1\dots i_d})\text{subject to}\mathrm{rank}(\mathcal{M})\le r.$$

$$\mathbb{E}[\tilde{x}]\equiv \sum _{i=1}^m{p}_i{v}_i.$$

$$n=\sqrt[d]{\prod _{k=1}^d{n}_k}\text{and}\bar{n}=\frac{1}{d}\sum _{k=1}^d{n}_k.$$

$${\mathbf{Y}}_{(k)}\underset{{\mathbf{Z}}_k}{\underbrace{({\mathbf{A}}_d\odot \cdots \odot {\mathbf{A}}_{k+1}\odot {\mathbf{A}}_{k-1}\odot \cdots \odot {\mathbf{A}}_1)}}.$$

$$y_i=\frac{\partial f}{\partial m}(x_i,m_i).$$

$${\mathbf{G}}_k={\mathbf{Y}}_{(k)}{\mathbf{Z}}_k,$$

$$y_i=\{\begin{array}{ll}\frac{\partial f}{\partial m}(x_i,m_i)& \text{if xi known},\\ 0& \text{if xi missing}.\end{array}$$

$${\tilde{\mathbf{G}}}_k={\tilde{\mathbf{Y}}}_{(k)}{\mathbf{Z}}_k\text{where}\mathbb{E}[\tilde{\mathcal{Y}}]=\mathcal{Y}\text{ and }\text{nnz}(\tilde{\mathcal{Y}})\le s\ll n^d.$$

$$\mathbb{E}[{\tilde{s}}_i]=\frac{s}{n^d}\text{for all}i\in \Omega .$$

$${\tilde{y}}_i={\tilde{s}}_i\frac{n^d}{s}{y}_i\text{for all}i\in \Omega .$$

$$i_k=\text{randi}(n_k)\text{for}k\in \{1,2,\dots ,d\}.$$

$$\mathbb{E}[{\tilde{s}}_i]=\frac{s_{\ell }}{|{\Omega }_{\ell }|}\text{where}i\in {\Omega }_{\ell }\text{for all}i\in \Omega .$$

$${\tilde{y}}_i={\tilde{s}}_i\frac{|{\Omega }_{\ell }|}{s_{\ell }}{y}_i\text{where}i\in {\Omega }_{\ell }\text{for all}i\in \Omega .$$

$$\frac{q}{1-\eta /n^d}=\frac{n^d}{\zeta }q\approx q.$$

$${\tilde{y}}_i={\tilde{p}}_i\frac{\eta }{p}(y_i-c_i)+{\tilde{q}}_i\frac{n^d}{q}{c}_i\text{where}{c}_i\equiv \frac{\partial f}{\partial m}(0,m_i).$$

$$\mathbb{E}[{\tilde{y}}_i]=\mathbb{E}[{\tilde{q}}_i]\frac{n^d}{q}{c}_i=\frac{q}{n^d}\frac{n^d}{q}{y}_i=y_i.$$

$$\mathbb{E}[{\tilde{y}}_i]=\mathbb{E}[{\tilde{p}}_i]\frac{\eta }{p}(y_i-c_i)+\mathbb{E}[{\tilde{q}}_i]\frac{n^d}{q}{c}_i=\frac{p}{\eta }\frac{\eta }{p}(y_i-c_i)+\frac{q}{n^d}\frac{n^d}{q}{c}_i=(y_i-c_i)+c_i=y_i.$$

$$\text{minimize }F(\mathcal{X},\mathcal{M})\equiv \sum _i{w}_{i}f(x_i,m_i)\text{subject to}\mathrm{rank}(\mathcal{M})\le r.$$

$$y_i=w_i\frac{\partial f}{\partial m}(x_i,m_i).$$

$$\hat{F}\equiv \sum _{i\in \Omega }{\tilde{s}}_i\frac{n^d}{s}f(x_i,m_i).$$

$$\hat{F}\equiv \sum _{x_i\ne 0}{\tilde{p}}_i\frac{\eta }{p}f(x_i,m_i)+\sum _{x_i=0}{\tilde{q}}_i\frac{\zeta }{q}f(x_i,m_i).$$

$$\text{empirical bias }=\Vert \hat{\mathbf{g}}-\mathbf{g}{\Vert }_2\text{where}\hat{\mathbf{g}}=\frac{1}{N}\sum _{\xi =1}^N{\tilde{\mathbf{g}}}_{\xi },\text{and}$$

$$\text{empirical variance }=\text{trace}\left(\frac{1}{N}\sum _{\xi =1}^N({\tilde{\mathbf{g}}}_{\xi }-\hat{\mathbf{g}})({\tilde{\mathbf{g}}}_{\xi }-\hat{\mathbf{g}}{)}^{⊺}\right)=\frac{1}{N}\sum _{\xi =1}^N\Vert {\tilde{\mathbf{g}}}_{\xi }-\hat{\mathbf{g}}{\Vert }_2^2.$$

$$v_i=\{\begin{array}{ll}{x}_i/m_i& \text{if }{x}_i\ne 0,\\ 0& \text{otherwise}.\end{array}$$

$$\rho \frac{n^d}{n^d-\text{nnz}(\mathcal{X})}{s}_{\text{zero}}=\rho \frac{1}{1-\text{nnz}(\mathcal{X})/n^d}{s}_{\text{zero}}=\rho \frac{s_{\text{zero}}}{p_{\text{zero}}}$$

$$s_{\text{reject}}=\text{icdf(‘Negative Binomial’, 0.999999, }{s}_{\text{zero}}\text{, }{p}_{\text{zero}}\text{)}.$$

$$\rho \ge (s_{\text{zero}}+s_{\text{reject}})\frac{p_{\text{zero}}}{s_{\text{zero}}}=\frac{s_{\text{zero}}+s_{\text{reject}}}{s_{\text{zero}}}{p}_{\text{zero}}.$$

$${\mathbf{A}}_k(i_k,j)\sim \text{uniform}(0,1)\text{for all}{i}_k=1,\dots ,n_k,j=1,\dots ,r,\text{ and }k=1,\dots ,d.$$

$$x_i\sim \text{gamma}(k,m_i)\text{with}k=1.$$

Full text

1 Introduction

Tensor decomposition is the higher-order analogue of matrix decomposition and is becoming an everyday tool for data analysis. It can be used for unsupervised learning, dimension reduction, tensor completion, feature extraction in supervised machine learning, data visualization, and more. For a given $d$-way data tensor $\bm{\mathscr{{X}}}$ of size $n_{1}\times n_{2}\times\cdots\times n_{d}$, the goal is to find a low-rank approximation $\bm{\mathscr{{M}}}$, i.e.,

[TABLE]

The low-rank structure of $\bm{\mathscr{{M}}}$ reveals patterns within the data, as defined by the factor matrices. The $k$th factor matrix of size $n_{k}\times r$ is denoted ${\bm{\mathbf{{A}}}}_{k}$. The $j$th factor (column) in mode $k$ is denoted ${\bm{\mathbf{{a}}}}_{k}(:,j)$. Each component of $\bm{\mathscr{{M}}}$ is a $d$-way outer product (denoted by $\circ$) of $d$ factors, forming a rank-one tensor. We say $\operatorname{rank}(\bm{\mathscr{{M}}})\leq r$ because $\bm{\mathscr{{M}}}$ can be written as the sum of $r$ rank-one tensors. We define $n\equiv(\prod_{k=1}^{d}n_{k})^{1/d}$ and assume $r\ll n^{d}$. The storage for $\bm{\mathscr{{X}}}$ is $n^{d}$ if it is dense and $O(\text{nnz}(\bm{\mathscr{{X}}}))$ if it is sparse. The storage for $\bm{\mathscr{{M}}}$ is $O(r\sum_{k=1}^{d}n_{k})$, which is usually much smaller than that for $\bm{\mathscr{{X}}}$.

The standard canonical polyadic or CANDECOMP/PARAFAC (CP) tensor decomposition seeks the best low-rank approximation with respect to sum of squared errors [9, 19]. The generalized CP (GCP) tensor decomposition is a novel approach that allows for an arbitrary elementwise loss function that is summed across all tensor entries [22]; i.e., the user provides a scalar loss function $f:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ resulting in the optimization problem

[TABLE]

GCP is useful in situations where non-standard choices of the scalar loss function $f$ may be appropriate. For instance, if the entries of $\bm{\mathscr{{X}}}$ are binary values in $\set{0,1}$, we may use the Bernoulli loss, i.e., $f(x,m)=\log(m+1)-x\log(m)$ where $m$ represents the odds of observing a one. For positive data that has a Gamma distribution, we may use $f(x,m)=x/m+\log(m)$. Standard CP tensor decomposition corresponds to $f(x,m)=(x{-}m)^{2}$, and Poisson tensor decomposition [54, 12] to $f(x,m)=m-x\log m$. See [22] for full details.

In this paper, we consider the problem of fitting GCP for large-scale tensors. The GCP gradient involves a sequence of $d$ matricized-tensor times Khatri-Rao products (MTTKRPs) with a dense tensor of size $n^{d}$ that costs $O(rn^{d})$ operations, even when $\bm{\mathscr{{X}}}$ is sparse.111Both standard and Poisson CP have special structure that allows the gradient to be formed implicitly, as described in Appendix A so that sparsity of $\bm{\mathscr{{X}}}$ can be exploited. Thus, the computational and storage costs of computing the gradient may be prohibitive. Our solution is a flexible framework for stochastic gradient computation for GCP, which can be used with stochastic gradient descent (SGD) or popular variants such as Adam [25]. Our stochastic framework replaces the dense tensor with a (stochastic) sparse tensor that equals the dense tensor in expectation. If the stochastic gradient samples $s\ll n^{d}$ tensor entries, the resulting sparse MTTKRP reduces the gradient cost to $O(sr)$ flops and $O(s)$ intermediate storage, albeit with a potentially significant sacrifice in accuracy due to stochasticity.

A distinguishing feature of our framework is stratified sampling. If $\bm{\mathscr{{X}}}$ is sparse, uniform sampling of the indices rarely samples nonzeros even though they are extremely important to the minimization. Stratified sampling randomly selects nonzeros disproportionately, which can reduce the variance of the stochastic gradient, accelerating convergence. Since stratified sampling can itself be expensive, we also introduce a semi-stratified sampling approach to further improve efficiency. Note that simply ignoring the zero entries, as is done in recommender systems where zeros indicate missing data, yields a fundamentally different problem. Section 2 surveys related work, Section 3 provides notation and background, and Section 4 discusses the framework.

We present computational experiments on both artificial and real-world problems in Section 5, demonstrating the reliability and efficiency of using stochastic gradients. Experimentally, we can set the number of gradient samples to be the sum of the tensor dimensions, i.e., $\sum_{k=1}^{d}n_{k}$. For dense problems, stochastic optimization can be an order of magnitude faster than non-stochastic optimization. For sparse problems, stochastic optimization enables computing GCP for sparse tensors that would otherwise be intractable. Additionally, stratified and semi-stratified sampling of zeros and nonzeros have a clear advantage over uniform sampling in this setting. We also compare to CP-APR [12, 18] on real-world Chicago crime data. In our experiments, the stochastic methods are faster and find equivalent solutions.

2 Related Work

This work is a follow-on to the introduction of GCP by Hong, Kolda, and Duersch [22]. They only considered small-scale problems that could be treated as dense. The goal of this work is to tackle larger scale problems via stochastic gradient methods.

A major motivation for the present work was that of Acar et al. [2] which focused on standard CP factorization for tensors with irretrievably missing data. They considered cases where the vast majority of the data is missing, i.e., scarce tensors in the terminology of [22], and observed that a scarce data tensor leads to a sparse MTTKRP in computing the gradient. That observation proves to be important here because our sampled data tensors are scarce.

In the specific case of third-order ($d=3$) and one sample per iteration ($s=1$), Beutel et al. [7] have proposed SGD for standard CP tensor decomposition. Since a single tensor element only updates one row in each of the three factor matrices, they are able to parallelize updates that do not operate on the same rows. They do not consider sampling strategies or the MTTKRP structure in the gradient.

Vervliet and De Lathauwer [51] pursue an alternative form of sampling for standard CP: rather than sampling elements, they sample indices from each mode of the tensor. The resulting random subtensor is used in either one outer iteration of alternating least squares (ALS) or one iteration of Gauss-Newton, with careful attention to step sizes as is needed in stochastic optimization. The updates only modify the rows of the factor matrices corresponding to the subsampled indices. This block sampling strategy is not readily amenable to sparse data tensors, which are a focus of our work.

For alternating least squares (CP-ALS), randomized matrix sketching can be used to solve the linear least squares problems in the dense [6] and sparse [11] cases; this approach does not directly apply to GCP because its subproblems are not least squares problems. Also in the vein of sketching, random projections can be used to build a sketch (or multiple sketches) of the complete tensor [38, 41, 59], including the special case of symmetric factorizations [53, 45].

In the streaming context, one dimension is typically treated as time with a new slice arriving at each time step and can be updated via SGD; see, e.g., [34]. Maehara, Hayashi, and Kawarabayashi [32] factorize tensors that arise as sums of streaming tensor samples. Similar to our approach, they apply SGD on the samples to update the factorization of the sum. In contrast to our focus on effective sampling, they primarily consider applications where a stream of samples is already given. There are also other works on updates for streaming tensors that are not based on SGD [37, 49, 31, 17].

In recommender systems and tensor completion, it is typically assumed that most data is missing, resulting in scarce tensors (sometimes incorrectly conflated with sparse tensors). Smith, Park, and Karypis [43] consider SGD for scarce tensors and only sample nonzeros as they correspond to observed entries. For symmetric tensors with missing data, Ge et al. [14] prove that SGD finds a local minimizer rather than a saddle point.

We close by highlighting that SGD has already had broad success in matrix factorization, primarily for recommender systems where zeros are treated as missing data (scarce matrices); see, e.g., [28] and references therein. However, we omit a full discussion of stochastic gradients in the matrix setting because the field is vast. To the best of our knowledge, nothing like the stratified and semi-stratified sampling proposed in this work has been proposed in the matrix regime. Gemulla et al. [15] propose a stratified version of SGD as a technique to partition the data into independent segments that can be processed in parallel, but their use of stratification is different than our proposal.

3 Notation and background

We provide some context for the remainder of the paper.

3.1 Probability notation and sampling background

In this paper, we use a tilde to indicate random variables and instances thereof, e.g., $\tilde{x}$. The expected value is denoted $\mathbb{E}[\tilde{x}]$. If $\tilde{x}\in\set{v_{1},v_{2},\dots v_{m}}$ is a discrete random variable that samples $v_{i}$ with probability $p_{i}$, then

[TABLE]

The random sample is uniform if each element has equal probability; i.e., $p_{i}=1/m$ for $i=1,2,\dots,m$. Sampling with replacement means that the same $v_{i}$ can be sampled more than once, i.e., that subsequent samples are independent identically distributed (i.i.d.) draws. If we let $\tilde{s}_{i}$ denote the number of times that $v_{i}$ is sampled from a uniform distribution over $s$ draws with replacement, then $\mathbb{E}[\tilde{s}_{i}]=s/m$ since there are $s$ independent draws and each draw selects sample $v_{i}$ with probability $1/m$.

3.2 Tensor notation

For a $d$-way data tensor $\bm{\mathscr{{X}}}$ of size $n_{1}\times n_{2}\times\cdots\times n_{d}$, we define

[TABLE]

These quantities provide convenient measures of how large the tensor is. We treat $d$ as constant for big-O notation. Throughout, we let $i=(i_{1},i_{2},\dots,i_{d})$ where $i_{k}\in\set{1,2,\dots,n_{k}}$ and refer to this as a multi-index or simply an index. Every multi-index has a corresponding linear index between 1 and $n^{d}$ [26]. We let $\Omega$ denote the set of all tensor indices, so necessarily $|\Omega|=n^{d}$. Excepting the discussion of missing data in Section 4.4, we assume all tensor entries are known. We let $\text{nnz}(\bm{\mathscr{{X}}})$ denote the number of nonzeros in $\bm{\mathscr{{X}}}$ and say that $\bm{\mathscr{{X}}}$ is sparse if $\text{nnz}(\bm{\mathscr{{X}}})\ll n^{d}$.

If $\bm{\mathscr{{X}}}$ is dense, then the storage is $n^{d}$. If $\bm{\mathscr{{X}}}$ is sparse, then only the nonzeros and their indices are stored, so the storage is $(d+1)\,\text{nnz}(\bm{\mathscr{{X}}})$ using a coordinate format [3]. The total storage for a rank-$r$ approximation $\bm{\mathscr{{M}}}$ as in Eq. 1 is the storage for the factor matrices, i.e., $\bar{n}dr$, and is usually significantly less than the storage for $\bm{\mathscr{{X}}}$

3.3 MTTKRP background

Given a tensor $\bm{\mathscr{{Y}}}$ of size $n_{1}\times n_{2}\times\cdots\times n_{d}$ and factor matrices ${\bm{\mathbf{{A}}}}_{k}$ (from a low-rank $\bm{\mathscr{{M}}}$) of size $n_{k}\times r$ for $k=1,\dots,d$, the MTTKRP in mode $k$ is

[TABLE]

The matrix ${\bm{\mathbf{{Z}}}}_{k}$ is of size $n^{d}/n_{k}\times r$ and is the Khatri-Rao product (denoted by $\odot$) of all the factor matrices except ${\bm{\mathbf{{A}}}}_{k}$, and the matrix ${\bm{\mathbf{{Y}}}}_{(k)}$ is of size $n_{k}\times n^{d}/n_{k}$ and denotes the mode-$k$ unfolding of $\bm{\mathscr{{Y}}}$.

Much work has gone into efficient computation of MTTKRP. Bader and Kolda [3] consider both dense and sparse $\bm{\mathscr{{Y}}}$ tensors, showing that the cost is $O(rn^{d})$ for dense $\bm{\mathscr{{Y}}}$ and $O(r\,\text{nnz}(\bm{\mathscr{{Y}}}))$ for sparse $\bm{\mathscr{{Y}}}$. Phan, Tichavsky, and Cichocki [39] propose methods to reuse partial computations when computing the MTTKRP for all $d$ modes in sequence. Much recent work has focused on more efficient representations of sparse tensors and parallel MTTKRP computations [44, 24, 29, 40]. There is also continued work on improving the efficiency of dense MTTKRP calculations [20, 5].

4 Stochastic GCP Gradient

GCP tensor decomposition generalizes CP tensor decomposition, minimizing the GCP loss function in Eq. 2. From [22, Theorem 3], the GCP gradient is calculated via a sequence of MTTKRP computations. Key to this calculation is the elementwise partial gradient tensor $\bm{\mathscr{{Y}}}$ that is the same size as $\bm{\mathscr{{X}}}$ and is defined as

[TABLE]

That is, $\bm{\mathscr{{Y}}}$ is the tensor of first derivatives of $f$ evaluated elementwise w.r.t. $\bm{\mathscr{{X}}}$ and $\bm{\mathscr{{M}}}$. Whether the data tensor $\bm{\mathscr{{X}}}$ is dense or sparse, the elementwise partial gradient tensor $\bm{\mathscr{{Y}}}$ is dense almost everywhere. The partial derivative, denoted ${\bm{\mathbf{{G}}}}_{k}$, of the objective $F$ in Eq. 2 w.r.t. ${\bm{\mathbf{{A}}}}_{k}$ is

[TABLE]

where ${\bm{\mathbf{{Z}}}}_{k}$ is the Khatri-Rao product defined in Eq. 3. Because $\bm{\mathscr{{Y}}}$ is dense (even when $\bm{\mathscr{{X}}}$ is sparse), $\bm{\mathscr{{Y}}}$ requires $n^{d}$ storage and the cost of computing MTTKRP for gradients ${\bm{\mathbf{{G}}}}_{k}$ is $O(rn^{d})$. Such costs may be prohibitive. For instance, if $n=1000$ and $d=4$, then $\bm{\mathscr{{Y}}}$ would require 8 TB of storage. There are a couple of special cases (see Appendix A) where dense computation can be avoided because $\bm{\mathscr{{Y}}}$ is formed implicitly, notably for standard CP [3], but this is not the case for general loss functions and so motivates stochastic approaches.

The stochastic approach is based on the following observation: $\bm{\mathscr{{Y}}}$ becomes sparse when the data tensor $\bm{\mathscr{{X}}}$ is scarce (i.e., most of the elements are missing). This is because, by [22, Theorem 3], the gradient in this case is the same as in Eq. 5 except that $\bm{\mathscr{{Y}}}$ is defined as

[TABLE]

If $\bm{\mathscr{{Y}}}$ is sparse, then the MTTKRP in Eq. 5 can be computed in time $O(r\,\text{nnz}(\bm{\mathscr{{Y}}}))$.

We develop a stochastic gradient whose general form is

[TABLE]

By linearity of expectation, $\mathbb{E}[\bm{\tilde{\mathscr{{Y}}}}]=\bm{\mathscr{{Y}}}$ implies $\mathbb{E}[{\bm{\tilde{\mathbf{{G}}}}}_{k}]={\bm{\mathbf{{G}}}}_{k}$. Making $\bm{\tilde{\mathscr{{Y}}}}$ sparse unlocks efficient sparse computation of the MTTKRP. Consequently, storage for $\bm{\tilde{\mathscr{{Y}}}}$ is $O(s)$ and the cost to compute the gradients is $O(rs)$, a reduction of roughly $n^{d}/s$ compared to computing the full gradient, albeit at the cost of lower accuracy.

In the remainder of this section, we discuss the pros and cons of several different choices for $\bm{\tilde{\mathscr{{Y}}}}$ based on uniform sampling, stratified sampling, and semi-stratified sampling.

4.1 Uniform Sampling

To create a random instance of $\bm{\tilde{\mathscr{{Y}}}}$, we sample $s$ indices uniformly with replacement. The number of times that $i$ is sampled is denoted as $\tilde{s}_{i}$, so $\sum_{i\in\Omega}\tilde{s}_{i}=s$. Using this sampling strategy, in expectation we have

[TABLE]

Note that the $\tilde{s}_{i}$ values are not stored as a dense object but rather as just a list of the samples or in some other sparse data structure. We then define the stochastic tensor $\bm{\tilde{\mathscr{{Y}}}}$ to be

[TABLE]

The stochastic $\bm{\tilde{\mathscr{{Y}}}}$ is sparse because at most $s$ entries are nonzero (since at most $s$ values of $\tilde{s}_{i}$ are nonzero). Clearly $\mathbb{E}[\bm{\tilde{\mathscr{{Y}}}}]=\bm{\mathscr{{Y}}}$ since $\mathbb{E}[\tilde{y}_{i}]=\mathbb{E}[\tilde{s}_{i}]\,\frac{n^{d}}{s}\,y_{i}=y_{i}$ for all $i\in\Omega$. A similar argument applies for sampling without replacement, but there is no practical difference when $s\ll n^{d}$.

The simplest way to generate a random index $i$ is to generate $d$ random mode indices $i_{k}$:

[TABLE]

Here, $\texttt{randi}(m)$ indicates selecting a random integer between 1 and $m$. This requires $d$ random numbers per index $i$. Alternatively, we can generate a random linear index via $\texttt{randi}(n^{d})$ and then convert it to a tensor multi-index.222Generating $d$ separate entries helps prevent overflow if $n^{d}>2^{64}$, i.e., the total number of entries in the tensor is more than the size of the largest unsigned $64$-bit integer (uint64 in MATLAB). Larger integers can be generated by instead using extended-precision arithmetic (e.g., int in Python). Either way is $O(1)$ work (we treat $d$ as a constant).

The procedure is presented in Algorithm 1. The function MTTKRP corresponds to Eq. 3, and this can be computed efficiently because specialized sparse implementations exist as discussed in Section 3.3. The model entries $m_{i}$ and elementwise partial gradient tensor entries $y_{i}$ are only computed for the $s$ randomly selected indices. The most expensive operations are computing the model entries $m_{i}$ and the corresponding MTTKRP calculations. In the implementation, we are able to share some intermediate computations across these two steps.

4.2 Stratified Sampling

Uniform sampling may not be appropriate for sparse tensors since nonzeros will rarely be sampled (each draw gets a nonzero with probability $\text{nnz}(\bm{\mathscr{{X}}})/n^{d}$). Intuitively, however, we expect nonzeros in a sparse tensor to be important to the factorization.

Non-uniform sampling has been used in other stochastic gradient contexts. Needell et al. [36] analyze SGD and argue that biased sampling toward functionals with larger Lipschitz constants can improve performance. Gopal [16] similarly biases sampling toward functionals with larger gradients to reduce the variance of the stochastic gradients, and similar ideas appear in [58, 57]. This idea has many potential applications in machine learning, and roughly translates in our case to up-sampling a tensor entry if the corresponding elementwise loss function $f_{i}\equiv f(x_{i},m_{i})$ has a higher Lipschitz constant for fixed $x_{i}$. Consider GCP with the Bernoulli odds loss [22]: $f(x,m)=\log(m+1)-x\log m$ where $x\in\set{0,1}$ and $m>0$ corresponds to the odds of $x=1$. For $x=0$, $|\frac{\partial f}{\partial m}(0,m)|=1/(m+1)\leq 1$; for $x=1$, $|\frac{\partial f}{\partial m}(1,m)|=1/(m^{2}+m)$ is unbounded as $m\downarrow 0$. Thus, the elementwise loss functions for nonzeros are not Lipschitz, their gradients can be very large, and sampling them more often by sampling zeros and nonzeros separately can reduce the variance of the stochastic gradients.

Consider a generic partition of $\Omega$ into $p$ partitions $\Omega_{1},\Omega_{2},\dots,\Omega_{p}$. If we partition into zeros and nonzeros, then $p=2$. Let $s_{\ell}>0$ be the number of samples from partition $\Omega_{\ell}$ so that the total number of samples is $s=\sum_{\ell}s_{\ell}$. Within each partition, we sample uniformly with replacement. Thus, the expected number of times index $i$ is sampled depends on its partition:

[TABLE]

We then define the stochastic tensor $\bm{\tilde{\mathscr{{Y}}}}$ to be

[TABLE]

Once again, at most $s$ entries of $\bm{\tilde{\mathscr{{Y}}}}$ are nonzero, and $\mathbb{E}[\bm{\tilde{\mathscr{{Y}}}}]=\bm{\mathscr{{Y}}}$ since $\mathbb{E}[\tilde{y}_{i}]=\mathbb{E}[\tilde{s}_{i}]\,\frac{|\Omega_{\ell}|}{s_{\ell}}\,y_{i}=y_{i}$ for all $i\in\Omega$. Interestingly, the weight (${|\Omega_{\ell}|}/{s_{\ell}}$) for index $i$ is completely independent of the overall sample or tensor size — it depends only on the partition size and the number of samples for that partition.

Algorithm 2 presents the stratified sampling procedure for sparse tensors, where we sample $p$ nonzeros and $q$ zeros. Sampling nonzeros is straightforward because they are stored as a list in coordinate format. However, sampling zeros requires sampling a random index as was done for uniform sampling and then rejecting the sample if it is actually a nonzero.

We can estimate how many random multi-indices are needed to obtain $q$ valid zero samples. Let $\eta=\text{nnz}(\bm{\mathscr{{X}}})$ and $\zeta=n^{d}-\eta$. We expect a very small proportion ($\eta/n^{d}$) of indices to be rejected. The number of samples needed to produce an average of $q$ zero indices is

[TABLE]

This is only on average, so we oversample to get sufficiently many zeros with high probability. Namely, we sample $\rho\,(n^{d}/\zeta)\,q$ indices where $\rho>1$ is the oversample rate. A default of $\rho=1.1$ is justified in Appendix B.

Since we have to check against the entire list of nonzeros for each sample, the cost of the rejection sampling for zeros can be significant. To achieve a speed that does not dominate the other costs, our MATLAB implementation uses various efficiencies such as conversion from multi-indices to linear indices, presorting the list of nonzero linear indices, and using the hidden builtin function _ismemberhelper. If $n^{d}\geq 2^{64}$, conversion to linear indices so that we can use these efficiencies is not possible. Beyond MATLAB, efficiency can be achieved using extended precision, hash tables, etc. As an alternative to specialized implementations for efficient rejection sampling (that are potentially language and architecture dependent), we propose a specialized algorithm that entirely avoids rejection sampling in the next subsection.

4.3 Semi-Stratified Sampling

To avoid rejection sampling entirely, we propose a variant that we call semi-stratified sampling. We sample zeros incorrectly but then correct for it when sampling nonzeros. Specifically, we sample without rejection to obtain “zeros”, knowing that a small proportion may actually be nonzeros. We still sample nonzeros explicitly but now add a correction to account for the possibility that they were also wrongly sampled as “zeros.”

Namely, we sample $p$ nonzeros and $q$ “zeros” from $\Omega$ (the entire set of indices). Let $\tilde{p}_{i}$ be the number of times that index $i$ was sampled as a nonzero where $\sum_{i}\tilde{p}_{i}=p$. Clearly, $\tilde{p}_{i}=0$ if $x_{i}=0$. Let $\tilde{q}_{i}$ be the number of times that index $i$ is sampled as a “zero” from the full set of possible indices where $\sum_{i}\tilde{q}_{i}=q$. It is possible that some nonzeros are sampled, i.e., we can have $\tilde{q}_{i}>0$ when $x_{i}\neq 0$. Using these counts and recalling $\eta=\text{nnz}(\bm{\mathscr{{X}}})$, we define $\bm{\tilde{\mathscr{{Y}}}}$ as

[TABLE]

This still satisfies $\mathbb{E}[\tilde{y}_{i}]=y_{i}$. For $i$ such that $x_{i}=0$, we have $\tilde{p}_{i}=0$ and $c_{i}=y_{i}$, so

[TABLE]

For $i$ such that $x_{i}\neq 0$, we have

[TABLE]

The procedure is implemented in Algorithm 3. The procedure for the nonzero samples is identical to that in Algorithm 2 except for the adjustment term $-g(0,m_{i})$ to ensure that the expectations are correct. The procedure for the “zeros” differs because it samples over the entire index space and does not reject nonzeros; it assumes that $x_{i}=0$ in computing $y_{i}$.

4.4 Adapting to Weighted Formulations

Any of the above sampling methods can be easily adapted to the weighted version of GCP with the weighted loss function

[TABLE]

In estimating the gradient, the only change in the methods is to replace $y_{i}$ in Eq. 4 with

[TABLE]

We do not explicitly study the weighted formulation in this work, but we briefly mention two scenarios where such methods may be useful.

Adapting to Weighted Formulations

In the context of recommender systems, zeros are generally treated as missing data. The data is assumed to be missing at random (MAR), meaning that the probability of a data item being present does not depend on its value. However, it has been argued in the matrix case that this may be a flawed assumption [35]. Instead, we may consider including the zero terms but down-weighting them by using a weighted scheme, e.g., $w_{i}=1$ for $x_{i}\neq 0$ and $w_{i}=0.1$ for $x_{i}=0$. Many large-scale tensor applications have a recommender system flavor, where down-weighting of the zero entries may be appropriate. More generally, down-weighting can be appropriate for entries that are less reliable or noisier, as has been done in various matrix factorization applications [13, 23, 47, 56]. See also [46, 48] for work on computing weighted matrix factorizations and [21] for a recent analysis of matrix factorization weights when columns have heterogeneous noise levels.

Weighted Formulations for Missing Data

Conversely, if some portion of data is missing, there are various strategies that can be used to avoid sampling missing elements. However, this can also be handled easily by setting the weights of missing entries to be zero. Ideally, elements with a weight of zero should be avoided during sampling.

5 Experimental Results

All experiments were run using MATLAB (Version 2018a) on a Dual Socket Intel E5-2683v3 2.00GHz CPU with 256 GB memory. The methods are implemented in gcp_opt in the Tensor Toolbox for MATLAB [4].

5.1 Stochastic optimization algorithm

The stochastic gradients can be used with any number of stochastic optimization methods. We use the popular Adam [25] method because it is less sensitive to the learning rate than standard SGD. The method is detailed in Algorithm 4. Based on the empirical results that follow, we recommend setting $s=d\bar{n}$. If the total number of gradient samples is $s$, then we use $s$ samples for uniform sampling in Algorithm 1 and $p=\lfloor s/2\rfloor$ and $q=\lceil s/2\rceil$ for the stratified sampling in Algorithm 2 and semi-stratified sampling in Algorithm 3. The parameter $\alpha$ is the learning rate and defaults to 0.01. The Adam parameters are set to the values recommended in the original paper: $\beta_{1}=0.9$, $\beta_{2}=0.999$, and $\epsilon=10^{-8}$. We employ a few standard modifications. We group the iterates into epochs, and the number of iterations per epoch defaults to $\tau=1000$. To track progress, we estimate the function value $F(\bm{\mathscr{{X}}},\bm{\mathscr{{M}}})$ once per epoch. Whenever the function value fails to decrease, we either decay the learning rate by $\nu=0.1$ (as motivated by [30]) or quit (after more than $\kappa=1$ failures). We enforce any lower bounds $\ell$ given for the parameters by simple projection. In all of our examples, we use $\ell=0$ (nonnegativity constraint).

We estimate $F(\bm{\mathscr{{X}}},\bm{\mathscr{{M}}})$ (via the function EstObj) in a way that is analogous to the stochastic gradient computation. There are two key differences. First, we use a much larger number of samples to ensure reasonable accuracy, which is less important for the gradient computation. Second, the set of samples used for function estimation are fixed across all iterations for consistent evaluation across epochs. For uniform sampling, let $\tilde{s}_{i}$ be the number of times that index $i$ is selected and then estimate

[TABLE]

For stratified sampling, let $\tilde{p}_{i}$ denote the number of times that nonzero $i$ is selected and $\tilde{q}_{i}$ be the same for zero $i$ and then estimate

[TABLE]

In either case, it is easy to show that $\mathbb{E}[\hat{F}]=F$. When we perform multiple runs of the same problem, we use the same set of samples for the function estimation across all runs so that they can be easily compared. The only exception is the non-stochastic method, which computes the full objective function.

5.2 Sample Size and Comparison to Full Method for Dense Tensors

We study the effect of sample size (also known as minibatch size) to understand the relevant trade-offs: larger sample sizes yield lower variance stochastic gradients but higher costs per iteration. We compute the GCP decomposition on an artificial four-way tensor of size $200\times 150\times 100\times 50$ and rank $r=5$ using the gamma loss function: $f(x,m)=x/m-\log m$ with a nonnegativity constraint on the factor matrices. Appendix C provides the details of the data generation.

Figure 1

shows the results of GCP-Adam with various sample sizes ranging from $s=125$ to $s=2000$. For comparison, we also include non-stochastic results based on the bound-constrained limited-memory BFGS method [8] using full gradients; this optimization approach is standard in MATLAB toolboxes such as Tensor Toolbox [1, 4] and TensorLab [52]. The same set of 25 initial guesses is used for every instance. The initial guesses comprise factor matrices with entries drawn uniformly from $(0,1)$. For GCP-Adam, we estimate the loss $\hat{F}$ using 100,000 uniformly sampled entries that are fixed across all epochs and trials. For the stochastic gradient, we use uniform sampling.

In the top two subfigures, we plot the function value versus time. In Fig. 1a, we consider just the stochastic methods and see the variation between them. There are $d\bar{n}r=2500$ free parameters and $\max_{k}n_{k}=200$. Common wisdom is to make one pass through the data per epoch, which would require $s=n^{d}/\tau=150,000$ samples per iteration. However, we see that two orders of magnitude fewer samples are needed in practice, arguably due to the low-rank structure in the data. Another option is to set $s$ large enough so that we get at least one sample per row of ${\bm{\tilde{\mathbf{{Y}}}}}_{(k)}$ and thus every row of ${\bm{\mathbf{{A}}}}_{k}$ is updated. At a minimum, therefore, we may want $s\geq\max_{k}n_{k}=200$.333This does not guarantee that every row is updated every time. To do so with high probability, one might use roughly $10\,\max_{k}n_{k}$ samples. From [27, Appendix A]: To sample $\rho t$ distinct “types” (i.e., row indices) from a set of $t$ types where $\rho\in(0,1)$, the expected number of draws is $t\log(1/(1-\rho))+O(1)$. Therefore, the expected number of samples needed to collect 99.99% of the members of a set with $t$ types is less than $10t$. Fewer samples generally corresponds to less progress per epoch; however, sometimes fewer samples is advantageous because it takes a different path to the solution. For instance, we see that $s=2000$ initially makes better progress, but $s=500$ and $s=1000$ find the minimum more quickly. At the other extreme, $s=125$ is the lowest cost per epoch, but its progress in reducing the loss is hindered by too few samples. In Fig. 1b we show a longer time range on the x-axis so that the non-stochastic method is visible. The non-stochastic method is approximately an order of magnitude slower.

The final function value may not tell the whole story, so we provide another measure of success. Each instance is run with 25 random starts. Figure 1c shows the number of times that each instance recovers the true solution, meaning that the cosine similarity score between the true solution and the recovered solution is at least 0.9. (See Appendix E for details of the cosine similarity score.) The only time the true solution is not recovered is one run for $s=125$ samples per gradient.

Figure 1d shows box plots of the cost per epoch, which is 1000 stochastic gradient evaluations (with $s$ specified on the x-axis) plus one function value estimation with 100,000 samples. In each box plot, the middle line indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually as red ’+’ symbols. Clearly, the cost per epoch grows linearly with the number of samples, but there are fixed costs that dominate the per iteration cost. Notably, $s=2000$ uses $16$ times as many samples as $s=125$ but is generally only around twice as slow.

5.3 Application to Dense Gas Measurements Tensor

This section illustrates the benefits of using stochastic gradients and the effect of sample size on GCP decomposition for dense real data. We consider a tensor based on chemo-sensing data collected by Vergara et al. [50].444Available at http://archive.ics.uci.edu/ml/datasets/Gas+sensor+arrays+in+open+sampling+settings. The dataset contains data for 11 gases; we used 7 that have more distinctive behaviors. It also has 6 sensor positions; we used position 3 (middle of the wind tunnel) where some of the interesting behaviors occur. The data is recorded at $\sim$100 Hz. We downsampled to 5 Hz, using the nearest measurement when needed, and skipped the first 9 seconds. We also removed sensor 33, which seemed to have erratic measurements. The dataset consists of measurements as a gas is blown over an array of conductometric metal-oxide sensors in a wind tunnel. The tensor modes correspond to 71 sensors, 1250 time points, 5 temperatures, and 140 trials (7 gases with 20 repeats each). This is a dense, relatively small 0.5 GB tensor. We note that Vervliet and De Lathauwer [51] and Battaglino, Ballard and Kolda [6] considered a similar 2 GB tensor derived from the same original dataset with more time points but only three gases; however, we do not compare directly to their approach because they focus on standard CP tensor decomposition. For the loss, we use $\beta$-divergence with $\beta=1/2$ yielding the loss function $f(x,m)=2m^{1/2}+2xm^{-1/2}$ with nonnegativity constraints. This loss function is not necessarily optimal, but it seemed to work well for this data in our experiments. Furthermore, it is an attractive choice since the tensor is nonnegative and may have some outliers. We use rank $r=4$ because this is the smallest rank that sufficiently distinguished the different gases in our experiments.

We run GCP-Adam using uniform sampling and the non-stochastic GCP under the same experimental conditions as in the previous subsection; initial guesses are scaled to match the magnitude of the data tensor.

Figure 2 shows the results of the best overall run in terms of the final objective value. The components are ordered by magnitude, the sensor mode is normalized to the magnitude of the component, and the other modes are normalized to unit length. The trial mode is color coded by gas. The first two components focus largely on sensor variations due to temperature. The final two components capture some temporal patterns for each gas that impact sensors near the middle of the array, where the gas is likely most concentrated. The factorization identifies generally smooth temporal profiles and tends to group the same gas (indicated by color) in the trial mode.

Figure 3 shows performance results.

Figure 3a shows boxplots of the runtimes for the varying sample sizes and the non-stochastic GCP (“full”). The smaller sample sizes take longer because they converge more slowly even though each epoch is cheaper. The non-stochastic instance is overall slowest. Figure 3b shows the range of objective function values for each instance, which is very small. For the stochastic instances, larger sample sizes tend to achieve a marginally better loss, as indicated by improved median and percentile losses. The stochastic instances with $s\geq 500$ samples per gradient perform at least as well as the non-stochastic instance, while being much faster. Figure 3c shows performance on a clustering task. For each run (i.e., a given random start and instance), the rows of the fourth factor matrix are clustered via k-means555We used 500 replicates with MATLAB’s built-in kmeans command to avoid local minima. and we measure what percentage of the 140 trials (7 gases with 20 trials each) get clustered correctly. Larger sample sizes have generally better clustering performance here. Non-stochastic GCP performs similarly to $s=2000$ samples per gradient.

5.4 Sample Size for Sparse Tensors

We consider a four-way sparse binary tensor of size $200\times 150\times 100\times 50$ and rank $r=5$ generated according to an odds model tensor, i.e., where $m_{i}$ gives the odds that $x_{i}=1$. The procedure to generate the data is described in detail in Appendix D. The factor matrices in the solution $\bm{\mathscr{{M}}}$ have $(r-1)$ sparse columns and one constant column. The result is a tensor that has 150,452 ‘structural’ nonzeros (from the sparse columns) and 374,435 ‘noise’ nonzeros (from the dense column), with an overall density of 0.35%. We use the loss corresponding to Bernoulli data with an odds link, i.e., $f(x,m)=\log(m+1)-x\log m$ and a nonnegativity constraint on the factor matrices.

The results of GCP-Adam with various sample sizes ranging from $s=125$ to $s=2000$ are shown in Fig. 4. We use stratified sampling to compute the stochastic gradient, and the gradient samples are evenly divided between zeros and nonzeros.

The same set of 25 initial guesses is used for every instance. The initial guesses comprise factor matrices with entries drawn uniformly from $(0,1)$ and then scaled to match the magnitude of the true solution tensor. For GCP-Adam, we estimate the loss $\hat{F}$ using 200,000 stratified sampled entries (evenly divided between zeros and nonzeros) that are fixed across all epochs and trials.

Figure 4a plots the individual runs. For the stochastic method, the overall run time is about four times more than for the (easier) Gamma case. We omit the non-stochastic method since it is again significantly slower. This is arguably a difficult test problem in terms of recovering the true factors, especially compared to the Gamma problem in Section 5.2. From Fig. 4b, observe that $s=125$ fails to find the true solution more often than it succeeds. For $s\geq 250$, the true solution is recovered in the majority of cases, and the recovery rate improves as the number of samples increases; all $25$ runs in this experiment succeeded for $s=2000$. Figure 4c shows boxplots of the time per epoch. As was the case with uniform sampling of a same-sized dense tensor in Section 5.2, the time per epoch mainly consists of costs that grow linearly with the number of samples and fixed costs that, at this scale, remain significant. As before, $s=2000$ uses $16$ times as many samples as $s=125$ but is generally only around twice as slow.

5.5 Comparison of Sampling Strategies for Sparse Tensors

This section shows the benefit of using stratified and semi-stratified sampling over uniform. We use the same procedure as Section 5.4 to create a sparse binary tensor (detailed in Appendix D). We generate a tensor of size $400\times 300\times 200\times 100$ that is 0.38% dense. This example has 4,402,374 ‘structural’ nonzeros and 4,788,052 ‘noise’ nonzeros, with a total of 9,181,549 nonzeros (less than the sum due to overlap). Storing it takes 0.37 GB as a sparse tensor, but would take 19 GB as a dense tensor.

The results of uniform, stratified, and semi-stratified sampling are shown in Fig. 5.

We calculate the estimated loss $\hat{F}$ once per epoch using 200,000 stratified sampled entries, evenly divided between zeros and nonzeros. We use $s=1000$ samples per stochastic gradient evaluation ($s=d\bar{n}$), evenly divided between nonzeros and zeros (or “zeros”) for the stratified samplers. We use 25 initial guesses with random positive values, scaled so that the norm of the initial guess is the same as the norm of the tensor. Note that uniform sampling on a sparse tensor is not as fast as it is on a dense tensor since only nonzeros are stored explicitly and every sampled index has to be checked against the list of nonzeros to determine its value.

The proportion of nonzeros is arguably higher than what is observed in many real-world datasets, which is a favorable condition for uniform sampling since every sample will more likely include nonzeros. Nevertheless, stratified and semi-stratified approaches clearly outperform uniform sampling. They converge faster and more often find the true solution, all at less cost per epoch. This result is expected because stratified sampling should reduce the variance. The time advantage of the semi-stratified approach is minimal for this small example. However, the speed of the MATLAB implementation of stratified sampling depends on the ability to use linearized indices, which means that the total size of the tensor must be less than $2^{64}$. For larger tensors where we cannot use this approach, the sampling efficiency can degrade by more than an order of magnitude. The semi-stratified approach has no such limitation.

Table 1 provides empirical validation that the variances of the stratified and semi-stratified sampling are lower than for uniform. We consider the variances with respect to the vectorized gradient, i.e., ${\bm{\mathbf{{g}}}}=[\textsc{vec}({\bm{\mathbf{{G}}}}_{1});\,\textsc{vec}({\bm{\mathbf{{G}}}}_{2});\,\dots;\,\textsc{vec}({\bm{\mathbf{{G}}}}_{d})]\in\mathbb{R}^{d\bar{n}r}$. Let ${\bm{\tilde{\mathbf{{g}}}}}$ denote the random variable representing the vectorized stochastic gradient, and recall that $\mathbb{E}[{\bm{\tilde{\mathbf{{g}}}}}]={\bm{\mathbf{{g}}}}$ because the sampling methods lead to unbiased estimators. Letting ${\bm{\tilde{\mathbf{{g}}}}}_{\xi}$ denote the $\xi$-th of $N$ realizations of the random variable, the quantities of interest are

[TABLE]

[TABLE]

Table 1 considers two distinct cases: an initial guess (far from the solution with a larger gradient) and the overall best final solution (close to the solution with a smaller gradient). In both cases, stratified and semi-stratified sampling had lower variances than uniform. Consequently, these have lower empirical biases as well. This helps to explain the superior convergence of the stratified and semi-stratified approaches.

5.6 Application to Sparse Count Crime Data and Comparison to CP-APR

We consider a real-world crime statistics dataset comprising more than 15 years of crime data from the city of Chicago. The data is available at www.cityofchicago.org, and we downloaded a 4-way tensor version from FROSTT [42]. The tensor modes correspond to 6,186 days from 2001 to 2017, 24 hours per day, 77 communities, and 32 types of crimes. Each $\bm{\mathscr{{X}}}(i,j,k,\ell)$ is the number of times that crime $\ell$ happened in neighborhood $k$ during hour $j$ on day $i$. The tensor has 5,330,673 nonzeros. Stored as a sparse tensor, it requires 0.21 GB of storage.

Since this is count data, we use GCP with a Poisson loss function, i.e., $f(x,m)=m-x\log m$ and nonnegativity constraints. We run GCP-Adam with both stratified and semi-stratified sampling using $s=d\bar{n}=6319$ samples for the stochastic gradient. We compare to the state-of-the-art for CP alternating Poisson Regression (CP-APR) [12], using the Quasi-Newton method described in [18]. We run each method with 20 different starting points. We compute rank $r=10$ factorizations.

Timing results are shown in Fig. 6. The CP-APR method has to do extensive pre-processing, which is why the first iteration does not complete until approximately 140 seconds. The GCP-Adam methods descend much more quickly but do not reduce the loss quite as much, though this failure to achieve the same final minima is likely an artifact of the function estimation and/or nuance of the Adam parameters.

Although the final loss functions values are slightly different, the factorizations computed by the three methods are similar. We show the results from the first random starting point in Fig. 7.

Each rank-1 component corresponds to a row in the figure. The first column is the day, shown as a line graph. Components in this column are scaled to capture the magnitude of the overall component, while all the other components are normalized. The second column is hour of the day, shown as a bar graph. The third column is the neighborhood, shown as a bar graph. (The third column is somewhat difficult to interpret visualized this way but we show it on a map in Appendix F.) The fourth column is the crime type, sorted by overall frequency, and only showing the 13 most prevalent crimes. When working with real-world data, we generally have to experiment. Nevertheless, certain trends emerge over and over again for different starting points, ranks, and methods. In this case, we see strong commonalities among the three methods, as follows. Component 1 for CP-APR is similar to component 2 for the GCP-Adam methods, with “theft” being the most prevalent crime and a similar pattern in time. Component 2 for CP-APR is similar to component 1 for the GCP-Adam methods. Component 3 for the semi-stratified GCP-Adam has a very strong seasonal signature, becoming most active in summer months. Component 4 for stratified GCP-Adam is similar, as is component 6 for CP-APR. Component 8 in CP-APR and component 6 in both GCP-Adam methods has a special pattern of a spike on the first of each year and again on the first of each month. There is also a spike at midnight in the hour mode. This is likely a feature of how the associated crimes were recorded in the dataset. To help further with interpretation of the rank-1 components, we zoom in on the components from the semi-stratified GCP-Adam solution in Appendix F, where we show each individual component, including drawing a heatmap of the neighborhoods on a map.

6 Conclusions

We propose a stochastic gradient for GCP tensor decomposition with general loss functions. The structure of the GCP gradient means that there is no general way to maintain sparsity in its computation even when the input tensor is sparse. Our investigation was prompted by the sparse case, but the stochastic approach applies equally well to dense tensors. A unique feature of our approach is the use of stratified and semi-stratified sampling in the gradient computation for sparse tensors.

We tested the stochastic gradient using Adam [25] for GCP tensor decomposition and made several findings. Stochastic gradient methods are effective in practice in terms of driving down the objective function and recovering the true solutions. Empirically, we find that the number of samples should be roughly equal to the sum of the dimensions, i.e., $s=\sum_{k=1}^{d}n_{k}=d\bar{n}$. This is much lower than would be required to cycle through the entire dataset each epoch, i.e., $n^{d}$ divided by the number of iterations per epoch. For dense problems, stochastic gradient methods can be much faster than the non-stochastic prior approach using L-BFGS-B [22]. For sparse problems, stochastic gradients enable us to circumvent formation of the dense tensor needed by the gradient, making it possible to solve much larger problems. Additionally, we propose stratified and semi-stratified sampling, which are typically superior to uniform sampling. We have not discussed how to determine the rank since that is a difficult problem even for standard CP.

Overall, stochastic methods have proved to be a promising approach for GCP tensor decomposition, especially for large-scale sparse tensors which have no viable alternative. However, many open questions remain. We can likely further improve the results by using more sophisticated stochastic optimization methods, e.g., weight decay strategies in Adam [30]. Likewise, more sophisticated sampling strategies such as leverage score sampling from matrix sketching [33, 55, 10] may further improve performance by reducing the variance. Another important line of investigation is developing appropriate theory to describe the improvement gains of the stratified approaches.

In terms of implementations, an interesting consequence of sampling in the context of parallel tensor decomposition [44, 24, 29, 40] is that we can reduce the computation and/or communication by sampling only a subset of the entries. Moreover, we may be able to stratify the samples in such a way that is amenable to more structured communications.

Appendix A Special Cases where Gradient Does Not Require Dense Calculations

Computing the gradient is not a major issue for standard CP due to its special structure. Specifically, the computation of ${\bm{\mathbf{{G}}}}_{k}$ can be simplified so that the primary work is computing ${\bm{\mathbf{{X}}}}_{(k)}{\bm{\mathbf{{Z}}}}_{k}$, which is a sparse MTTKRP whenever $\bm{\mathscr{{X}}}$ is sparse [22, Appendix A]. For Poisson CP [12], the primary work is computing ${\bm{\mathbf{{V}}}}_{\!(k)}{\bm{\mathbf{{Z}}}}_{k}$ where $\bm{\mathscr{{V}}}$ is the sparse tensor defined as

[TABLE]

In these cases, the gradient can be computed in $O(\text{nnz}(\bm{\mathscr{{X}}})\,rd)$ flops with $O(\text{nnz}(\bm{\mathscr{{X}}}))$ additional storage. Generally, however, we may not have such structure and we have to compute with a dense $\bm{\mathscr{{Y}}}$ tensor at a cost of $O(rn^{d})$ flops and $O(\text{nnz}(\bm{\mathscr{{X}}}))$ extra storage.

Appendix B Determining the Oversampling Rate

Section 4.2 mentions that we oversample to get sufficiently many zeros with high probability. Namely, we sample

[TABLE]

indices to get the desired $s_{\text{\tiny zero}}$ zeros, where $p_{\text{\tiny zero}}=1-\text{nnz}(\bm{\mathscr{{X}}})/n^{d}$ is the proportion of zeros in the tensor. Here we discuss the oversampling rate $\rho$.

We can use the inverse cumulative distribution function (CDF) of a negative binomial distribution to determine an appropriate $\rho$. The negative binomial distribution models the number of failures before a given number of successes with a specified success rate. In our case, we want $s_{\text{\tiny zero}}$ successes and the success rate is $p_{\text{\tiny zero}}$. We can use the inverse CDF to determine the number of rejections at the 99.9999% percentile. For instance, in MATLAB:

[TABLE]

This means that with probability 0.999999, no more than $s_{\text{\tiny reject}}$ nonzeros will be drawn before $s_{\text{\tiny zero}}$ zeros are obtained. So we want to choose the oversampling rate $\rho$ so that

[TABLE]

In Fig. 8, we plot the oversample rate that would be needed in different scenarios. The x-axis is the proportion of nonzeros. The y-axis is the oversample rate. We plot four lines corresponding to different values for $s_{\text{\tiny zero}}$. We observe two trends:

For fixed $p_{\text{\tiny zero}}$, $\rho$ decreases as $s_{\text{\tiny zero}}$ increases, and 2. 2.

For fixed $s_{\text{\tiny zero}}$, $\rho$ decreases as $p_{\text{\tiny zero}}$ increases.

For most real-world examples of sampling zeros from a sparse tensor, $p_{\text{\tiny zero}}\geq 0.99$ because the tensors are extremely sparse. Additionally, we usually use a sample size of at least $s_{\text{\tiny zero}}=1000$. Thus, oversampling by $\rho=1.1$ should be adequate for most scenarios we expect to encounter.

We could also determine the oversampling rate for any individual problem using this procedure, but the inverse CDF calculation can be expensive.

Appendix C Creating Gamma-Distributed Test Problems

To create the Gamma-distributed test problem used in Section 5.2, we generate factor matrices whose entries are drawn from the uniform distribution on $(0,1)$:

[TABLE]

Using these factor matrices, we create $\bm{\mathscr{{M}}}_{\rm true}$. The data tensor $\bm{\mathscr{{X}}}$ is generated as

[TABLE]

Appendix D Details of Creating Binary Test Problems

We assume an odds link with the data, so the factor matrices must be nonnegative. The probability of a one is given by $m/(1+m)$ where $m$ corresponds to the odds. For simplicity in the model and in generating the data tensor, we assume that factors 1 through $(r-1)$ are relatively sparse (i.e., sparsity specified by $\delta\in(0,1/2)$ and factor $r$ is dense. The idea here is the last dense component corresponds to noise in the model, i.e., random but infrequent observations of ones. Otherwise, the ones have a pattern as dictated by the sparse components.

We specify a density of factor matrix nonzeros and a probability of a one for nonzero values in the resulting model, denoted $\rho_{\rm high}$. To obtain that probability, the nonzero factor matrix entries should be $\sqrt[d]{\rho_{\rm high}/(1-\rho_{\rm high})}$. We modify that slightly by setting the nonzero factor matrix entries to be drawn from a normal distribution with the mean as the target value and a standard deviation of 0.5. Since only a few entries are nonzero, we can identify all the possible nonzeros corresponding to the first four factors and then compute the exact probability computed by the model and then generate an observation.

For the final dense component, we want the probability of a one, denoted $\rho_{\rm low}$, to be relatively low. This means that approximately $\rho_{\rm low}$ of the data tensor entries will correspond to this last “noise” column. The entries of the factor matrix are set to $\sqrt[d]{\rho_{\rm low}/(1-\rho_{\rm low})}$. We use this value exactly so that we can generate nonzero “noise” observations in bulk.

For the test problems in Section 5.2, we use $\delta=0.15$, $\rho_{\rm high}=0.9$, and $\rho_{\rm low}=0.0025$. For the test problems in Section 5.5, we use $\delta=0.15$, $\rho_{\rm high}=0.9$, and $\rho_{\rm low}=0.002$.

Appendix E Cosine Similarity Score

If the true factor matrices are known, we can compute a cosine similarity score between the true and recovered solutions. If the true solution is denoted by ${\bm{\mathbf{{A}}}}_{k}$ and the estimated solution is ${\bm{\hat{\mathbf{{A}}}}}_{k}$, then the cosine similarity score is

[TABLE]

where $\pi$ is a permutation that should yield the highest possible similarity. Recall that the cosine of two vectors ${\bm{\mathbf{{a}}}}$ and ${\bm{\mathbf{{b}}}}$ is ${\bm{\mathbf{{a}}}}^{\intercal}{\bm{\mathbf{{b}}}}/(\|{\bm{\mathbf{{a}}}}\|_{2}\|{\bm{\mathbf{{b}}}}\|_{2})$. We say that the true solution is recovered if the similarity score is at least 0.9. Assume that $\bm{\mathscr{{M}}}$ holds ${\bm{\hat{\mathbf{{A}}}}}_{k}$ and $\bm{\mathscr{{M}}}_{\rm true}$ holds ${\bm{\mathbf{{A}}}}_{k}$, the cosine similarity is computed using the Tensor Toolbox for MATLAB [4] via the following command:

{boxedverbatim}

score(M,Mtrue,’lambda_penalty’,false)

Appendix F Individual components of Chicago crime tensor factorization

In this appendix, we show the remaining 10 components for the factorization of the Chicago crime tensor discussed in Section 5.6 in Figs. 9, 10, 11, 12, 13, 14, 15, 16, 17 and 18. Here we have scaled the date to show the overall weight of the component, and the other components are normalized.

Appendix G Acknowledgments

We thank the referees for their helpful comments on an earlier draft of this paper, including bringing additional references to our attention. This paper was the result of a broader multi-year collaboration involving several institutions. We gratefully acknowledge Jed Duersch for preliminary investigations on handling large-scale tensors with GCP and many insightful discussions about the work presented in this paper. Thanks also to Cliff Anderson-Bergman for stimulating discussions about this project, including providing the information for Appendix B. We thank our colleague Eric Phipps for proposing an idea that ultimately led to semi-stratified sampling.