# Continuous Semi-Supervised Nonnegative Matrix Factorization

Michael R. Lindstrom The University of Texas Rio Grande Valley, School of
Mathematical and Statistical Sciences<EMAIL_ADDRESS>Xiaofu Ding
University of California Los Angeles Department of Mathematics Feng Liu
22footnotemark: 2 Anand Somayajula 22footnotemark: 2 Deanna Needell
22footnotemark: 2

###### 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 this paper, we show this factorization can be
combined with regression on a continuous response variable. In practice, the
method performs better than regression done after topics are identified and
retrains interpretability.

Keywords: topic modelling, regression, nonnegative matrix factorization,
optimization

##  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. More generally, it can be used to decompose
an image into identifiable patterns [1] and as a general-purpose
dimensionality reduction or preprocessing method before applying other machine
learning methods as has been done in studying various diseases [2, 3]. Like
Singular Value Decomposition (SVD) [4], NMF provides a low rank factorization.
In NMF, a nonnegative matrix $X\in\mathbb{R}^{n\times m}_{\geq 0}$
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}^{n\times r}_{\geq 0}$ and $H\in\mathbb{R}^{r\times m}_{\geq
0}.$ This non-negativity enforces that the data in $X$ is represented by a
non-negative combination of the dictionary atoms in the factorization, which
lends itself to human interpretability. For example, in the foundational work
[1], Lee and Seung show that NMF when applied to facial images decompose the
images into recognizable parts such as noses, eyes, and mouths.

When applied to a document-term matrix $X$ where row $i$ of $X$ represents
document $i$ and column $j$ represents the frequency of word $j$, the
classical NMF method amounts to

$(W,H)=\text{arg\,min}_{W\in\mathbb{R}^{n\times r}_{\geq
0},H\in\mathbb{R}^{r\times m}_{\geq 0}}||X-WH||_{F}^{2}$ (1)

where the $F$-subscript denotes the Frobenius norm with
$||A||_{F}^{2}=\mathrm{Tr}(A^{T}A)=\sum_{i=1}^{n}\sum_{j=1}^{m}|A_{ij}|^{2}$.
Other variations on the penalty norm exist including the Kulblack-Liebler
Divergence [5]. After computing $W$ and $H$, we interpret row $j$ of $H$ as
the $j$th topic, its components being the weight of each word in that topic,
and row $i$ of $W$ as the topic-encoding of document $j$, 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.

Prior authors have combined NMF with a linear regression procedure to maximize
the predictive power of a classifier [6, 7, 8, 9]. 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 [10].

In this paper, we combine this NMF with a linear regression model to predict
the value of a continuous response variable. We consider datasets that pair
written commentary with a real-valued observation. As our motivating example,
we consider Rate My Professor reviews [11] 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.

The rest of our paper is organized as follows: in section 2, we provide the
formulation of our method, its algorithmic implementation, and its theoretical
properties; proofs of the properties are given in section 3; in section 4, we
provide a proof of concept through synthetic data; in section 5, we test our
method on the Rate My Professor dataset; and finally we conclude our work in
section 6.

##  2 Model

We provide the framework for our proposed Continuous Semi-Supervised
Nonnegative Matrix Factorization method (CSSNMF).

###  2.1 Formulation

We consider having a document-term corpus $X\in\mathbb{R}_{\geq 0}^{n\times
m}$ for $n$ documents with their associated word frequencies in the $m$
columns. Each document has a corresponding value in $\mathbb{R}$ so that we
can associate with $X$ the vector $Y\in\mathbb{R}^{n}$. We choose
$r\in\mathbb{N}$ and $\lambda\geq 0$ as hyper-parameters where $r$ denotes the
number of topics and $\lambda$ is a regression weight. Given
$W\in\mathbb{R}_{\geq 0}^{n\times r}$, $H\in\mathbb{R}_{\geq 0}^{r\times m}$,
and $\theta\in\mathbb{R}^{r+1}$, we define the penalty function that combines
topic modelling with a linear regression based on the topic representation.

$\displaystyle F^{(\lambda)}(W,H,\theta;X,Y)$ $\displaystyle=N(W,H;X)+\lambda
R(W,\theta;Y),\text{ where }$ (2) $\displaystyle N(W,H;X)$
$\displaystyle:=||X-WH||_{F}^{2}$ (3) $\displaystyle R(W,\theta;Y)$
$\displaystyle:=||\tilde{W}\theta-Y||^{2}$ (4)

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

