Parallel Algorithms for Constrained Tensor Factorization via the Alternating Direction Method of Multipliers

I Introduction

Tensor factorizationIn the literature, the terms factorization and decomposition are often used interchangeably, even though the latter alludes to exact decomposition, whereas the former may include a residual term. has proven useful in a wide range of signal processing applications, such as direction of arrival estimation , communication signal intelligence , and speech and audio signal separation , as well as cross-disciplinary areas, such as community detection in social networks , and chemical signal analysis . More recently, there has been significant activity in applying tensor factorization theory and methods to problems in machine learning research - see .

There are two basic tensor factorization models: parallel factor analysis (PARAFAC) also known as canonical decomposition (CANDECOMP) , or CP (and CPD) for CANDECOMP-PARAFAC (Decomposition), or canonical polyadic decomposition (CPD, again); and the Tucker3 model . Both are sum-of-outer-products models which historically served as cornerstones for further developments, e.g., block term decomposition , and upon which the vast majority of tensor applications have been built. In this paper, we will primarily focus on the CP model.

Whereas for low-enoughE.g., relative to the sum of Kruskal-ranks of the latent factor matrices. Looser bounds can be guaranteed almost-surely. rank CP is already unique ‘on its own,’ any side information can (and should) be used to enhance identifiability and estimation performance in practice. Towards this end, we may exploit known properties of the sought latent factors, such as non-negativity, sparsity, monotonicity, or unimodality . Whereas many of these properties can be handled with existing tensor factorization software, they generally complicate and slow down model fitting.

Unconstrained tensor factorization is already a hard non-convex (multi-linear) problem; even rank-one least-squares tensor approximation is NP-hard . Many tensor factorization algorithms rely on alternating optimization, usually alternating least-squares (ALS), and imposing e.g., non-negativity and/or sparsity entails replacing linear least-squares conditional updates of the factor matrices with non-negative and/or sparse least-squares updates. In addition to ALS, many derivative-based methods have been developed that update all model parameters at once, see and references therein, and for recent work in this direction.

With few recent exceptions, all tensor factorization algorithms were originally developed for centralized, in-memory computation on a single machine. This model of computation is inadequate for emerging big data-enabled applications, where the tensors to be analyzed cannot be loaded on a single machine, the data is more likely to reside in cloud storage, and cloud computing, or some other kind of high performance parallel architecture, must be used for the actual computation.

A carefully optimized Hadoop/MapReduce implementation of the basic ALS CP-decomposition algorithm was developed in , which reported $100$-fold scaling improvements relative to the prior art. The jist of is to avoid the explicit computation of ‘blown-up’ intermediate matrix products in the ALS algorithm, particularly for sparse tensors, and parallelization is achieved by splitting the computation of outer products. On the other hand, is not designed for high performance computing (e.g., mesh) architectures, and it does not incorporate constraints on the factor matrices.

A random sampling approach was later proposed in , motivated by recent progress in randomized algorithms for matrix algebra. The idea of is to create and analyze multiple randomly sub-sampled parts of the tensor, then combine the results using a common piece of data to anchor the constituent decompositions. The downside of is that it only works for sparse tensors, and it offers no identifiability guarantees - although it usually works well for sparse tensors.

A different approach based on generalized random sampling was recently proposed in . The idea is to create multiple randomly compressed mixtures (instead of sub-sampled parts) of the original tensor, analyze them all in parallel, and then combine the results. The main advantages of over are that i) identifiability can be guaranteed, ii) no sparsity is needed, and iii) there are theoretical scalability guarantees.

Distributed CP decomposition based on the ALS algorithm has been considered in , and more recently in , which exploit the inherent parallelism in the matrix version of the linear least-squares subproblems to split the computation in different ways, assuming an essentially ‘flat’ architecture for the computing nodes. Regular (e.g., mesh) architectures and constraints on the latent factors are not considered in .

In this paper, we develop algorithms for constrained tensor factorization based on Alternating Direction Method of Multipliers (ADMoM). ADMoM has recently attracted renewed interest , primarily for solving certain types of convex optimization problems in a distributed fashion. However, it can also be used to tackle non-convex problems, such as non-negative matrix factorization , albeit its convergence properties are far less understood in this case. We focus on non-negative CP decompositions as a working problem, due to the importance of the CP model and non-negativity constraints; but our approach can be generalized to many other types of constraints on the latent factors, as well as other tensor factorizations, such as Tucker3, and tensor completion.

