A general method for regularizing tensor decomposition methods via pseudo-data

TL;DR

This paper introduces a versatile regularization technique for tensor decomposition methods that uses pseudo-data to improve model inference, applicable across various likelihood models and regularization goals.

Contribution

A novel, general approach to regularize tensor decomposition methods via pseudo-data, applicable to any differentiable regularizer and likelihood model.

Findings

Improves inference accuracy on synthetic and real data.

Enables regularization for transfer learning, sparsity, and interpretability.

Enhances orthogonality of learned parameters.

Abstract

Tensor decomposition methods allow us to learn the parameters of latent variable models through decomposition of low-order moments of data. A significant limitation of these algorithms is that there exists no general method to regularize them, and in the past regularization has mostly been performed using bespoke modifications to the algorithms, tailored for the particular form of the desired regularizer. We present a general method of regularizing tensor decomposition methods which can be used for any likelihood model that is learnable using tensor decomposition methods and any differentiable regularization function by supplementing the training data with pseudo-data. The pseudo-data is optimized to balance two terms: being as close as possible to the true data and enforcing the desired regularization. On synthetic, semi-synthetic and real data, we demonstrate that our method can improve inference accuracy and regularize for a broad range of goals including transfer learning, sparsity, interpretability, and orthogonality of the learned parameters.

Tables

Table 1. Table 1: Top headings for learned topics — Headings which are close to each other on the hierarchy tree are color coded with the same color (black signifies no neighbors within the topic)

Weight	Initial topics	Weight	Final topics
0.0402	Female	0.4451	Male
0.0393	Risk Factors	0.4281	Female
0.0367	Male	0.0349	Middle Aged
0.0361	Humans	0.0175	Adult
0.0359	Aged	0.0166	Aged
0.0330	Middle Aged	0.0045	Double-Blind Method
0.0277	Cardiovascular Diseases	0.0028	Aged, 80 and over
0.0272	Clinical Trials as Topic	0.0025	Follow-Up Studies
0.0717	Clinical Trials as Topic	0.4056	Anticholesteremic Agents
0.0597	Animals	0.3108	Hydroxymethylglutaryl-CoA …
0.0526	Hydroxymethylglutaryl-CoA …	0.2532	Hypolipidemic Agents
0.0496	Humans	0.0002	Fenofibrate
0.0376	Coronary Disease	0.0002	Fibrinolysis
0.0368	Anticholesteremic Agents	0.0001	Platelet Membrane Glycoproteins
0.0319	Hypercholesterolemia	0.0001	Erythrocytes
0.0284	Hypolipidemic Agents	0.0001	Mutagenicity Tests
0.0559	Male	0.0483	Male
0.0539	Female	0.0477	Female
0.0537	Middle Aged	0.0476	Anticholesteremic Agents
0.0485	Anticholesteremic Agents	0.0445	Hydroxymethylglutaryl-CoA …
0.0432	Hypercholesterolemia	0.0333	Humans
0.0384	Adult	0.0319	Middle Aged
0.0352	Aged	0.0288	Hypolipidemic Agents
0.0338	Humans	0.0281	Clinical Trials as Topic

Equations

$$M_2$$

$$M_3$$

$$L(X_T,X_P,\lambda )=-\log p(X_P|A_T)+\lambda R(A_{T\cup P}),$$

$$R(A_{T\cup P})=\log p(A_{T\cup P}|{\sigma }_m^2)=-\frac{DK}{2}\log (2\pi {\sigma }_m^2)-\frac{1}{2{\sigma }_m^2}\sum _{d,k}{a}_{d,k}^2,$$

$$R(A_{T\cup P})=\sum _{i\ne j}{a}_i\cdot a_j,$$

$$R(A_{T\cup P})=-\sum _{k=1}^K(\sum _{i\ne j}{a}_{i,k}{a}_{j,k}{O}_{i,j}^{-1}),$$

$$R(A_{T\cup P})=-\log p(A_{T\cup P}|{\alpha }_A).$$

$$M_2$$

$$M_3$$

$$b_n$$

$$z_{in}|b_n$$

$$w_{in}|z_{in},A$$

$${\hat{M}}_1$$

$${\hat{M}}_2$$

$${\hat{M}}_3$$

$$-\frac{K{\alpha }_B}{K{\alpha }_B+2}(\mathbb{E}[e_{w_1}\otimes e_{w_2}\otimes {\hat{M}}_1]$$

$$+\mathbb{E}[e_{w_1}\otimes {\hat{M}}_1\otimes e_{w_2}]+\mathbb{E}[{\hat{M}}_1\otimes e_{w_1}\otimes e_{w_2}])$$

$$+\frac{2(K{\alpha }_B{)}^2}{(K{\alpha }_B+2)(K{\alpha }_B+1)}{\hat{M}}_1\otimes {\hat{M}}_1\otimes {\hat{M}}_1,$$

$$h_n$$

$$x_n|h_n,A$$

$${\hat{M}}_2$$

$${\hat{M}}_3$$

$$+e_i\otimes \mathbb{E}[x]\otimes e_i+e_i\otimes e_i\otimes \mathbb{E}[x]),$$

