1 Introduction
Multiway data analysis have recently become increasingly popular supported by modern computational power [23, 37]
. Originally developed in the field of psychometrics and chemometrics, its applications can now also be found in signal processing (for example, for independent component analysis)
[14], neuroscience [29], and data mining [28]. Decomposition of multiway arrays (or tensors) into small number of factors have been one of the main concerns in multiway data analysis, because interpreting the original multiway data is often impossible. There are two popular models for tensor decomposition, namely the Tucker decomposition [44, 13] and the CANDECOMP/PARAFAC (CP) decomposition [11, 19]. In both cases, conventionally the estimation procedures have been formulated as nonconvex optimization problems, which are in general only guaranteed to converge locally and could potentially suffer from poor local minima. Moreover, a popular approach for Tucker decomposition known as the higher order orthogonal iteration (HOOI) may converge to a stationary point that is not even a local minimizer [15].Recently, convex formulations for the estimation of lowrank matrix, which is a special case of tensor, have been intensively studied. After the pioneering work of Fazel et al. [16], convex optimization has been used for collaborative filtering [38], multitask learning [3], and classification over matrices [40]. In addition, there are theoretical developments that (under some conditions) guarantee perfect reconstruction of a lowrank matrix from partial measurements via convex estimation [10, 32]. The key idea here is to replace the rank of a matrix (a nonconvex function) by the socalled trace norm (also known as the nuclear norm) of the matrix. One goal of this paper is to extend the tracenorm regularization for more than two dimensions. There have recently been related work by Liu et al. [27] and Signoretto et al. [36], which correspond to one of the proposed approaches in the current paper.
In this paper, we propose three formulations for the estimation of low rank tensors. The first approach is called “as a matrix” and estimates the lowrank matrix that is obtained by unfolding (or matricizing) the tensor to be estimated; thus this approach basically treats the unknown tensor as a matrix and only works if the tensor is lowrank in the mode used for the estimation. The second approach called “constraint” extends the first approach by incorporating the trace norm penalties with respect to all modes simultaneously. Therefore, there is no arbitrariness in choosing a single mode to work with. However, all modes being simultaneously lowrank might be a strong assumption. The third approach called “mixture” relaxes the assumption by using a mixture of tensors, where is the number of modes of the tensor. Each tensor is regularized to be lowrank in each mode.
We apply the above three approaches to the reconstruction of partially observed tensors. In both synthetic and realworld datasets, we show the superior predictive performance of the proposed approaches against conventional expectation maximization (EM) based estimation of Tucker decomposition model. We also demonstrate the effectiveness of a heuristic to improve the interpretability of the core tensor obtained by the proposed approaches on the amino acid fluorescence dataset.
This paper is structured as follows. In the next section, we first review the matrix rank and its relation to the trace norm. Then we review the definition of tensor mode rank, which suggests that a low rank tensor is a low rank matrix when appropriately unfolded. In Section 3, we propose three approaches to extend the tracenorm regularization for the estimation of lowrank tensors. In Section 4, we show that the optimization problems associated to the proposed extensions can be solved efficiently by the alternating direction method of multipliers [17]. In Section 5, we show through numerical experiments that one of the proposed approaches can recover a partly observed lowrank tensor almost perfectly from smaller fraction of observations compared to the conventional EMbased Tucker decomposition algorithm. The proposed algorithm shows a sharp threshold behaviour from a poor fit to a nearly perfect fit; we numerically show that the fraction of samples at the threshold is roughly proportional to the sum of the ranks of the underlying tensor when the tensor dimension is fixed. Finally we summarize the paper in Section 6
. Earlier version of this manuscript appeared in NIPS2010 workshop “Tensors, Kernels, and Machine Learning”.
2 Low rank matrix and tensor
In this section, we first discuss the connection between the rank of a matrix and the tracenorm regularization. Then we review the CP and the Tucker decomposition and the notions of tensor rank connected to them.
2.1 Rank of a matrix and the trace norm
The rank of an matrix
can be defined as the number of nonzero singular values of
. Here, the singularvalue decomposition (SVD) of
is written as follows:(1) 
where and are orthogonal matrices, and is the th largest singularvalue of . The matrix is called lowrank if the rank is less than . Unfortunately, the rank of a matrix is a nonconvex function, and the direct minimization of rank or solving a rankconstrained problem is an NPhard problem [32].
The trace norm is known to be the tightest convex lower bound of matrix rank [32] (see Fig. 1) and is defined as the linear sum of singular values as follows:
Intuitively, the trace norm plays the role of the norm in the subset selection problem [39], for the estimation of lowrank matrix^{1}^{1}1Note however that the absolute value is not taken here because singular value is defined to be positive.. The convexity of the above function follows from the fact that it is the dual norm of the spectral norm (see [6, Section A.1.6]). Since it is a norm, the trace norm is a convex function. The nondifferentiability of the trace norm at the origin promotes many singular values of to be zero when used as a regularization term. In fact, the following minimization problem has an analytic solution known as the spectral softthresholding operator (see [8]):
(2) 
where is the Frobenius norm, and is a regularization constant. The spectral softthresholding operation can be considered as a shrinkage operation on the singular values and is defined as follows:
(3) 
where is the SVD of the input matrix , and the operation is taken elementwise. We can see that the spectral softthresholding operation truncates the singularvalues of the input matrix smaller than to zero, thus the resulting matrix is usually lowrank. See also [42] for the derivation.
For the recovery of partially observed lowrank matrix, some theoretical guarantees have recently been developed. Candès and Recht [10] showed that in the noiseless case, samples are enough to perfectly recover the matrix under uniform sampling if the rank is not too large, where .
2.2 Rank of a tensor
For higher order tensors, there are several definitions of rank. Let be a way tensor. The rank of a tensor (see [24]) is defined as the minimum number of components required for rankone decomposition of a given tensor in analogy to SVD as follows:
(4) 
where denotes the outer product, denotes a way diagonal matrix whose th element is , and denotes the mode matrix product (see Kolda & Bader [23]); in addition, we define . The above decomposition model is called CANDECOMP [11] or PARAFAC [19]. It is worth noticing that finding the above decomposition with the minimum is a hard problem; thus there is no straightforward algorithm for computing the rank for higherorder tensors [24].
We consider instead the mode rank of tensors, which is the foundation of the Tucker decomposition [44, 13]. The mode rank of , denoted , is the dimensionality of the space spanned by the mode fibers of . In other words, the mode rank of is the rank of the mode unfolding of . The mode unfolding is the () matrix obtained by concatenating the mode fibers of
as column vectors. In MATLAB this can be obtained as follows:
X=permute(X,[k:K,1:k1]); X=X(:,:);
where the order of dimensions other than the first dimension is not important as long as we use a consistent definition. We say that a way tensor is rank if the mode rank of is (). Unlike the rank of the tensor, modek rank is clearly computable; the computation of the mode ranks of a tensor boils down to the computation the rank of matrices.
A rank tensor can be written as
(5) 
where is called a core tensor, and () are left singularvectors from the SVD of the mode unfolding of . The above decomposition is called the Tucker decomposition [44, 23].
The definition of a lowrank tensor (in the sense of Tucker decomposition) implies that a lowrank tensor is a lowrank matrix when unfolded appropriately. In order to see this, we recall that for the Tucker model (5), the mode unfolding of can be written as follows (see e.g., [23]):
Therefore, if the tensor is lowrank in the th mode (i.e., ), its unfolding is a lowrank matrix. Conversely, if is a lowrank matrix (i.e., ), we can set , , and other Tucker factors as identity matrices, and we obtain a Tucker decomposition (5).
Note that if a given tensor can be written in the form of CP decomposition (4) with rank , we can always find a rank Tucker decomposition (5) by orthonormalizing each factor in Equation (4). Therefore, the Tucker decomposition is more general than the CP decomposition.
However, since the core tensor that corresponds to singularvalues in the matrix case (see Equation (1)) is not diagonal in general, it is not straightforward to generalize the trace norm from matrices to tensors.
3 Three strategies to extend the tracenorm regularization to tensors
In this section, we first consider a given tensor as a matrix by unfolding it at a given mode and propose to minimize the trace norm of the unfolding . Next, we extend this to the minimization of the weighted sum of the trace norms of the unfoldings. Finally, relaxing the condition that the tensor is jointly lowrank in every mode in the second approach, we propose a mixture approach. For solving the optimization problems, we use the alternating direction method of multipliers (ADMM) [17] (also known as the split Bregman iteration [18]). The optimization algorithms are discussed in Section 4.
3.1 Tensor as a matrix
If we assume that the tensor we wish to estimate is (at least) lowrank in the th mode, we can convert the tensor estimation problem into a matrix estimation problem. Extending the minimization problem (2) to accommodate missing entries we have the following optimization problem for the reconstruction of partially observed tensor:
(6) 
where is the mode unfolding of , is the vector of observations, and is a linear operator that reshapes the prespecified elements of the input tensor into an dimensional vector; is the number of observations. In Equation (6), the regularization constant is moved to the denominator of the loss term from the numerator of the regularization term in Equation (2); this equivalent reformulation allows us to consider the noiseless case in the same framework. Note that
can also be interpreted as the variance of the Gaussian observation noise model.
Since the estimation procedure (6) is essentially an estimation of a lowrank matrix , we know that in the noiseless case samples are enough to perfectly recover the unknown true tensor , where and , if the rank is not too high [10]. This holds regardless of whether the unknown tensor is lowrank in other modes . Therefore, when we can estimate the mode unfolding of perfectly, we can also recover the whole perfectly, including the ranks of the modes we did not use during the estimation.
However, the success of the above procedure is conditioned on the choice of the mode to unfold the tensor. If we choose a mode with a large rank, even if there are other modes with smaller ranks, we cannot hope to recover the tensor from a small number of samples.
3.2 Constrained optimization of low rank tensors
In order to exploit the rank deficiency of more than one mode, it is natural to consider the following extension of the estimation procedure (6)
(7) 
This is a convex optimization problem, because it can be reformulated as follows:
(8)  
subject to  (9) 
where is the vectorization of (), is the matrix representation of mode unfolding (note that is a permutation matrix; thus ), is an auxiliary matrix of the same size as the mode unfolding of , and is the vectorization of . With a slight abuse of notation denotes the observation operator as a matrix.
This approach was considered earlier by Liu et al. [27] and Signoretto et al. [36]. Liu et al. relaxed the constraints (9) into penalty terms, therefore the factors obtained as the left singular vectors of does not equal the factors of the Tucker decomposition of . Signoretto et al. have discussed the general Shatten norms for tensors and the relationship between the regularization term in Equation (7) with (which corresponds to Shatten norm) and the function .
3.3 Mixture of lowrank tensors
The optimization problem (8) penalizes every mode of the tensor to be jointly lowrank, which might be too strict to be satisfied in practice. Thus we propose to predict instead with a mixture of tensors; each mixture component is regularized by the trace norm to be lowrank in each mode. More specifically, we solve the following minimization problem:
(10) 
Note that when for all , the problem (10) reduces to the problem (8) with .
3.4 Interpretation
All three proposed approaches inherit the lack of uniqueness of the factors from the conventional Tucker decomposition [23]. Some heuristics to improve the interpretability of the core tensor are proposed and implemented in the way toolbox [2]. However, these approaches are all restricted to orthogonal transformations. Here we present another simple heuristic, which is to apply PARAFAC decomposition on the core tensor . This approach has the following advantages over applying PARAFAC directly to the original data. First, the dimensionality of the core tensor is automatically obtained from the proposed algorithms. Therefore, the range of the number of PARAFAC components that we need to look for is much narrower than applying PARAFAC directly to the original data. Second, the PARAFAC problem does not need to take care of missing entries. In other words, we can separate the prediction problem and the interpretation problem, which are separately tackled by the proposed algorithms and PARAFAC, respectively. Finally, empirically the proposed heuristic seems to be more robust in recovering the underlying factors compared to applying PARAFAC directory when the rank is misspecified (see Section 5.2).
More precisely, let us consider the second “Constraint” approach. Let be the left singular vectors of the auxiliary variables . From Equation (5) we can obtain the core tensor as follows:
Let be the factors obtained by the PARAFAC decomposition of as follows:
Therefore, we have the following decomposition
(11) 
which gives the th factor as .
4 Optimization
In this section, we describe the optimization algorithms based on the alternating direction method of multipliers (ADMM) for the problems (6), (8), and (10).
4.1 Admm
The alternating direction method of multipliers [17] (see also [5]) can be considered as an approximation of the method of multipliers [31, 21] (see also [4, 30]). The method of multipliers generates a sequence of primal variables and multipliers by iteratively minimizing the so called augmented Lagrangian (AL) function with respect to the primal variables and updating the multiplier vector . Let us consider the following linear equality constrained minimization problem:
(12)  
subject to  (13) 
where and are both convex functions. The AL function of the above minimization problem is written as follows:
where is the Lagrangian multiplier vector. Note that when , the AL function reduces to the ordinary Lagrangian function. Intuitively, the additional penalty term enforces the equality constraint to be satisfied. However, different from the penalty method (which was used in [27]), there is no need to increase the penalty parameter very large, which usually makes the problem poorly conditioned.
The original method of multipliers performs minimization of the AL function with respect to and jointly followed by a multiplier update as follows:
(14)  
(15) 
Intuitively speaking, the multiplier is updated proportionally to the violation of the equality constraint (13). In this sense, can also be regarded as a stepsize parameter. Under fairly mild conditions, the above method converges superlinearly to a solution of the minimization problem (12); see [34, 41]. However, the joint minimization of the AL function (14) is often hard (see [41] for an exception).
The ADMM decouples the minimization with respect to and as follows:
(16)  
(17)  
(18) 
Note that the new value of obtained in the first line is used in the update of in the second line. The multiplier update step is identical to that of the ordinary method of multipliers (15). It can be shown that the above algorithm is an application of firmly nonexpansive mapping and that it converges to a solution of the original problem (12). Surprisingly, this is true for any positive penalty parameter [26]. This is in contrast to the fact that a related approach called forwardbackward splitting [26] (which was used in [36]) converges only when the stepsize parameter is chosen appropriately.
4.2 Stopping criterion
As a stopping criterion for terminating the above ADMM algorithm, we employ the relative duality gap criterion; that is, we stop the algorithm when the current primal objective value and the largest dual objective value obtained in the past satisfies the following equality
(19) 
Note that the multiplier vector computed in Equation (18) cannot be directly used in the computation of the dual objective value, because typically violates the dual constraints. See Appendix A for the details.
The reason we use the duality gap is that the criterion is invariant to the scale of the observed entries and the size of the problem .
4.3 ADMM for the “As a Matrix” approach
We consider the following constrained reformulation of problem (6)
(20) 
where is a vectorization of , is an auxiliary variable that corresponds to the mode unfolding of , and is the vectorization of . The AL function of the above constrained minimization problem can be written as follows:
(21) 
where is the Lagrangian multiplier vector that corresponds to the constraint . Note that we rescaled the Lagrangian multiplier vector by the factor for the sake of notational simplicity.
Starting from an initial point , we apply the ADMM explained in the previous section to the AL function (21). All the steps (16)–(18) can be implemented in closed forms. First, minimization with respect to yields,
(22) 
where is an dimensional vector that has one for observed elements and zero otherwise; is an dimensional vector filled with ones; denotes elementwise division. Note that when (no observational noise), the above expression can be simplified as follows:
(23) 
Here the observed entries of are overwritten by the observed values and the unobserved entries are filled with the mode tensorization of the current prediction . In the general case (22), the predicted values also affect the observed entries. The primal variable and the auxiliary variable becomes closer and closer as the optimization proceeds. This means that eventually the multiplier vector takes nonzero values only on the observed entries when .
Next, the minimization with respect to yields,
where is the spectral softthreshold operation (3) in which the argument is considered as a matrix.
The last step is the multiplier update (18), which can be written as follows:
(24) 
Note that the stepsize parameter does not appear in (24) due to the rescaling of in (21).
The speed of convergence of the algorithm mildly depends on the choice of the stepsize . Here as a guideline to choose , we require that the algorithm is invariant to scalar multiplication of the objective (20). More precisely, when the input and the regularization constant are both multiplied by a constant , the solution of the minimization (20) (or (6)) should remain essentially the same as the original problem, except that the solution is also multiplied by the constant . In order to make the algorithm (see (22)(24)) follow the same path (except that , , and are all multiplied by ), we need to scale inversely proportional to . We can also see this in the AL function (21); in fact, the first two terms scale linearly to , and also the last two terms scale linearly if scales inversely to . Therefore we choose as , where is a constant and
is the standard deviation of the observed values
.4.4 ADMM for the “Constraint” approach
The AL function of the constrained minimization problem (8)(9) can be written as follows:
Note that we rescaled the multiplier vector by the factor as in the previous subsection.
Starting from an initial point , we take similar steps as in (22)(24) except that the last two steps are performed for all . That is,
(25)  
(26)  
(27) 
By considering the scale invariance of the algorithm, we choose the stepsize as as in the previous subsection.
4.5 ADMM for the “Mixture” approach
We consider the following dual problem of the mixture formulation (10):
(28)  
subject to 
where is a dual vector; is an auxiliary variable that corresponds to the mode unfolding of , and is the vectorization of ; the indicator function is defined as , if , and , otherwise, where is the spectral norm (maximum singularvalue of a matrix).
The AL function for the problem (28) can be written as follows:
Similar to the previous two algorithms, we start from an initial point , and compute the following steps:
(29) 
The above steps can be computed in closed forms. In fact,
(30)  
(31)  
where the projection operator is the projection onto a radius spectralnorm ball, as follows:  
where is the SVD of the matricization of the input vector . Moreover, combining the two steps (31) and (29), we have (see [41])
(32) 
Note that we recover the spectral softthreshold operation by combining the two steps. Therefore, we can simply iterate steps (30) and (32) (note that the term in (30) can be computed from (29) as .)
In order to see that the multiplier vector obtained in the above steps converges to the primal solution of the mixture formulation (10), we take the derivative of the ordinary Lagrangian function with respect to and () and obtain the following optimality conditions:
where we used the relationship , and the fact that implies because the two functions and are conjugate to each other; see [33, Cor. 23.5.1]. By combining the above two equations, we obtain the optimality condition for the mixture formulation (10) as follows:
As in the previous two subsections, we require that the algorithm (29)(32) is invariant to scalar multiplication of the input and the regularization constant by the same constant . Since appears in the final solution, must scale linearly with respect to . Thus from (29), if and are constants with respect to , the stepsize must scale linearly. In fact, from (30) and (31), we can see that these two dual variables remain constant when , , and are multiplied by . Therefore, we choose .
5 Numerical experiments
In this section, we first present results on two synthetic datasets. Finally we apply the proposed methods to the Amino acid fluorescence data published by Bro and Andersson [7].
5.1 Synthetic experiments
We randomly generated a rank(7,8,9) tensor of dimensions (50,50,20) by drawing the core from the standard normal distribution and multiplying its each mode by an orthonormal factor randomly drawn from the Haar measure. We randomly selected some elements of the true tensor for training and kept the remaining elements for testing. We used the algorithms described in the previous section with the tolerance
. We choose for simplicity in the later two approaches. The stepsize is chosen as for the first two approaches and for the third approach with . For the first two approaches, (zero observation error) was used; see (23). For the last approach, we used . The Tucker decomposition algorithm tucker from the way toolbox [2] is also included as a baseline, for which we used the correct rank (“exact”) and the 20% higher rank (“large”). Note that all proposed approaches can find the rank automatically. The generalization error is defined as follows:where is the vectorization of the unobserved entries and is the prediction computed by the algorithms. For the “As a Matrix” strategy, error for each mode is reported. All algorithms were implemented in MATLAB and ran on a computer with two 3.5GHz Xeon processors and 32GB of RAM. The experiment was repeated 20 times and averaged.
Figure 2 shows the result of tensor completion using three strategies we proposed above, as well as the Tucker decomposition. At 35% observation, the proposed “Constraint” obtains nearly perfect generalization. Interestingly there is a sharp transition from a poor fit (generalization error) to an almost perfect fit (generalization error). The “As a Matrix” approach also show similar transition for mode 1 and mode 2 (around 40%), and mode 3 (around 80%), but even the first transition is slower than the “Constraint” approach. The “Mixture” approach shows a transition around 70% slightly faster than the mode 3 in the “As A Matrix” approach. Tucker shows early decrease in the generalization error, but when the rank is misspecified (“large”), the error remains almost constant; even when the correct rank is known (“exact”), the convergence is slower than the proposed “Constraint” approach.
The proposed convex approaches are not only accurate but also fast. Fig. 3 shows the computation time of the proposed approaches and EMbased Tucker decomposition against the fraction of observed entries. For the “As a Matrix” approach the total time for all modes is plotted. We can see that the “As a Matrix” and “Constraint” approaches are roughly 4–10 times faster than the conventional EMbased tucker decomposition.
We have further investigated the condition for the threshold behaviour using the proposed “Constraint” approach. Here we generated different problems of different core dimensions . The sum of mode ranks is defined as . For each problem, we apply the “Constraint” approach for increasingly large fraction of observations and determine when the generalization error falls below . Fig. 4 shows the fraction of observations required to obtain generalization error below (in other words, the fraction at the threshold) against the sum of mode ranks defined above. We can see that the fraction at the threshold is roughly proportional to the sum of the mode ranks of the underlying tensor. We do not have any theoretical argument to support this observation. Acar et al [1] also empirically discussed condition for successful recovery for the CP decomposition.
Figure 5 show another synthetic experiment. We randomly generated a rank tensor of the same dimensions as above. We chose the same parameter values , , , and . Here we can see that interestingly the “Constraint” approach perform poorly, whereas the “mode 3” and “Mixture” perform clearly better than other algorithms. It is natural that the “mode 3” approach works well because the true tensor is only lowrank in the third mode. In contrast, the “Mixture” approach can automatically detect the rankdeficient mode, because the regularization term in the formulation (10) is a linear sum of three () penalty terms. The linear sum structure enforces sparsity across . Therefore, in this case and were switched off, and “Mixture” approach yielded almost identical results to the “mode 3” approach.
5.2 Amino acid fluorescence data
The amino acid fluorescence data is a semirealistic data contributed by Bro and Andersson [7], in which they measured the fluorescence of five laboratorymade solutions that each contain different amounts of tyrosine, tryptophan and phenylalanine. Since the “factors” are known to be the three amino acids, this is a perfect data for testing whether the proposed method can automatically find those factors.
For the experiments in this subsection, we chose the same parameter setting , , , and . Setting corresponds to assuming no observational noise. This can be justified by the fact that the original data is already approximately lowrank (rank) in the sense of Tucker decomposition. The dimensionality of the original tensor is , which correspond to emission wavelength (250–450 nm), excitation wavelength (240–300 nm), and samples, respectively.
Fig. 6 show the generalization error obtained by the proposed approaches as well as EMbased Tucker and PARAFAC decompositions. Here PARAFAC is included because the dataset is originally designed for PARAFAC. We can see that the proposed “Constraint” approach show fast decrease in generalization error, which is comparable to the PARAFAC model knowing the correct dimension. Tucker decomposition of rank performs as good as PARAFAC models when more than half the entries are observed. However, a slightly larger rank Tucker decomposition could not decrease the error below .
Fig. 7 show the factors obtained by fitting directly threecomponent PARAFAC model, fourcomponent PARAFAC model, and applying a four component PARAFAC model to the core obtained by the proposed “Constraint” approach. The fraction of observed entries was . The two conventional approaches used EM iteration for the estimation of missing values. For the proposed model, the dimensionality of the core was ; this was obtained by keeping the singularvalues of the auxiliary variable that are larger than 1% of its largest singularvalue for each . Then we applied a fourcomponent (fullyobserved) PARAFAC model to this core and obtained the factors as in Equation (11). Interestingly, although the four componentPARAFAC model is redundant for this problem [7], the proposed approach seem to be more robust than applying fourcomponent PARAFAC directly to the data. We can see that the shape of the major three components (blue, green, red) obtained by the proposed approach (the right column) are more similar to the threecomponent PARAFAC model (the left column) than the fourcomponent PARAFAC model (the center column).
6 Summary
In this paper we have proposed three strategies to extend the framework of trace norm regularization to the estimation of partially observed lowrank tensors. The proposed approaches are formulated in convex optimization problems and the rank of the tensor decomposition is automatically determined through the optimization.
In the simulated experiment, tensor completion using the “Constraint” approach showed nearly perfect reconstruction from only 35% observations. The proposed approach shows a sharp threshold behaviour and we have empirically found that the fraction of samples at the threshold is roughly proportional to the sum of mode ranks of the underlying tensor.
We have also shown the weakness of the “Constraint” approach. When the unknown tensor is only lowrank in certain mode, the assumption that the tensor is lowrank in every mode, which underlies the “Constraint” approach, is too strong. We have demonstrated that the “Mixture” approach is more effective in this case. The “Mixture” approach can automatically detect the rankdeficient mode and lead to better performance.
In the amino acid fluorescence dataset, we have shown that the proposed “Constraint” approach outperforms conventional EMbased Tucker decomposition and is comparable to PARAFAC model with the correct number of components. Moreover, we have demonstrated a simple heuristic to obtain a PARAFACstyle decomposition from the decomposition obtained by the proposed method. Moreover, we have shown that the proposed heuristic can reliably recover the true factors even when the number of PARAFAC factors is misspecified.
The proposed approaches can be extended in many ways. For example, it would be important to handle nonGaussian noise model [12, 20]; for example a tensor version of robust PCA [9] would be highly desirable. For classification over tensors, extension of the approach in [40] would be meaningful in applications including braincomputer interface; see also [35] for another recent approach. It is also important to extend the proposed approach to handle large scales tensors that cannot be kept in the RAM. Combination of the firstorder optimization proposed by Acar et al. [1] with our approach is a promising direction. Moreover, in order to understand the threshold behaviour, further theoretical analysis is necessary.
Appendix A Computation of the dual objectives
In this Appendix, we show how we compute the dual objective values for the computation of the relative duality gap (19).
a.1 Computation of dual objective for the “As a Matrix” approach
The dual problem of the constrained minimization problem (20) can be written as follows:
(33)  
subject to  (34) 
Here is the linear operator that reshapes the elements of a given dimensional vector that correspond to the unobserved entries into an dimensional vector. In addition, is the matricization of and is the spectral norm (maximum singular value).
Note that the multiplier vector obtained through ADMM does not satisfy the above two constraints (34). Therefore, similar to the approach used in [45, 41], we apply the following transformations. First, we compute the projection by projecting to the equality constraint. This can be done easily by setting the elements of that correspond to unobserved entries to zero. Second, we compute the maximum singular value of the matricization of and shrink as follows:
Clearly this operation does not violate with the equality constraint. Finally we substitute into the dual objective (33) to compute the relative duality gap as in Equation (19).
a.2 Computation of dual objective for the “Constraint” approach
The dual problem of the constrained minimization problem (8) can be written as follows:
(35)  
subject to 
Here the antiobservation operator is defined as in the last subsection, and is the matricization of ().
In order to obtain a dual feasible point from the current multiplier vectors (), we apply similar transformations as in the last subsection. First, we compute the projection to the equality constraint. This can be done by computing the sum over for each unobserved entry. Then the sum divided by is subtracted from each corresponding entry for . Let us denote by the multiplier vectors after the projection. Next, we compute the largest singularvalues where is the matricization of the projected multiplier vector for . Now in order to enforce the inequality constraints, we define the shrinkage factor as follows:
(36) 
Using the above shrinkage factor, we obtain a dual feasible point as follows:
Finally, we substitute into the dual objective (35) to compute the relative duality gap as in Equation (19).
a.3 Computation of dual objective for the “Mixture” approach
The dual problem of the mixture formulation is already given in Equation (28). Making the implicit inequality constraints explicit, we can rewrite this as follows:
subject to 
Note that the norm in the second line should be interpreted as the spectral norm of the matricization of . Although the ADMM presented in Section 4.5 was designed to solve this dual formulation, we did not discuss how to evaluate the dual objective. Again the dual vector obtained through the ADMM does not satisfy the inequality constraints.
In order to obtain a dual feasible point, we compute the largest singularvalues for . From the singularvalues , we can compute the shrinkage factor as in Equation (36) in the previous subsection. Finally, a dual feasible point can be obtained as , which we use for the computation of the relative duality gap (19).
References
 [1] E. Acar, D.M. Dunlavy, T.G. Kolda, and M. Mørup, Scalable tensor factorizations with missing data, tech. report, arXiv:1005.2197v1 [math.NA], 2010.
 [2] C. A. Andersson and R. Bro, The way toolbox for MATLAB, Chemometrics & Intelligent Laboratory Systems, 52 (2000), pp. 1–4. http://www.models.life.ku.dk/source/nwaytoolbox/.
 [3] A. Argyriou, T. Evgeniou, and M. Pontil, Multitask feature learning, in Advances in Neural Information Processing Systems 19, B. Schölkopf, J. Platt, and T. Hoffman, eds., MIT Press, Cambridge, MA, 2007, pp. 41–48.
 [4] D. P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods, Academic Press, 1982.
 [5] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers, 2011. Unfinished working draft.
 [6] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.
 [7] R. Bro, PARAFAC. Tutorial and applications, Chemometrics and Intelligent Laboratory Systems, 38 (1997), pp. 149–171.
 [8] J.F. Cai, E. J. Candes, and Z. Shen, A singular value thresholding algorithm for matrix completion, tech. report, arXiv:0810.3286, 2008.
 [9] E. J. Candes, X. Li, Y. Ma, and J. Wright, Robust principal component analysis?, tech. report, arXiv:0912.3599, 2009.
 [10] E. J. Candes and B. Recht, Exact matrix completion via convex optimization, Foundations of Computational Mathematics, 9 (2009), pp. 717–772.
 [11] J.D. Carroll and J.J. Chang, Analysis of individual differences in multidimensional scaling via an nway generalization of “EckartYoung” decomposition, Psychometrika, 35 (1970), pp. 283–319.
 [12] E. C. Chi and T. G. Kolda, Making tensor factorizations robust to nongaussian noise, tech. report, arXiv: 1010.3043v1, 2010.
 [13] L. De Lathauwer, B. De Moor, and J. Vandewalle, On the best rank1 and rank() approximation of higherorder tensors, SIAM Journal on Matrix Analysis and Applications, 21 (2000), pp. 1324–1342.
 [14] L. De Lathauwer and J. Vandewalle, Dimensionality reduction in higherorder signal processing and rank() reduction in multilinear algebra, Linear Algebra and its Applications, 391 (2004), pp. 31–55.
 [15] L. Eldén and B. Savas, A Newton–Grassmann method for computing the best multilinear rank() approximation of a tensor,, SIAM J. Matrix Anal. Appl., 31 (2009), pp. 248–271.
 [16] M. Fazel, H. Hindi, and S. P. Boyd, A Rank Minimization Heuristic with Application to Minimum Order System Approximation, in Proc. of the American Control Conference, 2001.
 [17] D. Gabay and B. Mercier, A dual algorithm for the solution of nonlinear variational problems via finite element approximation, Comput. Math. Appl., 2 (1976), pp. 17–40.
 [18] T. Goldstein and S. Osher, The split Bregman method for L1 regularized problems, SIAM Journal on Imaging Sciences, 2 (2009), pp. 323–343.
 [19] R.A. Harshman, Foundations of the PARAFAC procedure: Models and conditions for an “explanatory” multimodal factor analysis, UCLA working papers in phonetics, 16 (1970), pp. 1–84.

