Analysis of the Stochastic Alternating Least Squares Method for the Decomposition of Random Tensors^††thanks: Submitted to the editors 04/01/2020.

Multi-modal data is represented by a tensor, or a multidimensional matrix. Tensor data is present in ares such as Natural Language Processing (NLP) [8], Blind Source Separation (BSS) [9, 28], and Phylogenetic Tree Reconstruction (PTR) [12, 2, 17]. In each of these areas, canonical decomposition (CANDECOMP)/parallel factors (PARAFAC), also known as CP tensor decomposition (representing a given tensor as a sum of rank-one tensors), provides important insights since the components of the rank-one terms, the factors, represent meaning in the data. For NLP given a tensor constructed from large text databases, the rank-one terms could represent topics, and the factors could represent words. For BSS given a tensor constructed from signals arriving at a single receiver from unknown sources, the rank-one terms could represent sources, and the factors could be used to locate the sources. For PTR given a tensor constructed as a contingency table for instances of neucleotide combinations from aligned DNA from several species, the factors represent model parameters for the phylogenetic tree.

In applications observations may represent samples of an underlying random tensor. For example, the word co-occurence frequencies used in NLP may come from a sample of a collection of documents. Here fast, reliable algorithms are desired for the CP decomposition of the average tensor. Because CP decomposition is known to be NP hard in general [15], and because the average tensor is usually only available through sampling, we investigate approximate decompositions based on stochastic optimization. In [22], Maehara et al. propose a stochastic version of the widely-used alternating least squares (ALS) method, the stochastic alternating least squares (SALS) method, and show the algorithm’s efficiency by means of numerical experiments. In this paper, we provide a detailed analysis of the algorithm, showing under mild regularization and a minimal set of verifiable assumptions that it is globally convergent to a local minimizer for any initial guess. In Section 4, we also include a detailed discussion the algorithm’s complexity and efficiency.

The alternating least squares (ALS) method, first proposed by Carroll and Chang [10], remains the most widely used algorithm for computing the CP decomposition of a tensor [13]. This block-nonlinear Gauss-Seidel method [5, 26] exploits the component-wise quadratic structure of the cost functional to compute iterates efficiently and with a low memory footprint. It has been modified [11] to exploit sparsity structure in data, which typically occurs in tensors representing co-occurence frequencies such as in the tree reconstruction above. Regularized versions of the ALS method, as well as the associated proximal minimization iterations considered in [21], help mitigate the potential ill-posedness of the CP problem to within sample noise. Although one may compute the tensor’s expectation a priori and decompose it by means of the standard ALS method (see Remark 2.1), such an approach results in a loss of sparsity and cannot seemlessly acommodate the arrival of new samples. In [22], Maehara et al. proposed the Stochastic Alternating Least Squares method, a block-stochastic minimization method that preserves the salient features of the ALS method, while also efficiently incorporating tensor samples in the optimization procedure. Recent work (see [3, 20]) has considered the related problem of using randomized (or stochastic) methods to decompose existing large-scale tensors. In particular, in [3], a variant of the ALS method is developed that approximates the component least squares problem efficiently by randomized/sketching methods. In [20], a dense tensor is expressed as the expectation of a sequence of sparse ones, thereby allowing the use of stochastic gradient descent (SGD) methods.

The convergence analysis of the SALS algorithm applied to the CP decomposition of a random tensor is complicated by the fact that the underlying cost functional is not convex. Moreover, the algorithm itself is a stochastic, block-iterative, second order method whose successive iterates do not have the Markovian dependence structure present in classical stochastic gradient descent (SGD) methods. The convergence of the ALS method was studied in [6], where a quotient-linear convergence rate was established. Block-coordinate techniques for the unconstrained optimization of general, possibly nonconvex, deterministic functionals were treated in [14] (see also [25, 4, 32] and references therein). Xu and Yin [31] discuss the convergence of block stochastic gradient methods for averages of convex and nonconvex functionals. They rely on assumptions (such as the uniformly bounded variance of gradient iterates) that, while standard in the literature (see e.g. [7]), are difficult to verify in practice since they pertain to the algorithm itself. The main contributions of this paper are to prove the convergence of the SALS algorithm, a block-stochastic Newton method, for the CP decomposition of the average of a random tensor, subject to a single verifiable assumption relating to the boundedness of the observed data.

