Efficient Learning of Mixed Membership Models
Zilong Tan Sayan Mukherjee






Efficient Learning of Mixed Membership Models
    Zilong Tan

Department of Computer Science

Duke University 

Email: ztan@cs.duke.edu
Sayan Mukherjee

Departments of Statistical Science

Computer Science, Mathematics, 

Biostatistics & Bioinformatics 

Duke University 

Email: sayan@stat.duke.edu
    December 30, 2023
=============================================================================================================================================================================================================================================================





We present an efficient algorithm for learning mixed membership models when the number of variables p is much larger than the number of hidden components k. This algorithm reduces the computational complexity of state-of-the-art tensor methods, which require decomposing an O(p^3) tensor, to factorizing O(p/k) sub-tensors each of size O(k^3). In addition, we address the issue of negative entries in the empirical method of moments based estimators. We provide sufficient conditions under which our approach has provable guarantees. Our approach obtains competitive empirical results on both simulated and real data.




§ INTRODUCTION


Mixed membership models <cit.> have been used extensively across applications ranging from modeling population structure in genetics  <cit.> to topic modeling of documents <cit.>. Mixed membership models use  Dirichlet latent variables to define cluster membership where samples can partially belong to each of k latent components. Parameter estimation for such latent variables models (LVMs) using maximum likelihood methods such as expectation maximization is computationally intensive for large data, for example, if number of samples n is large.

Parameter estimation using the method of moments for LVMs is an attractive scalable alternative that has been shown to have certain theoretical and computational advantages over maximum likelihood methods in the setting when n is large. For LVMs, method of moments approaches reduce to tensor methods—the moments of the model parameters are expressed as a function of statistics of the observations in a tensor form. Inference in this setting becomes a problem of tensor factorization. Computational advantages of using tensor methods have been observed for many popular models, including latent Dirichlet allocation <cit.>, spherical Gaussian mixture models <cit.>, hidden Markov models <cit.>, independent component analysis <cit.>, and multi-view models <cit.>. An appealing property of tensor methods is the guarantee of a unique decomposition under mild conditions <cit.>. 

There are two complications to using standard tensor decomposition methods <cit.> for LVMs. The first problem is computation and space complexity. Given p variables in the LVM, parameter inference requires factorizing typically a non-orthogonal estimator tensor of size O(p^3) <cit.>, which is prohibitive for large p. When the estimator is orthogonal and symmetric, this can be done in O(p^2 log p) <cit.>. Online tensor decomposition <cit.> uses dimension reduction to instead factorize a reduced k-by-k-by-k tensor. However, the dimension reduction can be slower than decomposing the estimator directly for large sample sizes, as well as suffer from high variance <cit.>. We introduce a simple factorization with improved complexity for the general case where the parameters are not required to be orthogonal.

The second problem arises from negative entries in the empirical moments tensor. LVMs for count data are constrained to have nonnegative parameters. However, the empirical moments tensor computed from the data may contain negative elements due to sampling variation and noise. Indeed, for small sample sizes or data with many small or zero counts, there will be many negative entries in the empirical moments tensor. General tensor decomposition algorithms <cit.>, including the tensor power method (TPM) <cit.>, do not guarantee the nonnegativity of model parameters. Approaches such as positive/nonnegative tensor factorization <cit.> also do not address this situation as they require all the elements of the tensor to be factorized to be nonnegative. With robust tensor methods <cit.>, sparse negative entries may potentially be treated as corrupted elements; however, these methods are not applicable in this setting since there can be many negative elements.

In this paper, we introduce a novel parameter inference algorithm called partitioned tensor parallel quadratic programming (PTPQP) that is efficient in the setting where the number of variables p is much larger than the number of latent components k. The algorithm is also robust to negative entries in the empirical moments tensor. There are two key innovations in the PTPQP algorithm. The first innovation is a partitioning technique which recovers the parameters through factorizing O(p/k) much smaller sub-tensors each of size O(k^3). This technique can also be combined with methods <cit.> to obtain further improved complexities. The second innovation is a parallel quadratic programming <cit.> based algorithm to factor tensors with negative entries under the constraint that the factors are all nonnegative. To the best of our knowledge, this is the first algorithm designed to address the problem of negative entries in empirical estimator tensors. We show that the proposed factorization algorithm converges linearly with respect to each factor matrix. We also provide sufficient conditions under which the partitioned factorization scheme is consistent, the parameter estimates converge to the true parameters.



§ PRELIMINARIES



  
Notations
We use bold lowercase letters to represent vectors and bold capital letters for matrices. Tensors are denoted by calligraphic capital letters. The subscript notation A_j refers to j-th column of matrix A. We denote the  j-th column
of the identity matrix as e_j and 1 is a vector of ones. We further write diag(x) for a diagonal matrix whose diagonal entries are x, and diag(A) to mean a vector of the diagonal entries of A.

Element-wise matrix operators include ≻ and ≽, e.g., A≽ 0 means that A has nonnegative entries. (·)_+ refers to element-wise max(·,0). * and ⊘ respectively represent element-wise multiplication and division. Moreover, × refers to the outer product and ⊙ denotes the Khatri-Rao product. ·_F and ·_2 represent the Frobenius norm and spectral norm, respectively. 



  
Tensor basics
This paper uses similar tensor notations as <cit.>. In particular, we are primarily concerned with Kruskal tensors in ℝ^d_1 × d_2 × d_3, which can be expressed in the form of

    𝒯 = ∑_j=1^r A_j ×B_j ×C_j,

where A, B, and C are respectively d_1-by-r, d_2-by-r, and d_3-by-r factor matrices. The rank of 𝒯 is defined as the smallest r that admits such a decomposition. The decomposition is known as the CP (CANDECOMP/PARAFAC) decomposition. The j-mode unfolding of 𝒯, denoted by 𝒯_(j), for j=1,2,3 is a d_j-by-(∏_t≠ j d_t) matrix whose rows are serializations of the tensor fixing the index of the j-th dimension. The unfoldings have the following well-known compact expressions:

    𝒯_(1) = A(C⊙B)^⊤,
      𝒯_(2) = B(C⊙A)^⊤,
      𝒯_(3) = C(B⊙A)^⊤.




§ LEARNING THROUGH METHOD OF MOMENTS




 §.§ Generalized Dirichlet latent variable models