The advantages of our approach are as follows. First, during each ADMoM iteration, we avoid the solution of constrained optimization problems, resulting in considerably smaller computational complexity per iteration compared to constrained least-squares based algorithms, such as alternating non-negative least-squares (NALS). Second, our approach leads naturally to distributed algorithms suitable for parallel implementation on regular high-performance computing (e.g., mesh) architectures. Finally, our approach can easily incorporate many other types of constraints on the latent factors, such as sparsity.

Numerical experiments are encouraging, indicating that ADMoM-based NTF has significant potential as an alternative to state-of-the-art approaches.

The rest of the manuscript is structured as follows. In Section II, we present the NTF problem and in Section III we present the general ADMoM framework. In Section IV, we develop ADMoM for NTF, while in Section V we develop distributed ADMoM for large NTF. In Section VI, we test the behavior of the developed schemes with numerical experiments. Finally, in Section VII, we conclude the paper.

II non-negative tensor factorization

where $f$ is a function measuring the quality of the factorization, ${\bf 0}$ is the zero matrix of appropriate dimensions, and the inequalities are element-wise. A common choice for $f_{\underline{\bf X}}$, motivated via maximum likelihood estimation for $\underline{\bf E}$ with Gaussian independent and identically distributed (i.i.d.) elements, is

Let $\underline{\bf W}=[{\bf A},{\bf B},{\bf C}]$ and ${\bf W}^{(1)}$, ${\bf W}^{(2)}$, and ${\bf W}^{(3)}$ be the matrix unfoldings of $\underline{\bf W}$, with respect to the first, second, and third dimension, respectively. Then,

and $f_{\underline{\bf X}}$ can be equivalently expressed as

These expressions are the basis for ALS-type CP optimization, because they enable simple linear least-squares updating of one matrix given the other two. Using NALS for each update step is a popular approach for the solution of (1), but non-negativity brings a significant computational burden relative to plain ALS and also complicates the development of parallel algorithms for NTF. It is worth noting that the above expressions will also prove useful during the development of the ADMoM-based NTF algorithm.

III ADMoM

ADMoM is a technique for the solution of optimization problems of the form

The augmented Lagrangian for problem (5) is

where $\rho>0$ is a penalty parameter. Assuming that at time instant $k$ we have computed ${\bf z}^{k}$ and ${\bf y}^{k}$, which comprise the state of the algorithm, the $(k+1)$-st iteration of ADMoM is Note that ${\bf y}^{kT}$ is shorthand notation for $\left({\bf y}^{k}\right)^{T}$.

It can be shown that, under certain conditions (among them convexity of $f$ and $g$), ADMoM converges in a certain sense (see for an excellent review of ADMoM, including some convergence analysis results). However, ADMoM can be used even when problem (5) is non-convex. In this case, we use ADMoM with the goal of reaching a good local minimum . Note that this is all we can realistically hope for anyway, irrespective of approach or algorithm used, since tensor factorization is NP-hard .

Let us consider the set-constrained optimization problem

where $g$ is the indicator function of set ${\cal X}$, that is,

Then it becomes clear that (7) can be solved via ADMoM. Assuming that at time instant $k$ we have computed ${\bf z}^{k}$ and ${\bf y}^{k}$, the $(k+1)$-st iteration of ADMoM is

where $\Pi_{\cal X}$ denotes projection (in the Euclidean norm) onto ${\cal X}$.

IV ADMoM for NTF

where, for any matrix argument ${\bf M}$,

The minimization problem in the first line of (11) is non-convex. Using the equivalent expressions for $f_{\underline{\bf X}}$ in (4), we propose the alternating optimization scheme of (12), at the top of the next page. The updates of (12) can be executed either for a predetermined number of iterations, or until convergence.In our implementations, we execute these updates once per ADMoM iteration. We observe that, during each ADMoM iteration, we avoid the solution of constrained optimization problems. This seems favorable, especially in the cases where the size of the problem is (very) large.

