---
title: "NNLM: A package For Fast And Versatile Nonnegative Matrix Factorization"
author: "Eric Xihui Lin"
date: "`r Sys.Date()`"
header-includes:
  - \usepackage{amsmath}
output:
  html_document:
    fig_width: 4.5
    fig_height: 4.5
    toc: true
vignette: >
  %\VignetteIndexEntry{Fast and versatile NMF}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: ref.bib
---

```{r setup, echo = FALSE, message = FALSE}
library(knitr);
opts_chunk$set(cache = TRUE, autodep = TRUE);
```


# Abstract

# Background

Non-negative matrix factorization (NMF or NNMF) has been widely used as a general method
for feature extraction on non-negative data. For example, meta-gene discovery from gene
expression profiles in @kim2007, @brunet2004.
There are currently a few algorithms for NMF decomposition, including the multiplicative algorithms
proposed by @lee1999, gradient decent and alternating non-negative least square (NNLS),
which is getting popular recently, due to its speed to convergence.

# Software and implementation

The most popular R implantation of NMF is the _NMF_ package (@nmf2010) which was first translated from a MATLAB package
and later optimized via C++ for some algorithms.
It implements various algorithms, such as the Lee's multiplicative updates based on square error
and Kullback-Leibler divergence distance, sparse alternating NNLS, etc. The sparse alternating NNLS (ALS) is supposed
to be a very fast algorithm, but in practice, when number of rows increased to tens of thousands,
this package gets very slow, probably due to its R implementation.

In NNLM, we adapt the ALS approach, but the NNLS problem is solved by
a coordinate-wise algorithm proposed by @voj2005, in which each unknown variable can be
solved sequentially and explicitly as simple quadratic optimization problems. Due to
this efficient NNLS algorithm and the usage of Rcpp, our NMF is fast. Choice of different
regularizations and many other unique features like integrating known profiles,
designable factorization, missing value handling are also available.


# Description

This package includes two main functions, `nnlm` and `nnmf`.

## Sequential coordinate-wise descent (SCD)

`nnlm` solves the following non-negative least square (NNLS)
$$ \min\limits_{\beta \ge 0} ||Y - X\beta||_F, $$

using the following sequentially coordinate-wise algorithm (@voj2005).

> 0. Let $V = X^T X.$
> 1. Initialize $\beta^{(0)} = 0$, $\mu^{(0)} = -X^T Y.$
> 2. Repeat until converge:
> for $k = 1$ to $p$,
> $$\beta_k^{(t+1)} =  \max\left(0,\,\beta_k^{(t)} - \dfrac{\mu_k^{(t)}}{v_{kk}}\right)$$
> $$\mu^{(t+1)} = \mu^{(t)} + \left( \beta_k^{(t+1)} - \beta_k^{(t)}\right) V_{\cdot k}$$
> where $V_{\cdot k}$ is the $k^{\mbox{th}}$ coloum of $V$ and $p$ is the length of $\beta$.

Note that in this problem, $Y$ and $X$ are not necessary non-negative.

## Alternating NNLS with regularization

A typical non-negative matrix factorization problem can be expressed as
$$ \min\limits_{W \ge 0, H \ge 0} \dfrac{1}{2}||A - WH||_F^2 + J_W(W) + J_H(H) $$
where
$$
\begin{aligned}
J_W(W) &= \alpha_1 J_1(W) + \alpha_2 J_2(W) + \alpha_3 J_3(W)\\
J_H(H) &= \beta_1 J_1(H^T) + \beta_2 J_2(H^T) + \beta_3 J_3(H^T)
\end{aligned}
$$
and
$$
\begin{aligned}
	J_1(X) &:= \dfrac12 ||X||_F^2 = \dfrac12 \mbox{tr}(XX^T)\\%  = \dfrac12 \sum\limits_{i, j} x_{ij}^2 \\
	J_2(X) &:= \sum\limits_{i < j} (X_{\cdot i})^TX_{\cdot j} = \dfrac12 \mbox{tr}(X(E-I)X^T) \\ %= \dfrac12 \sum\limits_{i\neq j, k} x_{ki} x_{kj}\\
	J_3(X) &:= \sum\limits_{i,j} |x_{ij}| = \mbox{tr}(XE) \\
	J_4(X) &:= \dfrac12 \sum\limits_k ||X_{k\cdot}||_1^2 = \dfrac12 \mbox{tr}(XEX^T) %= \dfrac12 \sum\limits_{i, j, k} x_{ki} x_{kj}.
\end{aligned}$$


In the above, $A \in \mathbb{R}^{n\times m}, \, W \in \mathbb{R}^{n\times K}, \, H \in \mathbb{R}^{K\times m}$,
$I$ is an identity matrix, $E$ is a matrix of proper dimension with all entries equal to 1, $w_{\cdot i}$ and $w_{i\cdot}$
are the $i^{\mbox{th}}$ column and row respectively, $h_{\cdot j}$ is the j-th column of $H$.  Obviously, $J_4 = J_1 + J_2$.

The above four types of regularizations can be used for different purposes. $J_1$ is a ridge penalty to control the magnitudes and smoothness.
$J_2(X)$ is used to minimize correlations among columns, i.e., to maximize independence or
the angle between $X_{\cdot i}$, $X_{\cdot j}$ (@zhang2008).
$J_3$ and $J_4$ (@kim2007) is a LASSO like penalty, which controls both magnitude and sparsity. However,
$J_3(X)$ tends control matrix-wise sparsity, but may also result in that some column of $X$ have all entries equal
to $0$, while $J_4(X)$ forces sparsity in a column-wise manner, which should seldom give zero column.