A generalized Dirichlet latent variable model (GDLM) was proposed in <cit.> for the joint distribution of n observations y_1,y_2,⋯,y_n. Each observation y_i consists of p variables y_i = (y_i1,y_i2,⋯,y_ip)^⊤. GDLM assumes a generative process involving k hidden components. For each observation, sample a random Dirichlet vector x_i = (x_i1,x_i2,⋯,x_ik)^⊤∈Δ^k-1 with concentration parameter α = (α_1,α_2,⋯,α_k)^⊤. The elements of x_i are the membership probabilities for y_i to belong to each of the k components. Specifically,

    y_ij∼∑_h=1^k x_ih g_j (θ_jh),

where g_j(θ_jh) is the density of the j-th variable specific to component h with parameter θ_j = (θ_j1,θ_j2,⋯, θ_jk). One advantage of GLDM is that y_ij can take categorical values. Let d_j denote the number of categories for the j-th variable (set d_j = 1 for scalar variables), θ_j becomes a d_j-by-k probability matrix where the c-th row corresponds to category c. We aim to accurately recover θ_j from independent copies of y_i involving variables of mixed data types, either categorical or non-categorical.



 §.§ Moment-based estimators


The moment estimators of latent variable models typically take the form of a tensor <cit.>. Consider the estimators of GDLM <cit.> for example. Let b_ij = e_y_ij if variable j is categorical; b_ij = y_ij otherwise. The second- and third- order parameter estimators for variable j, s, and t are written

    ℳ^js   =
    𝔼[b_ij×b_is] - α_0/α_0 + 1𝔼[b_ij] 𝔼[b_is]^⊤
    ℳ^jst   = 𝔼[b_ij×b_is×b_it]
    + 2α_0^2/(α_0 + 1)(α_0 + 2)𝔼[b_ij] ×𝔼[b_is] ×𝔼[b_it]
       = - α_0/α_0 + 2(
     𝔼[𝔼[b_ij] ×b_is×b_it]
    + 𝔼[b_ij×𝔼[b_is] ×b_it]
    + 𝔼[b_ij×b_is×𝔼[b_it]]
    ).

Alternatively, ℳ^js and ℳ^jst have the following CP decomposition into parameters θ_j:

    ℳ^js   = ∑_h≥ 1α_h/α_0 (α_0 + 1)θ_jh×θ_sh,   θ_uv∈ℝ^d_i
    ℳ^jst   = ∑_h≥ 12α_h/α_0(α_0 + 1)(α_0 + 2)θ_jh×θ_sh×θ_th,   θ_uv∈ℝ^d_i.

<ref> provides the derivation details of these estimators. For the special case of latent Dirichlet allocation, ℳ^js and ℳ^jst are scalar joint probabilities.

The parameters θ_j are typically obtained by factorizing the block tensor ℳ_2 whose (j,s)-th element is the empirical ℳ^js and/or ℳ_3 whose (j,s,t)-th element is the empirical ℳ^jst <cit.>. Note that θ_j are generally non-orthogonal, and thus preprocessing steps (see <ref>) are needed for orthogonal decomposition methods <cit.>. The preprocessing can be expensive and often leads to suboptimal performance <cit.>. Here, we highlight a few relevant observations:


  * ℳ^js alone does not yield unique parameters θ_j due to the well-known rotation problem. Suppose that θ_j^* and θ_s^* are the ground-truth parameters satisfying (<ref>) and any invertible R, there exists decomposition θ_j^' = θ_j^* R and θ_s^' = R^-1θ_s^* that also satisfy (<ref>) but are not ground-truth parameters.
 The ground-truth parameters are not uniquely identifiable through ℳ^js, this is true even when enforcing nonnegativity constraints on parameters <cit.>.

  * ℳ^jst is sufficient to uniquely recover the parameters under certain mild conditions <cit.>; for example, when any two of θ_j, θ_s, and θ_t have linearly independent columns and the columns of the third are pair-wise linearly independent <cit.>.

  * The empirical estimator ℳ^jst generally contains negative entries due to variance and noise. The fraction of negative entries can approach 50%, as we shall see in experiments. We address this issue in <ref>.

  * While the decomposition (<ref>) can be unique up to permutation and rescaling, the correspondence between each column of the factor matrix and each hidden component may not be consistent across multiple decompositions. Techniques for achieving consistency are developed in <ref>.




 §.§ Computational complexity


Tensor methods such as TPM typically decompose the O(p^3 d_max^3) full estimator tensor that includes all variables. More efficient algorithms have been developed for the case that parameters are orthogonal <cit.>, and when the sample size is small <cit.>. However, these methods do not apply in the general case where the parameters are non-orthogonal and the sample size can be potentially large. A key insight underlying our approach is that it is sufficient to recover the parameters by factorizing only O(p/k) much smaller sub-tensors each of size O(k^3). This technique can also be combined with the aforementioned methods to further improve the complexity in certain cases.



§ AN EFFICIENT ALGORITHM

In this section, we develop partitioned tensor parallel quadratic programming (PTPQP) an efficient approximate algorithm for learning mixed membership models. We first introduce a novel partitioning-and-matching scheme that reduces parameter estimation to factorizing a sequence of sub-tensors. Then, we develop a nonnegative factorization algorithm that can handle negative entries in the sub-tensors.



 §.§ Partitioned factorization


Factorizing the full tensor formed by all ℳ^jst is expensive while a three-variable tensor ℳ^jst in (<ref>) alone may not be sufficient to determine θ_j when k is large. In this section, we consider factorizing the sub-tensors corresponding to a cover of the set of variables [p] such that each sub-tensor admits an identifiable CP decomposition (<ref>), i.e. unique up to permutation and rescaling of columns. This gives the parameters for all variables. Suppose that p > k and the maximum number of categories d_max is a constant, the aggregated size of the sub-tensors can be much smaller, i.e., O(p k^2), than the size O(p^3) of the full estimator.

Let π^j, π^s, and π^t denote ordered subsets ⊆[p], with cardinality |π_j| = p_j, |π_s| = p_s, and |π_t| = p_t, respectively. Consider the p_j-by-p_s-by-p_t block tensor [For block tensor operations, see e.g., <cit.>.] ℳ^π^jπ^sπ^t whose (u,v,w)-th element is the tensor ℳ_uvw^π^jπ^sπ^t = ℳ^π_u^jπ_v^sπ_w^t. From (<ref>), we have that

    ℳ^π^jπ^sπ^t =
    ∑_h=1^k
    2α_h/α_0(α_0 + 1)(α_0 + 2)[   θ_π_1^jh;   θ_π_2^jh;          ⋮; θ_π_p_j^jh ]×[   θ_π_1^sh;   θ_π_2^sh;          ⋮; θ_π_p_s^sh ]×[   θ_π_1^th;   θ_π_2^th;          ⋮; θ_π_p_t^th ].

