Continuous Semi-Supervised Nonnegative Matrix Factorization

Abstract

Nonnegative matrix factorization can be used to automatically detect topics within a corpus in an unsupervised fashion. The technique amounts to an approximation of a nonnegative matrix as the product of two nonnegative matrices of lower rank. In certain applications it is desirable to extract topics and use them to predict quantitative outcomes. In this paper, we show Nonnegative Matrix Factorization can be combined with regression on a continuous response variable by minimizing a penalty function that adds a weighted regression error to a matrix factorization error. We show theoretically that as the weighting increases, the regression error in training decreases weakly. We test our method on synthetic data and real data coming from Rate My Professors reviews to predict an instructor’s rating from the text in their reviews. In practice, when used as a dimensionality reduction method (when the number of topics chosen in the model is fewer than the true number of topics), the method performs better than doing regression after topics are identified—both during training and testing—and it retrains interpretability.

1. Introduction

Nonnegative matrix factorization (NMF) is a highly versatile data science technique with far-reaching applications. It can identify thematic elements, i.e., groups of words that appear frequently together in a corpus, which together convey a common message [1]. More generally, it can be used to decompose an image into identifiable patterns [2] and as a general-purpose dimensionality reduction or preprocessing method before applying other machine learning methods, as has been done in studying various diseases [3,4]. Similar to singular value decomposition (SVD) [5], NMF provides a low rank factorization. In NMF, a nonnegative matrix $X\in {\mathbb{R}}_{\ge 0}^{n\times m}$ representing a corpus (or other nonnegative dataset) is factored into a low rank approximation $X\approx WH$ where the inner dimension, r, between W and H is such that $r\ll m$ and $r\ll n$; however, unlike SVD, there is an additional constraint that both W and H are nonnegative, i.e., $W\in {\mathbb{R}}_{\ge 0}^{n\times r}$ and $H\in {\mathbb{R}}_{\ge 0}^{r\times m}.$ This nonnegativity enforces that the data in X are represented by a nonnegative combination of the dictionary atoms in the factorization, which lends itself to human interpretability. For example, in the foundational work [2], Lee and Seung show that NMF, when applied to facial images, decomposes the images into recognizable parts such as noses, eyes, and mouths.

When applied to a document-term matrix [6] X, where row i of X represents document i and column j represents the frequencies of word j across the documents, the classical NMF method amounts to

$$(W,H)=\arg {\min }_{W\in {\mathbb{R}}_{\ge 0}^{n\times r},H\in {\mathbb{R}}_{\ge 0}^{r\times m}}\left|\right|X-{WH\left|\right|}_F^2$$

(1)

where, for $A\in {\mathbb{R}}^{n\times m}$, ${\left|\right|A\left|\right|}_F$ denotes the Frobenius norm of A with ${\left|\right|A\left|\right|}_F^2=\mathrm{tr}(A^{\mathrm{T}}A)={\sum }_{i=1}^n{\sum }_{j=1}^m{\left|A_{ij}\right|}^2$. Other variations on the penalty function exist, including the Kullback–Leibler divergence [7]. After computing W and H, we interpret row j of H as the jth topic, its components being the weight of each word in that topic, and row i of W as the topic-encoding of document i, i.e.,

$$X_{i,:}\approx \sum _{j=1}^r{W}_{i,j}{H}_{j,:}.$$

Throughout this manuscript we make use of “colon notation” where “:“ means the full range of indices for a row/column, “a:b” indicates a consecutive range of indices from a to b, etc.

The rest of our paper is organized as follows: the general context to our work and our main contributions are stated in Section 2; in Section 3, we provide the formulation of our method, its algorithmic implementation, its theoretical properties, and their proofs; in Section 4, we provide a proof of concept through synthetic data; in Section 5, we test our method on a real dataset from Rate My Professors; and finally we conclude our work in Section 6.

2. Relation to Current Work and Contributions

When designing algorithms to handle multiple objectives simultaneously, a weighted penalty function that combines the multiple objectives is often used [8]. For example, with LASSO regression [9], one seeks a linear model that also does model selection by zeroing out parameters that are less significant. This can be done by minimizing a least squares error with a weighted penalty for the ${\ell }^1$-norm of the parameters. Prior authors have combined NMF with a linear regression procedure to maximize the predictive power of a classifier [10,11,12]. This is accomplished through a penalty function that combines NMF with another objective function—a (semi) supervised approach. Semi-supervised NMF can also be applied to guide NMF to identify topics with desired keywords [13].

In this paper, we combine NMF with a linear regression model to predict the value of a continuous response variable. We consider synthetic data along with a real dataset that pairs written commentary with a real-valued observation. In particular, we use reviews from the Rate My Professors website [14,15] that include all student comments for a professor along with the mean rating in [1, 5]. Due to the averaging, the rating is effectively a continuous variable.