The alternative lease square (ALS) algorithm solves $W$ and $H$ iteratively. Due to the non-negative constraint,
the penalized NNLS is only slightly complicated. For example, When solving $H$ with $W$ fixed, i.e., minimizing for $H$ such that $H\ge0$, we have
$$
\begin{aligned}
&\frac{1}{2}||A - WH||_F^2 + J_H(H) \\
=\, & \mbox{tr}\left\{\frac12 H^T \left[W^TW + \beta_1 I + \beta_2 (E-I)\right] H - H^T\left[W^TA- \beta_3 E\right]\right\} + const.
\end{aligned}
$$

Since $E-I$ is semi-negative definite, therefore, to ensure the uniqueness and the convergence of the algorithm, one has to ensure that
$W^TW + \beta_1 I + \beta_2 (E-I)$ is positive definite. Indeed, we force a more stringent constraint that $\beta_1 \ge \beta_2$.
Similarly, $\alpha_1 \ge \alpha_2$, since the ma.  When the equality is reached, it is the $J_4$ constranint, i.e., row-wise LASSO for $W$ and
column-wise LASSO for $H$.

The following algorithm is used to solved penalized NNLS.

> 0. Let $V = W^T W + \beta_1 I + \beta_2 (E-I).$
> 1. Initialization. Set $H^{(0)} = 0$, $U^{(0)} = -W^TA + \beta_3 E.$
> 2. Repeat until converge:
> for $j = 1$ to $m$, $k = 1$ to $K$
> $$h_{kj}^{(t+1)} =  \max\left(0,\,h_{kj}^{(t)} - \dfrac{u_{kj}^{(t)}}{v_{kk}}\right)$$
> $$U_{\cdot j}^{(t+1)} = U_{\cdot j}^{(t)} + \left( h_{kj}^{(t+1)} - h_{kj}^{(t)}\right) V_{\cdot k}$$
> where $V_{\cdot k}$ is the $k^{\mbox{th}}$ coloum of $V$.


The alternating NNLS algorithm fixes $W$ and solve for $H$ using
NNLS, and then fixes $H$ and solve for $W$. This procedure is repeated until the change of
$A - W H$ is small enough.

Instead of initializing $H^{(0)} = 0$ for every iteration, we use _warm-start_, i.e., make use of the previous iteration result,
$$ H^{(0)} = H, \quad U^{(0)} = W^TWH - W^TA + \beta_3 E, $$
where $H$ is the solution from the previous iteration.


## Sequential quadratic approximation for Kullback-Leibler divergence loss

A well known problem about square error loss is that it is not robust to skewed distribution (e.g., count data) or outliers.
An alternative choice of loss function is the Kullback-Leibler divergence "distance" defined as
$$ L(A, \hat A) = \mbox{KL}(A | \hat A) = \sum_{i,j} a_{ij} \log\dfrac{a_{ij}}{\hat a_{ij}} - a_{ij} + \hat a_{ij}.$$
It can be proved that $\mbox{KL}(A | \hat A) \ge 0$ and equality is reached if and only if $A = \hat A$.  In NMF, $\hat A = W H$.

One explanation of the the $\mbox{KL}$ divergence is to assume that $a_{ij}$ is some count drawn independently as
$$ a_{ij} \,\sim\, \mbox{Poisson}(\lambda_{ij}). $$
The correpondent log likelihood is
$$ l(\Lambda) = \sum_{i, j} \left(a_{ij}\log\lambda_{ij} - \lambda_{ij}\right) + const.$$
Maximizing this log likelihood with respect to $\Lambda = \{\lambda_{ij}\}$ is equivalent minimizing $\mbox{KL}$ with respect to $\hat A$.
Assume that the observed $a_{ij}$ comes from a mixture of independent Poisson distributions with rates
$\lambda_{ij}^{(1)}, \dots, \lambda_{ij}^{(k)}$, then
$$a_{ij} \,\sim\, \mbox{Poisson}\left(\lambda_{ij} = \sum_{l=1}^K\lambda_{ij}^{(l)}\right).$$
If we force a structure $\lambda_{ij}^{(l)} = w_{il} h_{lj}$, we get NMF with $\mbox{KL}$ divergence loss.

With $\mbox{KL}$ divergence loss, a non-negative matrix factorization problem can be expressed as
$$ \min\limits_{W \ge 0, H \ge 0} \mbox{KL}(A | W H) + J_W(W) + J_H(H). $$

This problem can be solved using a similar alternating coordinate-wise algorithm, by approximating $\mbox{KL}(A | W H)$
with a quadratic function (Taylor expansion up to second order around its current value).

Assume $W$ is known and $H$ is to be solved. Let
$$ \begin{aligned}
b :=& \dfrac{\partial \mbox{KL}}{\partial h_{kj}} \left(H^0\right) =
\sum_{l}\left(w_{lk} - \frac{a_{lj}w_{lk}}{\sum_{q} w_{lq} h_{qj}^0}\right)\\
% + (\beta_1 - \beta_2) h_{kj}^0 + \beta_2 \sum_l h_{lj}^0 + \beta_3,\\
a :=& \dfrac{\partial^2\mbox{KL}}{\partial h_{kj}^2} \left(H^0\right) =
\sum_{l} a_{lj} \left(\dfrac{w_{lk}}{\sum_{q} w_{iq} h^0_{qj}}\right)^2 % + \beta_1,
\end{aligned}
$$
where $H^0$ is the current value of $H$ in the iteration procedure.