Clearly, the block tensor is identifiable if it has an identifiable sub-tensor. Suppose that a sub-tensor ℳ^π^uπ^vπ^w is identifiable, then one can construct an identifiable tensor ℳ^π^j'π^s'π^t' from  ℳ^π^jπ^sπ^t by setting 

    π^j' = π^j ∪π^u,   π^s' = π^s ∪π^v,   π^t' = π^t ∪π^w.

We further remark that a sub-tensor can be identifiable under mild conditions, for example, if the sum of the Kruskal rank of the three factor matrices is at least than 2 k + 2 <cit.>. 

Given an identifiable sub-tensor ℳ^π^uπ^vπ^w of anchor variables indexed by π^u, π^v, and π^w, the partitioning produces a set of sub-tensors (partitions) constructed through (<ref>), that includes all variables. Thus, ℳ^π^uπ^vπ^w is a common sub-tensor shared across all partitions. We choose anchor variables whose parameter matrices are of full column rank to obtain an identifiable ℳ^π^uπ^vπ^w. Finally, one can divide the rest of variables evenly and randomly into the partitions.



 §.§ Matching parameters with hidden components


Since the factorization of a partition (<ref>) can only be identifiable up to permutation and rescaling of the columns of constituent θ_j, the correspondence between the columns of θ_j and hidden components can differ across partitions. To enforce consistency, we associate a permutation operator ψ^j for each variable j such that (ψ^j θ_j)_h are the parameters specific to hidden component h across all variables j. Consider the following vector representation of ψ:

    ψ    = (ψ_1, ψ_2, ⋯, ψ_k),   ψ_i ∈[k]
    ψA   = [A_ψ_1, A_ψ_1, ⋯, A_ψ_k].

Observe that ψ^j=ψ^s=ψ^t within a factorization of ℳ^jst, and this also holds for the partitioned factorization (<ref>) of ℳ^π^jπ^sπ^t as well, i.e., ψ^x = ψ^y,   ∀ x, y ∈π^j∪π^s∪π^t.

Consider the factorizations of ℳ^π^jπ^sπ^t and ℳ^π^uπ^vπ^w and suppose that ∃ x∈(π^j∪π^s∪π^t) ∩(π^u∪π^v∪π^w). The permutation
operator for one factorization is determined given the other by column matching the parameters of variable x in both factorizations. Thus, an inductive way to achieve a consistent factorization is to start with one factorization, and let its permutation be the identity (1,2,⋯,k), then perform the factorization over new sets of variables with at least one variable in common with the initial factorization. Permutations for the sequential factorizations are determined via column matching parameter matrices of the common variables.

Given two factorized parameter matrices θ_j and θ_j^' of variable j, our goal is to find a consistent permutation ψ (of θ_j with respect to θ_j^') such that (ψθ_j)_h and θ_jh^' correspond to the same hidden component for all h ∈[k]. We now present an algorithm with provable guarantees to compute a consistent permutation.



  
Smallest angle matching

A simple matching algorithm is to match the two columns of the two parameter matrices that have the smallest angle between them. Consider the factorizations of ℳ^jst and ℳ^juv which yield respectively parameters θ_j and θ_j^' for the common variable j. Given the permutation ψ^j for ℳ^jst, the permutation ψ^u for ℳ^juv is computed by:

    ψ_s^u = max_t ( θ̅_j^'⊤ψ^j θ̅_j )_ts.

Here, θ̅_j and θ̅_j^' represent respectively the normalized θ_j and θ_j^' with each column having unit Euclidean norm. 

There are cases that ψ^u computed via (<ref>) is not consistent: 1) ψ^u contains duplicate entries and hence is ineligible; and 2) since θ_j and θ_j^' are the factorized parameter matrices which are generally perturbed from the ground-truth, the resulting ψ^u may differ from the consistent permutation. To cope with these cases, we establish in <ref> the sufficient conditions for ψ^u to be consistent.



  
Orthogonal Procrustes matching

One issue with the smallest angle matching is that each column is paired independently. It is easy for multiple columns to be paired with a common nearest neighbor. We describe a more robust algorithm based on the orthogonal Procrustes problem, and show improved guarantees. Since a consistent permutation is orthogonal, a natural relaxation is to only require the operator to be orthogonal. This is an orthogonal Procrustes problem, formulated in the same settings as <ref>

    min_Ψθ̅_j^'Ψ - ψ^j θ̅_j_F^2,
      s.t. Ψ^⊤Ψ = I.

Let θ̅_j^'⊤ψ^j θ̅_j = UΣV^⊤ be the singular value decomposition (SVD), the solution Ψ^* is given by the polar factor <cit.>

    Ψ^*    = UV^⊤.

Here, Ψ^* is orthogonal and does not immediately imply the desired permutation ψ^u. To compute ψ^u, one can additionally restrict Ψ to be a permutation matrix, and solve for ψ^u using linear programming <cit.>. Aside from efficiency, one fundamental question is that under what assumptions the objective (<ref>) yields the consistent permutation.

Given the solution Ψ^* to the Procrustes problem, we propose the following simple algorithm for computing ψ^u:

    ψ_s^u = max_t Ψ_ts^*
    .

We first establish through <ref> that if ψ^u obtained using (<ref>) is a valid permutation, i.e., no duplicate entries, then it is optimal in terms of the objective (<ref>).


The ψ^u obtained using (<ref>) satisfies

    ψ^uθ̅_j^' - ψ^j θ̅_j_F^2 ≤ψθ̅_j^' - ψ^j θ̅_j_F^2

for all permutations ψ.


First, rewrite the objective (<ref>) as follows

    θ̅_j^'Ψ - ψ^j θ̅_j_F^2
       = Ψ^⊤θ̅_j^'⊤θ̅_j^'Ψ
    + (ψ^jθ̅_j)^⊤ψ^jθ̅_j
    - 2Ψ^⊤θ̅_j^'⊤ψ^jθ̅_j 
       = θ̅_j^'_F^2 + θ̅_j_F^2
    - 2Ψ^⊤θ̅_j^'⊤ψ^jθ̅_j.

Recall the SVD θ̅_j^'⊤ψ^jθ̅_j = UΣV^⊤, and write Ψ = UV^⊤ + E. Keeping only terms that depend on E in (<ref>) to obtain -2E^⊤UΣV^⊤. Thus, the optimization (<ref>) is equivalent to

    max_EV^⊤E^⊤UΣ,
      s.t.  (UV^⊤ + E)^⊤(UV^⊤ + E) = I.

From the constraint, we obtain E^⊤E = -2V^⊤E^⊤U. The optimization now becomes

    min_EE^⊤EΣ = min_E∑_j (E^⊤E)_jjΣ_jj.

