# Estimation of low-rank tensors via convex optimization

In this paper, we propose three approaches for the estimation of the Tucker decomposition of multi-way arrays (tensors) from partial observations. All approaches are formulated as convex minimization problems. Therefore, the minimum is guaranteed to be unique. The proposed approaches can automatically estimate the number of factors (rank) through the optimization. Thus, there is no need to specify the rank beforehand. The key technique we employ is the trace norm regularization, which is a popular approach for the estimation of low-rank matrices. In addition, we propose a simple heuristic to improve the interpretability of the obtained factorization. The advantages and disadvantages of three proposed approaches are demonstrated through numerical experiments on both synthetic and real world datasets. We show that the proposed convex optimization based approaches are more accurate in predictive performance, faster, and more reliable in recovering a known multilinear structure than conventional approaches.

## Authors

• 26 publications
• 15 publications
• 12 publications
• ### Reexamining Low Rank Matrix Factorization for Trace Norm Regularization

Trace norm regularization is a widely used approach for learning low ran...

06/27/2017 ∙ by Carlo Ciliberto, et al. ∙ 0

read it

• ### Convex Coupled Matrix and Tensor Completion

We propose a set of convex low rank inducing norms for a coupled matrice...

05/15/2017 ∙ by Kishan Wimalawarne, et al. ∙ 0

read it

• ### Towards Faster Rates and Oracle Property for Low-Rank Matrix Estimation

We present a unified framework for low-rank matrix estimation with nonco...

05/18/2015 ∙ by Huan Gui, et al. ∙ 0

read it

• ### Convex recovery of tensors using nuclear norm penalization

The subdifferential of convex functions of the singular spectrum of real...

06/08/2015 ∙ by Stéphane Chrétien, et al. ∙ 0

read it

• ### On the Convergence of Projected-Gradient Methods with Low-Rank Projections for Smooth Convex Minimization over Trace-Norm Balls and Related Problems

Smooth convex minimization over the unit trace-norm ball is an important...

02/05/2019 ∙ by Dan Garber, et al. ∙ 0

read it

• ### Low rank methods for multiple network alignment

Multiple network alignment is the problem of identifying similar and rel...

09/21/2018 ∙ by Huda Nassar, et al. ∙ 0

read it

• ### Low-Rank Factorization of Determinantal Point Processes for Recommendation

Determinantal point processes (DPPs) have garnered attention as an elega...

02/17/2016 ∙ by Mike Gartrell, et al. ∙ 0

read it

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Multi-way 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 multi-way arrays (or tensors) into small number of factors have been one of the main concerns in multi-way data analysis, because interpreting the original multi-way 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 non-convex 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 low-rank 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], multi-task learning [3], and classification over matrices [40]. In addition, there are theoretical developments that (under some conditions) guarantee perfect reconstruction of a low-rank matrix from partial measurements via convex estimation [10, 32]. The key idea here is to replace the rank of a matrix (a non-convex function) by the so-called trace norm (also known as the nuclear norm) of the matrix. One goal of this paper is to extend the trace-norm 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 low-rank 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 low-rank 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 low-rank 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 low-rank in each mode.

We apply the above three approaches to the reconstruction of partially observed tensors. In both synthetic and real-world 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 trace-norm regularization for the estimation of low-rank 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 low-rank tensor almost perfectly from smaller fraction of observations compared to the conventional EM-based 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 trace-norm 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 singular-value decomposition (SVD) of

is written as follows:

 X =Udiag(σ1(X),σ2(X),…,σr(X))V⊤, (1)

where and are orthogonal matrices, and is the th largest singular-value of . The matrix is called low-rank 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 rank-constrained problem is an NP-hard 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:

 ∥X∥∗ =r∑j=1σj(X).

Intuitively, the trace norm plays the role of the -norm in the subset selection problem [39], for the estimation of low-rank matrix111Note 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 non-differentiability 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 soft-thresholding operator (see [8]):

 argminX 12∥X−Y∥2Fro+λ∥X∥∗, (2)