Each ADMoM iteration consists of simple matrix operations. Thus, rough estimates of its computational complexity can be easily derived (of course, accurate estimates can be derived after fixing the algorithms that implement the matrix operations).

A rough estimate for the computational complexity of the update of ${\bf A}^{k}$ (see the first update in (12)) can be derived as follows:

$O((K+J)F^{2})$ for the computation of the term $({\bf C}^{k}\odot{\bf B}^{k})^{T}({\bf C}^{k}\odot{\bf B}^{k})+\rho_{\bf A}{\bf I}_{F}$, and $O(F^{3})$ for its Cholesky decomposition. This is because $({\bf C}^{k}\odot{\bf B}^{k})^{T}({\bf C}^{k}\odot{\bf B}^{k})$ $=$ $\left(\left({\bf C}^{k}\right)^{T}{\bf C}^{k}\right)\circledast\left(\left({\bf B}^{k}\right)^{T}{\bf B}^{k}\right)$.

$O(F^{2}I)$ for the computation of the system solution that gives the updated value ${\bf A}^{k+1}$.

Analogous estimates can be derived for the updates of ${\bf B}^{k}$ and ${\bf C}^{k}$. Finally, the updates of the auxiliary and dual variables require, in total, $O\left((I+J+K)F\right)$ arithmetic operations.

IV-B Convergence

Let $\{\boldsymbol{Z}^{k}\}$ be a sequence generated by ADMoM for NTF that satisfies condition

Then, any accumulation point of $\{\boldsymbol{Z}^{k}\}$ is a KKT point of problem (8). Consequently, any accumulation point of $\{{\bf A}^{k},{\bf B}^{k},{\bf C}^{k}\}$ is a KKT point of problem (1).

Proof: The proof follows closely the steps of the proof of Proposition 2.1 of and is omitted.See report for a detailed proof. $\Box$

Proposition 1 implies that, whenever $\{\boldsymbol{Z}^{k}\}$ converges, it converges to a KKT point. We will further discuss ADMoM convergence from a practical point of view in the section with the numerical experiments.

IV-C Stopping criteria

The primal residual for variable ${\bf A}^{k}$ is defined as

Analogous conditions apply for the other residuals. Reasonable values for $\epsilon^{\rm rel}$ are $\epsilon^{\rm rel}\lessapprox 10^{-3}$, while the value of $\epsilon^{\rm abs}$ depends on the scale of the values of the latent factors.

We note that stopping criteria (17) and (18) involve quantities of the size of the latent factors which, in most cases, is small compared to the size of the tensor. Thus, their computation, even during every ADMoM iteration, is not computationally demanding.

IV-D Varying penalty parameters

We have found very useful in practice to vary the values of each one of the penalty parameters, $\rho_{\bf A}$, $\rho_{\bf B}$, and $\rho_{\bf C}$, depending on the size of the corresponding primal and dual residuals (see [26, Section 3.4]). More specifically, the penalty parameters $\rho_{\bf M}^{k}$, for ${\bf M}={\bf A},{\bf B},{\bf C}$, are updated as follows:

where $\mu>1$, $\tau^{\rm incr}>1$, and $\tau^{\rm decr}>1$ are the adaptation parameters. Large values of $\rho_{\bf M}$ place large penalty on violations of primal feasibility, leading to small primal residuals, while small values of $\rho_{\bf M}$ tend to reduce the dual residuals.

IV-E ADMoM for tensor factorization with structural constraints

Thorough study of ADMoM-based algorithms for tensor factorization and/or completion with more complicated structural constraints is a topic of future research.

V Distributed ADMoM for large NTF

In this section, we assume that all dimensions of tensor $\underline{\bf X}$ are large and derive an ADMoM-based NTF that is suitable for parallel implementation. Of course, our framework can handle the cases where only one or two of the dimensions of $\underline{\bf X}$ are large.

Let $\underline{\bf W}=[{\bf A},{\bf B},{\bf C}]$, and ${\bf A}$, ${\bf B}$, and ${\bf C}$ be partitioned as

We first derive partitionings of the matrix unfoldings of $\underline{\bf W}$ in terms of (the blocks of) matrices ${\bf A}$, ${\bf B}$, and ${\bf C}$. Towards this end, we write ${\bf W}^{(1)}$ as