$${\hat{M}}_{i,T\cup P}$$

$$={\hat{M}}_{i,T}+\frac{N_P}{N_T+N_P}({\hat{M}}_{i,P}-{\hat{M}}_{i,T}),$$

$$\underset{N_T\to \infty }{\lim }{\hat{M}}_{i,T\cup P}=\underset{N_T\to \infty }{\lim }\left[{\hat{M}}_{i,T}+\mathcal{O}\left(\frac{N_P}{N_T}\right)\right]={\hat{M}}_{i,T}.$$

$$\Vert {\hat{M}}_{i,P}-{\hat{M}}_{i,T}{\Vert }_2$$

$$\le \Vert {\hat{M}}_{i,T}-M_i{\Vert }_2+\Vert {\hat{M}}_{i,P}{\Vert }_2+\Vert M_i{\Vert }_2$$

$$\le ϵ+\gamma +\Vert M_i{\Vert }_2\stackrel{\text{def}}{=}{B}_{M_i}.$$

$$R(A)=-\sum _{k=1}^K(\sum _{i\ne j}{A}_{ik}{A}_{jk}{O}_{ij}^{-1}).$$

$$R(A)=-\frac{1}{2}\text{tr}(A^T{O}^{\ast }A-A^{T}A),$$

Full text

1 Introduction

Tensor decomposition methods (TDMs) have recently gained popularity as ways of performing inference for latent variable models (Anandkumar et al., 2014). The interest in these methods is motivated by the fact that they come with theoretical global convergence guarantees in the limit of infinite data (Anandkumar et al., 2012; Arora et al., 2013). However, a main limitation of these methods is that they lack natural methods for regularization or encouraging desired properties on the model parameters when the amount of data is limited.

Previous works attempted to alleviate this drawback by modifying existing tensor decomposition methods to incorporate specific constraints, such as sparsity (Sun et al., 2015), or incorporate modeling assumptions, such as the existence of anchor words (Arora et al., 2013; Nguyen et al., 2014). All of these works develop bespoke algorithms tailored to those constraints or assumptions. Furthermore, many of these methods impose hard constraints on the learned model, which may be detrimental as the size of the data grow—framed in the context of Bayesian intuition, when we have a lot of data, we want our methods to allow the evidence to overwhelm our priors.

We introduce an alternative approach which can be applied to encourage any (differentiable) desired structure or properties on the model parameters, and which will only encourage this “prior” information when the data is insufficient. Specifically, we adopt the common view of Bayesian priors as representing “pseudo-observations” of artificial data which bias our learned model parameters towards our prior belief (Bishop, 2006). We apply the tensor decomposition method of Anandkumar et al. (2014) to data sets comprised of both the actual data and an artificial pseudo-data. Gradient descent and automatic differentiation (Baydin et al., 2015; Maclaurin et al., 2015) can then be used to optimize our pseudo-data such that they maximize the desired properties on the inferred model parameters while still remaining as similar as possible to the actual data.

The resulting algorithm provides a method of imposing any regularizer on the standard TDMs and can be applied to any likelihood model which is learnable using TDMs. We provide theoretical analysis and prove that in the limit of infinite training data, the results of our algorithm converge to the results of standard TDMs which are known to be consistent. We empirically demonstrate our method can regularize for a wide range of properties—knowledge transfer, sparsity, interpretability, and orthogonality of the learned parameters—on two likelihood models—Gaussian mixtures and Latent Dirichlet Allocation (LDA).

2 Related Work

There exists a very large literature on algorithms for performing inference on latent variable models with various kinds of priors or regularizers. When placing this work in the context of previous and related work, we emphasize that while we draw inspiration from Bayesian statistics, where priors can be regarded as representing pseudo-observations, our formulation is not Bayesian—we do not translate between the pseudo-data and an actual probability distribution on the model parameters. Rather, the pseudo-data should be thought of as a regularization scheme particularly well-suited for tensor decomposition approaches.

Several works have adapted standard tensor decomposition methods to incorporate specific regularizers and constraints: Sun et al. (2015) develop methods designed to find sparse decompositions; Nguyen et al. (2014) produce more robust inference for topic models via the notion of anchor words (Arora et al., 2013). (While not strictly a regularization, the anchor words assumption is a form of prior knowledge imposed on the model parameters.) Cohen et al. (2013); Duchi et al. (2008) detail how to ensure that learned parameters do not have invalid values (e.g. topic models must be valid probability distributions). All of these approaches are specific to the kind of regularization or constraint, and often specific to a certain generative model (e.g. LDA).

Finally, closely related to tensor decomposition methods are methods based on the generalized method of moments (GMM) (Hansen, 1982): both learn the structure of a distribution from low order moments of the data. Recently, several methods have been proposed to regularize GMMs (Tran et al., 2016; Lewis and Syrgkanis, 2018; Yin et al., 2009). In general, however, most GMM algorithms do not leverage the moments eigenstructure to perform the provably optimal inference which is one of the main draws of TDMs.

3 Background and Notation