Therefore for $h_{kj}$, with all other entries fixed,
$$\mbox{KL}(A | W H) + J_H(H) \approx \frac12 a (h_{kj} - h_{kj}^0)^2 + b (h_{kj} - h_{kj}^0) +
\frac12 \beta_1 h_{kj}^2 + \beta_2 h_{kj}\sum_{l\neq k}h^0_{lj} + \beta_3 h_{kj} + const,$$
which can be minimized explicitly with constraint $h_{kj} \ge 0$, as
$$ h_{kj} = \max\left(0, \, \dfrac{a h^0_{kj} - b - \beta_2 \sum_{l\neq k} h^0_{lj} - \beta_3}{a + \beta_1} \right).$$

Obviously, this update can be parallelled in a column-wise manner.

Note that when an entry of $\hat A$ is 0, the $\mbox{KL}$ divergence is infinity. To avoid this,
we add a very small number both to $A$ and $\hat A$.

## Multiplicative updates with regularization

Two multiplicative updating algorithms are proposed in @lee1999 for square loss and KL divergence loss.
We modify these algorithms to integrate all the above regularizations as the followings.

With square loss
$$
\begin{aligned}
w_{ik} &= w_{ik} \dfrac{(AH^T)_{ik}}{\left(W\left[HH^T + \alpha_1 I + \alpha_2 (E-I)\right] + \alpha_3E \right)_{ik}},\\
h_{kj} &= h_{kj} \dfrac{(W^TA)_{kj}}{\left(\left[W^TW + \beta_1 I + \beta_2 (E-I)\right]H + \beta_3E \right)_{kj}}.
\end{aligned}
$$

With Kullback-Leibler divergence distance,
$$
\begin{aligned}
w_{ik} &= w_{ik} \dfrac{\sum_{l=1} h_{kl}a_{il}/\sum_{q}w_{iq}h_{ql}}
{\left(\sum_{l} h_{kl} + (\alpha_1-\alpha_2) w_{ik} + \alpha_2 \sum_{l} w_{il} + \alpha_3 \right)},\\
h_{kj} &= h_{kj} \dfrac{\sum_{l} w_{lk}a_{lj}/\sum_{q}w_{lq} h_{qj}}
{ \left( \sum_{l} w_{lk} + (\beta_1-\beta_2) h_{kj} + \beta_2 \sum_{l} h_{lj} +
\beta_3 \right) }.
\end{aligned}
$$

When $\alpha_i = 0, \beta_i = 0, i = 1, 2, 3$, these are the original multiplicative algorithms in @lee1999.

This multiplicative algorithm is straight forward to implemented, but it has a drawback that when $W$ or $H$ is initialized with a zero entry or positive,
it remains $0$ or positive through iterations. Therefore, true sparsity can not be achieved generally, unless a hard-thresholding is forced, as
many of entries should be small enough to be thresholded as 0. 

## Complexity and convergence speed

One can easily see that both the sequential coordinate descent (SCD) and Lee's multiplictive algorithms using MSE have complexity of
$\mathcal{O}\left(\left((m+n)k^2 N_i + 2nmK\right)N_o\right)$, while their KL conterparts
have complexity of $\mathcal{O}\left(nmk^2N_iN_o\right)$. Here $N_i$ is the number of inner iterations to solve the non-negative linear model and $N_o$ is
the number of outer iteration to alternate $W$ and $H$. $N_i x N_o$ is the total number of epochs, i.e., swaps over $W$ and $H$ matrices. Obviously
algorithm with MSE are faster then the KL based ones (by a factor of $k$) in terms of complexity, and can benefit from multiple
inner iterations $N_i$ (reducing the expensive computation of $W^TA$ and $AH^T$) as typically $k \ll m, n$, which generally should reduce $N_o$.  On the contrast, algorithm with KL
has not benefit from $N_i$ due to the re-calculation of $WH$ on each inner iteration.
Though the SCD and Lee's algorithm are close in terms of complexity, one can stil expect that SCD will converge much faster in practice.
This is because Lee's multiplicative algorithm is indeed a gradient descent with a special step size (@lee1999) which is first order method,
while the sequential coordinate-wise is a second order approach like the Newton-Raphson algorithm.


## Usage

One can install NNLM like

```{r install, echo = TRUE, eval = FALSE}
install.packages("NNLM");

# or get a dev-version from Github
library(devtools);
install_github('linxihui/NNLM')
```

### Interface

```{r interface, eval = FALSE}
nnlm(x, y, alpha = rep(0,3), method = c('scd', 'lee'), loss = c('mse', 'mkl'), init = NULL,
	mask = NULL, check.x = TRUE, max.iter = 10000L, rel.tol = 1e-12, n.threads = 1L)

nnmf(A, k = 1L, alpha = rep(0,3), beta = rep(0,3), method = c("scd", "lee"), loss = c("mse", "mkl"),
	init = NULL, mask = NULL, W.norm = -1L, check.k = TRUE, max.iter = 500L, rel.tol = 1e-4,
	n.threads = 1L, trace = 10L, verbose = 1L, show.warning = TRUE,
	inner.max.iter = ifelse('mse' == loss, 50L, 1L), inner.rel.tol = 1e-09)
```