$\bar{W}:=\begin{pmatrix}1&W_{1,1}&\ldots&W_{1,m}\\\
1&W_{2,1}&\ldots&W_{2,m}\\\ \vdots&\vdots&\vdots&\vdots\\\
1&W_{n,1}&\ldots&W_{n,m}\end{pmatrix}.$ (5)

The matrix $\bar{W}$ with its column of 1’s allows for an intercept: given a
topic representation $w\in\mathbb{R}^{r}$, we predict a value
$\hat{y}=\theta_{1}+\theta_{2}w_{1}+...+\theta_{r+1}w_{r}.$

When $\lambda>0$, we seek

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

And when $\lambda=0$, we define

$\displaystyle(W^{(0)},H^{(0)})$ $\displaystyle=\text{arg\,min}_{W,H}N(W,H;X)$
(7) $\displaystyle\theta^{(0)}$
$\displaystyle=\text{arg\,min}_{\theta}R(W^{(0)},\theta;Y).$ (8)

We also impose a normalization constraint, that

$\forall i,\quad\sum_{j=1}^{m}H_{ij}=1$ (9)

so that the topics have unit length in $\ell_{1}.$ If $X\approx WH$ is
normalized so its rows sum to $1$ then it is also the case that $\forall
i,\quad\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 through $H$ to be less
meaningful. In section 2.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}^{r}$
for the document and using that encoding in the linear model — see section
2.3.

######  Remark 1 (Uniqueness).

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

###  2.2 Theoretical Results

We present two important behaviours of CSSNMF with regards to increasing
$\lambda$ and its effect upon predicting the response variable. The proofs are
contained in section 3.

######  Proposition 1 (Regression Error with Nonzero $\lambda$).