Many common generative models have parameters that can be expressed as a matrix $A\in\mathbb{R}^{D\times K}$, where $D$ is the dimensionality of the data and $K$ is the number of latent variables. For example, the columns of $A$ could represent the means of Gaussian mixtures or the topic-word probabilities in LDA (Blei et al., 2003). Tensor decomposition methods leverage the relationship between the empirical moments of the data and the latent parameters of the model. Specifically, $A$ is learned by matching the low order moments of the model parameters,

[TABLE]

with empirical estimates, $\hat{M}_{2}$ and $\hat{M}_{3}$, that can be computed from data using expressions which are specific to the particular likelihood model and are presented in Appendix A. Here, $a_{k}$ is the $k^{th}$ column of $A$, and $\beta_{k}$ and $\gamma_{k}$ are constants that come from the model parameters (depend on each generative model). Throughout this paper, we use $a_{k}$ for the column vectors of $A$ and $a_{d,k}$ for individual elements of $A$.

The decomposition to solve for $\{a_{k}\}_{k=1}^{K}$ in Equation 1 given estimates of $M_{2}$ and $M_{3}$ is performed in two stages. We outline the idea here and refer the reader to Anandkumar et al. (2014) for details. First, one computes a whitening matrix $W$ such that $W^{T}\hat{M}_{2}W=I$, and use $W$ to project $\hat{M}_{3}$ to a $\mathbb{R}^{K\times K\times K}$ tensor which has an orthogonal decomposition. Then, the tensor power method can be used to decompose the reduced third order tensor. This process comes with theoretical guarantees on consistency and convergence for recovering $A$ with sufficiently large data sets.

4 Pseudo-data for Regularization

While tensor decomposition methods come with consistency and convergence guarantees, it is not obvious how to incorporate regularization appropriately. The core issue is that tensor decomposition methods provide an algorithm that produces a point estimate of the latent parameters $A$. There exists no log-likelihood and log-prior, as one would encounter in MAP estimation, nor is there an objective where one can simply add an arbitrary regularizer with some desired strength. If we modify the latent parameters, $A$, that are output from a tensor decomposition method, we have no notion of how much predictive quality we sacrifice.

To solve this problem we draw on the intuition which is often used when describing Bayesian priors, of viewing the priors as encoding pseudo-data which match our expectation of the data. Given a prior or a regularizer, we can choose our pseudo-data in a way which will drive the inferred latent parameters towards a form which is more in line with our expectations. We balance that goal with the requirement that our pseudo-data must also be likely under the model parameters we learn using only the real data, which ensures that the pseudo-data are not too different from the real data.

Objective

Formally, let $X_{T}\in\mathbb{R}^{D\times N_{T}}$ be a set of $N_{T}$ observations of true data. Let $X_{P}\in\mathbb{R}^{D\times N_{P}}$ be a set of $N_{P}$ observations of artificial pseudo-data. We seek a pseudo-dataset $X_{P}$ that maximizes the following objective

[TABLE]

where $A_{T}$ and $A_{T\cup P}$ are the parameter matrices learned by the TDM using either only the real training data, $X_{T}$, or a combination of both the real and pseudo-data, $X_{T\cup P}$, respectively. The first term in Equation 2 is the conditional probability of the pseudo-data, $X_{P}$, given the generative model for the data (ex. LDA or Gaussian mixtures) with the parameters $A$ learned using only the training data. This term encourages the pseudo-data, $X_{P}$, to be similar to the true data, $X_{T}$. The regularizer, $R(A_{T\cup P})$, could be any regularizer encoding our desired characteristics for the parameters $A$. The weight $\lambda$ controls the relative strengths of the regularizer and the likelihood of the pseudo-data.

In addition to the weight $\lambda$, the objective in Equation 2 requires us to choose the number of pseudo-data points, $N_{P}$. The number of pseudo-data points $N_{P}$ represents how much confidence we put in our prior knowledge of the parameters structure. An advantage of specifying the strength of the regularizer via the number of pseudo-data is that it limits how much the pseudo-data can influence the learned parameters: as the number of training samples $N_{T}$ increases, for any finite $\lambda$, the maximum effect of the pseudo-data on the inferred parameters diminishes. The dominance of the training data $X_{T}$ as $N_{T}\gg N_{P}$ represents the tendency to put less weight on the pseudo-data as more real data is collected. This notion of convergence is formalized in Theorem 1.

Theorem 1**.**

Fix a likelihood model, $p(X|A)$, and a regularizer, $R(A)$, that is bounded below by a constant $B_{R}$. For any fixed nonnegative $\lambda$ and $N_{P}$, as $N_{T}\rightarrow\infty$, minimizing the loss function, $L$, in Equation 2 with respect to $X_{P}$ results in $A_{T\cup P}\rightarrow A_{T}$.

The proof of Theorem 1 is presented in Appendix B. Furthermore, in Appendix B we also show that if we use plug-in estimators for the bounds of the standard TDMs (e.g. like those in Anandkumar et al. (2012); Hsu and Kakade (2013)), generally the convergence to the true parameters is of the form $A_{T\cup P}=A+\mathcal{O}(\frac{N_{P}}{N_{T}})$ as a consequence of Theorem 1.

