\documentclass{article}[11pt] \usepackage{graphicx} \usepackage{url} \usepackage{fancyvrb} \usepackage{relsize} \usepackage{listings} \setlength{\parindent}{0in} \setlength{\parskip}{0.1in} \lstset{language=R} \title{Quick Introduction to Nonnegative Matrix Factorization} \author{Norm Matloff \\ University of California at Davis } \date{} \begin{document} \maketitle \section{The Goal} Given an $u \times v$ matrix $A$ with nonnegative elements, we wish to find nonnegative, rank-$k$ matrices $W$ ($u \times k$) and $H$ ($k \times v$) such that \begin{equation} \label{awh} A \approx W H \end{equation} We typically hope that a good approximation can be achieved with \begin{equation} k \ll \textrm{ rank}(A) \end{equation} The benefits range from compression ($W$ and $H$ are much smaller than $A$, which in some applications can be huge) to avoidance of ``overfitting'' in a prediction context. \section{Notation} We'll use the following notation for a matrix $Q$ \begin{itemize} \item $Q_{ij}$: element in row $i$, column $j$ \item $Q_{i \cdot}$: row $i$ \item $Q_{\cdot j}$: column $j$ \end{itemize} Note the key relation \begin{equation} (WH)_{.j} = \sum_{i=1}^k H_{ij} W_{.i} \end{equation} In other words, in (\ref{awh}), we have that: \begin{itemize} \item Column $j$ of $A$ is approximately a linear combination of the columns of $W$, with the coefficients in that linear combination being the elements of column $j$ of $H$. \item Thus the columns $W_1,...,W_k$ then (approximately) form a basis for the vector space spanned by the columns of $A$.\footnote{Or, at least for the nonnegative {\it orthant} in that space, meaning all the vectors in that space with nonnegative components.} Since the relation is only approximate, let's use the term {\it pseudo-basis}. \item Similarly, we could view everything in terms of rows of $A$ and $H$: Row $i$ of $WH$ is a linear combination of the rows of $H$, with the coefficients being row $i$ of $W$. The rows of $H$ would then be a pseudo-basis for the rows of $A$. \item If (\ref{awh}) is a good approximation, then most of the information carried by $A$ is contained in $W$ and $H$. We can think of the columns of $W$ as ``typifying'' those of $A$, with a similar statement holding for the rows of $H$ and $A$. Here ``information'' could mean the appearance of an image, the predictive content of data to be used for machine text classification, and so on. \end{itemize} \section{Applications} This notion of {\it nonnegative matrix factorization} has become widely used in a variety of applications, such as: \begin{itemize} \item Image recognition: Say we have $n$ image files, each of which has brightness data for $r$ rows and $c$ columns of pixels. We also know the class, i.e.\ the subject, of each image, say car, building, person and so on. From this data, we wish to predict the classses of new images. Denote the class of image $j$ in our original data by $C_j$. We form a matrix $A$ with $rc$ rows and $n$ columns, where the $i^{th}$ column, $A_{\cdot i}$, stores the data for the $i^{th}$ image, say in column-major order. In the sense stated above, the columns of $W$ typify images. What about $H$? This is a little more complex, though probably more important: Each image contains information about $rc$ pixels. The $m^{th}$ of these, say, varies from image to image, and thus can be thought of as a random variable in the statistical sense. Row $m$ in $A$ can be thought of as a random sample from the population distribution of that random variable. So, the $k$ rows of $H$ typify these $rc$ random variables, and since $k$ is hopefully much smaller than $rc$, we have received what is called {\it dimension reduction} among those variables. It should be intuitively clear that this is possible, because there should be a high correlation between neighboring pixels, thus a lot of ``redundant'' information. This is very important with respect to the {\it overfitting} problem in statistics/machine learning. This is beyond the scope of this tutorial, but the essence of fhe concept is that having too complex a model can lead to ``noise fitting,'' and actually harm our goal of predicting the class of future images. By having only $k$ rows in $H$, we are reducing the dimensionality of the problem from $rc$ to $k$, thus reducing the chance of overfitting. Our prediction of new images works as follows. Given the $rc$-element vector $q$ of the new image of unknown class, we find its representation $q_p$ in the pseudo-basis for the rows of $A$, and then find $g$ such $G_{g \cdot}$ best matches $q_p$. Our predicted class is then $C_g$. \item Text classification: Here $A$ consists of, say, word counts. We have a list of $d$ key words, and $m$ documents of known classes (politics, finance, sports etc.). $A_{ij}$ is the count of the number of times word $i$ appears in document $j$. Otherwise, the situation is the same as for image recognition above. We find the NMF, and then given a new document to classify, with word counts $q$, we find its coordinates $q_p$, and then predict class by matching to the rows of $W$. \item Recommender systems: Here most of $A$ is unknown. For instance, we have users, who have rated movies they've seen, with $A_{ij}$ being the rating of movie $j$ by user $i$. Our goal is to fill in estimates of the unknown entries. One approach is to initially insert 0s for those entries, then perform NMF, producing $W$ and $H$.\footnote{A popular alternative uses {\it stochastic gradient} methods to find the best fit directly, without creating the matrix $A$.} We then compute $WH$ as our estimate of $A$, and now have estimates for the missing entries. \end{itemize} \section{The R Package NMF} The R package {\bf NMF} is quite extensive, with many, many options. In its simplest form, though, it is quite easy to use. For a matrix {\bf a} and desired rank {\bf k}, we simply run \begin{lstlisting} > nout <- nmf(a,k) \end{lstlisting} Here the returned value {\bf nout} is an object of class {\bf "NMF"} defined in the package. It uses R's S4 class structure, with \lstinline{@} as the delimiter denoting class members. As is the case in many R packages, {\bf "NMF"} objects contain classes within classes. The computed factors are in {\bf nout@fit@W} and {\bf nout@fit@H}. Let's illustrate it in an image context, using the following: \includegraphics[width=3.6in]{MtRush.png} Here we have only one image, and we'll use NMF to compress it, not do classification. First obtain $A$: \begin{lstlisting} > library(pixmap) # read file > mtr <- read.pnm('MtRush.pgm') > class(mtr) [1] "pixmapGrey" attr(,"package") [1] "pixmap" # mtr is an R S4 object of class "pixmapGrey" # extract the pixels matrix > a <- mtr@grey \end{lstlisting} Now, perform NMF, find the approximation to $A$, and display it: \begin{lstlisting} > aout <- nmf(a,50) > w <- aout@fit@W > h <- aout@fit@H > approxa <- w %*% h # brightness values must be in [0,1] > approxa <- pmin(approxa,1) > mtrnew <- mtr > mtrnew@grey <- approxa > plot(mtrnew) \end{lstlisting} Here is the result: \includegraphics[width=3.6in]{MtRush50.png} This is somewhat blurry. The original matrix has dimension $194 \times 259$, and thus presumably has rank 194.\footnote{This is confirmed by running the function {\bf rankMatrix()} in the {\bf Matrix} package.} We've approximated the matrix by one of rank only 50, a storage savings. Not important for one small picture, but possibly worthwhile if we have many. The approximation is not bad in that light, and may be good enough for image recognition or other applications. \section{Algorithms} How are the NMF solutions found? What is {\bf nmf()} doing internally? Needless to say, the methods are all iterative, with one approach being that of the Alternating Least Squares algorithm. By the way, this is not the default for {\bf nmf()}; to select it, set {\bf method = 'snmf/r'}. \subsection{Objective Function} We need an {\it objective function}, a criterion to optimize, in this case a criterion for goodness of approximation. Here we will take that to be the {\it Frobenius} norm, which is just the Euclidean ($L_2$) norm with the matrix treated as a vector:\footnote{The $L_p$ norm of a vector $v = (v_1,...,v_r)$ is defined to be $$ \left (\sum_i |v_i|^p \right )^{1/p} $$ } \begin{equation} \|Q\|_2 = \sqrt{ \sum_{i,j} Q_{ij}^2 } \end{equation} So our criterion for error of approximation will be \begin{equation} \label{errawh} \|A - WH\|_2 \end{equation} This measure is specified in {\bf nmf()} by setting {\bf objective = 'euclidean'}. \subsection{Alternating Least Squares} So, how does it work? It's actually quite simple. Suppose just for a moment that we know the exact value of $W$, with $H$ unknown. Then for each $j$ we could minimize \begin{equation} \label{errajwhj} \|A_{\cdot j} - W H_{\cdot j}\|_2 \end{equation} We are seeking to find $H_{\cdot j}$ that minimizes (\ref{errajwhj}), with $A_{\cdot j}$ and $W$ known. But since the Frobenius norm is just a sum of squares, that minimization is just a least-squares problem, i.e.\ {\it linear regression}; we are ``predicting'' $A_{\cdot j}$ from $W$. The solution is well-known to be\footnote{If you have background in regression analysis, you might notice there is no constant term, $\beta_0$, here.} \begin{equation} H_{\cdot j} = (W'W)^{-1} W' A_{\cdot j} \end{equation} R's {\bf lm()} function does this for us \begin{lstlisting} > h[,j] <- lm(a[,j] ~ w - 1) \end{lstlisting} for each $j$.\footnote{The -1 specifies that we do not want a constant term in the model.} On the other hand, suppose we know $H$ but not $W$. We could take transposes, \begin{equation} A' = H' W' \end{equation} and then just interchange the roles of $W$ and $H$ above. Here a call to {\bf lm()} gives us a row of $W$, and we do this for all rows. Putting all this together, we first choose initial guesses, say random numbers, for $W$ and $H$; {\bf nmf()} gives us various choices as to how to do this. Then we alternate: Compute the new guess for $W$ assuming $H$ is correct, then choose the new guess for $H$ based on that new $W$, and so on. During the above process, we may generate some negative values. If so, we simply truncate to 0. \subsection{Multiplicative Update} Alternating Least Squares is appealing in several senses. At each iteration, it is minimizing a {\it convex} function, meaning in essence that there is a unique local and global minimum; it is easy to implement, since there are many least-squares routines publicly available, such as {\bf lm()} here; and above all, it has a nice interpretation, predicting columns of $A$. Another popular algorithm is {\it multiplicative update}, due to Lee and Seung. Here are the update formulas for $W$ given $H$ and {\it vice versa}: \begin{equation} W \leftarrow W \circ \frac {AH'} {WHH'} \end{equation} \begin{equation} H \leftarrow H \circ \frac {W'A} {W'WH} \end{equation} where $Q \circ R$ and $\frac{Q}{R}$ represent \underline{elementwise} multiplication and division with conformable matrices $Q$ and $R$, and the juxtaposition $QR$ means ordinary matrix multiplication. \section{Predicting New Cases} Once we have settled on $W$ and $H$, what can we do with them? In the recommender system application mentioned earlier, we simply multiply them, and then retrieve the predicted values at the matrix positions at which we did not have user ratings. But recall the description given above for the image and text processing examples: \begin{quote} Given the $rc$-element vector $q$ of the new image, we find its representation $q_p$ in the pseudo-basis, and then find $g$ such $W_{g \cdot}$ best matches $q_p$. Our predicted class is then $C_g$. \end{quote} How do we find $q_p$? The answer is that again we can use the least-squares above. \begin{lstlisting} > qp <- lm(q ~ h - 1) \end{lstlisting} \section{Convergence and Uniqueness Issues} There are no panaceas for applications considered here. Every solution has potential problems. With NMF, an issue may be uniqueness --- there might not be a unique pair $(W,H)$ that minimizes (\ref{errawh}).\footnote{See Donoho and Stodden,{\it When Does Non-Negative Matrix Factorization Give a Correct Decomposition into Parts?}, \url{https://web.stanford.edu/~vcs/papers/NMFCDP.pdf}. } In turn, this may result in convergence problems. The NMF documentation recommends running {\bf nmf()} multiple times; it will use a different seed for the random initial values each time. The Alternating Least Squares method used here is considered by some to have better convergence properties, since the solution at each iteration is unique. \section{How Do We Choose the Rank?} This is not an easy question. One approach is to first decompose our matrix $A$ via {\it singular value decomposition}. Roughly speaking, SVD finds an orthonormal basis for A, but then treats the vectors in that basis as random variables. The ones with relatively small variance may be considered unimportant. \begin{lstlisting} > asvd <- svd(a) > asvd$d [1] 1.188935e+02 2.337674e+01 1.685734e+01 1.353372e+01 1.292724e+01 [6] 1.039371e+01 9.517687e+00 9.154770e+00 8.551464e+00 7.776239e+00 [11] 6.984436e+00 6.505657e+00 5.987315e+00 5.059873e+00 5.032140e+00 [16] 4.826665e+00 4.576616e+00 4.558703e+00 4.107498e+00 3.983899e+00 [21] 3.923439e+00 3.606591e+00 3.415380e+00 3.279176e+00 3.093363e+00 [26] 3.019918e+00 2.892923e+00 2.814265e+00 2.721734e+00 2.593321e+00 \end{lstlisting} So, the first standard deviation is about 119, the next about 23 and so on. The seventh is already down into the single-digit range. So we might take $k$ to be, say, 10. The columns of $A$ lie in a subspace of $R^{194}$ of dimension of about 10. \section{Why Nonnegative?} In the applications we've mentioned here, we always have $A \geq 0$. However, that doesn't necesarily mean that we need $W$ and $H$ to be nonnegative, since we could always truncate. Indeed, we could consider using SVD instead. There are a couple of reasons NMF may be preferable. First, truncation may be difficult if we have a lot of negative values. But the second reason is rather philosophical, as follows: In the image recognition case, there is hope that the vectors $W_{\cdot j}$ will be {\it sparse}, i.e.\ mostly 0s. Then we might have, say, the nonzero elements of $W_{\cdot 1}$ correspond to eyes, $W_{\cdot 2}$ correspond to nose and so on with other parts of the face. We are then ``summing'' to form a complete face. This may enable effective {\it parts-based recognition}. Sparsity (and uniqueness) might be achieved by using {\it regularization} methods, in which we minimize something like \begin{equation} \|A - WH\| + \lambda (\|W\|_1 + \|H\|_1) \end{equation} where the subscript 1 (``$L_1$ norm'') means the norm involves sums of absolute values rather than sums of squares. This guards against one of the factors becoming ``too large,'' and it turns out that this also can lead to sparse solutions. We can try this under {\bf nmf()} by setting {\bf method = 'pe-nmf'}. \end{document}