For $\lambda\geq 0$, let $W^{(\lambda)},H^{(\lambda)},\theta^{(\lambda)}$ be a
unique (as per Remark 1) global minimum to Eqs. (6)-(9) . Then
$R(W^{(\lambda)},\theta^{(\lambda)})\leq R(W^{(0)},\theta^{(0})$.

######  Theorem 1 (Weakly Decreasing Regression Error).

Let $0\leq\lambda_{1}<\lambda_{2}$ be given where
$W^{(\lambda_{i})},H^{(\lambda_{i})},\theta^{(\lambda_{i})}$ are the unique
(as per Remark 1) global minimizers of Eqs. (6)-(9) for $i=1,2$. Then
$R(W^{(\lambda_{2})},\theta^{(\lambda_{2})};Y)\leq
R(W^{(\lambda_{1})},\theta^{(\lambda_{1})};Y).$

######  Remark 2.

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

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 in fact
strictly monotonically decreases.

###  2.3 Algorithm

Our minimization approach is iterative and based on the alternating
nonnegative least squares [12] approach. Due to the coupling of NMF and
Regression errors, other approaches such as multiplicative or additive updates
[13] are less natural. Each iteration consists of: (1) holding $H$ and
$\theta$ fixed while optimizing each row of $W$ separately (non-negative least
squares); (2) holding $W$ and $\theta$ fixed while optimizing each column of
$H$ separately (non-negative 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.

####  $W$ and $H$ fixed.

If $W$ and $H$ are given and only $\theta$ can vary then Eq. (2) is minimized
when $||\bar{W}\theta-Y||^{2}$ is minimized. This happens when the error,
$\bar{W}\theta-Y$, is orthogonal to the column span of $W$ or that

$\theta=(\bar{W}^{T}\bar{W})^{-1}\bar{W}^{T}Y.$ (10)

See Algorithm 4. When $\bar{W}$ does not have full rank, we interpret
$(\bar{W}^{T}\bar{W})^{-1}$ as a pseudo-inverse.

####  $W$ and $\theta$ fixed.

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

$\displaystyle N(W,H;X)$ $\displaystyle=||X-WH||_{F}^{2}$
$\displaystyle=\sum_{j=1}^{m}||(X-WH)_{:,j}||^{2}$
$\displaystyle=\sum_{j=1}^{m}||X_{:,j}-WH_{:,j}||^{2}.$

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

$\displaystyle H_{:,j}=\text{arg\,min}_{h\in\mathbb{R}_{\geq 0}^{r\times
1}}||\bar{X}_{:,j}-Wh||^{2},\quad j=1,2,...,m,$ (11)

as given in Algorithm 3.

####  $H$ and $\theta$ fixed.

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

$\displaystyle F^{(\lambda)}$
$\displaystyle=||X-WH||_{F}^{2}+\lambda||\bar{W}\theta-Y||^{2}$
$\displaystyle=\sum_{i=1}^{n}||(X-WH)_{i,:}||^{2}+\lambda\sum_{i=1}^{n}(\bar{W}\theta-Y)_{i}^{2}$
$\displaystyle=\sum_{i=1}^{n}||X_{i,:}-W_{i,:}H||^{2}+\sum_{i=1}^{n}\left(\sqrt{\lambda}(\theta_{0}e+W_{i,:}\bar{\theta}-Y)\right)_{i}^{2}$
(12)

where $e=(1,1,...,1)^{T}\in\mathbb{R}^{n}$ and
$\bar{\theta}=(\theta_{2},...,\theta_{r+1})^{T}\in\mathbb{R}^{r}$. Defining
matrices

$\displaystyle\bar{X}$
$\displaystyle=\begin{bmatrix}X&|&\sqrt{\lambda}(\theta_{1}e-Y)\end{bmatrix}$
(13) $\displaystyle\bar{H}$
$\displaystyle=\begin{bmatrix}H&|&\sqrt{\lambda}\bar{\theta}\end{bmatrix}$
(14)

we can rewrite (12) as

$\displaystyle F^{(\lambda)}$
$\displaystyle=\sum_{i=1}^{n}||\bar{X}_{i,:}-W_{i,:}\bar{H}||^{2},$

which can be minimized through

$\displaystyle W_{i,:}=\text{arg\,min}_{w\in\mathbb{R}_{\geq 0}^{1\times
r}}||\bar{X}_{i,:}-w\bar{H}||^{2},\quad i=1,2,...,n.$ (15)

This is precisely Algorithm 2.

Input : A matrix $X\in\mathbb{R}_{\geq 0}^{n\times m}$,

a vector $Y\in\mathbb{R}^{n}$,

a positive integer $r\in\mathbb{N}$,

a scalar $\lambda\geq 0$,

a relative error tolerance $\tau>0$, and

a maximum number of iterations $maxIter$.

Output : Minimizers of Eqs. (6)-(9) : nonnegative matrix $W\in\mathbb{R}_{\geq
0}^{n\times r}$,

nonnegative matrix $H\in\mathbb{R}_{\geq 0}^{r\times m}$, and

vector $\theta\in\mathbb{R}^{r+1}.$

1 $relErr=\infty,err=\infty$

2 Elementwise, $W\sim Unif([0,||X||_{\infty}))$, $H\sim
Unif([0,||X||_{\infty}))$, $\theta\sim Unif([0,||X||_{\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

return _$W,H,\theta$_

Algorithm 1 Overall CSSNMF algorithm.

Input : A matrix $X\in\mathbb{R}_{\geq 0}^{n\times m}$,

a vector $Y\in\mathbb{R}^{n}$,

a matrix $W\in\mathbb{R}_{\geq 0}^{n\times r}$,

a matrix $H\in\mathbb{R}_{\geq 0}^{r\times m}$,

and a scalar $\lambda\geq 0.$

Output : A new value for $W$.

1 $\bar{\theta}=(\theta_{2},...,\theta_{r+1})^{T}$

2 $\bar{X}=\begin{bmatrix}X&|&\sqrt{\lambda}(\theta_{1}-Y)\end{bmatrix}$

3 $\bar{W}=\begin{bmatrix}H&|&\sqrt{\lambda}\bar{\theta}\end{bmatrix}$

4 for  _$i\leftarrow 1...n$_ do

5 $W_{i,:}\leftarrow\text{arg\,min}_{w\in\mathbb{R}_{\geq 0}^{1\times
r}}||\bar{X}_{i,:}-w\bar{H}||^{2}$

6 end for

return _$W$_

Algorithm 2 Updating $W$.

Input : A matrix $X\in\mathbb{R}_{\geq 0}^{n\times m}$,

a matrix $W\in\mathbb{R}_{\geq 0}^{n\times r}$, and

a matrix $H\in\mathbb{R}_{\geq 0}^{r\times m}.$

Output : A new value for $H$.

1 for  _$j\leftarrow 1...m$_ do

2 $H_{:,j}\leftarrow\text{arg\,min}_{h\in\mathbb{R}_{\geq 0}^{r\times
m}}||\bar{X}_{:,j}-Wh||^{2}$

3 end for

return _$H$_

Algorithm 3 Updating $H$.

Input : A vector $Y\in\mathbb{R}^{n}$, and

a matrix $W\in\mathbb{R}_{\geq 0}^{n\times r}$

Output : A new value for $\theta$.

1 $e=(1,1,...,1)^{T}\in\mathbb{R}^{n\times 1}$

2 $\bar{W}=\begin{bmatrix}e&|&W\end{bmatrix}$

return _$(\bar{W}^{T}\bar{W})^{-1}\bar{W}^{T}Y$_

Algorithm 4 Updating $\theta$.

Input : A matrix $W\in\mathbb{R}_{\geq 0}^{n\times r}$,

a matrix $H\in\mathbb{R}_{\geq 0}^{r\times m}$, and

a vector $\theta\in\mathbb{R}^{r+1}.$

Output : New values for $W$, $H$, and $\theta$.

1 $S\in\mathbb{R}_{\geq 0}^{r}$ 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}.$

return _$W$ , $H$, and $\theta$._

Algorithm 5 Normalization process.

Input : A matrix $H\in\mathbb{R}_{\geq 0}^{r\times m}$,

a vector $\theta\in\mathbb{R}^{r+1}$,

and a vector $x\in\mathbb{R}^{1\times m}$.

Output : Model prediction for response variable, $\hat{y}$.

1 Compute $w=\text{arg\,min}_{w\in\mathbb{R}_{\geq 0}^{1\times
r}}||wH-x||^{2}$.

2 Compute $\hat{y}=\theta_{1}+w\theta_{2:r+1}$.

return _$\hat{y}$_

Algorithm 6 Prediction process.

For our optimizations and linear algebra, we used Numpy [14] and SciPy [15].
Besides the steps outlined within these algorithms, we employed two additional
modifications: (1) we defined $\epsilon=10^{-10}$ and any entries in the $H$
less than $\epsilon$ were replaced by $\epsilon$ ( 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 [16].
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}^{m}$, into the topic basis, projecting to $r-$dimensions. With the
representation in topic-coordinates, we then use the linear model.

##  3 Proofs

Before proceeding to practical applications, we prove Proposition 1 and
Theorem 1.

###### Proof of Proposition 1.

If $\lambda>0$ then

$\displaystyle
F^{(\lambda)}(W^{(\lambda)},H^{(\lambda)},\theta^{(\lambda)};X,Y)$
$\displaystyle\leq F^{(\lambda)}(W^{(0)},H^{(0)},\theta^{(0)};X,Y)\implies$
$\displaystyle N(W^{(\lambda)},H^{(\lambda)};X)+\lambda
R(W^{(\lambda)},\theta^{(\lambda)};Y)$ $\displaystyle\leq
N(W^{(0)},H^{(0)};X)+\lambda R(W^{(0)},\theta^{(0)};Y)$
$\displaystyle\implies$
$\displaystyle\lambda(R(W^{(\lambda)},\theta^{(\lambda)};Y)-R(W^{(0)},\theta^{(0)};Y))$
$\displaystyle\leq N(W^{(0)},H^{(0)};X)-N(W^{(\lambda)},H^{(\lambda)};X)$
$\displaystyle\leq 0.$

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

Since we first assumed $\lambda>0$, we obtain
$R(W^{(\lambda)},\theta^{(\lambda)};Y)\leq 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

$\displaystyle
F^{(\lambda_{1})}(W^{(\lambda_{1})},H^{(\lambda_{1})},\theta^{(\lambda_{1})};X,Y)\leq
F^{(\lambda_{1})}(W^{(\lambda_{2})},H^{(\lambda_{2})},\theta^{(\lambda_{2})};X,Y)\implies$
$\displaystyle\lambda_{1}\left(R(W^{(\lambda_{1})},\theta^{(\lambda_{1})};Y)-R(W^{(\lambda_{2})},\theta^{(\lambda_{2})};Y)\right)\leq$
$\displaystyle
N(W^{(\lambda_{2})},H^{(\lambda_{2})};X)-N(W^{(\lambda_{1})},H^{(\lambda_{1})};X).$
(16)

We also have

$\lambda_{2}\left(R(W^{(\lambda_{2})},\theta^{(\lambda_{2})};Y)-R(W^{(\lambda_{1})},\theta^{(\lambda_{1})};Y)\right)\leq\\\
N(W^{(\lambda_{1})},H^{(\lambda_{1})};X)-N(W^{(\lambda_{2})},H^{(\lambda_{2})};X).$
(17)

Adding Eqs. (16) and (17) 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)\leq
0$

which, upon dividing by $\lambda_{2}-\lambda_{1}>0$, directly gives

$R(W^{(\lambda_{2})},\theta^{(\lambda_{2})};Y)\leq
R(W^{(\lambda_{1})},\theta^{(\lambda_{1})};Y).$

∎

##  4 Synthetic Datasets

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

###  4.1 Generating Synthetic Data

Our synthetic data generation can be summarized as follows:

  1. 1.

We fix values of $n=100,m=40$, $M=20$, and $r=4.$

  2. 2.

We then define $\eta_{x}=\eta_{y}=4.$

  3. 3.

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}.$

  4. 4.

We set $X=WH$.

  5. 5.

We pick $\theta\in\mathbb{R}^{r+1}$ such that each element is $\sim
Unif([-M/2,M/2)).$

  6. 6.

We set $Y=\bar{W}\theta.$

  7. 7.

We perturb $X$ with noise $\mathcal{D}_{X}$ and $Y$ with noise
$\sim\mathcal{D}_{Y}$.

  8. 8.

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 since noise is added to the matrix $X$ as well). When the training
errors are smaller than this estimate it suggests overfitting. Since our model
works with a sum of squared errors, it is expected when the errors are not
Gaussian, that the model will not perform as well. Indeed, there is some
degradation in testing errors in comparing the Gaussian with Uniform noise.

Figure 1: Illustration of decreasing objective function at $r=3$ topics for
$\lambda\in\\{0\\}\cup\\{10^{i/2}|i\in\mathbb{Z}\cap[-2,2]$ }.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

Figure 2: Regression errors with varying regression weight $\lambda$ with
different numbers of topics $r$ for Gaussian noise (a)-(g) and Uniform noise
(h)-(n). The $\lambda$ values used are the set
$\\{0\\}\cup\\{10^{i/2}|i\in\mathbb{Z}\cap[-8,8]\\}.$ For each $\lambda$ and
$r$, fifty trials were run and the regression errors corresponding to the best
overall objective function $F^{(\lambda)}$ were recorded. Points with
$\lambda>0$ for which the regression error exceeds $1.5$ times the regression
error at $\lambda=0$ are not displayed. The dashed horizontal line is the
estimated minimal mean regression error. The dashed vertical line is the
transition point between a linear and logarithmic $x-$scale.

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

###  5.1 Pre-processing

The corpus was first processed via TFIDF [17] with the TfidfVectorizer class
in Python’s scikit learn package [18]. 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 which gave the best results. See Figure 3. 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$.

.

Figure 3: Errors in training and validation on Rate My Professor dataset with
$r=11$ topics. Points with $\lambda>0$ for which the regression error exceeds
$1.5$ times the regression error at $\lambda=0$ are not displayed. The dashed
vertical line is the transition point between a linear and logarithmic
$x-$scale.

###  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 4 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 5.

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 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.

(a)

(b)

(c)

(d)

Figure 4: Histograms of the predicted rating for various ranges of true
ratings. The vertical dashed lines represent the mean values. The predicted
and true means are as follows: $2.206$ and $1.543$ for ratings in $[1,2]$,
$3.233$ and $2.529$ for ratings in $[2,3]$, $3.594$ and $3.593$ for ratings in
$[3,4]$ (the lines are indistinguishable), and $4.576$ and $4.494$ for ratings
in $[4,5]$. Figure 5: Scatterplot of Rate My Professor ratings.

###  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

$\theta=(2.39909812,2.82948873,-2.21028471,1.83876976,-4.77504984,\\\
-3.86467795,3.46353642,0.03914383,3.26619842,-5.51595505,\\\
4.15317532,3.90733652)^{T}.$

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

In Figures 6 and 7, we plot the words in the topics associated with positive
and negative ratings. The topics are interpretable. For positive topics, we
find Topic 10 (extra credit) and Topic 11 (being nice/enjoyable class) and
words like “recommend” in a couple of them. A few words seem out of place like
“hate” in Topic 11, but that can be explained by some positive reviews having
phrases like “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, both rants and sarcasm
(suggested by Topic 2) and being harder and failing students (suggested by
Topic 4) contribute negatively to the score, but being a harder teacher seems
to contribute more negatively to the rating than ranting.

Figure 6: Topics with positive $\theta-$weights. The $\theta-$ weight is given
as the topic weight. The strength of each word is given numerically beside
each of the top 10 words. Figure 7: Topics with negative $\theta-$weights. The
$\theta-$ weight is given as the topic weight. The strength of each word is
given numerically beside each of the top 10 words.

##  6 Conclusion 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 is weakly decreasing with the regression error
weight and that, in practical applications, the error in fact strictly
decreases. The topics identified can outperform the quantitative accuracy of
topics formed through NMF alone while retaining a high degree of
interpretability.

While 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.

## References

  * [1] Lee DD, Seung HS.  Learning the parts of objects by non-negative matrix factorization.  Nature. 1999;401(6755):788-91.
  * [2] Lao H, Zhang X.  Regression and Classification of Alzheimer’s Disease Diagnosis Using NMF-TDNet Features From 3D Brain MR Image.  IEEE Journal of Biomedical and Health Informatics. 2021;26(3):1103-15.
  * [3] Lai Y, Hayashida M, Akutsu T.  Survival analysis by penalized regression and matrix factorization.  The Scientific World Journal. 2013;2013.
  * [4] Stewart GW.  On the early history of the singular value decomposition.  SIAM review. 1993;35(4):551-66.
  * [5] Joyce JM.  Kullback-leibler divergence.  In: International encyclopedia of statistical science. Springer; 2011\. p. 720-2.
  * [6] Austin W, Anderson D, Ghosh J.  Fully supervised non-negative matrix factorization for feature extraction.  In: IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium. IEEE; 2018. p. 5772-5.
  * [7] Zhu W, Yan Y.  Joint linear regression and nonnegative matrix factorization based on self-organized graph for image clustering and classification.  IEEE Access. 2018;6:38820-34.
  * [8] Haddock J, L Kassab SL, Kryshchenko A, Grotheer R, Sizikova E, Wang C, et al.  Semi-supervised Nonnegative Matrix Factorization for Document Classification.  Proc 53rd Asilomar Conf on Signals, Systems and Computers. 2021.
  * [9] Haddock J, Kassab L, Li S, Kryshchenko A, Grotheer R, Sizikova E, et al.  Semi-supervised NMF Models for Topic Modeling in Learning Tasks.  arxiv. 2020.  Submitted.
  * [10] Li P, Tseng C, Zheng Y, J A Chew LH, Jarman B, Needell D.  Guided Semi-Supervised Non-negative Matrix Factorization on Legal Documents.  Algorithms. 2022;15(5):136.
  * [11] He J. Big Data Set from RateMyProfessor.com for Professors’ Teaching Evaluation. Mendeley Data; 2020.
  * [12] Kim H, Park H.  Nonnegative matrix factorization based on alternating nonnegativity constrained least squares and active set method.  SIAM journal on matrix analysis and applications. 2008;30(2):713-30.
  * [13] Lee D, Seung HS.  Algorithms for non-negative matrix factorization.  Advances in neural information processing systems. 2000;13.
  * [14] Harris CR, Millman KJ, van der Walt SJ, Gommers R, Virtanen P, Cournapeau D, et al.  Array programming with NumPy.  Nature. 2020 Sep;585(7825):357-62.
  * [15] Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods. 2020;17:261-72.
  * [16] Berry MW, Browne M, Langville AN, Pauca VP, Plemmons RJ.  Algorithms and applications for approximate nonnegative matrix factorization.  Computational statistics & data analysis. 2007;52(1):155-73.
  * [17] Joachims T.  A Probabilistic Analysis of the Rocchio Algorithm with TFIDF for Text Categorization.  Carnegie-mellon univ pittsburgh pa dept of computer science; 1996.
  * [18] Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al.  Scikit-learn: Machine Learning in Python.  Journal of Machine Learning Research. 2011;12:2825-30.