Optimization Procedure

Algorithm 1 describes the learning procedure for the regularized tensor decomposition method (RTDM). The function TDM() denotes the model parameters learned using the tensor decomposition algorithm, and ADAM() is the gradient descent step based on the ADAM algorithm (Kingma and Ba, 2014).

Computing the gradient of the loss $L(X_{T},X_{P},\lambda)$ requires taking the gradient of the two terms in Equation 2. The first depends directly on the pseudo-data $X_{P}$; most likelihoods are straight-forward to compute via standard auto-differentiation packages, e.g. the Python Autograd package (Maclaurin et al., 2015). The second term depends on $X_{P}$ implicitly via $A_{T\cup P}=\text{TDM}(X_{T\cup P})$. Fortunately, the standard TDM algorithm consists entirely of linear algebra computations which also renders the second term amendable to automatic differentiation. We initialize our pseudo-data by drawing samples from $p(X_{P}|A_{T})$ so that our pseudo-data start close to the training data, but our results are insensitive to this particular choice of initialization.

Computational Cost

We provide a full discussion of the computational cost of RTDM in Appendix C. We note that the limiting step in the algorithm is computing the whitening matrix, $W$, which involves performing SVD on the moment estimate $\hat{M}_{2}$, resulting in a computational cost of $\mathcal{O}(D^{2})$, and that the moments of the training data only need to be computed once and then cached.

5 Experiments

In the following, we provide a series of demonstrations on the properties and versatility of our approach on synthetic, semi-synthetic and real data under two generative models: Gaussian mixtures and LDA. First, we demonstrate that given prior knowledge of the latent variables structure, our method can improve inference when data is sparse and ignore the prior knowledge when data is abundant. Next, we demonstrate our approach on a range of different regularizers, including those targeting transfer, sparsity, orthogonality, and interpretability. Finally, we explore the effect of different choices for the parameters $\lambda$ and $N_{P}$ and discuss guidelines for choosing them.

Basic Demonstration: Encouraging Prior Properties (with Gaussian mixtures).

We first demonstrate our method for regularizing Gaussian mixture models. The generative process involves first selecting a mixture component $h_{n}$ out of $K$ choices with probability $w_{k}$ $(\sum_{k}w_{k}=1)$, and then sampling the data $x_{n}$ from a multivariate normal $\mathcal{N}(a_{h_{n}},\sigma^{2})$. We generate the matrix $A$ of mixture components means by sampling from a normal distribution $a_{k,d}\sim\mathcal{N}(0,\sigma_{m}^{2})$. The distribution from which $A$ is generated can be used as a prior when learning $A$, and therefore the regularization function is the log-likelihood of the learned $A_{T\cup P}$ under that prior:

[TABLE]

In Figure 1 we demonstrate the RTDM on a synthetic example with $D=10$ dimensions and $K=4$ mixture components drawn from a prior with $\sigma_{m}^{2}=1$ and data variance of $\sigma^{2}=100$. We generate $N_{T}=200$ training data points and use the standard TDM to learn the mixture means. Because of the large variance in the data, the true means are difficult to learn as can be seen by the large errors of the learned means (cyan circles) compared to the true means (green diamonds). Applying the RTDM to regularize the data with $N_{P}=50$ pseudo-observations and $\lambda=0.5$, we see that we are able to use our knowledge of the prior to learn a more accurate estimate of the means (yellow dots). The right plot of Figure 1 demonstrates the optimization where optimizing the loss function (Equations 2 and 3 - top) correlates with increasing the log-likelihood on a test data set of 1000 data points (bottom).

Regularizing for Transfer Learning: Demonstration on LDA with semi-synthetic autism spectrum disorder patient data.

We now demonstrate the RTDM on an LDA generative model, and show that regularization can be used for transfer learning — i.e. use prior knowledge about the parameter structure of our data to improve inference when data is limited, and ignore it when the data is abundant. We use semi-synthetic autism spectrum disorder (ASD) data — the data is a simulated dataset, but simulated from topics learned using real data, and we therefore expect these topics to include the sparsity and correlations which are representative of the true data (Arora et al., 2013). We use electronic health records of $D=64$ common diagnoses of children with autism (Doshi-Velez et al., 2014b). We use the real data to learn two topic matrices ($K=4$) representing common symptoms for two age groups, 6 to 7 and 8 to 9. We make the assumption that the symptoms of the two age groups share some similarities, and that we can transfer knowledge about one age group to better learn the characteristics of the other. In other words, we expect the topics learned for one age group to be an informative prior for inference on the other group.