where is the Frobenius norm, and is a regularization constant. The spectral soft-thresholding operation can be considered as a shrinkage operation on the singular values and is defined as follows:

 proxtrλ(Y)=Umax(S−λ,0)V⊤, (3)

where is the SVD of the input matrix , and the operation is taken element-wise. We can see that the spectral soft-thresholding operation truncates the singular-values of the input matrix smaller than to zero, thus the resulting matrix is usually low-rank. See also [42] for the derivation.

For the recovery of partially observed low-rank 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 rank-one decomposition of a given tensor in analogy to SVD as follows:

 X =r∑j=1λja(1)j∘a(2)j∘⋯∘a(K)j, =Λ×1A(1)×2A(2)⋯×KA(K) (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 higher-order 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:k-1]); 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, mode-k 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

 X =G×1U1×2U2⋯×KUK, (5)

where is called a core tensor, and () are left singular-vectors from the SVD of the mode- unfolding of . The above decomposition is called the Tucker decomposition [44, 23].

The definition of a low-rank tensor (in the sense of Tucker decomposition) implies that a low-rank tensor is a low-rank 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]):

 X(k) =UkG(k)(Uk−1⊗⋯⊗U1⊗UK⊗⋯⊗Uk+1)⊤.

Therefore, if the tensor is low-rank in the th mode (i.e., ), its unfolding is a low-rank matrix. Conversely, if is a low-rank 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 ortho-normalizing each factor in Equation (4). Therefore, the Tucker decomposition is more general than the CP decomposition.