3. Model

We provide the framework for our proposed continuous semi-supervised nonnegative matrix factorization method (CSSNMF).

3.1. Formulation

We consider having a document-term matrix [6] $X\in {\mathbb{R}}_{\ge 0}^{n\times m}$ for n documents with their associated word frequencies in the m columns—a “bag of words” where each document is represented only by the frequencies of its words. Each document has a corresponding value in $\mathbb{R}$ so that we can associate with X the vector $Y\in {\mathbb{R}}^{n\times 1}$. Put another way: each document is represented as a row of X, call it $x\in {\mathbb{R}}_{\ge 0}^{1\times m}$, which stores the frequencies of each of the m words within the corpus; then, to each such x there is an observation $y\in \mathbb{R}$ (and over n documents, this generates $Y\in {\mathbb{R}}^{n\times 1}$). We choose $r\in \mathbb{N}$ and $\lambda \ge 0$ as hyper-parameters where r denotes the number of topics and $\lambda$ is weight put on the regression error.

In the real data that we look at, each row of X will represent the reviews written for one university instructor, with frequencies of the words in the columns. For each instructor in the dataset, the mean value of their respective student ratings will be a single component of Y. From a predictive standpoint, we would like to predict the mean rating of an instructor based only on the words in their review, i.e., take a vector $x\in {\mathbb{R}}_{\ge 0}^{1\times m}$ of the word frequencies and make a prediction $\widehat{y}$ of their mean rating (the hat indicates a prediction). We want the prediction $\widehat{y}$ to be as close to the true mean rating as possible. The topic modelling aspect of this is that instead of using the full x vector of dimension m, we approximate x as a nonnegative linear combination of r topic vectors (interpretable vectors of word frequencies). We effectively compress x to a vector $w\in {\mathbb{R}}_{\ge 0}^{1\times r}$ of dimension $r,$ and we model the rating as a linear combination of the components of w.

Given $W\in {\mathbb{R}}_{\ge 0}^{n\times r}$, $H\in {\mathbb{R}}_{\ge 0}^{r\times m}$, and $\theta \in {\mathbb{R}}^{(r+1)\times 1}$, we define a penalty function that combines topic modelling with a linear regression based on the topic representations. The intuition with the weighting is that as the weight $\lambda$ increases, topic modelling is still done, but more and more emphasis is put on producing an accurate regression on Y. We define

$$\begin{array}{cc}\hfill F^{(\lambda )}(W,H,\theta ;X,Y)& =N(W,H;X)+\lambda R(W,\theta ;Y),\mathrm{where}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill N(W,H;X)& :=\left|\right|X-{WH\left|\right|}_F^2\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill R(W,\theta ;Y)& :=\left|\right|\tilde{W}\theta -Y{\left|\right|}^2\hfill \end{array}$$

(4)

and where $\overline{W}\in {\mathbb{R}}^{n\times (r+1)}$ is given by

$$\overline{W}:=\left(\begin{array}{cccc}1& W_{1,1}& \cdots & W_{1,m}\\ 1& W_{2,1}& \cdots & W_{2,m}\\ ⋮& ⋮& ⋮& ⋮\\ 1& W_{n,1}& \cdots & W_{n,m}\end{array}\right).$$

(5)

The matrix $\overline{W}$ with its column of 1s allows for an intercept: given a topic representation $w\in {\mathbb{R}}^{1\times r}$, we predict a value $\widehat{y}={\theta }_1+{\theta }_2{w}_1+\dots +{\theta }_{r+1}{w}_r.$

When $\lambda >0$, we seek

$$(W^{(\lambda )},H^{(\lambda )},{\theta }^{(\lambda )})=\arg {\min }_{W,H,\theta }{F}^{(\lambda )}(W,H,\theta ;X,Y).$$

(6)

When $\lambda =0$, we define

$$\begin{array}{cc}\hfill (W^{(0)},H^{(0)})& =\arg {\min }_{W,H}N(W,H;X)\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\theta }^{(0)}& =\arg {\min }_{\theta }R(W^{(0)},\theta ;Y).\hfill \end{array}$$

(8)

We also impose a normalization constraint, that

$$\forall i,\sum _{j=1}^m{H}_{ij}=1$$

(9)

so that the topics have unit length in ${\ell }^1$.

Remark 1

(Sum of Topic Representations). If $X\approx WH$ is normalized so its rows sum to 1 then it is also the case that $\forall i,{\sum }_{j=1}^r{W}_{ij}\approx 1$ by noting that $X_{ij}={\sum }_{k=1}^r{W}_{ik}{H}_{kj}$ and summing over j.