Thus, ${\bf W}^{(1)}$ can be partitioned as

where the $(n_{A},n_{C})$-th block of ${\bf W}^{(1)}$ is equal to the $I_{n_{A}}\times(JK_{n_{C}})$ matrix ${\bf W}^{(1)}_{n_{A},n_{C}}={\bf A}_{n_{A}}({\bf C}_{n_{C}}\odot{\bf B})^{T}$, for $n_{A}=1,\ldots,N_{A}$ and $n_{C}=1,\ldots,N_{C}$.

Similarly, it can be shown that ${\bf W}^{(2)}$ can be partitioned into blocks ${\bf W}^{(2)}_{n_{B},n_{C}}={\bf B}_{n_{B}}({\bf C}_{n_{C}}\odot{\bf A})^{T}$, of dimensions $J_{n_{B}}\times(IK_{n_{C}})$, for $n_{B}=1,\ldots,N_{B}$ and $n_{C}=1,\ldots,N_{C}$, and ${\bf W}^{(3)}$ can be partitioned into blocks ${\bf W}^{(3)}_{n_{C},n_{B}}={\bf C}_{n_{C}}({\bf B}_{n_{B}}\odot{\bf A})^{T}$, of dimensions $K_{n_{C}}\times(IJ_{n_{B}})$, for $n_{C}=1,\ldots,N_{C}$ and $n_{B}=1,\ldots,N_{B}$.An extension of the above partitioning scheme to higher order tensors appears in Appendix A.

If we partition ${\bf X}^{(1)}$, ${\bf X}^{(2)}$, and ${\bf X}^{(3)}$ accordingly, then we can write

These expressions will be fundamental for the development of the distributed ADMoM for large NTF.

V-B Distributed ADMoM for large NTF

The ADMoM for this problem is as follows:

The minimization problem in the first line of (25) is non-convex. Based on (22), we propose the alternating optimization scheme given in (26) at the next page.

Again, during each ADMoM iteration, we avoid the solution of constrained optimization problems. Furthermore, and more importantly, having computed all algorithm quantities at iteration $k$, the updates of ${\bf A}_{n_{A}}^{k}$, for $n_{A}=1,\ldots,N_{A}$, are independent and can be computed in parallel. Then, we can compute in parallel the updates of ${\bf B}_{n_{B}}^{k}$, for $n_{B}=1,\ldots,N_{B}$, and, finally, the updates of ${\bf C}_{n_{C}}^{k}$, for $n_{C}=1,\ldots,N_{C}$.

We note that we can solve problem (23) using the centralized ADMoM of Section IV. In fact, if we initialize the corresponding quantities of the two algorithms with the same values, then the two algorithms evolve in exactly the same way. As a result, the study (for example, convergence analysis and/or numerical behavior) of one of them is sufficient for the characterization of both.

Thus, via the distributed ADMoM, we simply uncover the inherent parallelism in the updates of the blocks of ${\bf A}^{k}$, ${\bf B}^{k}$, and ${\bf C}^{k}$. In Appendix B, we present a detailed proof of the equivalence of these two forms of ADMoM NTF.

V-C A parallel implementation of ADMoM for large NTF

In the sequel, we briefly describe a simple implementation of ADMoM for large NTF on a mesh-type architecture. In order to keep the presentation simple, we assume that (1) $N_{A}=N_{B}=N_{C}=N$ and (2) each of the matrix unfoldings ${\bf X}^{(1)}$, ${\bf X}^{(2)}$, and ${\bf X}^{(3)}$ has been split into $N^{2}$ blocks, with their $(i,j)$-th blocks stored at the $(i,j)$-th processing element, for $i,j=1,\ldots,N$ (for related results in the matrix factorization context see ).