- `x` : design matrix.
- `y` : a vector or matrix of responses.
- `A` : matrix to decompose.
- `k` : rank of NMF.
- `method`: sequential coordinate descent (`"scd"`) or Lee's multiplicative algorithm (`"lee"`).
- `metric`: loss function, mean square error or KL-divergence.
-  `init`:
    * `nnlm`: a non-negative matrix of initials.
    * `nnmf`: a list of named initial matrices for $W$ and $H$. One can also supply known matrices $W_0,\, H_0$,
       and initialize their correspondent matrices $H_1$ and $W_1$.
- `mask`:
    * `nnlm`: a logical matrix indicating if entries are fixed to 0 or initials if available.
    * `nnmf`: a list of named mask matrices for $W,\, H,\, H_1$ (if `init$W0` supplied), $W_1$ (if `init$H0` supplied),
    which should have the same shapes as $W,\, H,\, H_1$ and $W_1$ if specified. If initial matrices not specified, masked entries are fixed to 0.
- `alpha` and `beta`: a vector of length equal to or less than 3, indicating $\alpha_1, \alpha_2, \alpha_3$ and
$\beta_1, \beta_2, \beta_3$. If latter entries are missing, 0 is assumed.
- `trace`:  an integer $n$ indicating how frequent the error should be checked, i.e., $n$-iterations.
- `verbose`: either `0`/`FALSE`, `1`/`TRUE` or 2, indicating no printing, progression bar or iteration details.
- `n.threads`: number of openMP threads.
- `max.iter` and `rel.tol`:  number of outer alternating iterations and the relative tolerance.
- `inner.max.iter` and `inner.rel.tol`: number of inner iterations for solving NNLM and the relative tolerance. Default to 50 if `loss == 'mse'`
and 1 if `loss == 'mkl'`.


## Compare different algorithms

```{r alg-compare,  warning = FALSE}
library(NNLM);
set.seed(123);

k <- 15;
init <- list(W = matrix(runif(nrow(nsclc)*k), ncol = k),
	H = matrix(runif(ncol(nsclc)*k), nrow = k));
scd.mse  <- nnmf(nsclc, k, init = init, max.iter = 100, rel.tol = -1);
lee.mse  <- nnmf(nsclc, k, init = init, max.iter = 100, rel.tol = -1, method = 'lee');
scd.mkl  <- nnmf(nsclc, k, init = init, max.iter = 5000, rel.tol = -1, loss = 'mkl');
lee.mkl  <- nnmf(nsclc, k, init = init, max.iter = 5000, rel.tol = -1, loss = 'mkl',
	method = 'lee');
lee.mse1 <- nnmf(nsclc, k, init = init, max.iter = 5000, rel.tol = -1, method = 'lee',
	inner.max.iter = 1);
```

```{r alg-comp-plot, fig.show = 'hold'}
plot(NULL, xlim = c(1, 3000), ylim = c(0.15, 0.45), xlab = 'Epochs', ylab = 'MSE');
lines(cumsum(scd.mse$average.epochs), scd.mse$mse, col = 'firebrick1');
lines(cumsum(lee.mse$average.epochs), lee.mse$mse, col = 'orange');
lines(cumsum(scd.mkl$average.epochs), scd.mkl$mse, col = 'chartreuse3');
lines(cumsum(lee.mkl$average.epochs), lee.mkl$mse, col = 'deepskyblue4');
lines(cumsum(lee.mse1$average.epochs), lee.mse1$mse, col = 'orange', lty = 2);
legend('topright', bty = 'n', lwd = 1, lty = c(1,1,2,1,1),
	legend = c('SCD-MSE', 'LEE-MSE', 'LEE-MSE-1', 'SCD-MKL', 'LEE-MKL'),
	col = c('firebrick1', 'orange', 'orange', 'chartreuse3', 'deepskyblue4'));

plot(NULL, xlim = c(1, 3000), ylim = c(0.01, 0.034), xlab = 'Epochs', ylab = 'MKL');
lines(cumsum(scd.mse$average.epochs), scd.mse$mkl, col = 'firebrick1');
lines(cumsum(lee.mse$average.epochs), lee.mse$mkl, col = 'orange');
lines(cumsum(scd.mkl$average.epochs), scd.mkl$mkl, col = 'chartreuse3');
lines(cumsum(lee.mkl$average.epochs), lee.mkl$mkl, col = 'deepskyblue4');
lines(cumsum(lee.mse1$average.epochs), lee.mse1$mkl, col = 'orange', lty = 2);
legend('topright', bty = 'n', lwd = 1, lty = c(1,1,2,1,1),
	legend = c('SCD-MSE', 'LEE-MSE', 'LEE-MSE-1', 'SCD-MKL', 'LEE-MKL'),
	col = c('firebrick1', 'orange', 'orange', 'chartreuse3', 'deepskyblue4'));

summary.nnmf <- function(x) {
	if (x$n.iteration < 2) {
		rel.tol <- NA_real_;
	} else {
		err <- tail(x$target.loss, 2);
		rel.tol <- diff(err)/mean(err); 
		}
	return(c(
		'MSE' = tail(x$mse, 1), 'MKL' = tail(x$mkl, 1), 'Target' = tail(x$target.loss, 1),
		'Rel. tol.' = abs(rel.tol), 'Total epochs' = sum(x$average.epochs),
		'# Interation' = x$n.iteration, x$run.time[1:3]));
	};

library(knitr);
kable(
	sapply(
		X = list(
			'SCD-MSE' = scd.mse,
			'LEE-MSE' = lee.mse,
			'LEE-MSE-1' = lee.mse1,
			'SCD-MKL' = scd.mkl,
			'LEE-MKL' = lee.mkl
			),
		FUN = function(x) {z <- summary(x); sapply(z, sprintf, fmt = '%.4g')}
		),
	align = rep('r', 5),
	caption = 'Compare performance of different algorithms'
	);
```