When $\lambda =0$, $\theta$ has no effect upon $F^{(\lambda )}$ and we first perform regular NMF over W and H and, as a final step, we choose $\theta$ to minimize the regression error. In other words, if $\lambda =0$, we do NMF first and then find the best $\theta$ given the already determined weights for each document. It seems intuitive, however, that the regression could be improved if $\theta$ and W both were being influenced by the regression to Y, which is what our method aims to do when $\lambda >0$. From a practical perspective, if $\lambda \uparrow \infty$, then the regression error becomes dominant and we may expect the topics as found in H to be less meaningful. In Section 3.2, we state some theoretical properties of our method as it is being trained.

Once $H^{(\lambda )}$ and ${\theta }^{(\lambda )}$ are known, we can make predictions for the response variable corresponding to a document. This amounts to finding the best nonnegative topic encoding $w\in {\mathbb{R}}_{\ge 0}^{1\times r}$ for the document and using that encoding in the linear model—see Section 3.3.

Remark 2

(Uniqueness). Using our established notation, we remark that if $X^{\mathrm{*}}=WH$ and $Y^{\mathrm{*}}={\theta }_1+W{\theta }_{2:(r+1)}$, then $X^{\mathrm{*}}=\tilde{W}\tilde{H}$ and $Y^{\mathrm{*}}={\theta }_1+\tilde{W}\tilde{\theta }$, where $\tilde{W}=WS$, $\tilde{H}=S^{-1}H$, and $\tilde{\theta }=S^{-1}{\theta }_{2:(r+1)}$ for any invertible $S\in {\mathbb{R}}^{r\times r}$ with $SW$ and $S^{-1}H$ both having all entries nonnegative. Thus, uniqueness of an optima, if it exists, can only be unique up to matrix multiplications.

3.2. Theoretical Results

We present two important behaviors of CSSNMF with regard to increasing $\lambda$ and its effect upon predicting the response variable.

Proposition 1