In Figure 1, we depict the data flow for the computation of the blocks of ${\bf A}^{k+1}$. The inputs to the $N$ top processing elements are ${\bf C}_{n}^{k}$, for $n=1,\ldots,N$, as well as ${\bf B}^{k}$, which is common input to all top processing elements. Each processing element uses its inputs and memory contents and computes certain partial matrix sums. The communications between the processing elements are local and involve either the forwarding of the terms ${\bf C}_{n}^{k}$, for $n=1,\ldots,N$, and ${\bf B}^{k}$ (top-down communication), or the forwarding of the partial sums $\sum_{l=1}^{j}{\bf X}^{(1)}_{n,l}({\bf C}_{l}^{k}\odot{\bf B}^{k})$ and $\sum_{l=1}^{j}({\bf C}_{l}^{k}\odot{\bf B}^{k})^{T}({\bf C}_{l}^{k}\odot{\bf B}^{k})$ (left-right communication), of dimensions $\frac{I}{N}\times F$ and $F\times F$, respectively. The computation of ${\bf A}_{n}^{k+1}$, for $n=1,\ldots,N$, amounts to solution of $\rho$ systems of linear equations with common coefficient matrix and takes place at the rightmost computing elements.

Then, using a similar strategy, we can compute the blocks of ${\bf B}^{k+1}$ and, finally, the blocks of ${\bf C}^{k+1}$. The updates of the auxiliary and dual variables are very simple and can be performed locally (see at the rightmost computing elements of Figure 1).

As we see in Figure 1, in order to compute the blocks of the ${\bf A}^{k+1}$, we use the appropriate blocks of ${\bf C}^{k}_{n}$, for $n=1,\ldots,N$, as well as the whole matrix ${\bf B}^{k}$. When the size of ${\bf B}^{k}$ is not very large, the communication cost is not prohibitive (analogous arguments holds for the computation of the blocks of ${\bf B}^{k+1}$ and ${\bf C}^{k+1}$). Of course, if one or more latent factors are very large, the communication cost significantly increases.

Actual implementation of the distributed ADMoM for large NTF will depend on the specific parallel architecture and programming environment used. Since our aim in this paper is to introduce the basic methodology and computational framework, we leave those customizations and performance tune-ups, which are further away from the signal processing core, for follow-up work to be reported in the high-performance computing literature.

VI Numerical Experiments

In extensive numerical experiments, we have observed that the relative performance of the algorithms depends on the size and rank of the tensor as well as the additive noise power. Thus, we consider $12$ different scenarios, corresponding to the combinations of the following cases:

one, two, or three tensor dimensions are large;

For each scenario, we generate $R=50$ realizations of tensor $\underline{\bf X}$ as follows. We generate random matrices ${\bf A}^{o}$, ${\bf B}^{o}$, and ${\bf C}^{o}$ with i.i.d. ${\cal U}$ elements (using the ${\tt rand}$ command of Matlab) and construct $\underline{\bf X}=[{\bf A}^{o},{\bf B}^{o},{\bf C}^{o}]+\underline{\bf N}$, where $\underline{\bf N}$ consists of i.i.d. ${\cal N}(0,\sigma_{N}^{2})$ elements. For each realization, we solve the NTF problem with (1) NALS (${\tt parafac}$), (2) NLS (${\tt sdf\_nls}$), and (3) ADMoM.

We designed our experiments so that, upon convergence, all algorithms achieve practically the same relative factorization error. Towards this end, we set the values of the stopping parameters as follows: the parameter ${\tt Options(1)}$ of ${\tt parafac}$ is set to ${\tt Options(1)}=10^{-5}$, the parameter ${\tt TolFun}$ of ${\tt sdf\_nls}$ is set to ${\tt TolFun}=10^{-8}$, and the ADMoM stopping parameters are set to $\epsilon^{\rm abs}=10^{-4}$ and $\epsilon^{\rm rel}=10^{-4}$.

In all cases, the initial values of the ADMoM penalty terms are $\rho_{\bf M}=1$, for ${\bf M}={\bf A},{\bf B},{\bf C}$, while the ADMoM penalty term adaptation parameters are $\mu=8$, $\tau^{\rm incr}=4$, $\tau^{\rm decr}=2$.

In practice, convergence properties of ADMoM NTF depend on the (random) initialization point. In some cases, convergence may be quite fast while, in others, it may be quite slow. As we shall see in the sequel, this phenomenon seems more prominent in the cases where rank $F$ is large. In order to overcome the slow convergence properties associated with bad initial points, we adopted the following strategy. We execute ADMoM NTF for up to $n_{\max}=400$ iterations (we have observed that, in the great majority of the cases in the scenarios we examined, this number of iterations is sufficient for convergence when we start from a good initial point). If ADMoM does not converge within $n_{\max}$ iterations, then we restart it from another random initial point; we repeat this procedure until ADMoM converges.Of course, one may think of more elaborate strategies such as, for example, running in parallel more than one versions of the algorithm, with different initializations.