To test our method over a range of $N_{T}$, we sample observations from the LDA model using the topic matrix of ages 6 to 7. We refer to this topic matrix as $A_{true}$, as it represents the true topics that we wish to learn. As a regularizer, we choose $R(A_{T\cup P})=||A_{T\cup P}-A_{prior}||_{2}$, where $A_{prior}$ is the topic matrix representing ages 8 to 9. In other words, we use RTDM to generate pseudo-data which encourages the learned topics from data on ages 6 to 7 to be as similar as possible to the topics learned from ages 8 to 9. In Figure 2 (left) we demonstrate the results of our method for $N_{T}=100$ and $N_{P}=30$. We see that as we optimize our pseudo-data to decrease $R(A_{T\cup P})$ (solid blue line), the $L_{2}$ distance between the learned and true topics, $\varepsilon=||A_{T\cup P}-A_{true}||_{2}$, also decreases (dashed orange line). In Figure 2 (right), we demonstrate the results for a similar experiment, only this time we choose $N_{T}=10000$. For this much larger dataset $\varepsilon$ is very small, as we have enough data to estimate $A_{true}$ properly. In this case, the $30$ pseudo points are overwhelmed by the true data, and our algorithm has little effect on the final learned topics (dashed orange line value changes very little).

Regularizing for Anti-Correlation: Demonstration with LDA; Exploration of the effect of model parameters.

We now explore the effect of the choice of model parameters, $N_{P}$ and $\lambda$, on the performance of the RTDM in a setting where we use it on a synthetic LDA dataset to learn anti-correlated topics. Such topics structure leads to more diverse topics and can be useful, for example, in classification of text documents for which we wish each learned group of documents to have words which are unique to that specific group. Specifically, the regularizer we use is

[TABLE]

where $a_{i}\cdot a_{j}$ denotes the dot product. This function is zero when all topics are orthogonal. We choose this particular regularizer for its simplicity, but note that other diversity promoting regularizers have been proposed in the literature (Zou and Adams, 2012), and our method can easily be generalized to any other regularization function.

We generate synthetic data from an LDA model with $D=100$ dimensions and $K=4$ topics. In all experiments we use the same $X_{T}$ with $N_{T}=100$, but initialize a new set of $X_{P}$. In Figure 3 we show the topics correlation of $A_{T\cup P}$ for different choices of $\lambda$ and $N_{P}$. Because the full loss function in Equation (2) balances the regularizer and the log-likelihood of the pseudo-data, which is proportional to $N_{P}$, we plot the topics correlation vs. $\lambda/N_{P}$ rather than $\lambda$.

As either $\lambda$ or $N_{P}$ is increased, the pseudo-data pushes the learned topics towards higher diversity. However, for a given $N_{P}$, the influence of the pseudo-data is limited and at some point the topics correlation no longer decreases as $\lambda$ is increased. This implies that an alternative to tuning both parameters of our method could be to set $\lambda$ to a very high value, and adjust the weight we put on our regularizer by only tuning $N_{P}$.

Regularizing for Interpretability via Hierarchy Matching: Demonstration with real data of Medical Subject Headings hierarchy.

We demonstrate the performance of the RTDM with a complex regularization function used on real data. The National Library of Medicine (NLM) uses a hierarchically-organized terminology of medical subject headings (MeSH) for indexing medical articles111https://www.nlm.nih.gov/mesh/. Every article is labeled using several headings, which are given an assignment on a hierarchical tree, in which the root represents a general topic, and headings become more specific farther down the tree. An example of three generations of headings in the tree is “Adult [M01.060.116]”, “Aged [M01.060.116.100]” and “Aged, 80 and over [M01.060.116.100.080]”. Formally, each three digit number in the full heading represents a node, and the periods separating them represent edges.

Topic modeling on subject headings of papers can help in identifying publication and research trends by finding headings which occur together frequently. Because articles are hand labeled, there could be significant inconsistency in labeling—e.g., a particular paper could be given each of the three headings presented in the example, depending on the particular person who labeled it (Doshi-Velez et al., 2014a). A useful property of the topics which could help avoid missing information due to inconsistency in labeling is pushing for topics with headings which are close on the tree.

To encode tree-structured information into the model, we choose a regularizer that is much more complex than simply regularizing for sparsity or diversity, as we want to use our knowledge about the hierarchical indexing structure of the data. We use the following regularizer to achieve this property:

[TABLE]

where $O_{ij}$ is the distance on the tree between the $i^{th}$ and $j^{th}$ headings. In Appendix D we motivate the choice for this regularizer and explain why it promotes the desired property of learning topics with closely related headings.

We perform experiments on a labeled dataset of research articles on statins—a group of drugs for treating cardiovascular disease (Cohen et al., 2006). Our training data consists of $N_{T}=500$ documents using the $D=300$ most common headings, and we learn a model with $K=3$ topics.

In Figure 4 (left) we plot the value of the regularizer for the learned topics with different values of $N_{P}$ and $\lambda$, and observe that similarly to the last example, a given number of pseudo-observations is limited in the effect it can have on the topics, no matter how large $\lambda$ is. In Figure 4 (right) we plot the log-likelihood of a held out test set of 500 documents, using the learned topics. We see that in most cases, the interpretability of the topics comes at a relatively low cost in terms of prediction accuracy, which grows as $N_{P}$ is increased.