Let us now restrict each column of Ψ to be in {e_j  |  j=1,2,⋯,k }, but not necessarily distinct. Suppose that Ψ_j = e_y. We have that E_j = e_y - (UV^⊤)_j. Clearly, (<ref>) and hence (<ref>) are minimized with y = max_t (UV^⊤)_tj. 



In section <ref> we state sufficient conditions under which the objective (<ref>) yields a consistent permutation.



 §.§ Approximate nonnegative factorization


In previous sections, we reduced the inference problem to factorizing partitioned sub-tensors. We now present a factorization algorithm for the sub-tensors that contain negative entries. Our goal is to approximate a sub-tensor ℳ by a sub-tensor ℳ = ∑_j A_j ×B_j ×C_j where the factors A, B, and C are nonnegative. The Frobenius norm is used to quantify the approximation

    min_A, B, C≽ 0ℳ - ℳ_F.

Note that we do not assume that ℳ≽ 0 in (<ref>) which distinguishes our optimization problem 
from other approximate factorization algorithms <cit.>. In <ref>, we provide some details as to why negative entries are problematic for standard approximate factorization algorithms.
We can rewrite (<ref>) using the 1-mode unfolding as

    min_A, B, C≽ 0   ℳ_(1) - A(C⊙B)^⊤_F .

Equivalent formulations with respect to the 2-mode and 3-mode unfoldings can be readily obtained from (<ref>). 

We point out that another widely-used error measure — the I-divergence <cit.> — may not be suitable for our learning problem. The optimization using I-divergence is given by

    min_A, B, C≽ 0∑_u,v,w   [ℳ_uvwlogℳ_uvw/ℳ_uvw - ℳ_uvw + ℳ_uvw].

This optimization is useful for nonnegative ℳ when each entry follows a Poisson distribution. In this case, the objective is equivalent to the sum of Kullback-Leibler divergence across all entries of ℳ:

    ∑_u,v,w D_KL(Pois(x; ℳ_uvw)  Pois(x;ℳ_uvw)).

However, the Poisson assumption does not generally hold for the estimator tensor (<ref>).



 §.§ Handling negative entries in empirical estimators


We first illustrate that factorizing a tensor with negative entries using either positive tensor factorization <cit.> or nonnegative tensor factorization <cit.> will either result in factors that violate 
the the nonnegativity constraint or the result of the algorithm diverges. In addition, we show that general tensor decompositions cannot enforce the factor nonnegativity even after rounding the negative entries to zero.

We then present a simple method based on weighted nonnegative matrix factorization (WNMF) <cit.> that enforce the factor nonnegativity constraint. We further generalize this method using parallel quadratic programming (PQP) <cit.> to obtain a method with a provable convergence rate.



  
Issue of negative entries

If the tensor is strictly nonnegative, the optimization specified in (<ref>) can be reduced to nonnegative matrix factorization (NMF). Solvers abound for NMF including the celebrated Lee-Seung's multiplicative updates <cit.>. The reduction is done by viewing (<ref>) as Y - WH_F^2 with Y = ℳ_(1)^jst, W = A, and H = (C⊙B)^⊤, and alternating

    W_st← W_st(YH^⊤)_st/(WHH^⊤)_st,

over each unfolding and factor matrix W. Obviously, the updates may yield negative entries in W when the unfolding contains negative entries. In addition, convergence relies on the nonnegativity of the unfolding <cit.>. This issue extends to their tensor factorization variants <cit.> known as the positive tensor factorization and nonnegative tensor factorization. For these approaches, a naive resolution is to round negative entries of ℳ^jst to 0, this however lacks theoretical guarantees.

It is important to note that the rounding does not help general tensor decompositions like TPM. The following example illustrates that the unique decomposition (up to permutation and rescaling) of a positive tensor can contain negative entries. Consider a 2-by-2-by-2 positive tensor, whose 1-mode unfolding is given by

    [ [c c | c c]
    1             3             2             2;             2             2             2             2 ],

where the vertical bar separates two frontal slices. It has the following decomposition, written in the form of (<ref>):

    A = C = [ 1 1; 1 0 ],   B = [  2 -1;  2  1 ].

Since all factors are of full-rank, the decomposition is unique up to permutation and rescaling of columns <cit.>. Thus, a general tensor decomposition yields a B with negative entries regardless of rescaling.



 §.§ Factorization via WNMF

Since the ground-truth ℳ^jst are nonnegative, we may “ignore" the negative entries of ℳ^jst by treating them as missing values. This idea leads to the following modified objective:

    min_W,H≽ 0Ω * (Y - WH)_F^2

where Y, W, H are chosen identically as (<ref>), and we define

    Ω_uv = 
    1,    Y_uv≥ 0
    
    0,    Y_uv < 0
    .

The optimization can be carried out using WNMF. Here, we modify the original updates by introducing a positive constant ϵ to ensure that the updates are well-defined:

    W_uv← W_uv[(Ω*Y) H^⊤]_uv + ϵ/[((WH)*Ω) H^⊤]_uv + ϵ,
      
    H_uv← H_uv[W^⊤(Ω*Y) ]_uv + ϵ/[W^⊤(Ω*(WH)) ]_uv + ϵ.

<ref> states the correctness of the modified updates (<ref>).



The objective (<ref>) is non-increasing under the multiplicative updates (<ref>).


We prove the update for H, and the update for W follows by applying the update to Ω^⊤ * (v^⊤ - H^⊤W^⊤)_F. First, consider the error Frobenius norm for a column h of H, and the corresponding columns ω of Ω and v of V,

    F(h) = ω * (v - Wh)_F^2.

The following G(·,·) is an auxiliary function of F(·):

    G(h,h^t) = F(h) + (h - h^t)^⊤∇ F(h^t) + 1/2(h - h^t)^⊤K(h - h^⊤),

where we define

    K = diag([
    W^⊤diag(ω) Wh^t + ϵ1] ⊘h^t
    ).

Clearly, G(h,h) = F(h), and one can show that G(h,h^t) ≥ F(h) by rewriting

    F(h) = F(h) + (h - h^t)^⊤∇ F(h^t) + 1/2(h - h^t)^⊤W^⊤diag(ω) W(h - h^⊤),

where we note that ω*ω = ω from the Boolean definition of ω. Comparing (<ref>) with (<ref>), it is sufficient to show that K - W^⊤diag(ω) W is positive semi-definite. Now consider the scaled matrix

    U   = diag(h^t) Kdiag(h^t) - 
    diag(h^t) W^⊤diag(ω) Wdiag(h^t)
       = diag(W^⊤diag(ω) Wh^t + ϵ1) diag(h^t) - 
    diag(h^t) W^⊤diag(ω) Wdiag(h^t).

