{'fraction_non_alphanumeric': 0.036873238650631995, 'fraction_numerical': 0.01403498981556405, 'mean_word_length': 4.52141426590664, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Nonlinear Matrix Decomposition for Mining the Zeros of Sparse Data ∗

. We describe a simple iterative solution to a widely recurring problem in multivariate data analysis: given a sparse nonnegative matrix X , how to estimate a low-rank matrix Θ such that X ≈ f ( Θ ), where f is an elementwise nonlinearity? We develop a latent variable model for this problem and consider those sparsifying nonlinearities, popular in neural networks, that map all negative values to zero. The model seeks to explain the variability of sparse high-dimensional data in terms of a smaller number of degrees of freedom. We show that exact inference in this model is tractable and derive an expectation-maximization (EM) algorithm to estimate the low-rank matrix Θ . Notably, we do not parameterize Θ as a product of smaller matrices to be alternately optimized; instead, we estimate Θ directly via the singular value decomposition of matrices that are repeatedly inferred (at each iteration of the EM algorithm) from the model’s posterior distribution. We use the model to analyze large sparse matrices that arise from data sets of binary, grayscale, and color images. In all of these cases, we find that the model discovers much lower-rank decompositions than purely linear approaches.

1. Introduction. Many empirical disciplines depend increasingly on principled and transparent methods for high-dimensional data analysis [29,123]. The simplest methods arise from basic tools of linear algebra, for example, when data is stored in a large matrix X, we can use singular value decomposition (SVD) to approximate X by another matrix Θ of lower rank [30]. Very often such decompositions can help to analyze high-dimensional data in terms of a much smaller number of degrees of freedom [26,118].
Many types of data, however, have low-dimensional structure that is not revealed by such methods. This is especially true for nonnegative or binary data that is represented by a large sparse matrix (i.e., a matrix whose elements are mostly zero). Such data is better analyzed by models that introduce some additional constraint or nonlinearity to prohibit values that the data cannot realize [24,41,72,103,119]. Typically, these models parameterize their low-rank matrices as the product of smaller matrices-for instance, writing Θ = WH, where W is tall and thin and H is short and wide-and then optimize the factors W and H in an alternating fashion.
These types of alternating approaches are common for problems in high-dimensional data analysis that do not have closed-form solutions. In such problems, it is natural to identify the largest subsets of model parameters that can be efficiently optimized while holding the others fixed. The conceptual goal here is simple-to encapsulate the model updates as wellunderstood optimizations, and then to rely on the most powerful black-box solvers for those optimizations that we have at our disposal. In a nutshell, that is why many generalized lowrank models perform alternating least-squares or convex optimizations over factors W and H of a low-rank matrix Θ.
One might ask, though, whether this parameterization is always necessary for generalized low-rank models-whether, instead, it might be possible to optimize the low-rank matrix Θ directly, in a way that provides an even higher-level encapsulation of the model updates? There are few tractable optimizations over nonconvex sets of low-rank matrices, but as mentioned above, one notable exception is to compute (via SVD) an optimal low-rank approximation to some other matrix. This suggests using truncated SVDs as the basis of an iterative strategy for generalized low-rank modeling. Of course, such a strategy presumes that there exist generalized low-rank models of interest that can be estimated in this way. The motivation to find such models is compelling: if they do exist, then their fitting procedures can piggyback for the rest of time on better and faster algorithms for SVD. Note that this strategy takes a completely agnostic view of how the truncated SVD is computed. There are many choices for this purpose (e.g., bidiagonalization and column pivoting [22,38], alternating least-squares [46], other nonconvex optimizations [19,57], randomized methods [34,44,80,116,117]), but in this framework the choice becomes a lower-level operational decision that can be based on the specific problem instance. The goal of higher-level encapsulation is achieved no matter which routine is ultimately chosen.
In this paper we develop a nonlinear matrix decomposition (NMD) for sparse nonnegative data that achieves this goal. By this we mean that the decomposition is not computed by an alternating optimization over the factors W and H of a low-rank matrix Θ = WH. Instead we describe an iterative algorithm that leverages the full power of SVD to optimize Θ directly; in particular, at each iteration, the algorithm uses SVD to reestimate Θ as the optimal low-rank approximation to another matrix of inferred values. In addition, as in many previous studies, the decomposition exploits a single but essential nonlinearity [1,6,5,25,36,43,87,95,97,130,131]: given a sparse nonnegative matrix X, it attempts to estimate a low-rank matrix Θ such that X ≈ f (Θ), where f is an elementwise nonlinearity that does not take on negative values. As we shall see, the form of this nonlinearity can be purposefully tailored to reveal low-dimensional structure in sparse data.
Another contribution of our work lies in its appeal to latent variable modeling [3,9]. There are well-known latent variable models for many canonical problems in high-dimensional data analysis (e.g., Gaussian mixture distributions for clustering [85], factor analysis for dimensionality reduction [98]). These models have been widely applied, in large part because there are simple, provably convergent algorithms for estimating their parameters. In this paper, we derive an equally tractable latent variable model for NMD-one in which inference does not require sampling-based or variational approximations, and one in which parameter estimation does not require line searches, learning rates, or projected gradients. This model should be of broad interest because NMD can also be viewed as a framework for unsupervised learning in two-layer neural networks-specifically, those networks [74,83] with one hidden layer of binary threshold or rectified linear units (ReLU). This is perhaps the paper's main conceptual contribution: it elevates yet another core problem of unsupervised learning into the canon that can be studied by especially tractable latent variable models.
The paper is organized as follows. In section 2, we present several motivations for this work, and in section 3, we contrast our goals to those of previous approaches. In section 4, we show how to formulate NMD as a problem in latent variable modeling; here, in particular, we show how repeated SVDs can be used to estimate a low-rank matrix Θ such that X ≈ f (Θ). In section 5, we present experimental results on several illustrative data sets and examine the low-dimensional representations discovered by NMD. Finally, in section 6, we offer our main conclusions and discuss important directions for future work.
2. Motivation. Sparse nonnegative matrices arise in many areas of application. The elements of these matrices can record, for example, the edges of objects in grayscale images [15], the gene expression levels in cells [11], the presence of links in a social network [51], the word counts in a large corpus of documents [26], and the ratings or purchases of users on the internet [66]. In addition, it is widely believed that natural images and sounds are encoded by the brain as sparse distributed patterns of neural activity [28,33,49,91], so that any collection of these patterns (e.g., over time) can also be visualized as a sparse nonnegative matrix.
From many disparate applications, then, there arises the same question: how to discover low-dimensional structure in sparse high-dimensional data? The possibility of such structure is precisely what motivates the search for low-rank models [24,41,72,88,103,107,119]. The rank of a matrix is equal to the number of its columns (or rows) that are linearly independent. If the columns of such a matrix store the patterns of a data set, then the rank provides one way to quantify the data's underlying degrees of freedom. In practice, we can settle for a lowrank matrix that does not exactly match the data, but still provides a close approximation. Moreover, for certain types of data, we can hypothesize a nonlinear relationship between the data and its underlying degrees of freedom. The next sections motivate our low-rank model with these ideas in mind.

Strengths and limitations of SVD.
As is well known, the rank of a real-valued d × n matrix X can be immediately ascertained from its SVD. To do so, we write X = UΣV ⊤ , where U and V are orthogonal matrices (of size d × d and n × n, respectively) and Σ is a rectangular diagonal matrix that stores the singular values Σ ii ≥ 0. The rank of X is given by the number of its nonzero singular values.
The SVD of a matrix can also be used to find its optimal low-rank approximation, where the error of the approximation is measured by the Frobenius norm [30]. In particular, the SVD is used to solve the following optimization, min ∥X − Θ∥ F such that rank(Θ) = r, where r is less than the rank of X. In this case, the solution is given by the truncated SVD, with Θ = UΣ r V ⊤ , where Σ r has the same dimensions as Σ but retains only its top r singular values and replaces the others by zero.
The SVD is widely used in this way 1 for high-dimensional data analysis. Suppose, for instance, that the columns of X store the patterns of a data set, and also that X is well approximated by a matrix of lower rank. Then we can regard the rank r of the approximation as a measure of the data's most consequential degrees of freedom, and we can regard the error in the approximation as a measure of imprecision or noise.
It is worth pausing to appreciate the closed-form solution to (2.1) provided by the truncated SVD. The domain of the optimization is the set of d × n matrices with rank r, where r < min(d, n). Note that these matrices do not form a convex set, and as a result, this overall optimization is not convex. Nevertheless, the truncated SVD yields the globally optimal solution. This is the strength of the SVD, and naturally, one might hope to leverage this strength in the search for low-rank models that are specialized, in some way, to sparse nonnegative matrices. We develop this idea further in subsection 2.5, and then again more fully in section 4.

2.2.
Mining the zeros of sparse data. The goal of NMD is easily stated: given a sparse nonnegative matrix X, it attempts to estimate a low-rank matrix Θ such that X ≈ f (Θ), where f is an elementwise nonlinearity. In general, for nonnegative matrices, we will consider approximations of the form X ≈ max(0, Θ); this is the same hinge nonlinearity that appears in neural networks with ReLUs. However, for binary matrices in particular--where the elements are equal to either zero or one-we will consider approximations of the form X ≈ 1 2 [1+sign(Θ)] with a threshold nonlinearity. Thus, in both cases, the nonlinearity serves to restrict the approximated values of X to their underlying domain-an idea that has been widely explored for models of binary, ordinal, and mixed data [1,43,51,52,87,97,96,131]. This is not, however, the only purpose of the nonlinearity in our model. For nonnegative matrices, as we shall see, it is especially important that the nonlinearity maps all negative values of Θ into zeros of X. As a result, when X is sparse-with mostly zero elements-NMD has much more flexibility than SVD to discover low-rank decompositions.
This idea is illustrated in Figure 2.1. Here, X is a sparse nonnegative matrix of full rank. However, it is possible, by replacing the zeros of X with strategically chosen negative values, to find a matrix Θ of lower rank such that X = max(0, Θ). More generally, the goal is to find a lower-rank matrix Θ such that the reconstruction X ≈ max(0, Θ) is accurate to an acceptable degree of approximation.
As an aside, we note that a similar idea has been very extensively explored for the problem of matrix completion [14,46]. In this problem, one considers a matrix (not necessarily sparse or nonnegative) whose elements are only partially specified. Then, to complete the matrix, one fills in its missing elements by assuming that the matrix has the lowest possible rank. Thus the larger the number of missing elements, the more flexibility one has to complete the matrix with a low-rank model.
In our problem, this flexibility is derived from the zeros of the sparse nonnegative matrix X. In particular, as shown in Figure 2.1, the sparser the matrix X, the more zeros we can mine to lower the rank of Θ. As we shall see, NMD fully exploits this flexibility to replace the zeros of X by strategically chosen negative values. Finally, we note that an even greater flexibility is obtained when the matrix X is binary. In this case, NMD searches for approximations of the form X ≈ 1 2 [1 + sign(Θ)]; thus, in constructing Θ, not only can it replace the zeros of X by arbitrary negative values, but it can also replace the ones of X by arbitrary positive values. Top: X is a sparse nonnegative matrix containing mainly zeros (shown in red). Bottom: Θ is a lower-rank matrix satisfying X = max(0, Θ); it is obtained by replacing the zeros of X, as needed, with strategically chosen negative values (shown in yellow). In general, the sparser the matrix X, the more zeros we can mine to lower the rank of Θ.
While the example in Figure 2.1 illustrates the basic idea behind NMD, it does not particularly convey why such a decomposition might be useful. The matrix X in Figure 2.1 was randomly generated, and for that reason, it is not especially deserving of further study. Next we consider a more interesting example where the low-rank model provides exactly this type of insight.
2.3. Low-rank models of pattern manifolds. We adapt an example from [74] showing how an elementwise nonlinearity can model large disparities in rank between the matrices Θ and X. To this end, let α > 0, and consider the n × n circulant matrices with elements For these matrices, it is straightforward to show that rank(Θ) = 3 for all n ≥ 3, while X is of full rank for a range of values for α. Note that when α is very large, the matrix X in (2.3) reduces to the identity matrix, which is obviously of full rank. However, we will be focused on values of α for which X also has an interesting interpretation as a data set of sparse images; see Figure 2.2.
Let us quickly prove the above claims. First we show that rank(Θ) = 3. This is most easily done by defining the orthogonal column vectors c, s, u ∈ R n with elements for 1 ≤ i ≤ n. These vectors correspond in fact to the only eigenvectors of Θ with nonzero eigenvalues, as can be seen by writing Θ = (1−α) uu ⊤ + α(cc ⊤ + ss ⊤ ). This decomposition establishes that rank(Θ) = 3 for n ≥ 3. Note that this result holds for all values of α ̸ ∈ {0, 1}.
Next we prove that X has full rank. The idea of the proof is straightforward. First we show that the diagonal elements of X are all equal to one; then we show that the remaining (off-diagonal) elements of X are either zero or sufficiently small that none of its eigenvalues can deviate too far from unity. To simplify the analysis, we consider the fixed value α = 1 2 sin π n sin π 2n , (2.5) but it is not difficult to show that the results also hold for all larger values of α. We start by substituting (2.5) into (2.2), which after some simplification yields Immediately we see that Θ ii = 1, and hence from (2.3) we also have X ii = 1 for all terms on the diagonal. In addition, among the off-diagonal terms, we see that Θ ij < 0 (and hence X ij = 0) unless i − j ≡ 1 mod n. Next let us evaluate (2.6) for terms just above or below the diagonal. Then we see (for n ≥ 3) that each row of X has exactly two positive off-diagonal terms bounded by From this last result, we have shown not only that X ii = 1 everywhere on the diagonal, but also that j̸ =i |X ij | < 1. It follows from the Gershgorin circle theorem that X has no zero eigenvalues and hence must be of full rank. The particular choice of α in (2.5) was convenient for this proof, but note that any larger choice pushes X even closer to the identity matrix (and its eigenvalues closer to unity).
The example in (2.2) and (2.3) may seem contrived, but as shown in [74], it arises naturally from a data set of sparse one-dimensional images (i.e., each image consisting of a single row of n pixels) that are related by simple translations. In particular, consider a data set of n such images in which each image contains a symmetric blip-three pixels wide, and darkest in its center pixel-against a light background. Here we also allow the blip to wrap around the edge of the image if its center lies on a boundary; see Figure 2.2. If we order the images by the center of these blips, then we see that each successive image is formed by shifting the pixels of the previous one; moreover, the data set as a whole is invariant under this operation. This invariance is, in turn, reflected in the eigenvectors of Θ, which correspond to simple Fourier modes. In sum, the data set for X in (2.3) can be viewed as arising from a pattern manifold of one-dimensional translations, and the essential structure of this manifold is reflected in the much lower rank matrix Θ of (2.2). This transformation is illustrated in Figure 2 It is worth emphasizing two points about this transformation. First, it elicits a manifoldlike structure that is not immediately apparent in the original space of images. In particular, note that most pairs of images in this data set have zero overlap (i.e., they are orthogonal) Figure 2.2. Illustration of the matrices Θ and X in (2.2) and (2.3). The matrix X stores a data set of sparse, one-dimensional images in which a darkened blip (three pixels wide) is translated across a light background. The matrix Θ, of much lower rank, satisfies X = max(0, Θ). See text for further details.
when viewed as vectors of pixel values. Second, the transformation does not precisely reveal the number of degrees of freedom in the data, but it does provide, through the rank of Θ, an upper bound on this number. To be precise, the rank specifies a dimensionality in which the data's underlying manifold can be embedded. In Figure 2.2, for instance, the data have one essential degree of freedom, and the underlying one-dimensional manifold-the circle-is embedded in the three-dimensional space spanned by the columns of Θ.
It is widely believed that many data sets-and perhaps even the neural responses that underlie brain activity-can be understood or analyzed in terms of these pattern manifolds [102]. The degrees of freedom in these manifolds may arise from spatial transformations (e.g., translation, rotation) in the physical world; they may also reflect the continuous variabilities that are inherent to classes of diverse objects (e.g., the different shapes of faces). In section 5, we consider several sparse data sets, of increasing complexity, that are motivated by the example of this section.
2.4. Interpretation as a neural network. Low-rank models with elementwise nonlinearities have also been studied as a paradigm for unsupervised learning in neural networks [49,50,74,83]. Consider the network shown in Figure 2.3, with a bottom layer of r hidden units, a top layer of d visible units, and an r×d weight matrix W connecting the units in these layers. We assume that r < d, so that the network is attempting to explain a larger pattern of activities in the visible layer by a smaller pattern of activities in the hidden one. Such a network can be trained, in an unsupervised fashion, from a data set (x 1 , x 2 , . . . , x n ) of n visible patterns, each of dimensionality d. In this case, the goal of learning is to estimate a weight matrix W and infer a corresponding set (h 1 , h 2 , . . . , h n ) of hidden patterns, each of dimensionality r, such that where f is an elementwise nonlinearity of the sort we have described above. Equation (2.8) makes clear the connection to the low-rank models that we are considering in this paper. In particular, suppose that we represent the data set by a d × n matrix X and that we represent the network's inferred patterns by an r × n matrix H.
Unsupervised neural networks are commonly trained by alternating procedures for inference and learning. Roughly speaking, the goal of inference is to reveal target activities for the hidden units, and the goal of learning is to stabilize these target activities by adapting the network's weights. For the network in Figure 2.3, this can be implemented as an alternating optimization over the factors H and W of the low-rank matrix Θ = WH. As already mentioned, this paper explores a different approach, optimizing directly over the set of low-rank matrices from which Θ must be chosen. This approach can be viewed as a strategy for jointly estimating the network's weights along with the hidden unit activities for all of the training patterns. We develop this approach further in the next section.

The idea in a nutshell.
Suppose that X is a sparse nonnegative matrix. It is possible to convey the main idea behind our approach without fully developing it as a Gaussian latent variable model. We do so here for the problem of estimating a low-rank matrix Θ such that X ≈ max(0, Θ). (The binary case is similar.) Consider the following optimization over the matrix Θ and an additional auxiliary matrix Z, which is not required to be of low rank: We make two simple observations: first, the bottom constraint in (2.9) enforces that the elements of X can be perfectly recovered from those of Z, and second, the objective function in (2.9) is bounded below by zero and only obtains this minimum value when X = max(0, Θ).
This objective function can be minimized by a simple alternating optimization over the matrix variables Z and Θ. Suppose that Θ is fixed. Then Z is optimized by choosing Likewise, suppose that Z is fixed. Then Θ is optimized by computing (via a truncated SVD) the best matrix approximation of rank r to Z. These updates for Z and Θ can be alternated to compute a (possibly local) minimum of the objective function. Intuitively, if this minimum value is close to zero, then it is likely that X ≈ max(0, Θ).
We develop this approach more fully as a latent variable model in section 4. There is a useful analogy here to the k-means algorithm for clustering [79], which also involves an alternating optimization. The k-means algorithm may be viewed as a limiting case of the EM algorithm [27] for parameter estimation in Gaussian mixture models [85]. The optimization in (2.9) plays an analogous role for the latent variable model at the heart of this paper.
3. Relation to previous work. There is a large literature on low-rank models for multivariate data analysis, and our work builds on many previous approaches. In this section we survey these approaches, focusing on those for nonnegative or binary data, while highlighting the essential similarities and differences with our own work. We also return to some of these comparisons in section 6.
3.1. Nonnegative matrix factorization. Interest in low-rank models for nonnegative data [92] exploded after the work of Lee and Seung [72,73] on nonnegative matrix factorization (NMF). Given a nonnegative matrix X, the goal of NMF is to discover a low-rank factorization X ≈ WH where the factors W and H are also constrained to be nonnegative. In NMF these factors are estimated to minimize the error of the approximation, as measured by either the Frobenius norm or generalized Kullback-Liebler divergence. The popularity of NMF is owed in large part to closed-form, multiplicative updates that can be derived for updating W and H. These updates are not only simple to implement, but can also be shown to monotonically decrease the error of the approximation. In fact, the updates are widely used even though more traditional solvers (e.g., based on quasi-Newton methods [12]) are often much faster. It also remains an active area of research [2,20,31,55,63] to develop even faster approaches.
The popularity of NMF is owed in addition to the interpretability of its low-rank models. By constraining the factors W and H to be nonnegative, NMF is able to discover parts-based representations of the objects that it is used to model. For example, NMF can model images of faces as the compositions of different facial features, such as eyes, noses, and lips. These representations are notably different than those discovered by SVD or vector quantization, and they have been explored for many different applications [21,35,37,94,111]. Models based on similar representations have also been explored for factorizations of binary matrices [8] (whose elements are either zero or one) and stochastic matrices [54,100] (whose rows or columns sum to one).
Our work differs in motivation from NMF; whereas NMF's low-rank models seek partsbased representations of the data, ours are predicated on the search for pattern manifolds. These are two different types of low-dimensional structure that can exist (or coexist) in highdimensional data. The inventors of NMF were, in fact, quite aware that not all pattern manifolds can be described by purely linear models [74,102]. The reason is simple: some form of nonlinearity may be required to express the data in terms of its essential degrees of freedom. Note also that due to the nonnegativity constraints on W and H, NMF generally requires factorizations of higher rank than SVD to achieve the same degree of approximation. Thus, while NMF excels at discovering parts-based representations, it cannot reveal many other types of low-dimensional structure (e.g., the example in Figure 2.2).

Exponential-family Principal Component Analysis.
A different low-rank model for nonnegative or binary data was developed by Collins, Dasgupta, and Schapire [24]. Their approach can be viewed as an extension of generalized linear modeling [84] to the problem of matrix factorization. As usual, let X denote the matrix of data, and let Θ denote a low-rank matrix parameter of the same size as X. In this model, the type of data is paired with its appropriate distribution from the exponential family. Thus, for binary matrices, the model maximizes the log-likelihood whose form is derived from the Bernoulli distribution for binary random variables, while for matrices of whole numbers, it maximizes the log-likelihood whose form is derived from the Poisson distribution over nonnegative counts. (In the latter case, the matrix Θ is also constrained to be nonnegative.) The model is known as exponentialfamily principal component analysis (efPCA) because it contains, as a special case, the standard method of principal component analysis (PCA): this is the case when the model employs a Gaussian distribution over real-valued matrix elements. The model for binary data in (3.1) is known as logistic PCA.
As in NMF, these models are estimated by parameterizing Θ = WH as the product of smaller matrices and alternately optimizing over the factors W and H. In efPCA, each of these alternating optimizations is convex (although it is not generally possible to derive closed-form updates). For binary data, in particular, these optimizations take the form of logistic regressions. It is also possible, by introducing an auxiliary function, to optimize (3.1) by an alternating least-squares method [101].
Our model for NMD uses elementwise nonlinearities in a similar way to efPCA. But whereas efPCA is rooted in distributions from the exponential family, NMD is formulated as a Gaussian latent variable model-even for data that is binary or nonnegative. (In NMD, the values of observed data arise by quantizing or clipping the model's Gaussian latent variables, an idea that has been explored in a wide variety of contexts [1,6,5,25,36,43,87,97,131].) As we shall see, it is due to this formulation that NMD can optimize Θ directly without resorting to an alternating minimization over its smaller factors. Naturally, it is also possible to develop models that combine intuitions from both NMF and efPCA. This has been done mainly for binary data: such models maximize the loglikelihood in (3.1) while in addition constraining one or both of the factors W and H to be nonnegative [70,114].
3.3. Probabilistic models. Many researchers have more fully developed probabilistic models of matrix factorization [3,10,17,39,48,86,88,93,112,113]. These models often incorporate a prior distribution over possible factorizations, and some also seek to learn a proper distribution that describes how the data were generated. The estimation of these models can involve extra complications, such as additional iterative procedures-based on Markov chain Monte Carlo (MCMC) [89] or variational methods [61]-for approximate inference. On the other hand, these models are often better regularized and easier to interpret. In this section, we survey some especially prominent applications of these models to nonnegative and binary matrices.
One such application has been topic modeling. While NMF can be used to factor large word document matrices, it does not provide a fully probabilistic model of how words come to appear in text. Such a model was developed by Blei, Ng, and Jordan [10] to discover prevalent topics that run through a large corpus of documents. The model is known as latent Dirichlet allocation (LDA) because it places a Dirichlet prior over the distribution of topics in a document; the number of topics in LDA plays a role analogous to the rank of the factors in NMF. It is not possible to perform an exact inference over the latent variables in LDA, but it is possible to scale variational methods for approximate inference to very large corpora [53]. The model is trained to estimate a distribution over words for each topic (analogous to the nonnegative matrix W in NMF), and in the course of this training, it infers a posterior distribution over topics for each document (analogous to the nonnegative matrix H). Thus LDA discovers similar parts-based representations of documents as NMF, but the topics in LDA are even easier to interpret by virtue of the model's explicitly probabilistic semantics. As mentioned earlier, NMD was not conceived to learn parts-based representations of data, nor is it specialized to count-based data. It is possible, however, that it could serve as a complement to LDA for topic modeling; we discuss this possibility further in section 6.
Probabilistic models of matrix factorization have also been widely applied to recommender systems [66]; most relevant to NMD are those that explicitly model the user-item matrix as nonnegative or binary. One such framework was introduced by Mnih and Salakhutdinov [88], who used a Gaussian model with a sigmoid nonlinearity to bound matrix elements between minimum and maximum ratings. A related approach-known as logistic matrix factorization-was taken by Johnson [60] to model the probability of a user choosing an item; it extends the binary model from efPCA [24] by adding bias terms and Gaussian priors on the factors. Large-scale models for binary data have also been investigated in a Bayesian framework using variational methods for approximate inference [48]. While some user-item matrices record ratings on a binary or Likert scale, others record the number of times that users have purchased items; this type of data is most naturally treated in models of Poisson matrix factorization [39,40]. Like LDA, Poisson matrix factorization relies on variational methods for approximate inference of its latent variables. For very sparse matrices of implicit feedback data, it also has two key advantages: first, it can incorporate domain knowledge through informative priors, and second, its computations scale linearly in the number of nonzero elements of the user-item matrix.
Probabilistic models of matrix factorization have also been widely applied to linkage analysis in social and biological networks. Typically, these networks are represented by binary matrices where zeros and ones indicate the absence and presence of links. A low-rank factorization for such matrices was explored in an influential line of work by Hoff [51,52]. Hoff's initial eigenmodels were estimated using Bayesian methods, but later approaches by Wu, Levina, and Zhu [124] and Ma, Ma, and Yuan [78] led to more efficient algorithms based on projected gradient methods. These later approaches also substituted a logistic link function for the probit link function in Hoff's model. The factorizations in these models are more general than what we explore here, as they incorporate auxiliary information about individual nodes in the network. But even ignoring these terms, it is fair to say that these models differ considerably in both their mechanics and motivation from our approach. In particular, unlike Hoff's model, we do not attempt a Bayesian treatment (thereby avoiding the expense of MCMC methods), and unlike the projected gradient methods, the optimizations for NMD do not involve the tuning of variable or adaptive step sizes. Finally, NMD was motivated more generally for sparse rectangular matrices with nonnegative elements as opposed to square (and typically, symmetric) matrices with binary elements.
In addition to the above applications, there has also been a great deal of methodological and theoretical interest in Gaussian latent variable models with elementwise nonlinearities. These models have been widely developed for matrices of binary, ordinal, and mixed data [3,5,6,25,36,43,87,97], most recently through the use of copula methods [1,130,131]. All of these models use an elementwise nonlinearity to relate an observed (or partially observed) matrix of values to an unobserved matrix of Gaussian latent variables; this is also the starting point of our work. To the best of our knowledge, though, such models have not been empirically or systematically investigated for sparse nonnegative data using an elementwise ReLU nonlinearity. We develop these ideas further in section 4.

Theoretical guarantees.
Suppose that X ≈ f (Θ), where f is an elementwise nonlinearity and Θ is a low-rank matrix. Then each observed element of X provides some information about the corresponding element of Θ. It is known that under certain assumptions, a low-rank matrix Θ can be recovered (or completed) given the values of only a small fraction of its elements [14,18,62]. It is natural to ask whether these results can be extended to the nonlinear setting where Θ must be additionally inferred from observed elements of the higher-rank matrix X ≈ f (Θ).
This question has been investigated by a number of authors. Davenport et al. [25] studied the problem of 1-bit matrix completion, where X is a partially observed binary matrix, and presented a convex program to recover Θ with high accuracy. Bhaskar [5] studied the more general problem, where the elements of X are quantized, but not necessarily binary, and presented a globally convergent algorithm with similar guarantees of recovery; her paper also compares a number of related approaches [13,16,42,68,67,69,106] for these problems. More recently, Mazumdar and Rawat [83] studied this problem where f (z) = max(0, z) is a ReLU operation and proved that a maximum likelihood estimate could recover Θ with high probability and small error. (They did not, however, describe or empirically investigate an algorithm for maximum likelihood estimation.) Ganti Balzano, and Willett [36] described an even more general approach to jointly estimate both the low-rank matrix Θ and the elementwise nonlinearity f ; they assume only that the function f is monotonic, and within this framework, they also provide bounds on the mean squared error of the recovered matrix. Finally, more specialized guarantees are also available for the recovery of rank-one matrices [6,7], and nonrigorous results have been obtained in the thermodynamic limit of very large matrices [75].
At the moment, we cannot provide similar theoretical guarantees for the method in this paper. Our approach differs from the above lines of work in two main respects: first, we assume that the matrix X of data is fully observed, and second, we are motivated by the possibility that rank(Θ) ≪ rank(X), thus revealing the low-dimensional structure of some underlying manifold. The latter hypothesis may be at odds with typical assumptions that the matrix Θ is incoherent or not very spiky. In particular, many of the above analyses require the elements of Θ to be of order unity or similarly bounded; note, however, that this intuition does not hold for the example of subsection 2.3, where from (2.6) we see that max ij |Θ ij | ∼ O(n 2 ). In addition, many of the above guarantees are obtained by minimizing a convex surrogate for the rank of a matrix (e.g., the nuclear norm). But other work in matrix completion has shown that the manifold structure of data is sometimes best recovered by pursuing a nearly opposite objective (e.g., maximizing the trace of a positive semidefinite matrix [110,120,121]).
3.5. Generalized low-rank models. There have been many efforts to treat both probabilistic and nonprobabilistic models of matrix factorization in a unified framework [41,81,103,107,119]. These efforts have been able to identify core optimizations at the heart of many different models. As above, let us denote the factors in these models by W and H. Then it has been observed, across a very wide range of low-rank models, that their optimizations are biconvex in these factors. This observation has led to an extensive study of alternating minimization algorithms for low-rank models. Such algorithms reflect a common wisdom, namely, that in any nonconvex problem, it is often expeditious to identify and solve the largest subproblems that are convex or otherwise tractable. The latent variable model in this paper presents an intriguing special case where this wisdom can be even more powerfully applied. The model has a matrix parameter Θ of fixed rank, and its overall optimization is nonconvex. Nevertheless it is not necessary to parameterize the matrix Θ as the product of smaller factors; instead it is possible to update the low-rank matrix Θ directly, without line searches or learning rates. In addition, these updates are guaranteed to converge monotonically in the model's likelihood. Thus our approach demonstrates the potential for larger substructures to be exploited in the optimization of some generalized low-rank models.
3.6. Exceptions to the above. It is worth noting some exceptions to the above, where generalized low-rank models were estimated by a novel application of SVD rather than an alternating optimization over factors. There are two studies, in particular, that strongly motivated our work.
The first is the continuous latent variable model for binary data in Lee and Sompolinksy [74]. Their model was fit by a moment-matching method that takes an especially simple form for unbiased data (where each bit has equal probability to be zero or one). Our work extends their approach in two ways-first, by showing that similarly motivated models can be fit by maximum likelihood estimation, and second, by considering the case of sparse nonnegative data, to which their moment-matching method is less easily generalized. We also adapted the example in subsection 2.3 from their paper.
The second is the work of Srebro and Jaakkola [108] on weighted low-rank approximations. As discussed in subsection 2.1, the errors of low-rank approximations are conventionally measured by a Frobenius norm that treats all matrix elements with equal weight. Suppose, however, that the approximation is chosen to minimize a weighted sum of elementwise errors. For this case, Srebro and Jaakkola derived an EM algorithm where a truncated SVD is used to improve the approximation at each iteration. Our work is based on an analogous application of EM, but to low-rank models with an elementwise nonlinearity.
Finally, on a related note, we mention a recent study on nonnegative low-rank matrix approximations [104], in which a low-rank nonnegative matrix is estimated without assuming that it can be written as the product of two smaller nonnegative matrices. This is a novel generalization of NMF that can discover significantly better low-rank approximations-as measured by the Frobenius norm-to a nonnegative matrix of higher rank. We note, however, that the quality of these approximations is still (necessarily) worse than what is provided by a truncated SVD, which is not hampered by nonnegativity constraints. Our method differs from this approach in two important ways: first, it focuses exclusively on sparse matrices, and second, for these matrices, it exploits an elementwise nonlinearity to obtain decompositions of even lower rank. [64,65]. They studied a semidiscrete decomposition (SDD) whose goal is to express a matrix as a weighted sum of constrained outer products. Unlike the models considered earlier, however, this decomposition imposes a discrete constraint: the elements of vectors in these outer products are chosen from the set {−1, 0, 1}. This decomposition can provide comparable approximations as SVD while requiring much less storage.

Semidiscrete decompositions. Yet another type of matrix decomposition was explored in a series of papers by O'Leary and Peleg [90] and Kolda and O'Leary
Our motivation in this paper differs from that of SDD. In particular, with NMD we are not seeking a low-rank decomposition that requires the least amount of storage, but one that suggests or hints at the number of underlying degrees of freedom in the data. Generally speaking, while the quantized elements of SDD are ideally suited to compress dense matrices, they are less apt to model continuous modes of variability. On the other hand, while NMD is ideally suited to analyze sparse matrices, it is not designed to provide further levels of compression.

Model.
In this section we formally describe our approach. We start by formulating the latent variable model for NMD and deriving its log-likelihood. Next we discuss the twin problems of inference and parameter estimation; in NMD, as in most latent variable models, these problems are closely intertwined. Finally we describe how we initialize the model and test for convergence.

4.1.
Formulation. NMD shares many aspects of previous low-rank models, particularly in its interplay of Gaussian latent variables and elementwise nonlinearities. NMD is used to analyze a d × n matrix X, where each d-dimensional column of X stores a single instance of some data set and n denotes the number of such examples. Next we use a similarly sized matrix Θ of rank r < min(d, n) to parameterize a Gaussian latent variable model for the data [51,108,112]. In particular, given the matrix Θ, we generate a distribution over nonnegative or binary matrices in the following way. First, for each element Θ ij , we sample a Gaussian latent variable where the variance σ 2 is a single additional parameter of the model that we will estimate. Then, we obtain the matrix X deterministically from the elementwise nonlinear mapping where f (z) = 1 2 [1 + sign(z)] for binary data and f (z) = max(0, z) for nonnegative data. Thus the nonlinearity ensures that X ij ∈ {0, 1} in the former case and X ij ∈ [0, ∞) in the latter. Our main interest is in the case of sparse nonnegative data, but we present both models for completeness.
The model is estimated by maximizing the likelihood of the data X in terms of the matrix Θ and variance σ 2 . To obtain this likelihood, however, we must first compute the marginal distribution P (X ij |Θ ij , σ 2 ) for an individual observed element; this is done by integrating over those values of its corresponding latent variable Z ij that are consistent with the observation X ij = f (Z ij ). For zero elements of X, we note that f (z) = 0 if and only if z ≤ 0. Thus we have where Φ(·) denotes the cumulative distribution function for a normal distribution with zero mean and unit variance. For nonzero elements of X, the form of the marginal distribution depends on the type of data. For binary data, where all nonzero elements of X are equal to one, we can simply use the complement of the previous result, On the other hand, for nonnegative data, we note that X ij = Z ij if and only if Z ij > 0. Hence for positive values x > 0, we have where Z ij is normally distributed by (4.1). From the above results we easily obtain the overall log-likelihood log P (X|Θ, σ 2 ) of the data under this model. First we observe that in NMD, the elements of X are conditionally independent given those of Θ. Then it follows that log P (X|Θ, σ 2 ) = ij log P (X ij |Θ ij , σ 2 ). (4.6) The parameters Θ and σ 2 in our model are estimated by maximizing this sum, thus accounting for all the observed elements in X.
Before proceeding down this path, it is instructive to compare the form of the log-likelihood in (4.6) to the log-likelihoods for efPCA in subsection 3.2. For binary data, the log-likelihoods in (3.1) and (4.6) are qualitatively similar, except that the logistic model in efPCA is replaced by a probit model in NMD (as in many Gaussian latent variable models for binary-valued matrices [3,51]). On the other hand, for nonnegative data, the log-likelihood in (4.6) is qualitatively different than common forms of efPCA. In particular, unlike efPCA based on the Poisson distribution, NMD is not specialized to count-based data, nor does it constrain the mean and variance of the distribution P (X ij |Θ ij ) to have equal values. Also, unlike efPCA based on the exponential distribution, NMD employs a distribution P (X ij |Θ ij , σ 2 ) that can be peaked at or away from zero.
Finally, we observe that the model for binary data in (4.1) and (4.2) has a degeneracy in its parameterization: in particular, the model's predictions in (4.3) and (4.4) are invariant under the change of parameters Θ → λΘ and σ 2 → λ 2 σ 2 , where λ > 0. The model can be made identifiable in this case by fixing σ 2 = 1. We note, however, that a similar degeneracy arises in models of probit regression, and there it has been observed that the EM algorithm converges more quickly when a variance σ 2 is estimated alongside the model's other parameters [76]. We observed a similar phenomenon in NMD when the value of σ 2 was fixed to unity; we therefore retained and estimated the variance parameter σ 2 for all the experiments in this paper. More recent work also suggests that the EM algorithm may discover better solutions with overparameterized models [125].

Inference.
With the above formulation, we can follow a well-traveled path for latent variable modeling of high-dimensional data. Given a matrix of observed data X, we estimate the parameters Θ and σ 2 of the model by maximizing the likelihood in (4.6). To do this, we avail ourselves of the well-known EM procedure [27] for latent variable models. This procedure alternates between two steps: the E-step computes the posterior means and variances of the model's latent variables, and the M-step uses these posterior statistics to reestimate the model's parameters. The EM algorithm can also be viewed as a way of solving difficult, nonconvex optimizations via a sequence of simpler and better understood procedures [57] (e.g., fitting a Gaussian mixture model via a sequence of least-squares problems). In this section, we focus on the problem of probabilistic inference; these calculations are necessary to perform the E-step of the model's EM algorithm.
We begin by noting that exact inference in this model is tractable. In particular, all the required integrals can be expressed in terms of simple functions that arise from the model's elementwise nonlinearity. Also, by design, the model does not attempt a Bayesian treatment; i.e., it does not incorporate a prior distribution over its parameters Θ and σ 2 . As a result, it is not necessary to resort to approximate procedures based on MCMC methods [89] or variational inference [9,61]. Indeed, the integral in (4.3) is typical of the calculations required for inference.
The problem of inference in NMD is to compute the posterior statistics of the model's latent variables. First let us see why the posterior distributions in NMD take an especially simple form. Recall that for each observed matrix element X ij , there is a corresponding latent variable Z ij . Its posterior distribution is given by Bayes rule: Note that the first term in the numerator equals one if X ij = f (Z ij ) and zero otherwise. Thus there are three possible cases for the posterior distribution in (4.7): (i) for the case X ij = 0, it reduces to a right-truncated Gaussian, with no probability mass for positive values of Z ij ; (ii) for binary data with X ij = 1, it reduces to a left-truncated Gaussian, with no probability mass for negative values of Z ij ; (iii) finally, for nonnegative data with X ij > 0, it reduces to a Dirac delta function centered at Z ij = X ij . The last of these cases is trivial, and the first two give rise to truncated Gaussian distributions, of the sort shown in Figure 4.1.
As we shall see, the EM algorithm for NMD relies on repeatedly computing the posterior means and variances of the latent variables Z ij . As shorthand, we denote these statistics by The nontrivial statistics that we need are illustrated in Figure 4.1. For each zero element of X, we must compute the posterior mean and variance of a distribution that has been truncated on the right (shown in blue). Likewise, for each non-zero element of X (when X is binary), we must compute the posterior mean and variance of a distribution that has been truncated on the left (shown in red).
These statistics are most easily expressed in terms of some elementary functions, which we define next. To start, let φ(z) = 1 √ 2π e −z 2 /2 and Φ(z) = 1 2 1 + erf z/ √ 2 denote, respectively, the probability density and cumulative distribution function for a normal distribution with zero mean and unit variance. Also, as further shorthand, let ψ(z) = φ(z)/Φ(z) denote Table 4.2 Inference in NMD for nonnegative data, along with the equations in which these results are used. Results are expressed in terms of the dimensionless parameter γij = σ −1 Θij and the functions φ(z), Φ(z), and ψ(z) from subsection 4.2.
NMD for nonnegative data (x > 0) likelihood (4.6) P (Xij = 0|Θij, σ 2 ) = Φ(−γij)  Tables 4.1 and 4.2, respectively, for NMD with binary and nonnegative data. The above are not the only inferences that we can make from the model. We can also compute the expected value of an observed matrix element: The results of this calculation are also shown in Table 4.1 when f is a threshold nonlinearity (for binary data) and in Table 4.2 when f is a ReLU nonlinearity (for nonnegative data). As evident from (4.10), these expected values behave very simply in the limit of vanishing variance with E[X ij |Θ ij , σ 2 ] → f (Θ ij ) as σ 2 → 0. In practice, of course, the estimated value of σ 2 will never be exactly equal to zero. Instead, when σ 2 is small but finite, the expected value E[X ij |Θ ij , σ 2 ] approaches this limiting behavior as shown in Figure 4.2.

Learning.
We derive parameter updates for NMD using the EM procedure for maximum likelihood estimation in latent variable models [27]. In particular, given the model's current parameter estimates Θ and σ 2 , this procedure yields updated estimatesΘ andσ 2 that increase the log-likelihood in (4.6). Since the derivation may be of less interest than the final result, we begin by presenting the EM updates as a fait accompli: Note how simply in these expressions the model parameters are reestimated from the posterior statistics in (4.8) and (4.9) of the model's latent variables. We make a number of further observations. First, the update in (4.11) is obtained by performing a truncated SVD of the matrix Z of posterior means. Thus we see that the update leverages the full power of SVD to optimize (at each iteration) over the nonconvex set of low-rank matrices in its parameter space. Second, at each iteration, the updates should be performed in the order shown; this is necessary because the result forσ 2 in (4.12) depends on the result forΘ in (4.11). Finally, we observe that these updates have the desirable guarantee of increasing the likelihood in (4.6) at each iteration (except at stationary points). In practice, we iterate these steps to a desired level of convergence. The final result is a nonlinear low-rank decomposition X ≈ f (Θ), where the error of the approximation is modeled by the magnitude of σ 2 .
Let us now show, in more detail, how these updates are derived. The EM algorithm works, at each iteration, by calculating a surrogate for the log-likelihood that is easier to optimize [27]. For NMD, this surrogate is given by where Z and δZ 2 denote the posterior means and variances of the model's latent variables, as defined in (4.8) and (4.9). The M-step of the EM algorithm reestimates the model's parameters by maximizing the expression in (4.13). This maximization is straightforward and leads immediately to the reestimation formulas in (4.11) and (4.12).

Initialization and convergence.
The optimization of the log-likelihood in (4.6) is not convex and cannot be guaranteed to reach a global maximum. In practice, this means that the results from the EM algorithm can depend on how the model parameters are initialized. We adopted the following simple heuristic for initializing these parameters. First we computed the mean and variance of all the elements in X, given by Then we used these statistics to initialize the model parameters Θ and σ 2 . Specifically, for nonnegative data, we set Θ ij = x and σ 2 = δx 2 , while for binary data, we set Θ ij = Φ −1 (x) and σ 2 = 1. Note that in both cases Θ was initialized by a rank-one matrix in which every element is equal.
We also adopted a simple heuristic for evaluating the convergence of the EM updates. Specifically, we terminated the updates if either (a) they did not improve the log-likelihood (normalized per matrix element) by a minimal increment of 10 −5 , or (b) they had already been applied for a maximum of 512 iterations. In general, we observed that the lower-rank models for NMD were more likely to meet the first of these criteria for early stopping.
For the binary model of NMD, we supplemented these updates with one further optimization for the variance σ 2 . In preliminary experiments, we observed that the estimates for σ 2 converged rather slowly at the end of learning. (In this regime, it appears that the surrogate function in (4.13) does not closely track the effect of σ 2 on the log-likelihood.) To account for this, we performed one additional (golden-section) search of values for the variance σ 2 while holding Θ fixed. In practice this final search led to smaller estimates for the variance σ 2 and larger values for the log-likelihood in (4.6).

Experiments.
We investigated the performance of NMD on five large matrices derived from data sets of binary, grayscale, and color images. We describe these matrices and data sets in subsection 5.1, and we present our empirical results on them in subsection 5.2. Finally, in subsection 5.3, we compare the results from NMD on binary data to the closely related model of logistic PCA.

Data sets.
We describe the data sets in order from least to most complex. Our first data set (dots) consisted of 64 × 64 binary images, with each image containing a single black dot of fixed radius against a white background. Some representative images from this data set are shown in the top left plot of Figure 5.1. The data set had n = 1024 images, each containing d = 4096 pixels, of which a very small fraction (ρ = 0.011) were shaded. This data set can be viewed as a higher-dimensional analog of the example in Figure 2.2; in this case, the images lie on a pattern manifold with two translational degrees of freedom. One might expect these data to be modeled very well by an NMD of rank r = 5; intuitively, NMD should require two additional eigenvectors beyond those in Figure 2.2 to account for the extra translational degree of freedom. We will confirm this hypothesis in the next section.
Our second data set (sprites) also consisted of 64 × 64 binary images, but in this case, each image contained a single sprite of varying shape, size, location, and orientation. Some representative images from each class of shapes-squares, ovals, and hearts-are shown in the top right plot of Figure 5.1. The images in our experiments were evenly subsampled from a larger publicly available data set of sprite images [82]; our subsample had n = 15246 images, each containing d = 4096 pixels, of which a small fraction (ρ = 0.048) on average were shaded. This data set has more variability than the previous one, but the basic degrees of freedom remain simple to characterize.
Our third data set (digits) consisted of 28 × 28 grayscale images of MNIST handwritten digits [71]. A representative image from each digit class is shown in the center plot of Figure 5.1. We experimented on the n = 70000 combined images of the training and test set. These images are relatively sparse, with a small fraction (ρ = 0.191) of pixels on average having nonzero values. The variability in this data set arises not only from different classes of digits, but also from different templates within each class as well as variable aspects of penmanship (e.g., slant, thickness).
Our fourth and fifth data sets were derived from 96 × 96 color images in the STL-10 data set [23]. For these experiments, we used the subset of n = 13000 labeled images, each containing an object from one of ten classes. A representative image from each digit class is shown in the bottom plots of Figure 5.1. It is evident that these images are not sparse-and thus not suitable in their own right for models of NMD. However, there are many ways to derive sparse high-dimensional descriptors from dense color images, and we investigated two ways in particular, which we describe next.
For our fourth data set, we computed color histograms for the labeled subset of STL-10 images. To compute these histograms, we first requantized each color channel from 8 bits to 4 bits, so that each RGB pixel was represented by a 12-bit number. Then, for each image, we compiled a color histogram with 2 12 = 4096 bins. In this way we obtained a data set of n = 13000 sparse color histograms in which on average only a small fraction (ρ = 0.066) of the d = 4096 bins had nonzero counts. Note that each such histogram (conveniently, with 4096 = 64 2 bins) can be also visualized as 64 × 64 grayscale image in its own right. For illustration, we have accompanied each color image in the bottom left plot of Figure 5.1 by its corresponding histogram, as visualized in this way.
For our fifth data set, we computed approximate edge maps for the labeled subset of STL-10 images. This was done by computing the image gradient magnitude (across all three color channels) at each pixel in these images and then thresholding the values of these magnitudes, zeroing out all but those in the top 20%. In this way, we obtained the grayscale edge maps shown in the bottom right plot of Figure 5.1. Thus the resulting data set consisted of n = 13000 edge maps, each with d = 9216 grayscale pixels of which a fixed fraction (ρ = 0.2) per image were nonzero.

Empirical results.
The main goal of our experiments was to compare the accuracy of the low-rank approximations obtained by NMD versus SVD. To this end, we estimated models of widely varying sizes on all five of the data sets described in the previous section. Specifically, for the matrix X of each data set, we used the updates in (4.11) and (4.12) to learn low-rank matrices Θ NMD and variances σ 2 NMD maximizing the log-likelihood in (4.6). In addition, as a baseline, we performed truncated SVDs to obtain low-rank matrices Θ SVD minimizing ∥X−Θ SVD ∥ 2 F . We compared the quality of these approximations by attempting to reconstruct X from either Θ SVD or Θ NMD and then calculating the root mean squared error (RMSE) of these reconstructions. The reconstructions were computed aŝ X SVD = max(0, min(u, Θ SVD )), where the upper limit in (5.1) was set to u = 1 for binary data and u = ∞ for nonnegative data, and the RMSE was computed as Note that in (5.1), we clipped the predicted values from SVD to lie within the bounds of the data. This always lowers the RMSE for the model obtained by SVD, thus yielding a stronger baseline. Also, in (5.2), we computed the predictions for NMD from the expected values in (4.10), thus incorporating the model's variance parameter, and not merely by settinĝ X ij = f (Θ ij ). In practice, the former predictions were considerably better for small model sizes (where the estimated variance is higher).
The results of these experiments are shown in Figure 5.1. Each panel compares the RMSE from SVD and NMD on a particular data set and for different values of the rank r. The results show a clear pattern. On one hand, the reconstructions from SVD and NMD behave similarly for small values of the rank (where they both capture only the coarsest properties of the data) and also for sufficiently large values (where they both reconstruct the data with high accuracy). On the other hand, between these two extremes, SVD requires factorizations of much higher rank to achieve the same overall RMSE. Figure 5.1 provides quantitative evidence that NMD can learn more accurate low-rank models than SVD. But it is also interesting to examine the qualitative differences between these models. Some of these differences are illustrated in Figure 5.2, which shows the different reconstructions from NMD and SVD as images in their own right. The former are much sharper and contain many fewer artifacts, suggesting that NMD has better modeled the degrees of freedom at play in this data.

5.
3. Comparison to logistic PCA. As mentioned in subsection 3.2, NMD decomposes binary matrices in a similar manner to logistic PCA (henceforth, σPCA). These methods can be viewed as generalizations of probit and logistic regression to the problem of matrix factorization. It is well known that probit and logistic regression produce generally comparable models for binary classification, so we might also expect NMD and σPCA to obtain comparable low-rank decompositions.
To test this hypothesis, we trained a family of comparable models for σPCA on the dot and sprite data sets of binary images. For these models, we parameterized the matrix Θ in (3.1) as the product of two smaller matrices W and H, and we estimated these smaller matrices using an alternating least-squares method [101]. We initialized the matrices W and H with small random values and adopted the same convergence criteria as described in subsection 4.4. We also performed a final exhaustive sweep to fine-tune the overall scale of the matrix Θ.
The results of these experiments are shown in Figure 5.3. The figure compares NMD and σPCA over a wide range of model sizes. The main plots show that NMD and σPCA yield similar reconstruction accuracies, confirming our expectations based on the similarity of probit and logistic regression. The insets of these plots, however, compare the computation time per iteration of these algorithms, and here we see a marked difference. It is clear that the update for Θ in NMD scales better with the model size than the updates for the individual factors W and H in σPCA. We explain this next. To begin, we examine the time complexity of each EM iteration for NMD. The time complexity in these experiments was dominated by the truncated SVD that maximizes (4.11). The complete SVD for a d × n matrix requires O(min(dn 2 , d 2 n)) operations, and for matrices that are not too large, it is actually faster in MATLAB to perform a complete SVD than a truncated one. This is the case for the data set of dot images (d = 4096, n = 1024), and so in this case the time per iteration of the EM algorithm does not depend on the desired rank r of Θ. This analysis is consistent with the flat curve for NMD in the left inset plot of Figure 5.3. On the other hand, for larger matrices-such as arise from the data set of sprite images (d = 4096, n = 15426)-it is faster to perform a truncated SVD; this overall complexity is harder to estimate for the iterative method employed by MATLAB (based on Lanczos bidiagonalization 2 ), but from the right inset plot of Figure 5.3, this routine appears to scale as O(dnr). Next we discuss the time complexity for σPCA; specifically, we consider the approach based on alternating least-squares updates [101]. In each iteration of this approach, it is necessary to construct and solve multiple systems of linear equations; in particular, a different system is needed for each row of W (which is d × r) and each column of H (which is r × n). In total, it requires O(dnr 2 ) operations to construct these systems of linear equations and O((d+ n)r 2 ) operations to solve them. This analysis is consistent with the poorer scaling of σPCA in both inset plots of Figure 5.3.

SVD NMD
In general, we observed that σPCA and NMD took comparable numbers of iterations to converge. Thus even though these models lead to similar results in reconstruction accuracy, it appears that NMD has computational advantages for all but the lowest-rank decompositions.
6. Conclusion. In this paper we have investigated a generalized low-rank model for sparse nonnegative matrices. The crux of our model is a sparsifying elementwise nonlinearity; with this nonlinearity, the model can reveal decompositions of much lower rank than linear approaches. We found this to be the case in all the data sets that we analyzed. Our results also illuminate the role of these nonlinearities in unsupervised neural networks with ReLU hidden units. Though not explored here, one hopeful idea is that NMD might serve as a primitive for layerwise learning in deep neural networks [32,45,50,105,115,126,127,129]. This idea remains an important area for future work.
The starting point of our work was to introduce a low-rank matrix as the parameter of a Gaussian latent variable model. This matrix could then be reestimated, at each iteration of the model's EM algorithm, from the truncated SVD of an inferred matrix with the same number of rows and columns. It should be possible to apply NMD to much larger matrices than we considered here by computing this truncated SVD via randomized methods [44]. These methods are both simple to use and highly accurate. This is the benefit of encapsulating the model updates as truncated SVDs: the implementation of EM can seamlessly leverage other ongoing work in low-rank matrix approximation [116,117].
There may be other ways to optimize the parameters of the latent variable model for NMD. In subsection 5.3, we observed that the EM algorithm scaled better with the rank of the approximation than an earlier implementation of σPCA. In general, however, we have not shown that the EM updates (based on truncated SVDs) converge more quickly or discover better solutions than other approaches [19] to nonconvex optimization (e.g., conjugate gradient ascent, alternating least squares). Arguably, the main virtue of the EM algorithm is conceptual: it exposes and highlights large tractable substructures of the underlying intractable optimization [57]. It seems difficult to ignore these tractable substructures once they have been revealed; that is why, perhaps, it is more usual to accelerate the fitting of latent variable models within the framework of EM [47,58,59,76,99,128] than outside of it.
We noted in section 3 that NMD builds on many previous approaches for high-dimensional data analysis. Some benefits of NMD are most likely to be realized in conjunction with these other approaches. For instance, NMD seems ideally suited to further decompose the sparse parts-based representations discovered by algorithms such as NMF [72] and LDA [10]. NMF in particular is known to discover sparse basis vectors (e.g., pen strokes, facial features) that exhibit continuous modes of variability (e.g., translations, rotations) akin to those in the data sets we studied. Likewise, recall that for every document, LDA infers a high-dimensional vector of nonnegative topic proportions, and these vectors tend to be sparse for large corpora with fine-grained topics. It seems likely that NMD could detect low-rank structure in these vectors as well.
There remains a lingering debate over the merits of linear (parts-based) decompositions versus nonlinear (manifold-based) decompositions of nonnegative matrices [77]. On one hand, the former are observed to be more interpretable; on the other hand, the latter are claimed to be more compact. Suffice it to say, many data sets exhibit both parts-based and manifold-like structure. In general, it seems wiser to combine these intuitions than to insist on one at the expense of the other.
Another application of NMD may arise from problems in extreme multilabel classification [4,56]. Specifically, consider the problem of learning a large number of binary classifiers in parallel, where all the classifiers use the same features as input, and where the labels from parallel tasks are known to be correlated. Let Y denote the ℓ × n binary label matrix, where ℓ is the number of tasks and n is the number of training examples. For this problem, NMD could be used to learn a matrix Θ of rank r ≪ ℓ such that Y ≈ 1 2 (1 + sign(Θ)). Writing Θ = WH, we see that each r-dimensional column of H provides a low-dimensional latent encoding of the ℓ observed labels for its corresponding example. Thus we can reformulate the problem of ℓ-way multilabel classification as a single nonlinear regression from the shared feature space of training examples to a shared latent space of dimensionality r ≪ ℓ.
In this paper we have not addressed the problem of missing data, namely, when some elements of the matrix X are not observed. In fact, the EM updates of subsection 4.3 can incorporate missing elements simply by equating the prior and posterior statistics of their corresponding latent variables, i.e., setting Z ij = Θ ij and δZ 2 ij = σ 2 whenever X ij does not have an observed value. With this minor accommodation, the updates retain exactly the same form 3 as (4.11) and (4.12) for fully observed matrices. These updates may not be especially efficient, however, when the number of missing elements far exceeds the number of observed ones. A better solution may be to leverage more scalable methods for matrix completione.g., methods based on alternating least squares [46] or nuclear norm minimization [14]-that have been expressly developed for very large matrices with many missing entries. Likewise, in place of (4.11), we can also consider regularized decompositions that yield more interpretable solutions [122] or better performance in certain domains [66,88,109]. All of these methods can benefit in turn from the elementwise nonlinearity in NMD, which may yield lower-rank matrices for them to decompose. These possibilities suggest many interesting directions for future work.