Before proceeding, we mention that all the algorithms converged in all the realizations we run.

Since an accurate statement about the computational complexity per iteration of ${\tt parafac}$ is not easy, the metric we used for comparison of the algorithms is the ${\tt cputime}$ of Matlab. Despite the fact that ${\tt cputime}$ is strongly dependent on the computer hardware and the actual algorithm implementation, we feel that it is a useful metric for the assessment of the relative efficiency of the algorithms.For our experiments, we run Matlab 2014a on a MacBook Pro with a $2.5$ GHz Intel Core i7 Intel processor and $16$ GB RAM. The reason is that we used carefully developed, publicly available Matlab toolbox implementations of the baseline algorithms, and we carefully coded our ADMoM NTF implementation.

In Table I, we present the mean and standard deviation of ${\tt cputime}$, in seconds, denoted as ${\tt mean(t)}$ and ${\tt std(t)}$, respectively, for NALS, NLS, and ADMoM. We also present the mean relative factorization error (which is common to all algorithms up to four decimal digits), defined as

where $\underline{\bf X}_{k}$ is the $k$-th noisy tensor realization and ${\bf A}_{k}$, ${\bf B}_{k}$, and ${\bf C}_{k}$ are the factors returned by a factorization algorithm. Our observations are as follows:

There is no clear winner. Certainly, for high ranks, NLS has very good behavior.

In general, both NALS and NLS have more predictable behavior than ADMoM. Especially for high ranks, the ${\tt cputime}$ of our implementation of ADMoM has large variance.

For small ranks, ADMoM looks more competitive and, in the cases where one dimension is much larger than the other two, it behaves very well (we shall say more on this later).

In order to get a better feeling of the behavior of the three algorithms, we plot their ${\tt cputime}$, along the $50$ realizations we used to obtain the averages of Table I, for two different scenarios. In Figure 2, we consider the case for $I=3000$, $J=K=50$, $F=3$ and $\sigma_{N}^{2}=10^{-2}$. We observe that the behavior of the algorithms is stable, in the sense that there is a clear ordering among the three algorithms, with no large variations. In Figure 3, we keep the dimensions and the noise power the same as before and increase the rank to $F=30$. We observe that the variance of ADMoM ${\tt cputime}$ has significantly increased, while both NALS and NLS show stable behavior. When ADMoM starts from a good initial point, it converges faster than NALS and NLS while, when it starts from bad initial points, it needs one or more restarts.

In order to check if ADMoM maintains its advantage over NALS and NLS in the cases where one dimension is very large, compared with the other two, and the rank is relatively small, we performed an experiment with $I=10^{4}$, $J=K=50$, $F=10$, and $\sigma_{N}^{2}=10^{-2}$. However, in this case, we used somewhat relaxed stopping conditions for all algorithms; more specifically, we used ${\tt Options(1)}=10^{-3}$, ${\tt TolFun}=10^{-6}$, $\epsilon^{\rm abs}=10^{-3}$, and $\epsilon^{\rm rel}=10^{-3}$. In Table II, we present the mean relative factorization errors and the mean and standard deviation of ${\tt cputime}$. As we can see, both NALS and NLS are slightly less accurate than ADMoM, in terms of relative factorization error, which means that their stopping criteria are more relaxed. In terms of ${\tt cputime}$, we see that ADMoM is much faster than both NALS and NLS. In Figure 4, we plot the ${\tt cputime}$ of the three algorithms for the $50$ realizations of the experiment. Again, we see the significant difference between ADMoM and both NALS and NLS. We note that if we had used as values of the stopping parameters those of our initial experiments, then the gain of ADMoM, compared with NALS and NLS, would have been much greater. However, we believe that we have made clear that, in this case, ADMoM has a clear advantage. We have made analogous observations for larger $I$.

VI-B A closer look at ADMoM