One can see from the above summary, the SCD and Lee's algorithms have roughly the same
run time for each epoch, i.e., updating $W$ and $H$ entries once. However, SCD generally
converges much faster than Lee's, i.e., achieving the same accuracy in less epochs/iterations
and thus shorter time. Obviously, algorithms with mean KL loss are slower than those with
MSE for each epoch, but reducing error a bit more in each epoch. Lee's multiplicative algorith
with MSE is faster when multiple epochs used in each outer alternating iteration (LEE-MSE vs LEE-MSE-1).


# Applications


## Pattern extraction

In @lee1999, NMF is shown to be able to learn sparse basis images (column of $W$) that represent facial parts ,
like mouses, noses and eyes of different shapes.  Each face (column of $A$) is a positive combinations of a few
of these basis images.  In topic discovery, NMF can learn to group documents into topics, where
$A$ is a bag-of-word representation (count of a word in a document) of a batch of documents,
with _rows = vocabulary_, and _columns = documents_.  In NMF, the columns of $W$ (after normalization) represent topics (a topic
is represented as the distribution over vocabulary), and a column of $H$ shows the topics a document covered.

In bioinformatics, NMF can be used to discover 'meta-genes' from expression profile, which are linear combinations of genes that
may or may not related to some biology path ways.  In this case, $A$ is usually arranged as _rows = genes_ and _columns = patients_.
The columns of $W$ can then be interpreted as meta-genes, and $H$ are said to be the expression profile of meta-genes (@brunet2004, @kim2007).
@ludmil2013 used NMF to extract trinucleotide mutational signatures (basis) from next generation sequencing data (NGS) for human genes and discovered
that each cancer type is combination of these basis.


```{r nsclc, message = TRUE}
library(NNLM);
set.seed(123);

data(nsclc, package = 'NNLM')
str(nsclc)

decomp <- nnmf(nsclc[, 1:80], 3, rel.tol = 1e-5);
decomp
```

```{r nsclc1, fig.show = 'hold', message = TRUE}
heatmap(decomp$W, Colv = NA, xlab = 'Meta-gene', ylab = 'Gene', margins = c(2,2),
	labRow = '', labCol = '', scale = 'column', col = cm.colors(100));
heatmap(decomp$H, Rowv = NA, ylab = 'Meta-gene', xlab = 'Patient', margins = c(2,2),
	labRow = '', labCol = '', scale = 'row', col = cm.colors(100));
```

<!--

```{r nsclc2, fig.show = 'hold', fig.width = 5,fig.height = 5, fig.align = 'center'}
par(mar = c(5,5,4,5))
plot(decomp$mse, type = 'l', lwd = 2, xlab = 'Iteration', ylab = '',
	col = 'deepskyblue4', axes = FALSE, ylim = c(0.35, 0.45));
axis(1, ylim = seq_along(decomp$mse));
axis(2, ylim = c(0.35, 0.45), col = 'deepskyblue4', col.axis = 'deepskyblue4');
mtext("Mean Square Error", side = 2, col = 'deepskyblue4', line = 3)
par(new = TRUE);
plot(decomp$mkl, type = 'l', lwd = 2, xlab = '', ylab = '', col = 'firebrick4',
	axes = FALSE, ylim = c(0.025, 0.035));
axis(4, ylim = c(0.025, 0.035), col = 'firebrick4', col.axis = 'firebrick4');
mtext("Mean KL-divergence", side = 4, col = 'firebrick4', line = 3)
```

We see that the default alternating sequential coordinate-wise descent method coverage fairly quickly.
-->

We can also use the derived meta-genes and find their expressions in new patients, using `nnlm` which is
wrapped as by `predict` for convenience .

```{r nsclc3, fig.show = 'hold', message = TRUE}
# find the expressions of meta-genes for patient 81-100
newH <- predict(decomp, nsclc[, 81:100], which = 'H');
str(newH)
```

## Content deconvolution and designable factorization

There are several researches (@zhang2008, @alberto2006) shown the capacity of NMF in blind signal separation
similar to independent component analysis (ICA).  Here, we aim to focus on the application to
tumour content deconvolution in the area of bioinformatics, in which one of the signal is known.

Microarray is a popular technique for collecting mRNA expression. Indeed,
a mRNA profile (tumour profile) is typically a mixture of
cancer specific profile and healthy profile as the collected tumour tissues are 'contaminated' by healthy cells.
To exact the pure cancer profile for down-stream analysis, an NMF can be utilized. One can this NMF as
$$ A \approx W H + W_0 H_1,$$
where $W$ is unknown cancer profile, and $W_0$ is known healthy profile. The task here is
to deconvolute $W$, $H$ and $H_1$ from $A$ and $W_0$.


A more general deconvolution task can be expressed as
$$ A \approx W H + W_0 H_1 + W_1 H_0,$$
where $H_0$ is known coefficient matrix, e.g. a column matrix of 1. In this scenario,
$W_1$ can be interpreted as _homogeneous_ cancer profile within the specific cancer patients,
and $W$ is _heterogeneous_ cancer profile of interest for downstream analysis, such as
diagnostic or prognostic capacity, sub-type clustering.


