Title: MULTI-LABEL LEARNING FOR LARGE TEXT COR-PORA USING LATENT VARIABLE MODEL WITH PROV-ABLE GUARANTEES

Abstract: Here we study the problem of learning labels for large text corpora where each document can be assigned a variable number of labels. The problem is trivial when the label dimensionality is small and can be easily solved by a series of one-vs-all classifiers. However, as the label dimensionality increases, the parameter space of such one-vs-all classifiers becomes extremely large and outstrips the memory. Here we propose a latent variable model to reduce the size of the parameter space, but still efficiently learn the labels. We learn the model using spectral learning and show how to extract the parameters using only three passes through the training dataset. Further, we analyse the sample complexity of our model using PAC learning theory and then demonstrate the performance of our algorithm on several benchmark datasets in comparison with existing algorithms.
Multi-label learning for large text corpora is an upcoming problem in Large-Scale Machine Learning. Unlike the multi-class classification where a document is assigned only one label from a set of labels, here a document can have a variable number of labels. A basic approach to the problem is to use 1-vs-all classification by training a single binary classifier for every label. If the vocabulary size of a text corpus is D and the label dimensionality is L, then the 1-vs-all model require O(DL) parameters. Most of the text corpora have moderate to high vocabulary size (D), and 1-vs-all classification is feasible as long as D L. However, as the number of labels increases to a point when L ∼ D, the size of the parameter space increases to O(D 2 ), and we can no longer store the 1-vs-all models in the memory Yu et al. (2014).
A fairly intuitive approach to reducing the size of the parameter space is to use low-rank models. The existing low-rank approaches in the literature are mainly based on discriminative models that use a mapping Φ : R D → R L between the word space and the label space. If the rank of the model is restricted to K D, then such mappings have the form Z = HW , where W ∈ R D×K and H ∈ R L×K contain the low-rank features of the words and the labels respectively. The parameter space of such models has the size Θ ((D + L)K), which is much smaller than O(DL) for large values of L and can be conveniently stored the memory.
Amongst the recent literature on discriminative low-rank models, WSABIE Weston et al. ( 2011) defines weighted approximate pair-wise rank (WARP) loss on the mapping and optimises the loss on the training dataset. LEML Yu et al. (2014) generalises the loss function to squared-loss, sigmoid loss or hinge loss, which are typical of Linear Regression, Logistic Regression, and Linear SVM respectively. Both of LEML and WSABIE use gradient descent to optimise the loss on training datasets and therefore are susceptible to overfitting. Also, both of them run several iterations through the dataset and can slow down for large datasets containing millions of text instances.
Here we propose a latent variable model for multi-label learning. Our model can be viewed as a generative counterpart of the discriminative low-rank models like WSABIE or LEML. Unlike the usual cases where such latent variable models are trained using EM algorithm, we use Spectral Learning or Method of Moments Anandkumar et al. (2014) to extract the parameters. Method of Moments is an upcoming non-iterative approach for latent variable based models, and has successfully been applied in Spherical Gaussian Model Hsu & Kakade (2013), Hidden Markov Models

:

)). We first derive a closed form expression to estimate the model parameters using the Method of Moments. Further, we derive the convergence bounds on the estimated parameters, and then compare the performance of our algorithm with the discriminative low-rank models both from the perspective of accuracy and computation time.

There are other approaches beyond low-rank mappings, such as FastXML Prabhu & Varma (2014) that uses tree-based hierarchies for the labels, or SLEEC Bhatia et al. (2015) that uses kNN based embedding of the texts in the training dataset. However, these approaches use a combination of different algorithms and are geared towards end-to-end solutions. Whereas, we look into the problem from a modelling perspective and limit our discussion to the low-rank models since they have a similar model complexity to the latent variable model.

2. METHODOLOGY:

We build our model based on the "bag of words" assumption on the documents taking into account only the occurrence of the words, not their counts. Let us assume that there are N labeled documents in a corpus with vocabulary size D and label dimensionality L. Then using a latent variable h that assumes K states, the probability of a label l being assigned to a document d can be expressed as:

EQUATION

where the conditional probability of a document d given h can be modelled from the set of distinct words W d present in it as:

EQUATION

From here, P [h = k|d] can be expressed using Bayes Theorem as:

EQUATION

The latent variable model is described by the plate notation in Figure 1 , where the words are denoted by v, the labels by l and the latent variable by h. Our model parameters are P [h], P [v|h] and P [l|h], and the number of parameters is Θ ((D + L)K), same as the discriminative low-rank models.

Let us define π k ∈ [0, 1] as the probability of the latent variable h assuming the state k ∈ [K] .

EQUATION

Let us define µ k ∈ [0, 1] D as the probability vector of a word v conditional to the latent state k ∈ [K] .

EQUATION

Similarly, let us define γ k ∈ [0, 1] L as the probability vector of a label l conditional to the latent state k ∈ [K] .

EQUATION

One key difference of our parameterization from Anandkumar et al. (2012) and Wang & Zhu (2014) is that they define {µ k } K k=1 as the expected count of the words (E[v|h = k]) for each topic k, while we define {µ k } K k=1 as the probability of the words. A probabilistic formulation enables us to model the documents using Equation 2, which is not possible using expectations. We later use the document probabilities to predict the labels for a test document (Section 4). Therefore, unlike Anandkumar et al. (2012) or Wang & Zhu (2014) that uses moments on the word counts defined on the domain p R D (for pth order moment), we use moments based on the joint probability of the words defined on the domain p [0, 1] D . We outline the derivation for the probability moments in this section, and then show that our formulation achieves a tighter convergence bound on the parameters in Section 3.1. The label parameters {γ k } K k=1 are unique to our model.

2.1. MOMENT FORMULATION:

We first try to formulate the joint probability of a pair of words in terms of π and µ. Let us assume that we choose two words w 1 and w 2 from a document at random. The term

P [w 1 = v i , w 2 = v j ]

represents the probability by which any pair of words picked at random from a document turns out to be v i and v j , and it is nothing but

P [v i , v j ], where i, j ∈ [D]. Similarly, P [w 1 = v i ] is same as P [v i ],

and the same holds for

P [w 2 = v j ].