(Regression Error with Nonzero $\lambda$). For $\lambda \ge 0$, let $W^{(\lambda )},H^{(\lambda )},{\theta }^{(\lambda )}$ be a unique (as per Remark 2) global minimum to Equations (6)–(9). Then $R(W^{(\lambda )},{\theta }^{(\lambda )})\le R(W^{(0)},{\theta }^{(0})$.

Theorem 1

(Weakly Decreasing Regression Error). Let $0\le {\lambda }_1<{\lambda }_2$ be given where $W^{({\lambda }_i)},H^{({\lambda }_i)},{\theta }^{({\lambda }_i)}$ are the unique (as per Remark 2) global minimizers of Equations (6)–(9) for $i=1,2$. Then $R(W^{({\lambda }_2)},{\theta }^{({\lambda }_2)};Y)\le R(W^{({\lambda }_1)},{\theta }^{({\lambda }_1)};Y).$

Remark 3.

Proposition 1 and Theorem 1 are based on obtaining a global minimum. In practice, we may only find a local minimum.

Proposition 1 and Theorem 1 are statements pertaining to training the model. Assuming we have the optimal solutions, Proposition 1 tells us that the regression error for $\lambda >0$ is no worse than the regression error with $\lambda =0$ and could, in fact, be better. Thus, the intuition that selecting topics while paying attention to the regression error is practical. Then Theorem 1 says that the regression error is weakly monotonically decreasing as $\lambda$ increases. In practical application, we find the error strictly monotonically decreases.

Before proceeding to algorithmic procedures, we prove Proposition 1 and Theorem 1.

Proof of Proposition 1.

If $\lambda >0$ then

$$\begin{array}{cc}\hfill F^{(\lambda )}(W^{(\lambda )},H^{(\lambda )},{\theta }^{(\lambda )};X,Y)& \le F^{(\lambda )}(W^{(0)},H^{(0)},{\theta }^{(0)};X,Y)\Rightarrow \hfill \\ \hfill N(W^{(\lambda )},H^{(\lambda )};X)+\lambda R(W^{(\lambda )},{\theta }^{(\lambda )};Y)& \le N(W^{(0)},H^{(0)};X)+\lambda R(W^{(0)},{\theta }^{(0)};Y)\hfill \\ \hfill & \Rightarrow \hfill \\ \hfill \lambda (R(W^{(\lambda )},{\theta }^{(\lambda )};Y)-R(W^{(0)},{\theta }^{(0)};Y))& \le N(W^{(0)},H^{(0)};X)-N(W^{(\lambda )},H^{(\lambda )};X)\hfill \\ \hfill & \le 0.\hfill \end{array}$$

The first inequality comes from how $(W^{(\lambda )},H^{(\lambda )},{\theta }^{(\theta )})$ are defined by Equation (6). The final inequality comes from how $(W^{(0)},H^{(0)})$ are defined as minimizers in Equation  (7).

Since we first assumed $\lambda >0$, we obtain $R(W^{(\lambda )},{\theta }^{(\lambda )};Y)\le R(W^{(0)},{\theta }^{(0)};Y).$ Finally, if $\lambda =0$, then there is equality with $R(W^{(\lambda )},{\theta }^{(\lambda )};Y)=R(W^{(0)},{\theta }^{(0)};Y).$    □

Proof of Theorem 1.

Note that if ${\lambda }_1=0$, then Theorem 1 already applies, so we assume $0<{\lambda }_1<{\lambda }_2$. We have that

$$\begin{array}{ccc}\hfill & F^{({\lambda }_1)}(W^{({\lambda }_1)},H^{({\lambda }_1)},{\theta }^{({\lambda }_1)};X,Y)\le F^{({\lambda }_1)}(W^{({\lambda }_2)},H^{({\lambda }_2)},{\theta }^{({\lambda }_2)};X,Y)\Rightarrow \hfill & \hfill \\ \hfill & {\lambda }_1\left(R(W^{({\lambda }_1)},{\theta }^{({\lambda }_1)};Y)-R(W^{({\lambda }_2)},{\theta }^{({\lambda }_2)};Y)\right)\le \hfill & \hfill \\ \hfill & N(W^{({\lambda }_2)},H^{({\lambda }_2)};X)-N(W^{({\lambda }_1)},H^{({\lambda }_1)};X).\hfill & \hfill \end{array}$$

(10)

We also have

$$\begin{array}{cc}\hfill {\lambda }_2(R(W^{({\lambda }_2)},{\theta }^{({\lambda }_2)};Y)-R(W^{({\lambda }_1)},{\theta }^{({\lambda }_1)};Y))& \le \hfill \\ & \hfill N(W^{({\lambda }_1)},H^{({\lambda }_1)};X)-N(W^{({\lambda }_2)},H^{({\lambda }_2)};X).\end{array}$$

(11)

Adding Equations (10) and (11) together,

$$({\lambda }_1-{\lambda }_2)R(W^{({\lambda }_1)},{\theta }^{({\lambda }_1)};Y)+({\lambda }_2-{\lambda }_1)R(W^{({\lambda }_2)},{\theta }^{({\lambda }_2)};Y)\le 0$$

which, upon dividing by ${\lambda }_2-{\lambda }_1>0$, directly gives

$$R(W^{({\lambda }_2)},{\theta }^{({\lambda }_2)};Y)\le R(W^{({\lambda }_1)},{\theta }^{({\lambda }_1)};Y).$$

   □

3.3. Algorithm

Our minimization approach is iterative and based on the alternating nonnegative least squares [16] approach. Due to the coupling of NMF and regression errors, other approaches such as multiplicative or additive updates [17] are less natural. Each iteration consists of: (1) holding H and $\theta$ fixed while optimizing each row of W separately (nonnegative least squares); (2) holding W and $\theta$ fixed while optimizing each column of H separately (nonnegative least squares); and finally (3) holding W and H fixed while optimizing over $\theta$. The error is nondecreasing between iterations and from one optimization to the next. We now derive and justify this approach (Algorithm 1) in increasing complexity of cases.  

Algorithm 1: Overall CSSNMF algorithm.
	Input   :A matrix $X\in {\mathbb{R}}_{\ge 0}^{n\times m}$,
	a vector $Y\in {\mathbb{R}}^{n\times 1}$,
	a positive integer $r\in \mathbb{N}$,
	a scalar $\lambda \ge 0$,
	a relative error tolerance $\tau >0$, and
	a maximum number of iterations $maxIter$.
	Output:Minimizers of Equations (6)–(9): nonnegative matrix $W\in {\mathbb{R}}_{\ge 0}^{n\times r}$,
	nonnegative matrix $H\in {\mathbb{R}}_{\ge 0}^{r\times m}$, and
	vector $\theta \in {\mathbb{R}}^{(r+1)\times 1}.$
1	$relErr=\infty ,err=\infty$
2	Elementwise, $W\sim {Unif([0,\left|\right|X\left|\right|}_{\infty }))$, $H\sim {Unif([0,\left|\right|X\left|\right|}_{\infty }))$,
	$\theta \sim {Unif([0,\left|\right|X\left|\right|}_{\infty })).$
3	$iter=0$
4	while $relErr>\tau$ and $iter<maxIter$ do
5			$W\leftarrow newW$ as per Algorithm 2
6			$H\leftarrow newH$ as per Algorithm 3
7			$\theta \leftarrow new\theta$ as per Algorithm 4
8			Normalize W, H, and $\theta$ as per Algorithm 5
9			$errTemp=F^{(\lambda )}(W,H,\theta ;X,Y)$
10			if $err<\infty$ then
11					$relErr\leftarrow |err-errTemp|/err$
12			end if
13			$err\leftarrow errTemp$
14			$iter\leftarrow iter+1$
15	end while$\phantom{(}$
16	return $W,H,\theta$

Algorithm 2: Updating W.

Algorithm 3: Updating H.

Algorithm 4: Updating $\theta$.
	Input	:A vector $Y\in {\mathbb{R}}^{n\times 1}$, and
		a matrix $W\in {\mathbb{R}}_{\ge 0}^{n\times r}$
	Output	:A new value for $\theta$.
1	$e={(1,1,\dots .,1)}^{\mathrm{T}}\in {\mathbb{R}}^{n\times 1}$
2	$\overline{W}=\left[\begin{array}{ccc}e& |& W\end{array}\right]$
3	return ${\overline{W}}^{+}Y$

Algorithm 5: Normalization process.
	Input	:A matrix $W\in {\mathbb{R}}_{\ge 0}^{n\times r}$,
		a matrix $H\in {\mathbb{R}}_{\ge 0}^{r\times m}$, and
		a vector $\theta \in {\mathbb{R}}^{(r+1)\times 1}.$
	Output	:New values for W, H, and $\theta$.
1	$S\in {\mathbb{R}}_{\ge 0}^{r\times 1}$ a vector of row sums of H.
2	$S\leftarrow \mathrm{diag}(S)$
3	$W\leftarrow WS$.
4	$H\leftarrow S^{-1}H$.
5	${\theta }_{2:(r+1)}\leftarrow S^{-1}{\theta }_{2:(r+1)}.$
6	return W, H, and $\theta$.

• W and H fixed. 

If W and H are given and only $\theta$ can vary, then Equation (2) is minimized when $\left|\right|\overline{W}\theta -Y{\left|\right|}^2$ is minimized. If $\overline{W}$ has full rank or is overdetermined, this happens when the error, $\overline{W}\theta -Y$, is orthogonal to the column span of $\overline{W}$ or that

$$\theta ={({\overline{W}}^{\mathrm{T}}\overline{W})}^{-1}{\overline{W}}^{\mathrm{T}}Y.$$

(12)

When $\overline{W}$ is underdetermined or has full rank, we require $\overline{W}\theta =Y$ with $\left|\right|\theta \left|\right|$ minimized (for uniqueness) whereby

$$\theta ={\overline{W}}^{\mathrm{T}}{(\overline{W}{\overline{W}}^{\mathrm{T}})}^{-1}Y.$$

(13)

See Algorithm 4. Thus, when $\overline{W}$ does not have full rank, the solution is $\theta ={\overline{W}}^{+}Y$ where ${\overline{W}}^{+}$ is the pseudoinverse of $\overline{W}$. See [18] for a discussion of pseudoinverses. In applications, we only expect to see overdetermined systems because the number of topics r should be less than the number of documents $n.$ 

• $\mathbf{W}$ and $\theta$ fixed.  

If W and $\theta$ are given and only H can change, then minimizing Equation (2) requires minimizing $N(W,H;X)$. We can expand this error term out in the columns of H:

$$\begin{array}{cc}\hfill N(W,H;X)& =\left|\right|X-{WH\left|\right|}_F^2\hfill \\ \hfill & =\sum _{j=1}^m{\left|\right|(X-WH)}_{:,j}{\left|\right|}^2\hfill \\ \hfill & =\sum _{j=1}^m\left|\right|X_{:,j}-W{H}_{:,j}{\left|\right|}^2.\hfill \end{array}$$

Because columnwise the terms of the sum are independent, we can minimize each column $H_{:,j}$ of H separately to minimize the sum, i.e.,

$$\begin{array}{c}\hfill H_{:,j}=\arg {\min }_{h\in {\mathbb{R}}_{\ge 0}^{r\times 1}}\left|\right|X_{:,j}-Wh{\left|\right|}^2,j=1,2,\dots ,m,\end{array}$$

(14)

as given in Algorithm 3. 

• $\mathbf{H}$ and $\theta$ fixed.  

When H and $\theta$ are fixed, then Equation (2) can be written out as

$$\begin{array}{cc}\hfill F^{(\lambda )}& =\left|\right|X-{WH\left|\right|}_F^2+\lambda \left|\right|\overline{W}\theta -Y{\left|\right|}^2\hfill \\ \hfill & =\sum _{i=1}^n{\left|\right|(X-WH)}_{i,:}{\left|\right|}^2+\lambda \sum _{i=1}^n{(\overline{W}\theta -Y)}_i^2\hfill \\ \hfill & =\sum _{i=1}^n\left|\right|X_{i,:}-W_{i,:}H{\left|\right|}^2+\sum _{i=1}^n{\left(\sqrt{\lambda }({\theta }_{1}e+W_{i,:}\overline{\theta }-Y)\right)}_i^2\hfill \end{array}$$

(15)

where $e={(1,1,\dots ,1)}^{\mathrm{T}}\in {\mathbb{R}}^{n\times 1}$ and $\overline{\theta }={({\theta }_2,\dots ,{\theta }_{r+1})}^{\mathrm{T}}\in {\mathbb{R}}^{r\times 1}$. Defining matrices

$$\begin{array}{cc}\hfill \overline{X}& =\left[\begin{array}{ccc}X& |& \sqrt{\lambda }({\theta }_{1}e-Y)\end{array}\right]\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \overline{H}& =\left[\begin{array}{ccc}H& |& \sqrt{\lambda }\overline{\theta }\end{array}\right]\hfill \end{array}$$

(17)

we can rewrite Equation (15) as

$$\begin{array}{cc}\hfill F^{(\lambda )}& =\sum _{i=1}^n\left|\right|{\overline{X}}_{i,:}-W_{i,:}\overline{H}{\left|\right|}^2,\hfill \end{array}$$

which can be minimized through

$$\begin{array}{c}\hfill W_{i,:}=\arg {\min }_{w\in {\mathbb{R}}_{\ge 0}^{1\times r}}\left|\right|{\overline{X}}_{i,:}-w\overline{H}{\left|\right|}^2,i=1,2,\dots ,n.\end{array}$$

(18)

This is precisely Algorithm 2.

For our optimizations and linear algebra, we used Numpy [19] and SciPy [20]. The nonnegative least squares routine employs an active set method to solve the least squares minimization problem with inequality constraints [21]. The active set method amounts to gradually building up the set of active constraints (those for which their not being enforced in the unconstrained problem would result in constraint violation or equality, i.e., a regression variable with a nonnegativity constraint being less than or equal to 0) and then optimizing over the passive set (all variables not in the active set) once identified with equality constraints on the active set [22]. This method can also be parallelized [23].

In addition to the steps outlined within these algorithms, we employed two additional modifications: (1) we defined $ϵ={10}^{-10}$, and any entries in H less than $ϵ$ were replaced by $ϵ$ (otherwise on some occasions, the W update step would fail); and (2) the minimizations at times yielded worse objective errors than already obtained and, when this happened, we did not update to the worse value.

As noted with other NMF routines, we might not reach a global minimizer [24]. In practice, the minimization should be run repeatedly with different random initializations to find a more ideal local minimum.

From an application standpoint, we wish to run the model on documents it has not been trained on. Algorithm 6 stipulates how a prediction takes place. We first find the best nonnegative decomposition of the document, a vector in ${\mathbb{R}}^{1\times m}$, into the topic basis, projecting to $r-$dimensions. With the representation in topic-coordinates, we then use the linear model.   

Algorithm 6: Prediction process.
	Input	:A matrix $H\in {\mathbb{R}}_{\ge 0}^{r\times m}$,
		a vector $\theta \in {\mathbb{R}}^{(r+1)\times 1}$,
		and a vector $x\in {\mathbb{R}}^{1\times m}$.
	Output	:Model prediction for response variable, $\widehat{y}$.
1	Compute $w=\arg {\min }_{w\in {\mathbb{R}}_{\ge 0}^{1\times r}}\left|\right|wH-{x\left|\right|}^2$.
2	Compute $\widehat{y}={\theta }_1+w{\theta }_{2:(r+1)}$.
3	return $\widehat{y}$

An implementation of our algorithm can be found on our BitBucket repository, accessed on 21 February 2023.

4. Synthetic Datasets

In our synthetic data, we generate a matrix X that has nonnegative factors W and H, but we add noise. We also generate a response vector Y given as the matrix–vector product $\overline{W}\theta$ with noise. We investigate four items: (1) that the method does in fact work to decrease the objective function; (2) whether the regression errors decrease with increasing $\lambda$; (3) the effects of overfitting; and (4) the model robustness to noise.

4.1. Generating Synthetic Data

Our synthetic data generation can be summarized as follows:

• We fix values of $n=100,m=40$, $M=20$, and $r=4.$
• We then define ${\eta }_x={\eta }_y=4.$
• We pick $X\in {\mathbb{R}}^{n\times r}$ such that each entry is $\sim Unif([0,M))$. We likewise choose $H\in {\mathbb{R}}^{r\times m}.$
• We set $X=WH$.
• We pick $\theta \in {\mathbb{R}}^{(r+1)\times 1}$ such that each element is $\sim Unif([-M/2,M/2)).$
• We set $Y=\overline{W}\theta .$
• We perturb X with noise $\sim {\mathcal{D}}_X$ and Y with noise $\sim {\mathcal{D}}_Y$.
• Any negative X-entries are set to 0.

We consider two different forms for ${\mathcal{D}}_X$ and ${\mathcal{D}}_Y$:

• Being elementwise $\sim \mathcal{N}(0,{\eta }_x^2)$ and $\sim \mathcal{N}(0,{\eta }_y^2)$ or
• Being elementwise $\sim Unif([0,{\eta }_x))$ and $\sim Unif([0,{\eta }_y))$.

Note that in the synthetic data, the true number of topics is $r=4$. In testing our synthetic data, we run Algorithm 1 where $\tau ={10}^{-4}$ and $maxIter=100$. We use $70\%$ of the data for training and $30\%$ for testing.

4.2. Investigation

We confirm that the error in the objective function $F^{(\lambda )}$ decreases with each iteration of Algorithm 1 in Figure 1—done with Gaussian noise.

With the regression error being the mean squared prediction error, from Figure 2, we see the regression error in the training does tend to decrease with $\lambda$. (There are a few small exceptions, which we believe stem from randomizations leading to an assortment of different local optima.) The overall scale of the testing errors gets smaller as r goes from 1 to 4 and then stays steady or even gets slightly worse as r increases from 4. Indeed, $r=4$ is the “correct” synthetic value. Given the noise as either Gaussian or uniform, the variances of $\mathcal{N}(0,{\eta }_y^2)$, ${\eta }_y^2$, and $Unif([0,{\eta }_y))$, ${\eta }_y^2/12$, serve as loose estimates for the best possible testing loss (the loss could very well be higher as noise is added to the matrix X as well). When the training errors are smaller than this estimate, it suggests overfitting.

To study robustness with noise, we allow the level Gaussian noise to vary in the problem and evaluate the regression error on testing data. To ensure each noise level starts with the same ground truth, we start with unperturbed X and Y (as per Section 4.1 with ${\eta }_x={\eta }_y=0$) and compute matrices and vectors with the same size as X and Y with elementwise $\mathcal{N}(0,1)$ entries. Then, for each level of noise under investigation, we scale these unit normal distributions by a noise level $\eta \in \{0,4,8,12,16,20\}$ and examine the testing error as the $\lambda$ varies with $r=3$ (rank below true rank), $r=4$ (correct rank), and $r=5$ (number of chosen topics above correct rank). Figure 3 illustrates the results. The noise is handled well with $r=3$ in that a higher $\lambda$ does tend to improve the testing error; however, at higher noise levels $\eta$, it seems suitable minima are harder to find when $\lambda$ is large. With $r=4$ and $r=5$, the testing does not benefit with increasing $\lambda$ at any noise level.

Taken together, we anticipate that CSSNMF will perform well provided the number of topics chosen does not exceed the true number of topics in the dataset (difficult to assess). We expect that the optimal predictions on unseen data should occur at a $\lambda$ large enough that the testing errors have decreased and plateaued. In Figure 2, we see that for large $\lambda$, when overfitting is an issue, the testing performance is seldom better than where $\lambda =0$ (classical NMF and then regression) and, in fact, is often much worse.

5. Rate My Professors Dataset

Here we study our method on real data [15] coming from Rate My Professors where for each instructor in the dataset, all corresponding student comments are combined to generate a written narrative and we have the instructor’s rating (mean of all responses).

5.1. Pre-Processing

The corpus was first processed via term frequency–inverse document frequency (TF-IDF) [25] with the TfidfVectorizer class in Python’s scikit learn package [26]. We used arguments min_df=0.01, max_df=0.15, stop_words=’english’, norm=‘l1’, lowercase=True. We found the ratings were not balanced: there were 57 on the interval [1, 2), 235 on the interval $[2,3)$, 494 on the interval [3, 4), and 629 on the interval [4, 5]. To balance the dataset, we extracted only a random subset of 57 reviews in each interval (all ratings on $[1,2)$ were used). Overall, we obtained a corpus matrix X that was $228\times 1635$. The open right-end of the intervals ensures data are not duplicated.

5.2. Choice of Topic Number and Regression Weight

We did not know the true number of topics in the dataset and chose topics of $r=1,3,5,7,9$, and 11, with $\lambda \in \{0\}\cup \{{10}^{2i/3}|i\in [-12,0]\cap \mathbb{Z}\}$. We present the results for 11 topics that gave the best results. See Figure 4. We note that, for large enough $\lambda$, the testing error outperforms the testing error for $\lambda =0$. The optimal point was at $\lambda ={10}^{-2/3}\approx 0.215$.

We comment that it is generally difficult to know precisely where the testing error will be minimized, only that, based on observations of the synthetic data, the testing error is often better than the $\lambda =0$ case after the training error has dropped. We speculate that the level of noise in this dataset results in the testing errors not dropping below $\approx 0.75$.

5.3. Prediction

We examine the rating prediction by plotting histograms of predicted ratings where the true ratings were in $[1,2],[2,3],[3,4],$ and $[4,5]$—the closed intervals are used here. Figure 5 depicts these histograms along with the mean predicted rating and true rating. The predictions are often within range, and the mean predicted values are very close to the true means over each interval. We can also see the general predictive strength in the scatterplot of actual vs. predicted ratings in Figure 6.

These results suggest the model is able to identify topics and associated $\theta$-weights so as to generate predictions that are consistent with true ratings. For example, in the case where ratings are in $[1,2]$, we see the peak of the predictions is around 2, not exceeding 4, with some predictions as low as $-2$; then, in the case of ratings in $[4,5]$, the model peaks around $3.5$ and makes some predictions above 7. There is a clear capacity for the topics to shift the predictions.

5.4. Topics Identified

It is important that the method not only have predictive power, but also produce interpretable topics. We now look at the 11 topics with their associated $\theta$-weights. We find

$$\begin{array}{c}\theta =(2.39909812,2.82948873,-2.21028471,1.83876976,-4.77504984,\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} -3.86467795,3.46353642,0.03914383,3.26619842,-5.51595505,\hfill \\ \hfill {\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} 4.15317532,3.90733652)}^{\mathrm{T}}.\end{array}$$