[20]
K. Hayashi, T. Takenouchi, T. Shibata, Y. Kamiya, D. Kato, K. Kunieda,
K. Yamada, and K. Ikeda,
Exponential family tensor factorization for missingvalues prediction and anomaly detection
, in 2010 IEEE International Conference on Data Mining, 2010, pp. 216–225.  [21] M. R. Hestenes, Multiplier and gradient methods, J. Optim. Theory Appl., 4 (1969), pp. 303–320.
 [22] S. Ji and J. Ye, An accelerated gradient method for trace norm minimization, in Proceedings of the 26th International Conference on Machine Learning (ICML2009), New York, NY, 2009, ACM, pp. 457–464.
 [23] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM Review, 51 (2009), pp. 455–500.
 [24] J. B. Kruskal, Rank, decomposition, and uniqueness for 3way and nway arrays, in Multiway data analysis, R. Coppi and S. Bolasco, eds., Elsevier, NorthHolland, Amsterdam, 1989, pp. 7–18.
 [25] Z. Lin, M. Chen, L. Wu, and Y. Ma, The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted LowRank Matrices, Mathematical Programming, (2009). submitted.
 [26] P. L. Lions and B. Mercier, Splitting algorithms for the sum of two nonlinear operators, SIAM Journal on Numerical Analysis, 16 (1979), pp. 964–979.
 [27] J. Liu, P. Musialski, P. Wonka, and J. Ye, Tensor completion for estimating missing values in visual data, in Prof. ICCV, 2009.
 [28] M. Mørup, Applications of tensor (multiway array) factorizations and decompositions in data mining, Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 1 (2011), pp. 24–40.
 [29] M. Mørup, L.K. Hansen, C.S. Herrmann, J. Parnas, and S.M. Arnfred, Parallel factor analysis as an exploratory tool for wavelet transformed eventrelated EEG, NeuroImage, 29 (2006), pp. 938–947.
 [30] J. Nocedal and S. Wright, Numerical Optimization, Springer, 1999.
 [31] M. J. D. Powell, A method for nonlinear constraints in minimization problems, in Optimization, R. Fletcher, ed., Academic Press, London, New York, 1969, pp. 283–298.
 [32] B. Recht, M. Fazel, and P.A. Parrilo, Guaranteed minimumrank solutions of linear matrix equations via nuclear norm minimization, SIAM Review, 52 (2010), pp. 471–501.
 [33] R. T. Rockafellar, Convex Analysis, Princeton University Press, 1970.
 [34] R. T. Rockafellar, Augmented Lagrangians and applications of the proximal point algorithm in convex programming, Math. of Oper. Res., 1 (1976), pp. 97–116.
 [35] M. Signoretto, L. De Lathauwer, and J.A.K. Suykens, Convex multilinear estimation and operatorial representations, in NIPS2010 Workshop: Tensors, Kernels and Machine Learning (TKML), 2010.
 [36] , Nuclear norms for tensors and their use for convex multilinear estimation, Tech. Report 10186, ESATSISTA, K.U.Leuven, 2010.
 [37] A.K. Smilde, R. Bro, and P. Geladi, Multiway analysis with applications in the chemical sciences, Wiley, 2004.
 [38] N. Srebro, J. D. M. Rennie, and T. S. Jaakkola, Maximummargin matrix factorization, in Advances in Neural Information Processing Systems 17, Lawrence K. Saul, Yair Weiss, and Léon Bottou, eds., MIT Press, Cambridge, MA, 2005, pp. 1329–1336.
 [39] R. Tibshirani, Regression shrinkage and selection via the lasso, J. Roy. Stat. Soc. B, 58 (1996), pp. 267–288.
 [40] R. Tomioka and K. Aihara, Classifying matrices with a spectral regularization, in Proceedings of the 24th International Conference on Machine Learning (ICML2007), ACM Press, 2007, pp. 895–902.
 [41] R. Tomioka, T. Suzuki, and M. Sugiyama, Superlinear convergence of dual augmentedLagrangian algorithm for sparsity regularized estimation, tech. report, arXiv:0911.4046v2, 2009.
 [42] R. Tomioka, T. Suzuki, and M. Sugiyama, Augmented Lagrangian methods for learning, selecting, and combining features, in Optimization for Machine Learning, Suvrit Sra, Sebastian Nowozin, and Stephen J. Wright, eds., MIT Press, 2011.
 [43] R. Tomioka, T. Suzuki, M. Sugiyama, and H. Kashima, A fast augmented lagrangian algorithm for learning lowrank matrices, in Proceedings of the 27 th Annual International Conference on Machine Learning (ICML2010), Johannes Fürnkranz and Thorsten Joachims, eds., Omnipress, 2010.
 [44] L. R. Tucker, Some mathematical notes on threemode factor analysis, Psychometrika, 31 (1966), pp. 279–311.
 [45] S. J. Wright, R. D. Nowak, and M. A. T. Figueiredo, Sparse reconstruction by separable approximation, IEEE Trans. Signal Process., 57 (2009), pp. 2479–2493.