However, since the core tensor that corresponds to singular-values 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 trace-norm 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 low-rank 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) low-rank 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:

 minimizeX∈Rn1×⋯×nK 12λ∥Ω(X)−y∥2+∥X(k)∥∗, (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 low-rank 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 low-rank 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.

Various advanced methods [42, 43, 25, 22] for the estimation of low-rank matrices can be used for solving the minimization problem (6). Here we use ADMM to keep the presentation concise; see Section 4 for the details.

### 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)

 minimizeX∈Rn1×⋯×nK 12λ∥Ω(X)−y∥2+K∑k=1γk∥X(k)∥∗. (7)

This is a convex optimization problem, because it can be reformulated as follows:

 minimizex,Z1,…,ZK 12λ∥Ωx−y∥2+K∑k=1γk∥Zk∥∗, (8) subject to Pkx=zk(k=1,…,K), (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 low-rank tensors

The optimization problem (8) penalizes every mode of the tensor to be jointly low-rank, 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 low-rank in each mode. More specifically, we solve the following minimization problem:

 minimizeZ1,…,ZK 12λ∥∥Ω(∑Kk=1Pk⊤zk)−y∥∥2+K∑k=1γk∥Zk∥∗. (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:

 G =X×1U1⊤×2U2⊤⋯×KUK⊤.

Let be the factors obtained by the PARAFAC decomposition of as follows:

 G =Λ×1A(1)×2A(2)⋯×KA(K).

Therefore, we have the following decomposition

 X =Λ×1(U1A(1))×2(U2A(2))⋯×K(UKA(K)), (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:

 minimizex∈Rn,z∈Rm f(x)+g(z), (12) subject to Ax=z, (13)

where and are both convex functions. The AL function of the above minimization problem is written as follows:

 Lη(x,z,α) =f(x)+g(z)+α⊤(Ax−z)+η2∥Ax−z∥2,

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:

 (xt+1,zt+1) =argminx∈Rn,z∈RmLη(x,z,αt), (14) αt+1 =αt+η(Axt+1−zt+1). (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 step-size parameter. Under fairly mild conditions, the above method converges super-linearly 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:

 xt+1 =argminx∈RnLη(x,zt,αt), (16) zt+1 =argminz∈RmLη(xt+1,z,αt), (17) αt+1 =αt+η(Axt+1−zt+1). (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 forward-backward splitting [26] (which was used in [36]) converges only when the step-size 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

 (p(xt,zt)−maxt′=1,…,td(~αt′))/p(xt,zt) <ϵ. (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)

 minimizex∈RN,Z∈Rnk×¯n∖k12λ∥Ωx−y∥2+∥Z∥∗,\rm subject toPkx=z, (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:

 Lη(x,Z,α) =12λ∥Ωx−y∥2+∥Z∥∗+ηα⊤(Pkx−z)+η2∥Pkx−z∥2, (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,

 xt+1 =(Ω⊤y+ληPk⊤(zt−αt))./(1Ω+λη1N), (22)

where is an -dimensional vector that has one for observed elements and zero otherwise; is an -dimensional vector filled with ones; denotes element-wise division. Note that when (no observational noise), the above expression can be simplified as follows:

 xt+1i ={(Ω⊤y)i,i∈Ω,(Pk⊤(zt−αt))i,i∉Ω(i=1,…,N). (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 non-zero values only on the observed entries when .

Next, the minimization with respect to yields,

 Zt+1=proxtr1/η(Pkxt+1+αt),

where is the spectral soft-threshold 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:

 αt+1 =αt+(Pkxt+1−zt+1). (24)

Note that the step-size 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 step-size . 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:

 Lη(x,{Zk}Kk=1,{αk}Kk=1) =12λ∥Ωx−y∥2+K∑k=1γk∥Zk∥∗ +K∑k=1(ηαk⊤(Pkx−zk)+η2∥Pkx−zk∥2).

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,

 xt+1 =(Ω⊤y+λη∑Kk=1Pk⊤(ztk−αtk))./(1Ω+ληK1N), (25) Zt+1k =proxtrγk/η(Pkxt+1+αtk)(k=1,…,K), (26) αt+1k =αtk+(Pkxt+1−zt+1k)(k=1,…,K). (27)

By considering the scale invariance of the algorithm, we choose the step-size as as in the previous subsection.

### 4.5 ADMM for the “Mixture” approach

We consider the following dual problem of the mixture formulation (10):

 minimizeα∈RM,Wk∈Rnk×¯n∖k λ2∥α∥2−α⊤y+K∑k=1δγk(Wk), (28) subject to wk=PkΩ⊤α(k=1,…,K),

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 singular-value of a matrix).

The AL function for the problem (28) can be written as follows:

 Lη(α,{Wk}Kk=1,{zk}Kk=1) =λ2∥α∥2−α⊤y+K∑k=1δγk(Wk) +K∑k=1(zk⊤(PkΩ⊤α−wk)+η2∥PkΩ⊤α−wk∥2)

Similar to the previous two algorithms, we start from an initial point , and compute the following steps:

 αt+1 =argminαLη(α,{Wtk}Kk=1,{ztk}Kk=1) Wt+1k =argminWkLη(αt+1,{Wk}Kk=1,{ztk}Kk=1) zt+1k =ztk+η(PkΩ⊤αt+1−wt+1k). (29)

The above steps can be computed in closed forms. In fact,

 αt+1 =1λ+ηK(y−Ω∑Kk=1Pk⊤(ztk−ηwtk)), (30) Wt+1k =projγk(PkΩ⊤αt+1+ztk/η), (31) where the projection operator projλ is the projection onto a radius λ-spectral-norm ball, as follows: projλ(w) :=Umin(S,λ)V⊤,

where is the SVD of the matricization of the input vector . Moreover, combining the two steps (31) and (29), we have (see [41])

 zt+1k =proxtrγkη(ztk+ηPkΩ⊤αt+1). (32)

Note that we recover the spectral soft-threshold 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:

 α =1λ(y−Ω∑Kk=1Pk⊤zk), PkΩ⊤α ∈∂γk∥Zk∥∗(k=1,…,K),

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:

 −1λPkΩ⊤(y−Ω∑Kk=1Pk⊤zk)+∂γk∥Zk∥∗∋0(k=1,…,K).

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 step-size 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 step-size 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:

 error =∥ypred−ytest∥∥ytest∥,

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 EM-based 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 EM-based 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 low-rank in the third mode. In contrast, the “Mixture” approach can automatically detect the rank-deficient 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 semi-realistic data contributed by Bro and Andersson [7], in which they measured the fluorescence of five laboratory-made 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 low-rank (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 EM-based 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 three-component PARAFAC model, four-component 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 singular-values of the auxiliary variable that are larger than 1% of its largest singular-value for each . Then we applied a four-component (fully-observed) PARAFAC model to this core and obtained the factors as in Equation (11). Interestingly, although the four component-PARAFAC model is redundant for this problem [7], the proposed approach seem to be more robust than applying four-component 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 three-component PARAFAC model (the left column) than the four-component 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 low-rank 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 low-rank in certain mode, the assumption that the tensor is low-rank 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 rank-deficient mode and lead to better performance.

In the amino acid fluorescence dataset, we have shown that the proposed “Constraint” approach outperforms conventional EM-based Tucker decomposition and is comparable to PARAFAC model with the correct number of components. Moreover, we have demonstrated a simple heuristic to obtain a PARAFAC-style 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 non-Gaussian 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 brain-computer 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 first-order 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:

 maximizeα∈RN −λ2∥∥ΩPk⊤α∥∥2+y⊤ΩPk⊤α (33) subject to ¯ΩPk⊤α=0,∥A∥≤1. (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:

 ~αt=min(1,1/σ1)^αt.

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:

 maximize{αk}Kk=1 −λ2∥∥∥ΩK∑k=1Pk⊤αk∥∥∥2+y⊤ΩK∑k=1Pk⊤αk, (35) subject to ¯ΩK∑k=1Pk⊤αk=0,∥Ak∥≤γk(k=1,…,K).

Here the anti-observation 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 singular-values 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:

 c=min(1,γ1/σ1,1,γ2/σ2,1,…,γK/σK,1). (36)

Using the above shrinkage factor, we obtain a dual feasible point as follows:

 ~αtk =c^αtk(k=1,…,K).

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:

 maximizeα∈RM −λ2∥α∥2+α⊤y, subject to ∥PkΩ⊤α∥≤γk(k=1,…,K).

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 singular-values for . From the singular-values , 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.
• [3] A. Argyriou, T. Evgeniou, and M. Pontil, Multi-task 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 n-way generalization of “Eckart-Young” decomposition, Psychometrika, 35 (1970), pp. 283–319.
• [12] E. C. Chi and T. G. Kolda, Making tensor factorizations robust to non-gaussian noise, tech. report, arXiv: 1010.3043v1, 2010.
• [13] L. De Lathauwer, B. De Moor, and J. Vandewalle, On the best rank-1 and rank-() approximation of higher-order tensors, SIAM Journal on Matrix Analysis and Applications, 21 (2000), pp. 1324–1342.
• [14] L. De Lathauwer and J. Vandewalle, Dimensionality reduction in higher-order 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” multi-modal 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 missing-values 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 3-way and n-way arrays, in Multiway data analysis, R. Coppi and S. Bolasco, eds., Elsevier, North-Holland, 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 Low-Rank 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 event-related 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 minimum-rank 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 10-186, ESAT-SISTA, K.U.Leuven, 2010.
• [37] A.K. Smilde, R. Bro, and P. Geladi, Multi-way analysis with applications in the chemical sciences, Wiley, 2004.
• [38] N. Srebro, J. D. M. Rennie, and T. S. Jaakkola, Maximum-margin 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, Super-linear convergence of dual augmented-Lagrangian 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 low-rank 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 three-mode 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.