Following the generative model in Figure 1 , any two words w 1 and w 2 in a document are conditionally independent given h. Therefore,

P [w 1 = v i , w 2 = v j ] = K k=1 P [w 1 = v i , w 2 = v j |h = k]P [h = k] = K k=1 P [w 1 = v i |h = k]P [w 2 = v j |h = k]P [h = k] = K k=1 P [v i |h = k]P [v j |h = k]P [h = k] = K k=1 π k µ ki µ kj

where µ ki is the ith element of the vector µ k and so on.

EQUATION

If M 2 ∈ [0, 1] D×D is the matrix containing the joint probability of the pairs of words with

M 2i,j = P [v i , v j ] ∀i, j ∈ [D]

, then it can be expressed as:

M 2 = K k=1 π k µ k µ k (8)

Similarly, if the tensor M 3 ∈ [0, 1] D×D×D is defined as the moment containing the joint probability of a triplet of words with M 3i,j,κ = P

[v i , v j , v κ ], ∀i, j, κ ∈ [D], then M 3 = K k=1 π k µ k ⊗ µ k ⊗ µ k (9)

where ⊗ stands for the tensor product.

Finally, we define M L ∈ [0, 1] L×D×D as the probability moment of a label occurring with a pair of words such that M Lκ,i,j

= P [l κ , v i , v j ], ∀κ ∈ [L] and ∀i, j ∈ [D]. Then, M L = K k=1 π k γ k ⊗ µ k ⊗ µ k (10) 2.2 CLOSED FORM EXPRESSION OF THE PARAMETERS

We derive the equations for extracting π, µ and γ in this section. The first step is to whiten the matrix M 2 , where we try to find a matrix low rank W such that W M 2 W = I. This is similar to the whitening in ICA Hyvärinen et al. (2004) , with the covariance matrix being replaced by the cooccurrence probability matrix. The whitening is usually done through singular value decomposition of M 2 . If the K maximum singular values of M 2 are {ν k } K k=1 , and {ω k } K k=1 are the corresponding left singular vectors, then the whitening matrix of rank K is computed as

W = ΩΣ -1/2 , where Ω = ω 1 |ω 2 | . . . |ω K and Σ = diag (ν 1 , ν 2 , . . . , ν K ).

Upon whitening M 2 takes the form

W M 2 W = W K k=1 π k µ k µ k W = K k=1 √ π k W µ k √ π k W µ k = K k=1 μk μ k = I (11)

Hence μk = √ π k W µ k are orthonormal vectors. Multiplying M 3 along all three dimensions by W , we get

M3 = M 3 (W, W, W ) = K k=1 π k (W µ k ) ⊗ (W µ k ) ⊗ (W µ k ) = K k=1 1 √ π k μk ⊗ μk ⊗ μk (12)

M3 is a tensor in the domain R K×K×K . Upon the factorisation of M3 , if the eigenvalues and eigenvectors are {λ k } K k=1 and {u k } K k=1 respectively, then

λ k = 1 √ π k =⇒ π k = λ -2

k , and

u k = μk = √ π k W µ k = 1 λ k W µ k (13)

Hence, {µ k } K k=1 can be recovered as

µ k = λ k W † u k , where W † = W W W -1 is the pseudo- inverse of W . Since we compute P [v|h] by normalising µ k as P [v|h = k] = µ kv v µ kv , it suffices to compute µ k = W † u k as λ k will be cancelled during normali- sation.

It is possible to compute γ through the factorisation of second and third order moments of the labels in a similar way to µ. However, there is no guarantee that the label probabilities of γ k will correspond to the same topic as the word probabilities in µ k for every k when we compute them separately. Without the topic alignment between µ and γ, label prediction is not possible.

To overcome this limitation, we use the cross-moment M L to compute γ. If we multiply M L twice by W , then

M L (W, W ) = K k=1 π k γ k ⊗ W µ k ⊗ W µ k = K k=1 γ k ⊗ √ π k W µ k ⊗ √ π k W µ k = K k=1 γ k ⊗ μk ⊗ μk (14)

If u k is the kth eigenvector of M3 , then

u k M L (W, W )u k = μ k M L (W, W )μ k = μ k K k=1 γ k ⊗ μk ⊗ μk μk = γ k (15)

since μk s are orthonormal. Computing γ k by this method ensures that γ k and μk correspond to the same topic for every k, which in turn ensures that γ k and µ k correspond to the same topic.

3. PARAMETER ESTIMATION:

So far we have shown how to extract the parameters from the global moments M 2 , M 3 and M L defined on population. In practice, we cannot compute the population parameters; all we can do is to estimate them from the sample in hand and compute an error bound for the estimation. We denote the estimated values of M 2 , M 3 , M L , π, µ, γ, W and U from the sample as M2 , M3 , ML , π, μ, γ, Ŵ and Û , conforming with the notations used in the previous literature.

If we take into account only the occurrence of each distinct word in a document, then using "one hot encoding," we can represent the words as a binary vector x ∈ {0, 1} D , and the labels as another y ∈ {0, 1} L . If there are N samples, then we can represent the words in the entire corpus as X ∈ {0, 1} N ×D , and the associated labels as Y ∈ {0, 1} N ×L . The pairwise counts of the words can be estimated by X X, whose sum of all elements is,

D v1=1 D v2=1 (X X) v1,v2 = D v1=1 D v2=1 N i=1 x i,v1 x i,v2 = N i=1 D v1=1 D v2=1 x i,v1 x i,v2 = N i=1 nnz(x i ) 2

where x i is the row of X corresponding to the ith document, x i,v is the vth element in it (v ∈ [D]), and nnz(x i ) is the number of non-zero elements in x i , i.e., the number of distinct words in the ith document.

Therefore, the joint probability matrix of the words (M 2 ) can be estimated as,

EQUATION

Similarly, the triple-wise occurrences of the words can be estimated by the tensor X ⊗ X ⊗ X, and the sum of all of the elements of this tensor is

D v1=1 D v2=1 D v3=1 (X ⊗ X ⊗ X) v1,v2,v3 = N i=1 nnz(x i ) 3

Therefore, M 3 can be estimated as

EQUATION