Observe that U is strictly diagonally dominant as U1≻ 0 and the off-diagonal entries are negative. Also note that all diagonal entries of U are positive, it follows that U is positive semi-definite. We thereby conclude that K - W^⊤diag(ω) W is positive semi-definite.

Let h^t+1 = min_h G(h,h^t), we have that F(h^t) = G(h^t,h^t) ≥ G(h^t+1,h^t) ≥ F(h^t+1). The minimizer h^t+1 is obtained by setting ∇_h^t+1 G(h^t+1,h^t) = 0, which yields

    - ∇ F(h^t)    = K(h^t+1 - h^t)
    W^⊤[ω*(v - Wh^t)] 
       = Kh^t+1 - W^⊤diag(ω) Wh^t - ϵ1
    h^t+1   = h^t * [W^⊤(ω*v) + ϵ1] ⊘[W^⊤(ω * (Wh^t)) + ϵ1].

The particular choice of Ω guarantees that h^t+1 is always positive.




 §.§ Parallel quadratic programming

We now generalize the WNMF approach using parallel quadratic programming to obtain a convergence rate. Let 𝕊_++ denote the set of symmetric positive definite matrices, we consider the following optimization problem

    min_x1/2x^⊤Qx + z^⊤x  s.t.  x≥ 0,  Q∈𝕊_++,

which can be solved by iterating multiplicative updates <cit.>. We use the parallel quadratic programming (PQP) algorithm <cit.> to solve (<ref>), partly because it has a provable linear convergence rate. The PQP multiplicative update for (<ref>) takes the following simple form:

    x←x * (Q^- x + z^-) ⊘(Q^+ x + z^+),

with

    Q^+    = (Q)_+ + diag(γ),      Q^-    = (-Q)_+ + diag(γ)
    z^+    = (z)_+ + ϕ,      z^-    = (-z)_+ + ϕ.

Here γ and ϕ are arguments to PQP, we will discuss these arguments in section <ref>. The update maintains nonnegativity since all items are nonnegative. We make the following observation.



The multiplicative updates for Lee-Seung and WNMF are special cases of PQP.


Since the WNMF (<ref>) generalizes Lee-Seung, which is the case that Ω has all ones, we need only to prove for WNMF. Let Λ = Ω*Ω and γ = 0, some matrix algebra reveals the following PQP updates

    W_uv   ← W_uv[((Λ*Y) H^⊤)_+]_uv + Φ_uv/[((WH)*Λ) H^⊤ + ((-Λ*Y) H^⊤)_+ ]_uv + Φ_uv
    
    H_uv   ← H_uv[(W^⊤(Λ*Y))_+]_uv + Φ_uv^'/[W^⊤(Λ*(WH) ) + (-W^⊤(Λ*Y) )_+ ]_uv + Φ_uv^'.

Comparing (<ref>) to (<ref>), they are equivalent if Φ_uv = Φ_uv^' = ϵ.


We can now solve the approximate nonnegative factorization problem stated in (<ref>) using (<ref>). <ref> states the multiplicative updates. A more detailed discussion of Φ is included in <ref>. We present pseudo-code in 
<ref>.



For optimization (<ref>), the following update converges linearly to a local optimum

    A←A * [(-Z)_+ + Φ] ⊘[AQ + (Z)_+ + Φ]

with

    Q = (C^⊤C)*(B^⊤B), 
        Z = -ℳ_(1)(C⊙B)
    Φ≻1/2(
    √(diag(ZQ^-1Z^⊤)/λ_min(Q))diag(Q)^⊤
    - |Z|
    )_+,

where λ_min(·) is the smallest eigenvalue. Similar updates for B and C are obtained using (<ref>).


We apply PQP updates (<ref>) to each row of A. Let v_j: and A_j: be the j-th row of ℳ_(1) and A, respectively. Fixing the current factor estimates B and C, the optimization with respect to A_j: follows from (<ref>):

    min_A_j:≽ 0v_j:^⊤ - (C⊙B)A_j:^⊤_F 
    = min_A_j:≽ 01/2A_j:(C⊙B)^⊤(C⊙B) A_j:^⊤
    - v_j:(C⊙B) A_j:^⊤.

Now the updates (<ref>) can be applied immediately, where we set γ = 0 and Φ according to <ref> in <ref>. Using the identity (C⊙B)^⊤(C⊙B) =  (C^⊤C)*(B^⊤B) and performing the updates simultaneous for all rows of A give (<ref>).








 §.§ Proposed approach

To summarize, the proposed approach, referred to as PTPQP, consists of three steps. Given the indexes of anchor variables π^u∪π^v∪π^w, the variables [p]\(π^u∪π^v∪π^w) are first evenly divided into r partitions, and the anchor variables are added to each partition. The second step consists of forming and factorizing the sub-tensor of each partition using <ref>, this step can be parallelized. Third, normalize the anchor matrix [θ^π^u⊤,θ^π^v⊤,θ^π^w⊤]^⊤ formed by the anchor variable parameters to have unit column Euclidean norm, and then use either (<ref>) or (<ref>) to match over the anchor matrix.



  
Efficiency
Most of the computational cost is in the factorization. Consider one partition, and let ℳ^π_j π_s π_t be the corresponding sub-tensor, the sub-tensor size is ∏_π∈{π^j,π^s,π^t}∑_h∈π d_h. The maximum number of categories for a variable is generally a constant for the GDLM. Under smallest partitioning, this size is determined by the sub-tensor of anchor variables, i.e., O(k^3), which corresponds to  (p/k) partitions. One benefit of PTPQP is that the number of sub-tensor factorizations is linear in p due to the partitioned factorization, this results in significant efficiency gains when p ≫ k. Furthermore, PTPQP is easy to be parallelized across multiple CPUs and machines, since the computation as well as data are not distributed across partitions. 



§ PROVABLE GUARANTEES


In this section, we state the main theoretical results of the proposed partitioned factorization and tensor PQP factorization.



 §.§ Sufficient conditions for guaranteed matching

<ref> and <ref> state that when the anchor parameter matrices from two factorizations are “close", the proposed matching algorithms obtain a consistent permutation.



Suppose that θ_j is the ground-truth matrix for variable j. Solving (<ref>) results in a consistent permutation if for all factors θ_j of variable j

    θ_jh - θ_jh_2/θ_jh_2
    < 1 - √(1/2 + √(1/8(1+max_u<v(θ̅_j^⊤θ̅_j)_uv)))