This general deconvolution is implemented in `nnmf` via the alternating NNLS algorithm.
The known profile $W_0$ and $H_0$ can be passed via arguments `W0` and `H0`. $L_2$ and $L_1$
constrain for unknown matrices are also supported.


Once can add mask matrices for $W$ and $H$, where the masked entries are fixed to 0. Indeed, this
is a form of hard-regularization. The above known-profile feature is a special case of this mask
technique, in which masked entries are fixed to their initial values.
This feature is designed to incorporate domain knowledge, like gene sub-networks, path ways, etc.

Assume $\mathcal{S} = \{S_1, ..., S_L\}$, where $S_l, \, l = 1, ..., L$ is a set of genes in sub-network $S_l$. One can
design $W$ as a matrix of $L$ columns (or more), with $W_{i, l} = 0, \, i \not\in S_l$.  Then the NMF factorization
will learn the _weights_ of real genes in expression profile $W_{i, l}, \, i\in S_l$ from the data. This is implemented
in `nnmf` with a logical mask matrix $M_W$ = $\{\delta_{i\in S_l, l}\}$. Similar mask matrix $M_H$ with
the same shape of $H$ can be specified in the `nnmf` function.

This feature can be used for meta-analysis of different cancer types, to force some of the $k$ meta-genes
to be cancer specific (and others are shared). For example, assume $A_1, A_2$ are expressions of lung cancer and
prostate cancer microarray. By setting parts of the coefficient matrix $H$ to 0 like
$$\left(A_1\,A_2\right) = \left(W_0\,W_1\,W_2\right)
\left(\begin{array}{cc}
H_{01} & H_{02} \\
H_1 & 0\\
0 & H_2
\end{array}\right),
$$
we can expect that $W_1$ and $W_2$ are lung and prostate cancer specific profiles.


_ISOpureR_ (@isopure2015) is an R package for tumour content deconvolution based on a
graphic model introduced in @isopure2013. As argued in the above, NMF can also be used for this purpose.

```{r isopure-data, fig.align = 'center'}
set.seed(123);
if (!require(ISOpureR)) {
	install.packages('ISOpureR', repos = "http://cran.stat.sfu.ca/");
	library(ISOpureR);
	}

path.to.data <- file.path(system.file(package = 'ISOpureR'), 'extdata/Beer');
# normal profile
load(file.path(path.to.data, 'beer.normaldata.1000.transcripts.RData'));
# transcriptome of 30 patients (part of the Beer dataset)
load(file.path(path.to.data, 'beer.tumordata.1000.transcripts.30.patients.RData'));

# assume k = 3, beer.normaldata is the known healthy profile
beer.nmf <- nnmf(beer.tumordata, k = 3, init = list(W0 = beer.normaldata));

# compute proportion of tumour content
tm.frac.nmf <- with(beer.nmf, colSums(W[,1:3] %*% H[1:3,])/colSums(W %*% H));

# tumour content from ISOpureR using the full dataset
tm.frac.full <- read.delim(file.path(path.to.data, "alphapurities_full_dataset.txt"));
```

```{r isopure-plot, fig.align = 'center', fig.width = 5, fig.height = 5}
plot(tm.frac.full$alphapurities*100, tm.frac.nmf*100,
	xlim = c(20, 100), ylim = c(20, 100), col = 'firebrick1',
	xlab = "% tumour from ISOpureR full data",
	ylab = "% tumour from NMF partial data");
abline(a = 0, b = 1, lty = 2)
```

## Missing value imputation and application in recommendation system

Since matrix $A$ is assumed to have a low rank $K$, information in $A$ is redundant
for such a decomposition. Hence it is possible to allow some entries in $A$ to be absent.
Such an observation implies that on can used complete entries to impute the missing ones in $A$,
which is a classic task in
recommendation system. For example, on Netflix, each customer scores only a small proportion of
the movies and each movie is scored by a fraction of customers.
Once can expect that such a movie-customer score matrix is fairly sparse (lots of missings).
Using a NMF that allows missing values, one can predict a customer's scores on movies he/she has not watched.
A recommendation can be made simply based on the predicted scores.  In addition, the resulting $W$ and $H$
can be used to further cluster movies and customers. This feature is implemented in the `nnmf`
function and used whenever there is a missing in $A$ by using only complete entries.

The advantage of NMF imputation, compared to other model based methods, is that it takes into account
all the complete entries when imputing a single missing entry, which means it can capture complex dependency
among entries.  While a typical missing value imputation algorithm usually models missings in a feature-by-feature
(column-by-column or row-by-row) manner, and iterates over all features multiple times to capture complex dependency.


The following example shows the capacity of NMF for imputation. We use the NSCLC microarray data
as an example and set 30\% percent of the entries to `NA`. One can see that NMF imputation is faster
and better than tradition imputation algorithm like median substitution and
MICE (multivariate imputation by chained equations @mice), and close to result from  _missForest_ (@missForest)
which utilizes the black-box random forest algorithm that is highly nonlinear. NMF imputation is
the fastest which makes it suitable for big matrix.

```{r missing-data, fig.show = 'hold', message = FALSE}
set.seed(123);
nsclc2 <- nsclc;
index <- sample(length(nsclc2), length(nsclc2)*0.3);
nsclc2[index] <- NA;
```