In order to get a more detailed view of the convergence properties of ADMoM, we return to the scenario with $I=3000$, $J=K=50$, $F=30$, and $\sigma_{N}^{2}=10^{-2}$, whose ${\tt cputime}$ we plot in Figure 3. We recall that, in order to converge in this case, ADMoM needed often restarts. In Figure 5, we plot the total number of ADMoM iterations, denoted as ${\tt iters}$, and the number of ADMoM iterations during its final way to convergence, which is equal to ${\rm mod}({\tt iters},n_{\max})$. As expected, ${\tt iters}$ is compatible with the corresponding ${\tt cputime}$ (see the red line in Figure 3). Quantity ${\rm mod}({\tt iters},n_{\max})$ shows how many iterations are required for convergence if ADMoM always starts from good initial points. We observe that ${\rm mod}({\tt iters},n_{\max})$ is quite stable around its mean, which is approximately equal to $320$. This gives an estimate of the fastest possible ADMoM convergence in this case.

VI-C ADMoM NTF with under- and over-estimated rank

In the sequel, we consider ADMoM behavior in the cases where we under- or over-estimate the true rank, in both noisy and noiseless cases. Towards this end, we fix $I=J=K=100$ and $F=30$, and investigate ADMoM with exact rank as well as with rank under- and over-estimated by $1$. We expect that, in this case, all versions of ADMoM may need restarts. In the sequel, we examine the influence of under- and over-estimating the rank on (1) factorization accuracy and (2) number of restarts. Of course, in under-modeled cases, we expect that the relative factorization error will be higher than that of the true rank case. However, we know nothing in advance about ADMoM behavior in over-modeled cases. In order to get insight into these issues, we perform the following experiment. We set stopping parameters $\epsilon^{\rm abs}=\epsilon^{\rm rel}=10^{-4}$ and run each of the three versions of ADMoM for $n_{\rm max}=500$ iterations. For each version, we proceed as follows: if it converges within $n_{\rm max}$ iterations, we stop; otherwise, we restart, and repeat until convergence. Thus, finally, the number of iterations for each ADMoM version will be a multiple of $n_{\rm max}$. For the computation of the trajectory of the mean relative factorization error we use only the last $n_{\rm max}$ values; in this way, we avoid the influence of bad initial points. However, we keep count of the restarts of each version and, thus, can assess the time it needs to achieve convergence.

In Figure 6, we plot the average relative factorization errors (computed over $50$ realizations in the way we mentioned before), versus the iteration number, for $\sigma_{N}^{2}=10^{-2}$. As was expected, the ADMoM version with under-estimated rank converges to a higher relative factorization error. We observe that the average relative factorization errors for ADMoM with exact rank and rank over-estimated by $1$ follow almost the same trajectory. The average numbers of restarts for the three ADMoM versions are ${\tt mean}({\rm restarts}_{\tt exact})=1.08$, ${\tt mean}({\rm restarts}_{\tt under})=1.22$, and ${\tt mean}({\rm restarts}_{\tt over})=1.92$. Thus, in the cases of over-estimated rank, we finally achieve a relative factorization error trajectory as good as in the exact rank case, but we may need more restarts and, thus, more time. This implies that the probability of bad initial points may increase.

In Figure 7, we plot the same quantities for noiseless data. Again, the ADMoM behavior in the under-modeled case is as expected. Interestingly, we observe that there is no relative factorization error floor neither for the exact rank nor for the over-estimated by $1$ rank case. Reasonably, after a certain precision level, the over-modeled case converges slower. The average numbers of restarts for the three ADMoM versions are ${\tt mean}({\rm restarts}_{\tt exact})=1.14$, ${\tt mean}({\rm restarts}_{\tt under})=1.54$, and ${\tt mean}({\rm restarts}_{\tt over})=1.06$.

We have observed similar behavior for more drastic rank under- and over-estimation.

VI-D ADMoM for NTF with box-linear constraints

In our final experiment, we briefly consider NTF for the case where two of the latent factors, say ${\bf A}$ and ${\bf B}$, are non-negative while ${\bf C}$ is subject to box-linear constraints in the sense that each row of ${\bf C}$ is a probability mass function, that is, has non-negative elements with sum equal to 1.

The only difference between the ADMoM for this case and the ADMoM for NTF is that, instead of computing ${\bf C}^{k+1}$ as the solution of an unconstrained least-squares problem, we compute it as the solution of linearly constrained least-squares; note that both cases exhibit closed-form solutions.