for all h∈[k], where θ̅_jh = θ_jh / θ_jh_2.


Consider the smallest pair-wise angle α_min between the columns of θ̅_j, we have that

    cosα_min = max_u<v(θ̅_j^⊤θ̅_j)_uv.

Denote by α the maximum angle between the column of a factorized parameter matrix θ_j and the corresponding column of the ground-truth. It is sufficient to ensure that

    α < 1/4α_min.

Consider any two columns s≠ t of the ground-truth parameter matrix, and the corresponding perturbed columns {θ_js, θ_jt} and {θ_js^', θ_jt^'} from two factorizations. We have that

    ∠(θ_js, θ_js^')
       ≤ 2α
    ∠(θ_js, θ_jt^')
       ≥∠(θ̅_js, θ̅_jt)
    - 2α
       ≥α_min - 2α.

From (<ref>), we have that 

    ∠(θ_js, θ_js^')
    <
    ∠(θ_js, θ_jt^'),

as desired for (<ref>) to work correctly. Now consider the inner product of a perturbed column and the ground-truth, it holds that

    <θ_jh/θ_jh,
    θ_jh + ϵ/θ_jh + ϵ>    =
    θ_jh^2 - ϵ^2
    + θ_jh + ϵ^2/2θ_jhθ_jh+ϵ
       ≥θ_jh- ϵ/2θ_jh +
    θ_jh+ϵ/2θ_jh
       ≥ 1 - ϵ/θ_jh.

Thus, a sufficient condition for (<ref>) to yield the consistent permutation is

    1 - ϵ/θ_jh > cos(1/4α_min),

which written in analytic form proves the theorem.


<ref> states that one obtains a consistent permutation by solving (<ref>) in the columns of the ground-truth parameter matrix
are distinct from each other in angles and the factorized parameter matrix is near the ground-truth in Frobenius norm. Thus, a good anchor variable for the partitioned factorization (<ref>) is one whose parameter matrix has distant columns in angles.

The bound in <ref> can be made sharp for certain θ_j, and thus the smallest angle matching algorithm has general guarantees only when the perturbation is small, i.e., the relative error ratio is less than 1- √(2+√(2))/2 ≈ 1/13.



Suppose that θ and θ^' are two factorized parameter matrices for a variable. Solving (<ref>) results in a consistent permutation ψ, if

    E_2 < σ_k(θ^⊤θ)
      and  
    -E_2/ρlog(
    1 - ρ/ν)
    < 2-√(2)/4

with

    ρ = σ_1(E) + σ_2(E),   ν = σ_k(θ^⊤θ)
    + σ_k-1(θ^⊤θ)

where the error matrix is define as E = (ψθ)^⊤( θ^' -ψθ), and σ_j(·) denotes the j-th largest singular value.

The proof of Theorem <ref> follows from the following two Lemmas.


Suppose that ψ is the consistent permutation of θ with respect to θ^'. Formula (<ref>) is guaranteed to recover ψ, if

    diag(UV^⊤) ≻√(2)/21,

where U and V are the left and right singular matrices of (ψθ)^⊤θ^'.


We need to show that (<ref>) yields ψ for the orthogonal Procrustes problem min_Ψ^⊤Ψ = IθΨ - θ^'_F. From the solution (<ref>), it is easy to show that the minimizer Ψ^* of min_Ψ^⊤Ψ = Iθ^'Ψ - θ_F and the minimizer Ψ^' of min_Ψ^⊤Ψ = I(ψθ) Ψ - θ^'_F satisfy

    Ψ^'⊤ = ψΨ^*.

Note that Ψ^*⊤ is the desired minimizer of  min_Ψ^⊤Ψ = IθΨ - θ^'_F, and thus it remains to show that (<ref>) gives ψ when applied to Ψ^*⊤, or equivalently 

    max_t Ψ_st^* = ψ_s.


Since the row and column vectors of Ψ^* have unit Euclidean norm, the following dual statements imply each other

    max_t Ψ_ts^' = j   ⇔  max_t Ψ_jt^' = s     ∀ j,s ∈[k],

if condition (<ref>) holds. Under this condition, we also have that (<ref>) gives the identity permutation [1,2,⋯,k] for the orthogonal Procrustes problem min_Ψ^⊤Ψ = I(ψθ) Ψ - θ^'_F. Thus, applying (<ref>) to both sides of (<ref>) yields

    max_t Ψ_tψ_s^* = s,

which implies (<ref>) from (<ref>).


 (Mathias)

Suppose that A∈ℝ^n × n is nonsingular. Then for any E∈ℝ^n× n with σ_1(E) < σ_n(A) and any unitarily invariant norm ·, it holds that

    μ(A + E) - μ(A)≤
    -2E/E_2log(
    1 - E_2/σ_n(A) +
    σ_n-1(A)),

where μ(·) represents the unitary factor of the polar decomposition, and ·_k is the Ky Fan k-norm.


Let H = (ψθ)^⊤ψθ, which has the same singular values as θ^⊤θ. Denote by μ(·) the unitary factor of the polar decomposition. Using the fact that μ(H) = I, the sufficient condition of <ref> is restated as

    diag(μ(H) - μ(H + E)) ≺(1 - √(2)/2)1.

Also note that

    max_t|diag(μ(H) - μ(H + E))_t| ≤μ(H) - μ(H + E)_2.

Thus, it suffices to enforce the right term to be less than 1-√(2)/2. From <ref>, this can be achieved by letting

    -2E_2/E_2log(
    1 - E_2/σ_n(A) +
    σ_n-1(A))
    ≤
    1 - √(2)/2.




The first condition in <ref> requires that at least one of θ and θ^' must have full column rank. We may exchange θ and θ^' in <ref> to first obtain the consistent permutation of θ^' with respect to θ, ψ then follows immediately.

<ref> states that solving (<ref>) recovers a consistent permutation whenever the error spectral norm is small as compared to the smallest singular value of θ^⊤θ. This is especially useful for θ∈ℝ^d × k with the number of rows d much larger than the number of columns k. In particular, for θ with independent and identically distributed subgaussian entries, σ_k(θ^⊤θ) is at least of the order (√(d) - √(k-1))^2 <cit.>.



 §.§ Convergence


The following theorem states a sufficient condition for PQP to achieve linear convergence rate. The theorem statement and proof is an adaptation of results stated in <cit.>—the proof in <cit.> overlooks a required condition on ϕ and the condition
 γ≽diag(Q_jj) in the original proof is unnecessary.



The PQP algorithm given by (<ref>) monotonically decreases the objective (<ref>) and has linear convergence, if

    γ≽(-Q)_+ 1  and  ϕ≻1/2(
    √(z^⊤Q^-1z/λ_min(Q))diag(Q) - |z|
    )_+,

where λ_min(·) is the smallest eigenvalue.


First, the condition γ≥(-Q)_+ 1 suffices to ensure that the updates monotonically decrease (<ref>) <cit.>. Thus, it remains to show the condition on ϕ. Suppose that the i-th element of the optimum x^* is perturbed by a non-zero ϵ > - x_i^*. Let x = x^* + ϵe_i, and applying one update gives x^'. Denote the i-th row of Q^+, Q^-, and Q respectively by P_i, N_i, and Q_i, then it holds that P_i e_i = Q_ii + γ_i and N_i e_i = γ_i by definition. We now consider the ratio of errors between successive iterations:

    |x_i^' - x_i^*/x_i - x_i^*|    = |1/ϵ(N_i (x^* + ϵe_i) + z_i^-/P_i(x^* + ϵe_i) + z_i^+x_i - x_i^*)
    |
       = | N_ix^* + ϵγ_i + z_i^-/P_ix^* + ϵ Q_ii + ϵγ_i + z_i^+x_i^* + ϵ/ϵ - x_i^*/ϵ|
       = | N_ix^* + ϵγ_i + z_i^-/P_ix^* + ϵ Q_ii + ϵγ_i + z_i^+
    - x_i^*/ϵQ_ix^* + z_i + ϵ Q_ii/P_ix^* + ϵ Q_ii + ϵγ_i + z_i^+|.

From the KKT first-order optimality condition x_i^* (Q_i x^* + z_i) = 0, we simplify the ratio as

    |x_i^' - x_i^*/x_i - x_i^*|    = 
    | N_ix^* + ϵγ_i + z_i^- - x_i^* Q_ii/P_ix^* + ϵ Q_ii + ϵγ_i + z_i^+|.

Observe that the denominator is nonnegative. We also have that the denominator is greater than the numerator using the KKT optimality condition Q_i x^* + z_i ≥ 0:

    P_ix^* + ϵ Q_ii + ϵγ_i + z_i^+
    - (N_ix^* + ϵγ_i + z_i^- - x_i^* Q_ii)
    > Q_i x^* + z_i ≥ 0.

To achieve linear convergence rate, we may enforce the ratio to be less than one. Equivalently,

    P_ix^* + ϵ Q_ii + ϵγ_i + z_i^+
    + N_ix^* + ϵγ_i + z_i^- - x_i^* Q_ii > 0.

It suffices to set

    ϕ_i >
    1/2(Q_ii x_i^* - |z_i|)_+.


To get rid of x^* in (<ref>), we have the following inequality

    | 1/2x^*⊤Qx^* + z^⊤x^*| ≤1/2z^⊤Q^-1z,

where the right term is the negative of the minimum of the unconstrained problem, assuming that Q is non-singular. If Q is singular, then x^* can be unbounded. Further simplify the inequality using KKT optimality conditions as

    |x^*⊤Qx^* |    ≤ z^⊤Q^-1z
    λ_min(Q)x^*_2^2    ≤ z^⊤Q^-1z
    x^*_∞   ≤√(z^⊤Q^-1z/λ_min(Q)).

Combining with (<ref>) completes the proof.






§ RESULTS ON REAL AND SIMULATED DATA


We compare the proposed algorithm ptpqp with state-of-the-art approaches including: 1) the tensor power method tpm <cit.> and matrix simultaneous diagonalization, nojd0 and nojd1 <cit.>—two general tensor decomposition methods; 2) nonnegative tensor factorization hals <cit.>; and 3) generalized method of moments meld <cit.>. We use the online code provided by the corresponding authors.



 §.§ Learning GDLMs on simulated data

We adapt a simulation study from <cit.> to compare runtime and accuracy of parameter estimation. We consider a GDLM 
where each variable takes categorical values {0,1,2,3} and the parameters of the 
Dirichlet mixing distribution are {α_j = 0.1}_j=1^k. We initially consider 25 variables. The true parameters for each
hidden component h  are drawn from the Dirichlet distribution Dir(0.5,0.5,0.5,0.5). The resulting moment estimator is a 100-by-100-by-100 tensor. We vary the number of components k and add noise by replacing a fraction δ of the observations with draws from a discrete uniform distribution. We also vary the number of samples n=100,500,1000,5000, number of clusters  k=3,5,10,20, and contamination δ=0,0.05,0.1. Across these settings we found that the empirical third-order estimator typically exhibits between 20% and 50% negative entries.



  
Accuracy of inference
Accuracy is measured by root-mean-square error (RMSE) which we compare across algorithms as a function of the number of components 
for various sample sizes and levels of contamination, see <ref>. Both hals and ptpqp are consistently among the top estimators, and ptpqp outperforms hals as n grows. For small sample sizes and many hidden components meld achieves the smallest RMSE.  The RMSE of tpm is relatively large,  probably due to the whitening technique used to approximately transform the nonorthogonal factorization into an orthogonal one, see <cit.>. The most relevant observation is that  ptpqp outperforms other methods for large, noisy data. 





  
Computational cost We examined how runtime scales as a function of the number of partitions. For the same model we set p=1000 variables and n=1000 samples. The tensor is now 4000-by-4000-by-4000. We evaluated the runtime of ptpqp (without parallelization) with the 
number of partitions set to {30,40,50,100,200}. On a laptop with Intel i7-4702HQ@2.20GHz CPU and 8GB memory, ptpqp with 100 partitions completes within 3.5 min, 4 min, and 5 min for k=4,8,12, respectively. In addition, the runtime monotonically decreases with the number of partitions. Further speedups can be obtained by parallelizing the factorization of partitions across multiple CPUs or machines.



 §.§ Predicting crowdsourced labels

In <cit.>, a combination of EM and tensor decompositions was used to predict crowdsourcing annotations. The task is to predict the true label given incomplete and noisy observations from a set of workers, this is a mixed membership problem <cit.>. In <cit.> a third-order tensor estimator was proposed to obtain an initial estimate for the EM algorithm. We compare the predictive performance on five data sets of several tensor decomposition methods as well as the EM algorithm initialized with majority voting by the workers (MV+EM). The fraction of incorrect predictions and the size of each dataset are in the table below. Note that ptpqp matches or outperforms the other tensor methods on all but one dataset, and even outperforms MV+EM on two datasets.





§ CONCLUSIONS

We proposed an efficient algorithm for learning mixed mixture models based on the idea of partitioned factorizations. The key challenge is to consistently match the partitioned parameters with the hidden components. We provided sufficient conditions to ensure consistency. In addition, we have also developed a nonnegative approximation to handle the negative entries in the empirical method of moments estimators, a problem not addressed by several recent tensor methods. Results on synthetic and real data corroborate that the proposed approach achieves improved inference accuracy as well as computational efficiency than state-of-the-art methods.



§ CODE


Code for all the simulations is available from Zilong Tan's GitHub repository 

<https://github.com/ZilongTan/ptpqp>. 



§ ACKNOWLEDGEMENTS

Z.T. would like to thank Rong Ge for sharing helpful insights. S.M.  would like to thank Lek-Heng Lim for insights. Z.T. would like to acknowledge the support of grants NSF CNS-1423128,  NSF IIS-1423124, and NSF CNS-1218981. S.M. would like to acknowledge the support of grants NSF IIS-1546331, NSF DMS-1418261, NSF IIS-1320357, NSF DMS-1045153, and NSF DMS-1613261. 


IEEE



§ DIRICHLET MOMENTS


For a Dirichlet random vector x with concentration parameters α, the component moments can be easily shown by the integral

    𝔼[x_i]    = α_i/α_0     𝔼[x_i x_j]    =
    α_i α_j/α_0(α_0 + 1),    i≠ j
    𝔼[x_i^2]    =α_i(α_i+1)/α_0(α_0+1)     𝔼[x_j x_s x_t]    = α_j α_s α_t/α_0(α_0+1)(α_0+2),    j≠ s ≠ t
    𝔼[x_j^2 x_s]    = α_j(α_j + 1) α_s/α_0(α_0+1)(α_0+2),    j≠ s
         𝔼[x_j^3]    = α_j(α_j + 1) (α_j + 2)/α_0(α_0+1)(α_0+2),

where α_0 = ∑_i≥ 1α_i. Comparing the second and third order component moments, we arrive at the following cross-moments:

    𝔼[x]    = 1/α_0α
    𝔼[x×x]    = 1/α_0(α_0 + 1)αα^⊤ + 1/α_0(α_0 + 1)diag(α) = α_0/α_0+1𝔼[x] 𝔼[x]^⊤ + diag(α)/α_0(α_0+1)
    𝔼[x×x×x]    = 
    1/α_0 (α_0 + 1) (α_0 + 2) [
    α×α×α +
    ∑_i≥ 1α_i (α×e_i ×e_i) +
    ∑_i≥ 1α_i (e_i ×α×e_i ) . 
       = +.
    ∑_i≥ 1α_i (e_i ×e_i ×α) + 2∑_i≥ 1α_i
    (e_i ×e_i ×e_i)
    ].

To express the parameters as third-order cross-moments, first observe that the following holds for a Dirichlet random vector x:

    1/α_0^2 (α_0 + 1)∑_i≥ 1α_i (α×e_i×e_i + e_i×α×e_i + e_i×e_i×α)
    = 𝔼[𝔼[x] ×x×x]
    + 𝔼[x×𝔼[x] ×x]
    + 𝔼[x×x×𝔼[x]] - 3α_0/α_0 + 1𝔼[x] ×𝔼[x] ×𝔼[x].

This is an immediate result from the the second-order component moments. Combining with (<ref>) yields

    ∑_i≥ 12α_i e_i ×e_i ×e_i/α_0 (α_0+1) (α_0+2)   = 𝔼[x×x×x]
    + 2α_0^2/(α_0 + 1)(α_0 + 2)𝔼[x] ×𝔼[x] ×𝔼[x] 
    - α_0/α_0 + 2(
     𝔼[𝔼[x] ×x×x]
    + 𝔼[x×𝔼[x] ×x]
    + 𝔼[x×x×𝔼[x]]
    ).




 §.§ Derivation of moment estimators

Our goal is to derive the estimators of parameter vectors θ_j for each variable j using the first- and second- order empirical cross-moments of b_ij. In GDLM, the expectation of variable j conditioned on x is written

    𝔼[b_ij| x] = θ_j x.

Thus, the expected observation of variable j is given by

    𝔼[b_ij] = 𝔼[𝔼[b_ij| x]] = θ_j 𝔼[x] = θ_j α/α_0.

Now consider two variables b_ij and b_is which are generated with the same latent factors x. Combining (<ref>) and (<ref>) to obtain

    ∑_r≥ 1α_r/α_0 (α_0 + 1)θ_jr×θ_sr 
    =
    𝔼[b_ij×b_is] - α_0/α_0 + 1𝔼[b_ij] 𝔼[b_is]^⊤.

For three variables b_ij, b_is, and b_it, we can write  𝔼[b_ij×b_is×b_it] = 𝔼[x×x×x] ×_1 θ_j ×_2 θ_s ×_3 θ_t. Using (<ref>), we establish that

    ∑_i≥ 12α_i θ_j ×θ_s ×θ_t/α_0 (α_0+1) (α_0+2)   = 𝔼[b_ij×b_is×b_it]
    + 2α_0^2/(α_0 + 1)(α_0 + 2)𝔼[b_ij] ×𝔼[b_is] ×𝔼[b_it] 
    - α_0/α_0 + 2(
     𝔼[𝔼[b_ij] ×b_is×b_it]
    + 𝔼[b_ij×𝔼[b_is] ×b_it]
    + 𝔼[b_ij×b_is×𝔼[b_it]]
    ).




§ APPROXIMATE ORTHOGONALIZATION IN THE TENSOR POWER METHOD


TPM requires the tensor to be decomposed to be symmetric, and the factor matrices to be orthogonal. Specifically, it performs the following decomposition

    ℳ_3^' = ∑_i=1^r λ_i u_i ×u_i ×u_i,

where u_i are orthonormal vectors. Thus, TPM does not immediately apply to the general CP decomposition (<ref>).

The general resolution is to first use the symmetric tensor embedding <cit.>, forming a larger symmetric tensor ℳ_3 that contains the asymmetric tensor to be decomposed. The formed ℳ_3 is a sparse (∑_i=1^p d_i)-by-(∑_i=1^p d_i)-by-(∑_i=1^p d_i) tensor of which 7/9 entries are zero. The space and computation complexities rapidly become prohibitive when the number of variables p and the category counts d_j grow. 

Next, TPM requires an addition empirical second-order estimator ℳ_2 for orthogonalizing the factor matrices of ℳ_3 to obtain ℳ_3^' <cit.>. This is done by computing the whitening transformation from ℳ_2. However, the whitening technique based on empirical ℳ_2 is often a cause of suboptimal performance <cit.>.