This paper is structured as follows. In Section 2, we introduce the CP decomposition problem for random tensors and describe the stochastic alternating least squares algorithm (Algorithm 2). Section 3 contains our main theoretical contributions. In Section 3.1, we exploit the special multinomial structure of the cost functional to quantify the regularity of its component-wise gradient and Hessian in terms of the size of the algorithm’s iterate vectors (see Lemma 3.5 and its corollaries). In Section 3.2 we show that the iterates themselves can be bounded in terms of the random tensor to be decomposed (Lemma 3.15), which naturally leads to our single, verifiable assumption on the latter’s statistical properties (3.16). In Section 3.3, we show that the iterates generated by the SALS algorithm converge to a local minimizer. We validate the SALS method numerically via a few computational examples in Section 4 and offer concluding remarks in Section 5.

We to follow notational conventions that are closely aligned with the literature on the CP decomposition (see [19]), as well as on stochastic optimization (see [7, 31]). Use use lower case letters to denote scalars, bold letters for vectors, uppercase letters for matrices, and uppercase calligraphic letters for tensors. We use superscripts to indicate iteration (or sub-iteration) indices and subscripts for components, i.e. $x_{i}^{k}$ is the $i$-th component at the $k$-th iteration. Multi-indices are denoted by bold letters and sums, products, and integrals over these should be interpreted to be in iterated form. For the block coordinate minimization method in Algorithm 2, it is convenient to express the vector ${\bm{x}}=({\bm{x}}_{1},...,{\bm{x}}_{p})$ in terms of its $i$-th block component ${\bm{x}}_{i}$ and the remaining components, denoted ${\bm{x}}_{i^{*}}:=({\bm{x}}_{1},\ldots,{\bm{x}}_{i-1},{\bm{x}}_{i+1},\ldots,{\bm{x}}_{p})$. Using this notation, we write ${\bm{x}}=({\bm{x}}_{i},{\bm{x}}_{i^{*}})$ with the implicit understanding that the block components are in the correct order. The Frobenius norm of a vector, matrix, or tensor, i.e. the root mean square of its entries, is denoted by $\|\cdot\|$. For $1\leq p\leq\infty$, $\|\cdot\|_{p}$ denotes the standard Euclidean $p$-norm for vectors and the induced matrix $p$-norm for matrices. We use ‘$\circ$’ to denote the outer product, ‘$*$’ for the component-wise (or Hadamard) product, ‘$\odot$’ for the column-wise Khatri-Rao product, and ‘$\otimes$’ for the Kronecker product.

Let $(\Omega,\mathcal{F},d\mu)$ be a complete probability space and let $\mathcal{X}:\Omega\rightarrow\mathbb{R}^{n_{1}\times n_{2}\ldots\times n_{p}}$ be a measurable map, also known as a $p$-th order random tensor. In practice the underlying probability space is rarely known explicitly. Instead, its law can often be observed indirectly through sample realizations $\mathcal{X}^{1},\mathcal{X}^{2},\ldots,\mathcal{X}^{n}$ of $\mathcal{X}$ that are assumed to be independent and identically distributed (iid). We use $\mathbb{E}$ to denote expectation:

for any integrable function $f:\mathbb{R}^{n_{1}\times\ldots\times n_{p}}\rightarrow\mathbb{R}^{m}$.

The rank-$r$ decomposition problem for $\mathcal{X}$ (or its sample realizations) amounts to finding a set of $r$ rank-one deterministic tensors $\{\hat{\mathcal{X}}_{j}\}_{j=1}^{r}$, so that the quantity