In Figure 8, we illustrate the behavior of ADMoM in this case by plotting the trajectories of the norms of the average (over $50$ realizations) relative estimation errors of the latent factors, as computed by function ${\tt cpderr}$ of tensorlab, versus the iteration number, for a noiseless case with $I=J=100$, $K=50$ and $F=5$. We observe that ADMoM works to very high precision.

VI-E Discussion

Our numerical results are encouraging and suggest that, in many cases, ADMoM NTF can efficiently achieve close to state-of-the-art factorization accuracy. The fact that ADMoM is suitable for high-performance parallel implementation (the first NTF algorithm with this property, as far as we know) can only increase its potential. Thus, we believe that it will be a valuable tool in the NTF toolbox.

Obviously, in order to fully uncover the pros and cons of ADMoM NTF, more extensive experimentation is required. But our intention in this paper is to give the fundamental ideas and some basic performance metrics. Experiments with real-world data (using ADMoM for tensor completion and factorization) as well as constraints well beyond non-negativity are ongoing work.

A weak point of the version of ADMoM we developed in this manuscript is the high ${\tt cputime}$ variance in cases of high rank. The improvement of the behavior of ADMoM in these cases remains a very interesting problem. To achieve this goal, it might be possible to combine elements of NLS and ADMoM and derive a more efficient algorithm. However, more research efforts are needed in this direction.

As we mentioned, if the centralized and the distributed algorithms start from the same initial point, they evolve in exactly the same way. Thus, distributed ADMoM inherits the convergence properties of centralized ADMoM.

VII Conclusion

Motivated by emerging big data applications, involving multi-way tensor data, and the ensuing need for scalable tensor factorization tools, we developed a new constrained tensor factorization framework based on the ADMoM. We used non-negative factorization of third order tensors as an example to work out the main ideas, but our approach can be generalized to higher order tensors, many other types of constraints on the latent factors, as well as other tensor factorizations and tensor completion. Our numerical experiments were encouraging, indicating that, in many cases, the ADMoM-based NTF has high potential as an alternative to the state-of-the-art and, in some cases, it may become state-of-the-art. The fact that it is naturally amenable to parallel implementation can only increase its potential. The improvement of its behavior in the high rank cases remains a very interesting problem.

Appendix A Extension to higher order tensors

In this appendix, we highlight how our approach can be extended to higher order tensors. We focus on fourth-order tensors, with the general case being obvious. If $\underline{\bf W}=[{\bf A},{\bf B},{\bf C},{\bf D}]$, then its matrix unfoldings satisfy relations

Partitioning matrices ${\bf A}$, ${\bf B}$, ${\bf C}$, and ${\bf D}$ as in subsection V-A, we obtain that matrix ${\bf W}^{(1)}$ can be partitioned into $N_{A}\times N_{D}$ blocks, with the $(i,j)$-th block being equal to

Analogous partitionings apply to the other matrix unfoldings. Then, development of ADMoM NTF (centralized and distributed) is rather easy.

Appendix B On the equivalence of the centralized and the distributed ADMoM NTF

A simple proof of the equivalence of the centralized and the distributed ADMoM NTF is as follows. We focus on the update of ${\bf A}^{k}$ of the centralized algorithm and the updates of its blocks, ${\bf A}^{k}_{n_{A}}$, for $n_{A}=1,\ldots,N_{A}$, of the distributed algorithm, and prove that they are equivalent. We remind that

Using the partitionings of ${\bf C}^{k}\odot{\bf B}^{k}$ (see subsection V-A), it can be shown that

Rewriting the update of ${\bf A}^{k}$ in terms of partitioned matrices, we obtain (27) at the top of this page. If we focus on a certain block of ${\bf A}^{k+1}$ in (27), then we obtain the corresponding update of the distributed algorithm (see (26)). We observe that the matrix inverse in the second line of (27) is common to all blocks, and should be computed once.

Analogous statements hold for the updates of ${\bf B}^{k}$ and ${\bf C}^{k}$. The equivalence of the updates of the rest of the variables is trivial.

Thus, in fact, using the partitionings of subsection V-A, the distributed ADMoM simply uncovered the inherent parallelism of the centralized ADMoM.

References