```{r missing-nmf, results = 'hold', fig.show = 'hold', message = FALSE}
# NMF imputation
system.time(nsclc2.nmf <- nnmf(nsclc2, 2));
nsclc2.hat.nmf <- with(nsclc2.nmf, W %*% H);
```

```{r missing-mic, results = 'hold', fig.show = 'hold', message = FALSE}
# multivariate imputation by chained equations (MICE)
if(!require(mice)) {
	install.packages('mice', repos = "http://cran.stat.sfu.ca/");
	library(mice);
	}

# logarithm for positivity and the normality assumption
system.time(nsclc2.mice <- mice(log(nsclc2), m = 1, printFlag = FALSE)); # one imputation
nsclc2.hat.mice <- exp(as.matrix(complete(nsclc2.mice)));
```

```{r missing-missforest, results = 'hold', fig.show = 'hold', message = FALSE}
# imputation using random forest
if(!require(missForest)) {
	install.packages('missForest', repos = "http://cran.stat.sfu.ca/");
	library(missForest);
	}

system.time(capture.output(nsclc2.missForest <- missForest(log(nsclc2))));
nsclc2.hat.missForest <- exp(nsclc2.missForest$ximp);
```

```{r missing-plot, fig.show = 'hold'}
# simple imputation, fill-in median expression values of genes
nsclc2.hat.median <- matrix(
	apply(nsclc2, 1, median, na.rm = TRUE),
	nrow = nrow(nsclc2), ncol = ncol(nsclc2)
	);

# compare different imputations
library(knitr);
kable(
	sapply(
		X = list(
			Baseline = mean(nsclc2, na.rm = TRUE),
			Medians = nsclc2.hat.median[index],
			MICE = nsclc2.hat.mice[index],
			MissForest = nsclc2.hat.missForest[index],
			NMF = nsclc2.hat.nmf[index]
			),
		FUN = mse.mkl, # mean square error, mean KL-divergence
		obs = nsclc[index]
		),
	caption = "A comparison of different imputation methods",
	digits = 4
	);

plot(nsclc[index], nsclc2.hat.median[index], col = 'chartreuse3', cex = 0.3, pch = 20,
	ylim = c(2, 14), xlab = "Observed", ylab = "Medians")
abline(a = 0, b = 1, lwd = 2, lty = 2)

plot(nsclc[index], nsclc2.hat.mice[index], col = 'orange', cex = 0.3, pch = 20,
	ylim = c(2, 14), xlab = "Observed", ylab = "MICE imputed")
abline(a = 0, b = 1, lwd = 2, lty = 2)

plot(nsclc[index], nsclc2.hat.missForest[index], col = 'darkorchid', cex = 0.3, pch = 20,
	ylim = c(2, 14), xlab = "Observed", ylab = "MissForest imputed")
abline(a = 0, b = 1, lwd = 2, lty = 2)

plot(nsclc[index], nsclc2.hat.nmf[index], col = 'deepskyblue', cex = 0.3, pch = 20,
	ylim = c(2, 14), xlab = "Observed", ylab = "NMF imputed")
abline(a = 0, b = 1, lwd = 2, lty = 2)
```

## Determine rank $k$ via missing value imputation

Tuning hyper-parameter is a typical challenge for all unsupervised learning algorithms. The rank $k$ is
the only but very crucial parameter, which is unknown before hand. @brunet2004 suggests to run
multiple times of $k$ and uses a consensus matrix to determine $k$. This idea assumes that sample assignment
to clusters would vary little from run to run if a clustering into $k$ classes is strong. However, this
assumption is not validated and the purpose of NMF is not always for clustering. Another idea, brought from
denoising autoencoder (@dA), is to add noise to matrix $A$, factorize the noisy version and compare
the reconstructed matrix to the original $A$, and the $k$ that gives the smallest error is picked. This
could a general approach for many unsupervised learning algorithms, but in NMF, the choice of 'noise' is not
easy as the noisy version of $A$ has to be non-negative as well, which implies 'noise' may introduce bias.


Given the powerful missing value imputation in NMF, we come up with a novel idea, adapting
the well known training-validation split idea in supervised learning.
Some entries are randomly deleted from $A$ and then imputed by NMF with a set
of $k$'s.  These imputed entries are later compared to their observed values, and the $k$ that gives
the smallest error will be our choice, as only the correct $k$, if exists, has the right decomposition
that recovers the missing entries.  As contrast to training-validation split in
supervised learning, due the typically big number of entries in $A$, we have a much large sample size.
One can also easily adapt the idea of cross validation to this approach.  This idea should be applicable
to any unsupervised learning methods that supports missing value imputation.


To illustrate, we do a simulation study as follows.  As we could see, different runs give consistent
results. The mean square errors(MSEs) decrease as rank $k$ increase, but the increases rate slows down
when $k = 3$ (the true rank). Meanwhile, the MSEs for imputed values are minimized at $k = 3$ for all
runs.

```{r rank-sim, fig.width = 5, fig.height = 5, fig.align = 'center', warning = FALSE}
set.seed(678);

n <- 400;
m <- 50;
k <- 3;
W  <- matrix(runif(n*k),  n, k);
H  <- matrix(10*runif(k*m),  k, m);
noise <- matrix(rnorm(n*m), n, m);
A <- W %*% H + noise;
A[A < 0] <- 0;

plot(-1, xlim = c(1,6), ylim = c(0.5, 2.5), xlab = "Rank", ylab = "MSE")
cols <- c('deepskyblue', 'orange', 'firebrick1', 'chartreuse3');
for (col in cols) {
	ind <- sample(length(A), 0.3*length(A));
	A2 <- A;
	A2[ind] <- NA;
	err <- sapply(X = 1:6,
		FUN = function(k) {
			z <- nnmf(A2, k);
			c(mean((with(z, W %*% H)[ind] - A[ind])^2), tail(z$mse, 1));
			}
		);
	invisible(lines(err[1,], col = col, type = 'b'));
	invisible(lines(err[2,], col = col, type = 'b', lty = 2));
	}
```