Note that ${\theta }_1\approx 2.4$ suggests that for a set of reviews with no topics, the average rating would be around $2.4$—this suggests it is the presence of positive/negative topics that raise/lower the rating.

In Figure 7 and Figure 8, we plot the words in the topics associated with positive and negative ratings. The topics are interpretable. For positive topics, we find Topic 10 (suggestive of extra credit), Topic 6 (suggestive of being inspiring), Topic 8 (suggestive of being approachable/brilliant), and Topic 11 (suggestive of being nice/enjoyable class) and words such as “recommend” in a couple of them. A few words seem out of place, such as “hate” in Topic 11, but that can be explained by some positive reviews having phrases such as “i hated chemistry in high school and after taking her class i don t [sic] hate chem as much.” Among the negative topics, we see Topic 4 (being horrible) and Topic 5 (being unfair).

As a whole, the topics are consistent with intuitive notions of what would be associated with higher or lower ratings. It is also interesting to look at the $\theta$-topic weights quantitatively. For example, Topic 2 (mentioning rants and sarcasm) and Topic 4 (suggestive of being harder and failing students) contribute negatively to the score, but being a harder teacher seems to contribute more negatively to the rating than ranting. To elaborate more: because the row sums of the corpus matrix X and those of the topic matrix H sum to 1, we have that the sum of the topic weights for each document are approximately 1 as per Remark 1. All topics exist on the “same scale” within a document (roughly in $[0,1]$) and must roughly sum to 1; therefore, the sizes of the weights in $\theta$ for each topic can be ordered by their positive/negative contributions.