In Appendix E we illustrate the effect of regularization on the interpretability of the topics, by comparing the top eight headings of all three topics for the original and regularized topics with $N_{P}=100$ and $\lambda=1000$. We manually color headings on similar branches of the hierarchy tree for clarity. We see that the regularized topics include more headings from similar branches of the tree, and more noticeably, these headings have significantly higher weights than all other headings.

Regularizing for Sparsity: Finding Structure in Noisy Data.

We now consider an example in which the true structure is sparse but the data has been corrupted by additional noise (not part of the assumed generative model). Specifically, we generate synthetic data according to the LDA model, but add to each element in $X$ noise sampled from a Poisson distribution, $\varepsilon\sim\text{Poisson}(\lambda_{P})$. The topics themselves are sampled from a Dirichlet distribution ($a_{k}\sim\text{Dirichlet}(\alpha_{A}1_{D})$; with $\alpha_{A}=0.1$). The noise causes the learning algorithm to attribute some probability mass to all dimensions, obscuring the sparse structure of the true topics. To address this problem, we use the log-liklihood of the Dirichlet probability distribution to regularize the topics,

[TABLE]

In Figure 5 (left) we demonstrate that as the optimization learns topics that minimize $R(A_{T\cup P})$, the $L_{2}$ norm between the true and inferred topics also decreases. The way in which the algorithm improves inference can be understood by comparing the true topics with the topics learned before and after optimization. The right plot in Figure 5 represents one of the $K=4$ topics, where the x-axis represents the index out of the $D=100$ dimensions and the y-axis the value of that dimension on the simplex. The true—very sparse—topics are represented by green diamonds. Because of the noise, the initial topics learned by the standard tensor decomposition algorithm (cyan circles) are not sparse. The algorithm finds the dimensions which have a significant probability mass, but are just barely above the noise. The Dirichlet prior over the topics prefers sparser topics, which it generates by amplifying the value of the dimensions above noise level, resulting in learned topics which are much closer to the true topics (yellow squares).

We note that unlike all other examples in this paper, where our method’s convergence to the standard tensor decomposition method is desirable when data is abundant, in this case abundance of data will prevent our method from sparsifying the topics, as the noise-free model is misspecified. We may increase $N_{P}$ to allow our algorithm to sparsify topics learned from large noisy datasets—corresponding to using very sharply peaked priors.

6 Discussion

The proposed tensor decomposition regularization algorithm requires two parameters—$N_{P}$ and $\lambda$. We demonstrated throughout this paper that for a given $N_{P}$, at some point increasing $\lambda$ no longer influences the final regularized topics. The intuition behind this saturation is that there is a limit to how big of an effect a small fraction of the data can have on the learned topics. In practice, this means we can choose $\lambda$ to be in the saturated regime (high $\lambda$), and only tune $N_{P}$ to control the strength of our regularizer, reducing the number of parameter choices required.

More broadly, we provide a general approach for regularizing tensor decomposition methods via pseudo-data. The designer has full control over the properties they wish to induce via the regularization, and in this sense, our approach is different from, for example, simply bootstrapping the data (which might provide robustness, but not induce desired properties) or bespoke methods based on specific properties (e.g. anchor words). A strength of this method lies in its ability to use regularization and prior knowledge to improve learning when data is limited, and ignore our prior beliefs when data is abundant, much like the effect of priors in a Bayesian setting. Future work could look at alternate optimization techniques for the pseudo-data, including ways to encourage the collection of pseudo-data to have similar properties to the training data. Another open direction is drawing a tighter connection between our method and Bayesian inference.

Appendix A Empirical moments expressions for LDA and Gaussian mixtures

For convenience, we present here the expressions for the empirical estimates of the low order moments of the data required to perform TDM on the models for which we present experimental results — LDA and Gaussian mixture models. For a full description of TDM we refer the reader to Anandkumar et al. [2014].

As a reminder from Section 3 in the main text, tensor decomposition methods learn the parameters of a latent variable model given as a matrix $A\in\mathbb{R}^{D\times K}$, where $D$ is the dimensionality of the data and $K$ the number of latent variables. TDM works by leveraging the relationship between the empirical moments of the data and the latent parameters of the model. Specifically, $A$ is learned by matching the theoretical moments of the model,

[TABLE]

with their empirical estimates, which can be computed from data and are provided below. Here $a_{k}$ is the $k^{th}$ column of $A$, and $\beta_{k}$ and $\gamma_{k}$ are constants that come from the model parameters (depend on each generative model).

A.1 Empirical moments for LDA

The generative model for LDA is

[TABLE]

where $b_{n}$ is the distribution over topics of sample $n$. $w_{in}$ is the feature assignment of the $i^{th}$ feature in sample $n$, where $w_{in}\in 1,..,D$. $z_{in}$ is the topic assignment of the $i^{th}$ feature in sample $n$. We note that a common formulation of LDA also includes a prior on the topics, $a_{k}$. We omit this prior here as it is not used in TDM.

Given a dataset, we can learn the topics structure by computing the following empirical moments of the data:

[TABLE]