For imputation, as in the previous example for NSCLC, we chose $k = 2$. But why?
To determine the rank, we set a fraction of the complete entries to missing, and
impute them. MSEs are computed to determine the optimal rank. As from the resulting
graph below, one can tell that $k = 2$ is the optimal.

```{r rank-nsclc, fig.width = 5, fig.height = 5, fig.align = 'center', warning = FALSE}
set.seed(567);

plot(0, xlim = c(1,10), ylim = c(0.4, 1.4), xlab = "Rank", ylab = "MSE")
cols <- c('deepskyblue', 'orange', 'firebrick1', 'chartreuse3');
for (col in cols) {
	index2 <- sample(which(!is.na(nsclc2)), 2000);
	nsclc3 <- nsclc2;
	nsclc3[index2] <- NA;
	err <- sapply(X = 1:10,
		FUN = function(k, A) {
			z <- nnmf(A, k, verbose = FALSE);
			mean((with(z, W%*%H)[index2] - nsclc2[index2])^2)
			},
		A = nsclc3
		);
	invisible(lines(err, col = col, type='b', lwd = 2, cex = 1));
	}
```


## Noise reduction

This is obvious as the reconstruction from $W$ and $H$ lies on a smaller dimension, and should therefore
give a smoother reconstruction. This noise reduction is particularly useful when the noise is not Gaussian
which cannot be done using many other methods where Gaussian noise is usually assumed.


# NNLM vs NMF

The R package _NMF_ by @nmf2010 is an excellent and complete suite for basic NMF decomposition. It has a careful design,
various choices of different algorithms (some of which are
implemented in C for speed), different built-in initializations, a large amount of utility functions
to process resulting NMF object, such as visualization,
meta-gene thresholding. It supports multiple runs via the _foreach_ framework, which parallels
multiple runs on a multiple-cores machine and HPC cluster.

One the other hand, _NNLM_ does not provide any functionality other than NMF decomposition. It focuses
on speed and applications. The built-in coordinate NNLS solver is highly efficient to
solve NMF via the alternating scheme. The algorithm is parallelled within a single run via openMP. Thus it is fast
for a single NMF, which makes NMF suitable for big matrix. We believe our users can parallel multiple runs
via _foreach_, _doMC_ and  _doMPI_ without much effort, thanks to the simplicity of _foreach_.
The NMF algorithms in NNLM are given more choices of regularizations, both traditional penalty methods and hard
regularization like partially known $W$ and $H$, which can be powerful to integrate previous or domain
knowledge into NMF. All built-in NMF algorithms are capable for matrices with missing values,
which is desirable in some applications like recommendation system and imputation. Although only random
initialization is built-in, user specified initial matrices is accepted, which can come from an
SVD, ICA, etc.


\clearpage

# Reference

<!---
# Examples
## Non-small Cell Lung Cancer expression data


## Simulated deconvolution

```{r deconvol, message = TRUE}
n <- 200; m <- 50;
k <- 5; k1 <- 2; k2 <- 1;

set.seed(123);
W  <- matrix(runif(n*k),  n, k); # unknown heterogeneous cancer profile
H  <- matrix(runif(k*m),  k, m);
W0 <- matrix(runif(n*k1), n, k1); # known healthy profile
H1 <- matrix(runif(k1*m), k1, m);
W1 <- matrix(runif(n*k2), n, k2); # unknown common cancer profile
H0 <- matrix(1, k2, m);
noise <- 0.1*matrix(runif(n*m), n, m);

# A is the observed profile to be de-convoluted
A <- W %*% H + W0 %*% H1 + W1 %*% H0 + noise;

deconvol <- nnmf(A, k = 5, init = list(W0 = W0, H0 = H0), method = 'lee');
```

Check if $W$ and $H$, our main interest, are recovered.

```{r deconvol2, message = TRUE}
round(cor(W, deconvol$W[, seq_len(k)]), 2);
# round(cor(t(H), t(deconvol$H[seq_len(k), ])), 2);
```

We see that $W$, $H$ are just permuted. However, as we known that
the minimization problem for NMF usually has not unique solutions
for $W$ and $H$. Therefore, $W$ and $H$ cannot be guaranteed to
be recovered exactly(different only with a permutation and a scaling).

```{r deconvol3, eval = c(T, F, T), echo = c(TRUE, FALSE, TRUE),  message = TRUE}
permutation <- c(3, 1, 3, 5, 4);
round(cor(W, deconvol$W[, permutation]), 2);
#round(cor(t(H), t(deconvol$H[permutation, ])), 2);
```


As from the following result, $H_1$, coefficients of health profile and
$W_1$, common cancer profile, are recovered fairly well.

```{r deconvol4, message = TRUE}
round(cor(t(H1)), 2);
round(cor(t(H1), t(deconvol$H[6:7, ])), 2);
round(cor(W1, deconvol$W[, 8]), 2);
```

## Tumour deconvolution


## Missing value imputation
-->