is a good overall representation of $\mathcal{X}$. Each rank-one tensor $\hat{\mathcal{X}}_{j}$, $j=1,\ldots,r$, is formed by the outer product $\hat{\mathcal{X}}_{j}={\bm{a}}_{1,j}\circ\ldots\circ{\bm{a}}_{p,j}$ where ${\bm{a}}_{i,j}=(a_{i,j,1},\ldots,a_{i,j,n_{i}})\in\mathbb{R}^{n_{i}}$ for each $i=1,\ldots p$, i.e. $\hat{\mathcal{X}}_{j}\in\mathbb{R}^{n_{1}\times\ldots\times n_{p}}$ is defined component-wise by

We use the mean squared error $\mathbb{E}\left[\|\mathcal{X}-\hat{\mathcal{X}}\|^{2}\right]$ to measure the quality of the representation. Other risk measures may be more appropriate, depending on the application, possibly requiring a different analysis and approximation technique.

For analysis and computation, it is convenient to consider two other representations of the design variable. We define the $i$-th factor matrix $A_{i}=[{\bm{a}}_{i,1},\ldots,{\bm{a}}_{i,r}]$ in $\mathbb{R}^{n_{i}\times r}$ to consist of the $i$-th component vectors of each term in decomposition (1). Let ${\bm{x}}=\mathrm{vec}([A_{1},\ldots,A_{r}])\in\mathbb{R}^{nr}$, with $n=\sum_{i=1}^{p}n_{i}$, be the vectorization of the factor matrices, i.e. the concatenation of their column vectors. We also write ${\bm{x}}$ in block component form as ${\bm{x}}=({\bm{x}}_{1},\ldots,{\bm{x}}_{p})$, where ${\bm{x}}_{i}=\mathrm{vec}(A_{i})\in\mathbb{R}^{rn_{i}}$ for $i=1,\ldots,p$. To emphasize dependence of $\hat{\mathcal{X}}$ on ${\bm{x}}$, we write the following (which also defines $\left\llbracket\cdot\right\rrbracket$)

No efficient closed form solution of the factorization problem exists (even for deterministic tensors) [15]. So it is commonly reformulated as the optimization problem

Tensor decompositions have scale ambiguity. Indeed, for any direction $i=1,..,p$ and any scalars $\beta_{i,1},\ldots,\beta_{i,r-1}>0$, let $\beta_{i,r}=1/\prod_{j=1}^{r-1}\beta_{i,j}$. Then

The optimizer of (2) therefore lies on a simply connected manifold (see [30]), which can lead to difficulties in the convergence of optimization algorithms. To ensure isolated minima that can be readily located, it is common to enforce bounds on the size of the factors, either directly through an appropriate normalization (see [30]), or indirectly through the addition of a regularization term [22]. We pursue the latter strategy, leading to the problem

While the regularization term biases the minimizers of Problem (2) its inclusion is key to guaranteeing well-posedness in the presence of noise. It plays a pivotal role in (i) proving the existence of a minimizer (Lemma 2.2), (ii) ensuring the positive-definiteness of the Hessian (Corollary 3.11), and (iii) guaranteeing the boundedness of iterates in terms of random tensor $\mathcal{X}$ (Lemma 3.15). Heuristic methods are typically used to choose the value of $\lambda$ that balances the bias of the optimizers against stability considerations, the most well-known of which is the Morozov discrepancy principle [24].

Here we discuss the convergence of Algorithm 2. Since Problem (3) is nonconvex, there can in general be no unique global minimizer. Theorem 3.25 establishes the convergence of the SALS iterates ${\bm{x}}^{k,i}$ to a local minimizer in expectation. Yet, the special structure of the CP problem allows us to forego most of the standard assumptions made in the SGD framework. In fact, we only assume boundedness of the sampled data $\mathcal{X}$. Indeed, we show that the regularity of the component-wise problem, as the Lipschitz constant of the gradient, depends only on the norm of the current iterate ${\bm{x}}^{k,i}$, which is in turn bounded by $\|\mathcal{X}\|$ (Lemma 3.15).