where $e_{w_{i}}\in\{0,1\}^{D}$ is the vector whose only non-zero element corresponds to the feature assignment, $w_{i}$. The topics are found by decomposing the empirical moments in Equations 13 and A.1 and matching them with the moments in Equations 7 and 8, with $\beta_{k}=\beta=\alpha_{B}/(K\alpha_{B}+1)K\alpha_{B}$ and $\gamma_{k}=\gamma=2\alpha_{B}/(K\alpha_{B}+2)(K\alpha_{B}+1)K\alpha_{B}$.

A.2 Empirical moments for Gaussian mixtures

The generative model for (spherical) Gaussian mixture models is

[TABLE]

where $a_{h_{n}}$ is the $(h_{n})^{th}$ column $A\in\mathbb{R}^{D\times K}$ representing the mixture means, and $w\in\mathbb{R}^{K}$ represents the probability of data points to be drawn from each topic $(\sum_{k=1}^{K}w_{k}=1)$.

The empirical moments of the data are given by [Hsu and Kakade, 2013]

[TABLE]

where $\sigma^{2}$ is estimated by the smallest eigenvalue of the covariance matrix, $\mathbb{E}[x\otimes x]-\mathbb{E}[x]\otimes\mathbb{E}[x]$. The topics are found by decomposing the empirical moments in Equations 15 and A.2 and matching them with the moments in Equations 7 and 8, with $\beta_{k}=\gamma_{k}=w_{k}$.

Appendix B Convergence of RTDM to true model parameters

Theorem 2**.**

Fix a likelihood model, $p(X|A)$, and a regularizer, $R(A)$, that is bounded below by a constant $B_{R}$. For any fixed nonnegative $\lambda$ and $N_{P}$, as $N_{T}\rightarrow\infty$, minimizing the cost function, $L$, in Equation 2 with respect to $X_{P}$ results in $A_{T\cup P}\rightarrow A_{T}$.

Proof.

The model parameters matrix, $A_{T\cup P}$, is the result of applying a standard tensor decomposition method with the estimated $i^{th}$ order moments $\hat{M}_{i,T\cup P}$, for $(i\in\{2,3\})$, computed using the augmented data, $X_{T\cup P}$, containing the training and the pseudo-data. Since $\hat{M}_{i,T\cup P}$ is an estimated expectation of some function of the data (depending on the generative model being learned), we can write $\hat{M}_{i,T\cup P}$ as a weighted sum of the $i$-th moment estimated using the training data, $\hat{M}_{i,T}$, and the that estimated using the pseudo data, $\hat{M}_{i,T}$:

[TABLE]

Therefore, if $\|\hat{M}_{i,P}-\hat{M}_{i,T}\|_{2}$ can be bounded by some constant $B_{M_{i}}$ for all $N_{T}$ and all training data $X_{T}$ of size $N_{T}$, then we have that the training data dominates the pseudo data in the moment estimation as as $N_{T}$ increases. More formally, we can conclude that

[TABLE]

It then follows that as the training data size $N_{T}$ increases, $A_{T\cup P}$ converges to the parameter matrix computed by a standard TDM using the training data alone, i.e. $A_{T\cup P}$ converges to $A_{T}$.

It remains to show that $\|\hat{M}_{i,P}-\hat{M}_{i,T}\|_{2}\leq B_{M_{i}}$ for all $N_{T}$ and $|X_{T}|=N_{T}$. We first note that $\hat{M}_{i,T}$ converge to the true (finite) moment $M_{i}$ of the model parameters, that is, for any $\epsilon>0$ there is some integer $N_{\epsilon}$ such that $\|\hat{M}_{i,T}-M_{i}\|_{2}\leq\epsilon$ for all $|X_{T}|>N_{\epsilon}$. On the other hand, since the regularizer term, $R(A_{T\cup P})$, in Equation 2 is bounded below, for any training data $X_{T}$, there must exist some $\delta_{X_{T}}$-ball, $\mathcal{B}(\delta_{X_{T}})$, centered at the the mean of the theoretical distribution of $X_{T}$ such that $log(X_{P}|A(X_{T}))>B_{R}$ for all $X_{P}$ sampled outside $\mathcal{B}(\delta_{X_{T}})$. Thus, the pseudo data $X_{P}^{*}$ that minimizes the cost function $L$ in Equation 2 must lie within $\mathcal{B}(\delta_{X_{T}})$. As $N_{T}$ become sufficiently large, we can assume that $\delta_{X_{T}}$ is constant for all training data set $|X_{T}|=N_{T}$, and, hence, that $\|\hat{M}_{i,p}\|_{2}$ is bounded, say by some constant $\gamma$. Thus, for $N_{T}>N_{\epsilon}$, we have the following for any $|X_{T}|=N_{T}$:

[TABLE]

The above shows that $\|\hat{M}_{i,P}-\hat{M}_{i,T}\|$ is bounded by a fixed constant $B_{M_{i}}$ and thus completes the proof. ∎

The assumption of boundedness of the regularizing function $R(A_{T\cup P})$ can be relaxed: for $X_{P}$ sampled from $\delta$-ball centered at the origin, we assume that the regularization term grows more slowly than the log-likelihood term as $\delta\to\infty$. However, most reasonable regularizers have some notion of "optimal" model parameters, in the sense that there is a particular minimal value for the models parameters which would have been chosen if the training data would have been completely ignored, and therefore the lower boundedness of the regularizer is a very mild assumption from a practical standpoint.