The dimensions of M2 and M3 are D 2 and D 3 respectively. But in practice, these quantities are extremely sparse. We do not need to compute M3 in practice. If the whitening matrix computed through the Singular Value Decomposition of M2 is Ŵ , then M3 can be estimated straight away using Equation 12 as, M3 = 1

EQUATION

Since M3 has a dimension of K 3 and K D, it can be conveniently stored in the memory. The computation of M2 takes one pass through the dataset, and the computation of M3 takes another. Also, the counts of the labels occurring with the pairs of words can be estimated by the tensor Y ⊗ X ⊗ X, whose sum of all element is

L l=1 D v1=1 D v2=1 (Y ⊗ X ⊗ X) l,v1,v2 = N i=1 nnz(y i )nnz(x i ) 2

where y i represents the ith row of Y and nnz(y i ) is the number of labels associated with the ith document. Therefore, M L can be estimated as, 

ML = 1 N i=1 nnz(y i )nnz(x i ) 2 Y ⊗ X ⊗ X (19)

Ŵ = ΩΣ -1/2 ∈ R D×K 4. Estimate M3 = (X Ŵ ⊗ X Ŵ ⊗ X Ŵ )/ N i=1 nnz(x i ) 3 (pass # 2) 5. Compute eigenvalues { λk } K k=1 and eigenvectors {û k } K k=1 of M3 . Assign Û = [û 1 |û 2 | . . . |û K ] 6. Estimate μk = Ŵ † ûk , where Ŵ † = Ŵ ( Ŵ Ŵ ) -1 , and πk = λ-2 k , ∀k ∈ 1, 2 . . . K. 7. Estimate P [v|h = k] = μkv v μkv , ∀k ∈ 1 . . . K, v ∈ v 1 . . . v D 8. Estimate Γ = Y (X Ŵ Û ) • (X Ŵ Û )

, where • stands for the element-wise product (pass # 3) 9. Estimate

P [l|h = k] = Γlk l Γlk , ∀k ∈ 1 . . . K, l ∈ l 1 . . . l L If kth eigenvector of M3 is ûk , then from Equation 15 γk = û k ML ( Ŵ , Ŵ )û k = 1 N i=1 nnz(y i )nnz(x i ) 2 û k Y ⊗ X Ŵ ⊗ X Ŵ ûk = 1 Z(X, Y ) Y ⊗ X Ŵ ûk ⊗ X Ŵ ûk = 1 Z(X, Y ) Y (X Ŵ ûk ) • (X Ŵ ûk )

where • stands for the element-wise product, and

EQUATION

where

Û = [û 1 |û 2 | . . . |û K ]. Since we estimate P [l|h] by normalising the columns of Γ as P [l|h = k] = Γlk l Γlk , it is sufficient to compute Γ = Y (X Ŵ Û ) • (X Ŵ Û )

as the constant Z(X, Y ) will be cancelled out during normalisation. Instead of computing all the γk s separately, we can compute Γ in one step by just one pass through the entire dataset. The overall method takes three passes through the dataset as outlined in Algorithm 1. There is no need to compute ML explicitly since Γ can be computed straight away from X, Y and the intermediate parameters.

3.1. CONVERGENCE BOUNDS:

Theorem 1. Let us assume that we run Algorithm 1 on N i.i.d. samples with word vectors n 1 , n 2 , n 3 ), where   1 .

x 1 , x 2 . . . x N and label vectors y 1 , y 2 . . . y N . Let us define ε 1 = 1 + log(1/δ) 2 and ε 2 = 1 + log(2/δ) 2 for some δ ∈ (0, 1). If N ≥ max(

n 1 = c 2 log K + log log K c1 • πmax πmin 2. n 2 = Ω ε1 d2sσK (M2) 2 3. n 3 = Ω K 2 10 d2sσK (M2) 5/2 + 2 √ 2 d3sσK (M2) 3/2 2 ε 2 1

for some constants c 1 and c 2 , then following bounds on the estimated parameters hold with probability at least 1δ,

||µ k -μk || ≤ 160 σ 1 (M 2 ) d2s σ K (M 2 ) 5/2 + 32 2σ 1 (M 2 ) d3s σ K (M 2 ) 3/2 + 4 σ 1 (M 2 ) d2s σ K (M 2 ) ε 1 √ N ||γ k -γk || ≤ 160 d2s σ K (M 2 ) 7/2 + 32 √ 2 d3s σ K (M 2 ) 5/2 + 2 + √ 2 d2s σ K (M 2 ) 2 2ε 1 √ N + 8ε 2 dls σ K (M 2 ) √ N |π k -πk | ≤ 200 σ K (M 2 ) 5/2 + 40 √ 2 σ K (M 2 ) 3/2 ε 1 d3s √ N

The terms σ 1 (M 2 ) . . . σ K (M 2 ) are the K largest eigenvalues of the matrix M 2 , whereas d2s = E nnz(x) 2 , d3s = E nnz(x) 3 and dls = E nnz(y)nnz(x) 2 , with nnz(x) representing the number of distinct words and nnz(y) representing the number of labels present in an instance. The proof is included in the appendix.

The convergence bounds on π and µ are very similar to those in Anandkumar et al. (2012) and Wang & Zhu (2014) , except for the terms d2s , d3s and dls in the denominator. These terms arise in our convergence bounds since we use the probability moments, whereas Anandkumar et al. (2012) and Wang & Zhu (2014) use the moments of word counts. The parameter γ is unique to our model, although its bound contains those terms too.

We need at least one document with 3 distinct words or more to construct the third order moment M 3 , i.e., nnz(x) ≥ 3 for at least one document. Therefore, d2s = E nnz(x) 2 > 1 and d3s = E nnz(x) 3 > 1. Also, since every document has least one label, nnz(y) ≥ 1 for all the documents, and therefore, and dls > 1. In practice, these terms can be much larger than 1 for any real-life text corpus (The estimates of them for different experimental datasets are in Table 1 ). Therefore, our algorithm achieves a tighter convergence bound than Anandkumar et al. (2012) or Wang & Zhu (2014) . Further, since the terms d2s and d3s appear in the denominators of n 2 and n 3 too, using probability moments also lowers the minimum number of samples required.

3.2. COMPUTATIONAL COMPLEXITY:

The bottleneck of the algorithm is the whitening step of M2 , especially for the large datasets. M2 computed using Equation 16 is a symmetric p.s.d. matrix with very high sparsity. The number of elements of M2 is O N i=1 nnz(x i ) 2 , with the worst case occurring when no two documents have any word in common, and every element in X X is 1. The top K eigenvalues for large symmetric matrices are usually computed through Lanczos method. For a sparse matrix of size

D × D with O N i=1 nnz(x i ) 2 number of elements, this step has a complexity of O N i=1 nnz(x i ) 2 K + DK 2 Mahadevan (2008).

We do not compute the third order tensors M3 or ML explicitly. The step to compute M3 using Equation 18 has a complexity of O(N K 3 ), while the computation of ΓL using Equation 20 has a complexity of O N i=1 nnz(y i ) K 2 . We carry out the tensor factorisation of M3 ∈ R K×K×K using the Tensor Toolbox Bader et al. (2015) , and this step has a complexity of O(K 4 log ) to compute each of the K eigenvalues up to an accuracy of Kolda & Mayo (2011) . These steps contribute the most to the computational cost. The overall time complexity of Algorithm 1 is

O N i=1 nnz(x i ) 2 K + D + N i=1 nnz(y i ) K 2 + N K 3 + K 4 log .

As per as space requirement is concerned, the storage of M2 takes up O N i=1 nnz(x i ) 2 space, whereas the word parameters like Ŵ and μ1:K take up a space of Θ(DK), the label parameter Γ takes up a space of Θ(LK), and M3 takes up a space of Θ(K 3 ). The overall space complexity is 

O N i=1 nnz(x i ) 2 + (D + L)K + K 3 .

4. LABEL PREDICTION:

P [d|h = k] = v∈W d P [v|h = k] v / ∈W d 1 -P [v|h = k]

This step will take D multiplications involving floating point operations. However, since 

v P [v|h = k] = 1 for v ∈ {v 1 . . . v D }

P [d|h = k] = v∈W d P [v|h = k] v / ∈W d P [v|h=k]>η 1 -P [v|h = k]

Usually only a few of P [v|h = k] will be greater than the threshold, and this will significantly reduce the number of multiplications when v / ∈ W d without sacrificing accuracy. Also, since a document has only a few distinct words, i.e. |W d | D, the left multiplicand (when v ∈ W d ) accounts for only a few multiplications. From there, we can estimate the following using Bayes rule.

P [h = k|d] = P [h = k] P [d|h = k] K k=1 P [h = k] P [d|h = k] = πk v∈W d P [v|h = k] v / ∈W d P [v|h=k]>η 1 -P [v|h = k] K k=1 πk v∈W d P [v|h = k] v / ∈W d P [v|h=k]>η 1 -P [v|h = k]

And finally, the label probabilities for the document ( P [l|d]) can be estimated from P [l|h] and P [h|d] using Equation 1.

5. EXPERIMENTAL RESULTS:

We carry out our experiments on six datasets ranging from small datasets like Bibtex with 4, 880 training instances to large datasets like WikiLSHTC with more than 1.7M training instances. The description of the datasets are listed in Table 1 . We categorise the datasets into three groups:

1. Small: Bibtex and Delicious "airlines and airplanes" (0.34), "terrorism" (0.30), "united states international relations" (0.27), "elections" (0.22), "armament, defense and military forces" (0.18), "internationalrelations" (0.18), "bombs and explosives" (0.15), "murders and attempted murders " (0.13), "biographical information" (0.13), "islam" (0.12)

"terrorism" (0.12), "united states international relations" (0.08), "airlines and airplanes" (0.07), "world trade center (nyc)" (0.07), "hijacking" (0.07), "united states armament and defense" (0.07), "pentagon building" (0.03), "bombs and explosives" (0.03), "islam" (0.02), "missing persons " (0.02) "armament, defense and military forces", "civil war and guerrilla warfare", "politics and government"

"civil war and guerrilla warfare" (0.62), "united states international relations" (0.39), "united states armament and defense" (0.23), "armament, defense and military forces" (0.23), "internationalrelations" (0.17), "oil (petroleum) and gasoline" (0.11), "surveys and series" (0.10),"military action" (0.09), "foreign aid" (0.08), "independence movements" (0.08)

"civil war and guerrilla warfare" (0.09), "united states international relations" (0.09), "united states armament and defense" (0.06), "politics and government" (0.04), "armament, defense and military forces" (0.03), "internationalrelations" (0.02), "immigration and refugees" (0.02), "foreign aid" (0.02), "terrorism" (0.02), "economic conditions and trends" (0.02) Table 2 : Examples of label prediction from the NYTimes dataset. The numbers in parenthesis are the scores for the top 10 labels. The scores of LEML and MoM have different ranges.

Since LEML is shown to outperform WSABIE and other benchmark algorithms on various small and large-scale datasets in Yu et al. (2014) , we benchmark the performance of our method against LEML. Also, there are other methods proposed on topic based embedding of the labels, most recent of which is Rai et al. (2015) that extends Latent Dirichlet Allocation Blei et al. (2003) to Multi-label learning, and uses Bayesian Learning through Gibbs Sampling. However, use of Gibbs sampler limits the use of the algorithms only to the datasets of limited size. The largest dataset used in Rai et al. (2015) is EurLex that is similar to NYTimes dataset in size. Using any MCMC based sampling scheme is not viable for the large datasets such as AmazonCat or WikiLSHTC containing millions of training instances.

In our experiments, we chose to measure AUC (of Receiver Operator Characteristics) against the latent dimensionality (K). AUC is a versatile measure and has been used in a range of problems from binary classification to ranking. Also, it is shown that there exists a one-to-one relation between AUC and Precision-Recall curves in Davis & Goadrich (2006) , i.e., the same algorithm achieving a higher AUC will also produce a better Precision-Recall curve. We carried out our experiments on Unix Platform on a single machine with Intel i5 Processor (2.4GHz) and 16GB memory, and no multi-threading or any other performance enhancement method is used in the code. 1 For the label prediction step of MoM, we chose η = 0.001. For LEML, we ran ten iterations for the small datasets and five iterations for the medium and large datasets, since the authors of LEML chose a similar number of iterations for their experiments in Yu et al. (2014) . For AmazonCat and WikiLSHTC datasets, we ran LEML on an i2.4xlarge instance of Amazon EC2 with 122 GB of memory since LEML needs significantly larger memory for these two datasets (Figure 2 ). N . However, when we compute the AUC on test dataset, the AUC of LEML decreases with latent dimensionality(K) for some datasets, including the larger dataset of AmazonCat containing more than 1M training instance. This shows the possibility of over-fitting in LEML. MoM, on the other hand, is not an optimization algorithm, and the parameters are extracted from Moment Factorization rather than optimizing any target function. It is not susceptible to over-fitting, which is evident from its performance. On the other hand, MoM has the requirement N ⌦(K 2 ) on the number of documents in the training set, and it will not work if N < ⇥(K 2 ). However, for smaller text corpora where N < ⇥(K 2 ) hold, 1-vs-all classifiers are usually sufficient to predict the labels. We need dimensionality reduction techniques for large text corpora where 1-vs-all classifiers fail, and MoM provides a very competitive choice for such cases. We computed AUC for every test documents and performed a macro-averaging across the documents, and repeated the experiments for K = {50, 75, 100, 125, 150} (Figure 2 ). Both LEML and Method of Moments perform very similarly, but the memory footprint of MoM is significantly less than LEML. MoM takes longer to finish for the small datasets since tensor factorisation takes much more time compared to the LEML iterations. However, as the size of the datasets grows, the LEML iterations become more and more costly. For the medium and large datasets, MoM takes a fraction of the time taken by LEML. For WikiLSHTC, LEML takes more than two days to finish, while MoM finishes within a few hours. The training times of LEML and MoM for different datasets are listed in Table 3 .

6. CONCLUSION:

Here we propose a latent variable model for multi-label learning in large-scale text corpus and show how to extract the model parameters based on the spectral decomposition of the probability moments of the words and the labels. Our method gives similar performance in comparison with state-of-art algorithms like LEML while taking a fraction of time and memory for the medium and the large datasets. Our method takes only three passes through the training dataset to extract all the parameters, which contributes to its superior time performance. Also, the memory requirement of our method is nominal when compared to the existing algorithms, and it scales to large datasets like WikiLSHTC containing millions of Wikipedia articles just on a single node with 16GB of memory. Lastly, since our method consists of only linear algebraic operations, it is embarrassingly parallel and can easily be scaled up in any parallel ecosystem using linear algebra libraries. In our implementation, we used Matlab's linear algebra library based on LAPACK/ARPACK, although we did not incorporate any parallelisation.

LEML has a bound of O 1/ √ N on the training loss Yu et al. (2014) , whereas our method (MoM) has a convergence bound of O 1/ √ N on the model parameters w.r.t. the number of training instances N . However, when we compute AUC on the test dataset, the AUC of LEML decreases with the latent dimensionality(K) for three datasets including AmazonCat with more than a million of training instances. It shows the possibility of overfitting in LEML. MoM, on the other hand, is not an optimisation algorithm, and the parameters are extracted through spectral decomposition of matrices and tensors, rather than optimising any target function. It is not susceptible to over-fitting, which is evident from its performance. But MoM has the requirement of N ≥ Ω(K 2 ), and will not work when N < Θ(K 2 ). However, for smaller text corpora where N < Θ(K 2 ) holds, 1-vs-all classifiers are usually sufficient to predict the labels. We need low-rank models for large text corpora where 1-vs-all classifiers fail, and MoM provides a very competitive choice for such cases.

A MATRIX NORM INEQUALITIES:

Here we paraphrase the inequalities regarding the matrix norms from Anandkumar et al. (2014) , which we will use further in our final proof. Let the true pairwise probability matrix and the third order probability moment be M 2 = p(v, v) and M 3 = p (v, v, v) , where v stands for the words. Let us assume that we select N i.i.d. samples x 1 , . . . x N from the population, and the estimates of pairwise matrix and third order moment are M2 = p(y, y) and M3 = p(y, y, y). Let ε M2 = ||M 2 -M2 || 2 . We use the second order operator norm of the matrices here. Let us assume ε M2 ≤ σ K (M 2 )/2, where σ K is the Kth largest eigenvalue of M 2 . We will derive the conditions which satisfies this later.

If Σ = diag(σ 1 , σ 2 . . . σ K ) are the top-K eigenvalues of M 2 , and U are the corresponding eigenvectors, then the whitening matrix

EQUATION

Let Ŵ be the whitening matrix for M2 , i.e., Ŵ M2 Ŵ = I K×K . Then by Weyl's inequality,

σ k (M 2 ) -σ k ( M2 ) ≤ ||M 2 -M2 ||, ∀k = 1, 2 . . . K.

Therefore,

EQUATION

Also, by Weyl's Theorem,

EQUATION

Let U * be the eigenvectors of Ŵ M 2 Ŵ , and Λ be the corresponding eigenvalues. Then we can write,

EQUATION

EQUATION

EQUATION

B TENSOR NORM INEQUALITIES:

Let us define the second order operator norm of a tensor T ∈ R D ×D×D with D, D ∈ Z + as,

||T || 2 = sup u {||T (•, u, u)|| : u ∈ R D &||u|| = 1} (27) Lemma 1. For a tensor T ∈ R D ×D×D , ||T || 2 ≤ ||T || F , where ||T || F is the Frobenius norm defined as, ||T || F = i,j,k (T i,j,k ) 2

Proof. The tensor T can be unfolded as an array of D matrices each of size

EQUATION

Let us assume that the singular value decomposition of T d has the form (d ∈ [D ]),

T d = σ 1 ψ 1 φ 1 + σ 2 ψ 2 φ 2 + • • • + σ D ψ D φ D

where {ψ i } D i=1 are the left singular vectors, {φ i } D i=1 are the right singular vectors and {σ i } D i=1 are the singular values.

Each of [ψ 1 , ψ 2 . . . ψ D ] and [φ 1 , φ 2 . . . φ D ] forms a spanning set in R D , and any vector u ∈ R D can be spanned using both ψ 1:D and φ 1:D as the basis.

Let u = α 1 ψ 1 + α 2 ψ 2 + • • • + α D ψ D = β 1 φ 1 + β 2 φ 2 + • • • + β D φ D , where α 1:D and β 1:D are scalars. Since ||u|| = √ < u, u > = 1, ||α|| = α 2 1 + α 2 2 + . . . α 2 D = 1 and ||β|| = β 2 1 + β 2 2 + . . . β 2 D = 1 Therefore, |u T d u| = |α 1 σ 1 β 1 + α 2 σ 2 β 2 . . . α D σ D β D | ≤ σ 2 1 + σ 2 2 + • • • + σ 2 D α 2 1 β 2 1 + α 2 2 β 2 2 . . . α 2 D β 2 D (Cauchy-Schwartz Inequality) ≤ σ 2 1 + σ 2 2 + • • • + σ 2 D α 2 1 + α 2 2 • • • + α 2 D β 2 1 + β 2 2 • • • + β 2 D = σ 2 1 + σ 2 2 + • • • + σ 2 D = ||T d || F

Equality holds when T d is of rank 1. Since this holds for any vector u ∈ R D such that ||u|| = 1, sup

||u||=1 |u T d u| ≤ ||T d || F for d ∈ [D ]

. Therefore, from Equation 28,

||T || 2 2 ≤ ||T 1 || 2 F + |T 2 || 2 F + • • • + ||T D || 2 F = ||T || 2 F =⇒ ||T || 2 ≤ ||T || F Lemma 2.

(Robust Power Method from (Anandkumar et al., 2014 )

) If T = T + E ∈ R K×K×K ,

where T is an symmetric tensor with orthogonal decomposition

T = K k=1 λ k u k ⊗ u k ⊗ u k with each λ k > 0, and E has operator norm ||E|| 2 ≤ . Let λ min = min K k=1 {λ k } and λ max = max K k=1 {λ k }.

Let there exist constants c 1 , c 2 such that ≤ c 1 • (λ min /K), and N ≥ c 2 (log K + log log (λ max / )). Then if Algorithm 1 in (Anandkumar et al., 2014) is called for K times, with L = poly(K) log(1/η) restarts each time for some η ∈ (0, 1), then with probability at least 1η, there exists a permutation Π on [K], such that,

||u Π(k) -ûk || ≤ 8 λ Π(k) , |λ k -λ Π(k) | ≤ 5 ∀k ∈ [K] (29) Since ≤ c 1 • (λ min /K) and λ k = 1 √ π k , ∀k ∈ [K], we need N ≥ c 2 log K + log log Kλ max c 1 λ min = c 2 log K + log log K c 1 π max π min (30)

This contributes in the first lower bound (n 1 ) of N in Theorem 1.

C TAIL INEQUALITY:

Lemma 3. If we draw N i.i.d. sample documents x 1 , x 2 . . . x N , and probability mass function, pairwise probability and triplet probability of the words v estimated from these N samples are p(v), p(v, v) and p(v, v, v) respectively, whereas the true probabilities are p(v), p (v, v) and p(v, v, v) respectively with v ∈ {v 1 , v 2 . . . v D } , then with probability at least 1δ with δ ∈ (0, 1),

EQUATION

where, d1s = E [nnz(x)], d2s = E nnz(x) 2 , d3s = E nnz(x) 3 , and nnz(x) is the non-zero entries in the binary vector representing the words in the documents as described in Section 3 (main paper).

Proof. Since the number of distinct words in a document is always bounded, we can assume ||x|| ≤ 1 ∀x without loss of generality. Then from Lemma 7 of supplementary material of Wang & Zhu (2014) , with probability at least 1δ with δ ∈ (0, 1),

EQUATION

where E stands for true expectation, and Ê stands for the expectation estimated from the N samples, i.e.,

Ê[x] = 1 N N i=1 x i = 1 N X 1 Ê[xx ] = 1 N N i=1 x i ⊗ x i = 1 N X X Ê[x ⊗ x ⊗ x] = 1 N N i=1 x i ⊗ x i ⊗ x i = 1 N X ⊗ X ⊗ X

Now, since the samples contains binary data, probability of the items can be computed from as

EQUATION

Similarly,

EQUATION

Now, since x is a binary vector of dimension D, the sum across its dimensions is

D v=1 x v = nnz(x),

and therefore,

D v=1 E[x v ] = E[ D v=1 x v ] = E[nnz(x)] ≈ Ê[nnz(x)] = D v=1 Ê[x v ]. Assigning d1s = E[nnz(x)],

from Equation 37 and 38, we get

EQUATION

Similarly, summing across the respective dimensions, we get 

(x ⊗ x ⊗ x) v1,v2,v3 = nnz(x) 3 . Therefore, D v1=1 D v2=1 E[(xx ) v1,v2 ] = E[nnz(x) 2 ] = d2s ≈ D v1=1 D v2=1 Ê[(xx ) v1,v2 ],

and

D v1=1 D v2=1 D v3=1 E[(x ⊗ x ⊗ x) v1,v2,v3 ] = E[nnz(x) 3 ] = d3s ≈ D v1=1 D v2=1 D v3=1 Ê[(x ⊗ x ⊗ x) v1,v2,v3 ].

From here, we can show that,

EQUATION

Plugging these equations in Equation 34, 35 and 36, we can complete the proof.

Also, if y represents the label vector of the documents, since the number of labels per documents is limited, we can assume ||y|| ≤ 1 without a loss of generality. Then,

EQUATION

From Equation 34and 35,

P Ê[y] -E[y] F ≤ 2 √ N 1 + log(1/δ) 2 ≥ 1 -δ P Ê[xx ] -E[xx ] F ≤ 2 √ N 1 + log(1/δ) 2 ≥ 1 -δ

Reversing the tail inequalities,

P Ê[y] -E[y] F ≥ 2 √ N 1 + log(1/δ) 2 ≤ δ P Ê[xx ] -E[xx ] F ≥ 2 √ N 1 + log(1/δ) 2 ≤ δ

If E 1 and E 2 are two events,

P (E 1 ∩ E 2 ) = P (E 1 ) + P (E 2 ) -P (E 1 ∪ E 2 ) ≤ P (E 1 ) + P (E 2 )

Therefore,

EQUATION

Or,

P Ê[y] -E[y] F + Ê[xx ] -E[xx ] F ≤ 4 √ N 1 + log(1/δ) 2 ≥ 1 -2δ

From Equation 41,

EQUATION

where l stands for the labels (l ∈ {l 1 , l 2 . . . l L }). Therefore, with probability at least 1δ,

||p(l, v, v) -p(l, v, v)|| F ≤ 4 dls √ N 1 + log(2/δ) 2 (43) D COMPLETING THE PROOF Assigning ε 1 = 1 + log(1/δ) 2

in the inequalities of Lemma 3, we get

EQUATION

ε M3 = ||M 3 -M3 || 2 = ||p(v, v, v) -p(v, v, v)|| 2 ≤ ||p(v, v, v) -p(v, v, v)|| F ≤ 2ε 1 d3s √ N

since operator norm is smaller than Frobenius norm for both matrices and tensors (Lemma 1).

Therefore, to satisfy

ε M2 ≤ σ K (M 2 )/2, we need N ≥ Ω ε1 d2sσK (M2)

2 . This contributes in the second lower bound (n 2 ) of N in Theorem 1.

Similarly, assigning ε 2 = 1 + log(2/δ) 2 in Equation 43, we get

EQUATION

From Lemma 1, ||M 3 || 2 ≤ ||M 3 || F ≤ 1, because M 3 is a tensor with individual elements as probabilities, with the sum of all elements being 1.

From Lemma 2, ≤ c 1 • (λ min /K), and we can assign as the upper bound of ε tw . To satisfy this, we need Here, we will derive the convergence bounds for the parameters. Since µ k = W † u k (Algorithm 1 in main paper), with probability at least 1δ, 44, we can say that with probability at least 1δ,

10 d2s σ K (M 2 ) 5/2 + 2 √ 2 d3s σ K (M 2 ) 3/2 2ε √ N ≤ c 1 λ min K , or, 10 d2s σ K (M 2 ) 5/2 + 2 √ 2 d3s σ K (M 2 ) 3/2 2ε √ N ≤ c 1 1 K √ π max

EQUATION

EQUATION

since |λ k -λk | ≤ 5 from Lemma 2. Therefore, with probability at least 1δ, we get

EQUATION

where ε = 1 + log(1/δ) 2 .

Further, since γ = u k M L (W, W )u k ,

EQUATION

since u k = ûk = 1. Also, from Lemma 1, M L 2 ≤ M L F ≤ 1, since M L 2 is a tensor with the sum of its elements as 1.

Now,:

M L (W, W ) -ML ( Ŵ , Ŵ ) 2 = M L (W, W ) -M L ( Ŵ , Ŵ ) + M L ( Ŵ , Ŵ ) -ML ( Ŵ , Ŵ ) 2 ≤ M L (W, W ) -M L ( Ŵ , Ŵ ) 2 + Ŵ 2 M L -ML 2 = M L (W, W ) -M L (W, Ŵ ) + M L (W, Ŵ ) -M L ( Ŵ , Ŵ ) 2 + Ŵ 2 M L -ML 2 ≤ M L (W, W ) -M L (W, Ŵ ) 2 + M L (W, Ŵ ) -M L ( Ŵ , Ŵ ) 2 + Ŵ 2 M L -ML 2 ≤ W W -Ŵ M L 2 + Ŵ W -Ŵ M L 2 + Ŵ 2 M L -ML 2 = ( W + Ŵ ) M L 2 ε W + Ŵ 2 ε M L ≤ ( W + Ŵ )ε W + Ŵ 2 ε M L (50)

Therefore,

γ k -γk ≤ 2 u k -ûk W 2 + ( W + Ŵ )ε W + Ŵ 2 ε M L ≤ 16 λ k W 2 + ( W + Ŵ )ε W + Ŵ 2 ε M L ≤ 16 σ K (M 2 ) + ( W + Ŵ )ε W + Ŵ 2 ε M L Since 1/λ k = √ π k ≤ 1.

Assigning as the upper limit of ε tw in Equation 44, with probability

1 -δ, γ k -γk ≤ 16 σ K (M 2 ) 10 d2s σ K (M 2 ) 5/2 + 2 √ 2 d3s σ K (M 2 ) 3/2 2ε 1 √ N + 1 + √ 2 σ K (M 2 ) 2 σ K (M 2 ) 3/2 2ε 1 √ N d2s + 2 σ K (M 2 ) 4ε 2 √ N dls ≤ 80 d2s σ K (M 2 ) 7/2 + 16 √ 2 d3s σ K (M 2 ) 5/2 + 1 + √ 2 d2s σ K (M 2 ) 2 4ε 1 √ N + 8ε 2 dls σ K (M 2 ) √ N

where, ε 1 = 1 + log(1/δ) 2 and ε 2 = 1 + log(2/δ)

2:

. This completes the proof of Theorem 1.

To be shared later

Figures and Tables:

Figure 1: Latent Variable Model

and D can be very large, the values of P [v|h = k] are much less than one for the most, and 1 -P [v|h = k] ≈ 1. If we threshold P [v|h = k] when v / ∈ W d using some threshold η (say η = 0.01 or 0.001), then P [d|h = k] takes the form:

AUC and Memory (GB) of WikiLSHTC Dataset

Figure 2: AUC and Memory vs Latent State Dimensionality(K) for LEML and MoM for Different Datasets

Figure 2: AUC and Memory vs. Latent Dimensionality (K)

we need N ≥ Ω K This contributes to n 3 in Theorem 1.

45) where ||u kûk || ≤ 8 /λ k from Lemma 2, and ||û k || = 1. Since 1/λ k = √ π k ≤ 1, assigning as the upper bound of ε tw in Equation

Algorithm 1 Three-pass Algorithm for Parameter ExtractionInput: Sparse Binary Data X ∈ {0, 1} N ×D , Labels Y ∈ {0, 1} N ×L and K ∈ Z +

Table content: <table><tr><td>Output: π, P [v|h] and P [l|h]</td></tr><tr><td>1. Estimate M2 = (X X)/ 2. Compute K maximum singular values of M2 as {ν k } K N i=1 nnz(x i ) 2 k=1 , and corresponding left singu-(pass # 1) lar vectors as {ω k } K k=1 . Define Ω = ω 1 |ω 2 | . . . |ω K , and Σ = diag (ν 1 , ν 2 , . . . , ν K ) 3. Estimate the whitening matrix</td></tr></table>

Description of the Datasets

Table content: <table><tr><td>Name Bibtex Delicious Wiki-31K NYTimes AmazonCat 1,186,239 306,782 # of Train Points # of Test Points 4,880 2,515 12,920 3,185 14,146 6,616 14,669 15,989 WikiLSHTC 1,778,351 587,084</td><td>Feature Dimen-sion 1,836 500 101,938 24,670 203,882 1,617,899 325,056 Label Dimen-sion 159 983 30,938 4,185 13,330</td><td>Average # of Features 68.67 18.29 669.05 373.91 71.09 42.15</td><td>Median # of Features 69 6 513 354 45 30</td><td>Average # of Labels 2.4 19.02 18.64 5.40 5.04 3.19</td><td>Median # of Labels 2 20 19 5 4 2</td><td>d2s 5,957 1,892 559K 175K 6,301 2,210</td><td>d3s 597K 363K 545M 95M 651K 135K</td><td>dls 15K 24K 11M 755K 27K 7,594</td></tr><tr><td>True Labels "airlines and airplanes", "hijacking", "terrorism"</td><td>LEML</td><td/><td/><td/><td>MoM</td><td/><td/><td/></tr></table>



Table content: <table><tr><td>: Training Time and Speed-up</td></tr><tr><td>only three passes through the training dataset to extract all the parameters. Since MoM consists of only linear algebraic operations, it is embarrassingly parallel, and can easily be scaled up in any parallel eco-system using linear algebra libraries. In our implementation, we used Matlab's linear algebra library based on LAPACK/ARPACK, although we did not incorporate any parallelization. Both LEML and MoM have error bound of O(1/ p N ) on training performance w.r.t. the number of training samples</td></tr></table>

Training Time

Table content: <table/>

Peer review (iclr2018):
 The paper addresses the problem of multi-label learning for text corpora and proposes to tackle the problem using tensor factorization methods. Some analysis and experimental results for the proposed algorithm are presented. QUALITY: I find the quality of the results in this paper rather low. The proposed probabilistic model is defined ambiguously. The authors then look at joint probability distributions of co-occurence of two and three words, which gives a matrix and a tensor, respectively. They propose to match these matrix and tensor to their sample estimates and refer to such procedure as the moment matching method, which it is not. They then apply a standard two step technique from the moment matching literature consisting of whitening and orthogonal tensor factorization. However, in their case this does not have much statistical meaning. Indeed, whitening of the covariance matrix is usually justified by the scaling unidentifiability of the problem. In their case, the mathematics works because of the orthogonal unidentifiability of the square root of a matrix. Furthermore, the proposed sample estimators do not actually estimate densities they are dealing with (see, e.g., Eq. (16) and (17)). Their theoretical analysis seems like a straightforward extension of the analysis by Anandkumar, et al. (2012, 2014), however, I find it difficult to assess this analysis due to numerous ambiguities in the problem formulation and method development. This justifies my statement in the beginning of the paragraph. CLARITY: The paper is not well written and, therefore, is difficult to assess. Many important details are omitted, the formulation of the model is self contradicting, the standard concepts and notations are sometimes abused, some statements are wrong. I provide some examples in the detailed comment below. ORIGINALITY AND SIGNIFICANCE: The idea to apply tensor factorization approaches to the multi-label learning is novel up to my knowledge and is a pro of the paper. However, I have problems to find other pros in this submission because the clarity is quite low and in the present form there is no novelty in the proposed procedure. Moreover, the authors claim to work with densities, but end up estimating other quantities, which are not guaranteed to have the desirable form. They also emphasize the fact that there is the simplex constraint on the estimated parameters, but this constraint is completely ignored by the algorithm and, in general, won*t be satisfied in practice. If think the authors should do some more work before this paper can be published. DETAILED COMMENTS: Since I am quite critical about the paper, I point out some examples of drawbacks or flaws of this paper: - The proposed model (Section 2) is not well defined. In particular, the description in Section 2 is not sufficient to understand the proposed model; the plate diagram in Figure 2 is not consistent with the text. It is not mentioned how at least some conditional distributions behave (e.g., tokens given labels or states). The diagram in Fig. 1 does not help since it isn*t consistent with the text (e.g. the elements of labels or states are not conditionally / independent). The model is very close to latent Dirichlet allocation by Blei, et al. (2003), but differences are not discussed. - The standard terminology is often abused. For example, the proposed approach is referred to as the method of moments when it is not. In Section 2.1, the authors aim to match joint distributions (not the moments) to their empirical approximations (which are also wrong; see below). The usage of tokes and documents is interchanged without any explanations. - The use of the whitening approach is not justified in their setting working with joint distributions of couples and triples and it has no statistical meaning. No explanation is provided. I would definitely not call this whitening. - In Section 2.2, the notation is not defined and is different from what is usually used in the literature. For example, Eq. (15) does not make much sense as is. One could guess from the context that they are talking about the eigenvectors of an orthogonal tensor as defined in, e.g. Anandkumar, et al. (2014). - In Section 3, the authors emphasize the fact that their parameters are constrained to the probability simplex, but this constraint is not ensured in the proposed algorithm (Alg. 1). - Importantly, the estimators of the matrix M_2 and tensor M_3 do not make much sense to me. For example, for estimating M_2 it would be reasonable to average over all word pairs, i.e. something like [M_2]_{ij} = 1/L .sum_{w_k .not = w_l} P(w_k = v_i, w_l = v_j), where L is the number of pairs. This is different from the expression in Eq. (16), which is just a rescaled non-central second moment. Similar issue is true for the order-3 estimator. - The factorization procedure does not ensure non-negativity of the obtained parameters and, therefore, the rescaling is not guaranteed to belong to the probability simplex. I could not find any explanations of this issue. - I explain good plots in the experimental section, potentially, by the fact that the authors do algorithmically something different from what they aim to do, because the estimators do not estimate the desired entities (i.e. are not consistent). The procedure looks to me quite similar to the procedure for LDA, hence the reasonable results. However, the authors do not justify their proposed method.