The fact that instructor difficulty has a negative effect on rating and the easiness has a positive effect has been found in another study [27] that looked at Rate My Professors data on instructor easiness and overall quality ratings. It is also interesting that many of the “dimensions” detected through our study, such as being approachable/nice (niceness), being hard/easy (difficulty), and being inspiring, were dimensions naturally identified by other scholars [28] who analyzed Rate My Professors data by hand through reading comments and classifying key phrases within the comments. In this latter study, however, each dimension could be positive or negative, depending on how it was used.

6. Conclusions and Future Work

We have developed CSSNMF as a means to combine NMF with regression on a continuous response variable. We accomplished this by minimizing an objective function that combines an NMF error with a weighted regression error. We have shown that the regression error in training is weakly decreasing with the regression error weight and that, in practical applications, the error strictly decreases. The topics identified can outperform the quantitative accuracy of topics formed through NMF alone while retaining a high degree of interpretability (as found with real data). The method is robust to noise and tends to perform best on new data when there is a dimensionality reduction, i.e., when the number of topics fit for is fewer than the true number of topics.

We expect that our methods can be applied to very large datasets as our algorithmic steps involve hierarchical alternating nonnegative least squares [16], which could even be handled in parallel [29] (with or without parallel nonnegative least squares), and a least squares problem (Equation (12)) with a small $r\times r$ matrix to invert.

Although our analysis focused on the case of linear regression, incorporating nonlinearities would be of interest. We also noted the challenge in choosing the appropriate $\lambda$ given only training data. A more theoretical understanding of when testing errors drop substantially could be explored, but this may be dataset-specific.