Because TDM are provably converge to the true model parameters when there is no model mismatch [Anandkumar et al., 2014, 2012, Hsu and Kakade, 2013], a corollary of Theorem 1, is that the RTDM algorithm converges to the true model parameters, $A$. Furthermore, the convergence of the TDM is shown through perturbation analysis of the TDM to perturbation of the estimated tensors, $\hat{M}_{i}$. If these analyses are robust to small perturbations, the additional perturbation of $\hat{M}_{i}$ due to the pseudo-data, which is shown in the proof above to be of order $\mathcal{O}\left(\frac{N_{P}}{N_{T}}\right)$ for large $N_{T}$ can be plugged-in to these analyses to obtain the convergence rate of the RTDM to the true model parameters. In particular, for LDA [Anandkumar et al., 2012] and Gaussian mixtures [Hsu and Kakade, 2013], the parameters reconstruction error bound was shown to be linear in the moments perturbation, and therefore for these models the RTDM converges to the true topics at a rate of $\mathcal{O}\left(\frac{N_{P}}{N_{T}}\right)$.

Appendix C Computational Cost

The computational complexity of the tensor decomposition algorithm as it appears in Anandkumar et al. [2014] is $\mathcal{O}(D^{3})$, where the limiting step is in computing $\hat{M}_{3}\in\mathbb{R}^{D\times D\times D}$ and whitening it to the $\mathbb{R}^{K\times K\times K}$ tensor, $\hat{M}_{3,w}$. Zou et al. [2013] demonstrated that for sparse data, the tensor decomposition can be performed in $\mathcal{O}(DK+nnz(X))$ where $nnz(X)$ is the number of non-zero elements in $X$. Because our algorithm is based on differentiating the results of the tensor decomposition algorithm with respect to its input, if we wish for our algorithm to be flexible enough to impose any regularizer, $X_{P}$ will generally not be sparse, and we cannot use the method introduced in Zou et al. [2013]. Instead, we first compute $\hat{M}_{2}$, and use the whitening matrix, $W$ to whiten the data. We then compute $\hat{M}_{3,w}$ directly from the whitened data $X_{w}$, and never explicitly compute $\hat{M}_{3}$. This makes the limiting step in the algorithm the SVD computation of $\hat{M}_{2}$, and the computational complexity of the algorithm is $\mathcal{O}(D^{2})$. We note that even though for high dimensional data with a large number of samples, $N_{T}$, the cost of constructing $\hat{M}_{2}$ is of order $\mathcal{O}(N_{T}D^{2})$, this step must be preformed only once. For each subsequent optimization step we only need to recompute the pseudo-data contribution to $\hat{M}_{2}$ which has computational cost of $\mathcal{O}(N_{P}D^{2})$, and in such cases we expect $N_{P}$ to be much smaller than $N_{T}$.

Appendix D Choice of regularizer for MeSH dataset

For the task of regularizing the model we learn for the MeSH data we choose a regularizer which is much more complex than simply regularizing for sparsity or diversity, as we want to use our knowledge about the hierarchical indexing structure of the data. We use the following regularizer to achieve this property—

[TABLE]

where $O_{ij}$ is the distance on the tree between the $i^{th}$ and $j^{th}$ headings. Minimizing this regularizer rewards topics with several headings which are close to each other on the tree ($-O_{ij}^{-1}$ is more negative), while simultaneously pushing for topics with at least more than one highly weighted heading. The latter property is important in understanding why we chose this regularizer over the more straightforward regularizer of $\sum_{k=1}^{K}\sum_{i\neq j}A_{ik}A_{jk}O_{ij}$—the total distance of headings on the tree weighted by their probabilities. Minimizing the weighted total distance between headings can be artificially reduced by making sparser topics and is trivially zero for a topic which has all of its probability on one heading, whereas the regularizer we use prefers the probability mass to be spread across as many headings as possible.222To see this, consider the case where all elements in the vector $a$ are zero except for $D_{*}$ entries with the value $1/D_{*}$. Then $-\sum_{i\neq j}a_{i}a_{j}=-{{D_{*}}\choose{2}}\frac{1}{D_{*}^{2}}=-\frac{1}{2}(1-\frac{1}{D_{*}})\sim\frac{1}{D_{*}}$, implying that the $-\sum_{i\neq j}a_{i}a_{j}$ term in the regularizer scales in inverse proportion to the number of non-zero elements. While sparsity itself is often considered a desirable property for interpretability, in practice we observed that minimizing the weighted distance between headings tended to produce topics which were too sparse, leading us to choose the regularizer we present in Equation 5.

To make computation easier and more efficient, $R$ can be written in matrix notation as

[TABLE]

where $O^{*}_{ij}=O^{-1}_{ij}$ for $i\neq j$ and $O^{*}_{ij}=1$ for $i=j$.

Appendix E Learned topics for MeSH data