---
abstract: |
    Non-negative matrix factorization (NMF) has become a well-established class of methods for the analysis of non-negative data. In particular, a lot of effort has been devoted to probabilistic NMF, namely estimation or inference tasks in probabilistic models describing the data, based for example on Poisson or exponential likelihoods. When dealing with time series data, several works have proposed to model the evolution of the activation coefficients as a non-negative Markov chain, most of the time in relation with the Gamma distribution, giving rise to so-called temporal NMF models. In this paper, we review three Gamma Markov chains of the NMF literature, and show that they all share the same drawback: the absence of a well-defined stationary distribution. We then introduce a fourth process, an overlooked model of the time series literature named BGAR(1), which overcomes this limitation. These four temporal NMF models are then compared in a MAP framework on a prediction task, in the context of the Poisson likelihood.

    **Keywords:** Non-negative matrix factorization, Time series data, Gamma Markov chains, MAP estimation
author:
- |
    Louis Filstroff$^{1}$ Olivier Gouvert$^1$ Cédric Févotte$^1$ Olivier Cappé$^2$\
    $^1$ IRIT, Université de Toulouse, CNRS, France\
    $^2$ DI ENS, CNRS, INRIA, Université PSL
bibliography:
- 'articlebib.bib'
title: 'A Comparative Study of Temporal Non-Negative Matrix Factorization with Gamma Markov Chains'
---

Introduction {#sec-intro}
============

Non-negative matrix factorization
---------------------------------

Non-negative matrix factorization (NMF) [@paatero1994positive; @lee1999learning] has become a widely used class of methods for analyzing non-negative data. Let us consider $N$ samples in $\mathbb{R}^{F}_{+}$. We can store these samples column-wise in a matrix, which we denote by $\mathbf{V}$ (therefore of size $F \times N$). Broadly speaking, NMF aims at finding an approximation of $\mathbf{V}$ as the product of two non-negative matrices: $$\mathbf{V} \simeq \mathbf{WH},
\label{eq-nmf-approx}$$ where $\mathbf{W}$ is of size $F \times K$, and $\mathbf{H}$ is of size $K \times N$. $\mathbf{W}$ and $\mathbf{H}$ are referred to as the dictionary and the activation matrix, respectively. $K$ is usually chosen such that $K \ll \min(F,N)$, hence producing a low-rank approximation of $\mathbf{V}$. This factorization is often retrieved as the solution of an optimization problem, which we can write as: $$\min_{\mathbf{W} \geq 0,~\mathbf{H} \geq 0} D(\mathbf{V}|\mathbf{WH}),
\label{eq-nmf-min}$$ where $D$ is a measure of fit between $\mathbf{V}$ and its approximation $\mathbf{WH}$, and the notation $\mathbf{A} \geq 0$ denotes the non-negativity of the entries of the matrix $\mathbf{A}$. One of the key aspects to the success of NMF is that the non-negativity of the factors $\mathbf{W}$ and $\mathbf{H}$ yields an interpretable, part-based representation of each sample: $\mathbf{v}_n \simeq \mathbf{Wh}_n$ [@lee1999learning].

Various measures of fit have been considered in the literature, for instance the family of $\beta$-divergences [@fevotte2011algorithms], which includes some of the most popular cost functions in NMF, such as the squared Euclidian distance, the generalized Kullback-Leibler divergence, or the Itakura-Saito divergence. As it turns out, for many of these cost functions, the optimization problem described in Eq.  can be shown to be equivalent to the joint maximum likelihood estimation of the factors $\mathbf{W}$ and $\mathbf{H}$ in a statistical model, that is: $$\max_{\mathbf{W},\mathbf{H}} p(\mathbf{V}|\mathbf{W},\mathbf{H}).
\label{eq-nmf-max}$$ This leads the way to so-called *probabilistic* NMF, i.e., estimation or inference tasks in probabilistic models whose observation distribution may be written as: $$\mathbf{v}_{n} \sim p(~.~;\mathbf{Wh}_{n}, \boldsymbol{\Theta}), \quad \mathbf{W} \geq 0, \quad \mathbf{H} \geq 0,
\label{eq-prob-nmf}$$ that is to say that the distribution of $\mathbf{v}_n$ is parametrized by the dot product of the factors $\mathbf{W}$ and $\mathbf{h}_n$. Other potential parameters of the distribution are generically denoted by $\boldsymbol{\Theta}$. Most of the time these distributions are such that $\mathbb{E}(\mathbf{v}_n) = \mathbf{Wh}_n$.

This large family encompasses many well-known models of the literature, for example models based on the Gaussian likelihood [@schmidt2009bayesian] or the exponential likelihood [@fevotte2009nonnegative; @hoffman2010bayesian]. It also includes factorization models for count data, which are most of the time based on the Poisson distribution[^1] [@canny2004gap; @cemgil2009bayesian; @zhou2012beta; @gopalan2015scalable], but can also make use of distributions with a larger tail, e.g., the negative binomial distribution [@zhou2018nonparametric]. Finally, more complex models using the compound Poisson distribution have been considered [@simsekli2013learning; @basbug2016hierarchical; @gouvert2019recommendation], allowing to extend the use of the Poisson distribution to various supports $(\mathbb{N}, \mathbb{R}_+, \mathbb{R}, ...)$.

In the vast majority of the aforementioned works, prior distributions are assumed on the factors $\mathbf{W}$ and $\mathbf{H}$. This is sometimes referred to as *Bayesian* NMF. In this case, the columns of $\mathbf{H}$ are most of the time assumed to be independent: $$p(\mathbf{H}) = \prod_{n=1}^{N} p(\mathbf{h}_n).
\label{eq-h-idp}$$ The factors being non-negative, a standard choice is the Gamma distribution, which can be sparsity-inducing if the shape parameter is chosen to be lower than one. The inverse Gamma distribution has also been considered.

Temporal structure of the activation coefficients
-------------------------------------------------

In this work, we are interested in the analysis of specific matrices $\mathbf{V}$ whose columns cannot be treated as exchangeable, because the samples $\mathbf{v}_n$ are correlated. Such a scenario arises in particular when the columns of $\mathbf{V}$ describe the evolution of a process over time.

From a modeling perspective, this means that correlation should be introduced in the statistical model between successive columns of $\mathbf{V}$. This can be achieved by lifting the prior independence assumption of Eq. , thus introducing correlation between successive columns of $\mathbf{H}$. In this paper, we consider a Markov structure on the columns of $\mathbf{H}$: $$p(\mathbf{H}) = p(\mathbf{h}_1) \prod_{n \geq 2} p(\mathbf{h}_n | \mathbf{h}_{n-1}).
\label{eq-h-markov}$$ We will refer to such a model as a $\textit{dynamical}$ NMF model. Note that very recent works go beyond the Markovian assumption, i.e., assume dependency with multiple past time steps, and are labeled as “deep” [@gong2017deep; @guo2018deep].

Several works [@fevotte2013non; @schein2016poisson] assume that the transition distribution $p(\mathbf{h}_n | \mathbf{h}_{n-1})$ makes use of a transition matrix $\boldsymbol{\Pi}$ of size $K \times K$ to capture relationships between the different components. In this case, the distribution of $h_{kn}$ depends on a linear combination of all the components at the previous time step: $$p(\mathbf{h}_n | \mathbf{h}_{n-1}) = \prod_{k} p(h_{kn}|\sum_l \pi_{kl} h_{l(n-1)}).
\label{eq-trans-h-matrix}$$

In this work, we will restrict ourselves to $\boldsymbol{\Pi} = \mathbf{I}_K$. Equivalently, this amounts to assuming that the $K$ rows of $\mathbf{H}$ are a priori independent, and we have $$p(\mathbf{H}) = \prod_{k} p(h_{k1}) \prod_{n \geq 2}p(h_{kn}|h_{k(n-1)}).
\label{eq-h-markov-idp}$$ We will refer to such a model as a *temporal* NMF model.

A first way of dealing with the temporal evolution of a non-negative variable is to map it to $\mathbb{R}_+$. It is then commonly assumed that this variable evolves in Gaussian noise. This is for example exploited in the seminal work of @blei2006dynamic on the extension of latent Dirichlet allocation to allow for topic evolution[^2]. A similar assumption is made in @charlin2015dynamic, which introduces dynamics in the context of a Poisson likelihood (factorizing the user-item-time tensor). Gaussian assumptions allow to use well-known computational techniques, such as Kalman filtering, but result in loss of interpretability.

We will focus in this paper on naturally non-negative Markov chains. Various non-negative Markov chains have been proposed in the NMF literature [@cemgil2007conjugate; @fevotte2009nonnegative; @acharya15nonparametric]. They are all built in relation with the Gamma (or inverse Gamma) distribution. As a matter of fact, these models exhibit the same drawback: the chains all have a degenerate stationary distribution. This can lead to undesirable behaviors, such as the instability or the degeneracy of realizations of the chains. We emphasize that this is problematic from the probabilistic perspective only, since these prior distributions may still represent an appropriate regularization in a MAP setting.

Contributions and organization of the paper
-------------------------------------------

The contributions of this paper are 4-fold:

-   We review the existing non-negative Markov chains of the NMF literature and discuss some of their limitations. In particular we show that these chains all have a degenerate stationary distribution;

-   We present an overlooked non-negative Markov chain from the time series literature, the first-order autoregressive Beta-Gamma process, denoted as BGAR(1) @lewis1989gamma, whose stationary distribution is Gamma. To the best of our knowledge, this particular chain has never been considered to model temporal dependencies in matrix factorization problems;

-   We derive majorization-minimization-based algorithms for maximum a posteriori (MAP) estimation in the NMF models with all presented prior structures on $\mathbf{H}$, including BGAR(1);

-   We compare the performance of all these models on a prediction task on three real-world datasets.

The paper is organized as follows. Section \[sec-study\] introduces and compares non-negative Markov chains from the literature. Section \[sec-map\] presents MAP estimation in temporal NMF models. Experimental work is conducted in Section \[sec-exp\], before concluding in Section \[sec-ccl\].

Comparative study of Gamma Markov chains {#sec-study}
========================================

This section reviews existing models of Gamma Markov chains, i.e., Markov chains which evolve in $\mathbb{R}_{+}$ in relation with the Gamma distribution. We have identified three different models in the NMF literature:

1.  Chaining on the rate parameter of a Gamma distribution (Section \[sec-rate\]);

2.  Chaining on the rate parameter of a Gamma distribution with an auxiliary variable (Section \[sec-cd-rate\]);

3.  Chaining on the shape parameter of a Gamma distribution (Section \[sec-shape\]).

As shall be discussed in these subsections, these three models are all built around the assumption $\mathbb{E}(h_{kn}|h_{k(n-1)}) \propto h_{k(n-1)}$ (which roughly means that the chain should not drift too far away from its previous value), but lack a well-defined stationary distribution, which leads to the degeneracy of the realizations of the chains. A fourth model from the time series literature, called BGAR(1), is presented in Section \[sec-bgar\]. It is built to have a well-defined stationary distribution (it is marginally Gamma distributed), but does not share the property $\mathbb{E}(h_{kn}|h_{k(n-1)}) \propto h_{k(n-1)}$. The realizations of the chain are not degenerate and exhibit some interesting properties. To the best of our knowledge, this kind of process has never been used in a probabilistic NMF problem to model temporal evolution.

Throughout the section, $(h_n)_{n \geq 1}$ denotes the (scalar) Markov chain of interest, where the index $k$ as in Eq.  has been dropped for enhanced readability. It is further assumed that $h_1$ is set to a fixed, deterministic value.

Chaining on the rate parameter {#sec-rate}
------------------------------

### Model

Let us consider a general Gamma Markov chain model with a chaining on the rate parameter: $$h_{n}|h_{n-1} \sim \text{Gamma} \left( \alpha, \frac{\beta}{h_{n-1}} \right).
\label{eq-gmc-rate-def}$$ As it turns out, Eq.  can be rewritten as a multiplicative noise model: $$h_{n} = h_{n-1} \times \phi_n,
\label{eq-gmc-rate-noise}$$ where $\phi_n$ are i.i.d. Gamma random variables with parameters $(\alpha, \beta)$. We have $$\mathbb{E}(h_{n}|h_{n-1}) = \frac{\alpha}{\beta}h_{n-1}, \quad
\text{var}(h_{n}|h_{n-1}) = \frac{\alpha}{\beta^2}h^2_{n-1}.$$

This model was introduced in @fevotte2009nonnegative to add smoothness to the activation coefficients in the context of audio signal processing. The parameters were set to $\alpha > 1$ and $\beta = \alpha - 1$, such that the mode would be located at $h_{n} = h_{n-1}$. A similar inverse Gamma Markov chain was also considered in @fevotte2009nonnegative and in @fevotte2011majorization.

### Analysis

From Eq.  we can write: $$h_{n} = h_1 \prod_{i=2}^{n} \phi_i.$$

The independence of the $\phi_i$ yields: $$\begin{aligned}
\mathbb{E}(h_{n}) & = h_1 \left( \frac{\alpha}{\beta} \right)^{n-1}, \\
\text{var}(h_{n}) & = h_1^2 \left[ \left( \frac{\alpha^2}{\beta^2} + \frac{\alpha}{\beta^2} \right)^{n-1} - \left( \frac{\alpha^2}{\beta^2} \right)^{n-1} \right].
\label{eq-gmc-rate-moments}\end{aligned}$$

We enumerate all the possible regimes ($n \rightarrow +\infty)$, which all give rise to degenerate stationary distributions for different reasons:

-   $\beta > \sqrt{\alpha(\alpha+1)}$: both mean and variance go to zero;

-   $\beta = \sqrt{\alpha(\alpha+1)}$: variance converges to 1, however the mean goes to zero;

-   $\beta \in \left]\alpha;\sqrt{\alpha(\alpha+1)}\right[$: variance goes to infinity, mean goes to zero;

-   $\beta = \alpha$: mean is equal to 1, but the variance goes to infinity;

-   $\beta < \alpha$: both mean and variance go to infinity.

Each subplot of Figure \[fig-ch4-gmc-rate\] displays ten independent realizations of the chain, for a different set of parameters $(\alpha,\beta)$. As we can see, the realizations of the chain either collapse to 0, or diverge.

Hierarchical chaining with an auxiliary variable {#sec-cd-rate}
------------------------------------------------

### Model

Let us consider the following Gamma Markov chain model introduced in @cemgil2007conjugate: $$\begin{aligned}
z_{n}|h_{n-1} & \sim \text{Gamma}(\alpha_z, \beta_z h_{n-1}), \label{eq-gmc-cd-def1} \\
h_{n}|z_{n} & \sim \text{Gamma}(\alpha_h, \beta_h z_{n}). \label{eq-gmc-cd-def2}\end{aligned}$$ As it turns out, this model can also be rewritten as a multiplicative noise model: $$h_{n} = h_{n-1} \times \tilde{\phi}_n,
\label{eq-gmc-cd-noise}$$ where $\tilde{\phi}_n$ are i.i.d. random variables defined as the ratio of two independent Gamma random variables with parameters $(\alpha_h, \beta_h)$ and $(\alpha_z, \beta_z)$. The distribution of $\tilde{\phi}_n$ is actually known in closed form, namely $$\tilde{\phi}_n \sim \text{BetaPrime} \left( \alpha_h, \alpha_z, 1, \tilde{\beta} \right),
\label{eq-gmc-cd-betaprime}$$ with $\tilde{\beta} = \frac{\beta_z}{\beta_h}$ (see Appendix \[app-a\] for a definition). We have $$\begin{aligned}
\mathbb{E}(h_{n}|h_{n-1}) & = \tilde{\beta} \frac{\alpha_h}{\alpha_z-1} h_{n-1} & \text{for~} \alpha_z > 1, \\
\text{var}(h_{n}|h_{n-1}) & = \tilde{\beta}^2 \frac{\alpha_h ( \alpha_h + \alpha_z - 1)}{(\alpha_z-1)^2 (\alpha_z-2)} h^2_{n-1} & \text{for~} \alpha_z > 2.\end{aligned}$$

This model is less straightforward in its construction than the previous one, as it makes use of an auxiliary variable $z_n$ (note that a similar inverse Gamma construction was proposed as well in @cemgil2007conjugate). There are two motivations behind the introduction of this auxiliary variable:

1.  Firstly, it ensures what is referred to as “positive correlation” in @cemgil2007conjugate, i.e., $\mathbb{E}(h_n|h_{n-1}) \propto h_{n-1}$ (something the model described by Eq.  does as well).

2.  Secondly, it ensures the so-called conjugacy of the model, i.e., the conditional distributions $p(z_{n}|h_{n-1},h_{n})$ and $p(h_{n}|z_{n},z_{n+1})$ remain Gamma distributions. Indeed, these are the distributions of interest when considering Gibbs sampling or variational inference. This property is not satistfied by the model described by Eq.  (i.e., $p(h_n|h_{n-1},h_{n+1})$ is neither Gamma, nor a known distribution).

This particular chain has been used in the context of audio signal processing in @virtanen2008bayesian (under the assumption of a Poisson likelihood, which does not fit the nature of the data), and also to model the evolution of user and item preferences in the context of recommender systems [@jerfel17dynamic; @do2018gamma].

### Analysis

From Eq. , we can write: $$h_{n} = h_1 \prod_{i=2}^{n} \tilde{\phi}_{i}.$$

We have by independence of the $\tilde{\phi}_i$: $$\begin{aligned}
\mathbb{E}(h_{n}) & = h_{1} \left( \tilde{\beta} \frac{\alpha_h}{\alpha_z - 1} \right)^{n-1} \qquad \qquad \qquad \text{for~}\alpha_z > 1,  \\
\text{var}(h_{n}) & = h_1^2 \tilde{\beta}^{2(n-1)} \left[
\left( \frac{\alpha_h^2}{(\alpha_z - 1)^2} + \frac{\alpha_h ( \alpha_h + \alpha_z - 1)}{(\alpha_z-1)^2 (\alpha_z-2)} \right)^{n-1} \right. \notag \\
& \left. \qquad \qquad \qquad - \left( \frac{\alpha_h^2}{(\alpha_z - 1)^2} \right)^{n-1}
\right]~\text{for~} \alpha_z > 2.\end{aligned}$$

As in the previous model, we can show that either the expectation or the variance diverges or collapses as $n \rightarrow \infty$ for every possible choice of parameters, which means that they all give rise to a degenerate stationary distribution of the chain. Each subplot of Figure \[fig-ch4-gmc-cd\] displays ten independent realizations of the chain, for a different set of parameters $(\alpha_z,\beta_z,\alpha_h,\beta_h)$. As we can see, the realizations of the chain either collapse to 0, or diverge.

Chaining on the shape parameter {#sec-shape}
-------------------------------

### Model

Let us consider a general Gamma Markov chain model with a chaining on the shape parameter: $$h_{n}|h_{n-1} \sim \text{Gamma}(\alpha h_{n-1}, \beta).
\label{eq-gmc-shape-def}$$ We have $$\mathbb{E}(h_{n}|h_{n-1}) = \frac{\alpha}{\beta}h_{n-1}, \quad
\text{var}(h_{n}|h_{n-1}) = \frac{\alpha}{\beta^2}h_{n-1}.$$

In contrast with the two models presented previously, this model cannot be rewritten as a multiplicative noise model. This model is therefore more intricate to interpret. It was introduced in @acharya15nonparametric in the context of Poisson factorization. It is mainly motivated by a computational trick that can be used when working with a Poisson likelihood, hence making a Gibbs sampling feasible in the model. The authors set the value of $\alpha$ to 1 (though the same trick can be applied for any value of $\alpha$).

### Analysis

Using the law of total expectation and total variance, it can be shown that $$\mathbb{E}(h_{n}) = h_1 \left( \frac{\alpha}{\beta} \right)^{n-1}, \quad
\text{var}(h_{n}) = h_{1} \frac{1}{\beta} \left( \frac{\alpha}{\beta} \right)^{n-1} \sum_{i=1}^{n-1} \left( \frac{\alpha}{\beta} \right)^i.$$

The discussion is hence driven by the value of $r = \frac{\alpha}{\beta}$.

-   If $r < 1$, mean and variance go to zero;

-   If $r = 1$, mean is fixed but variance goes to infinity (linearly);

-   If $r > 1$, mean and variance go to infinity.

This chain only exhibits degenerate stationary distributions. Each subplot of Figure \[fig-ch4-gmc-shape\] displays ten independent realizations of the chain, for a different set of parameters $(\alpha,\beta)$. As we can see, the realizations of the chain either collapse to 0, or diverge.

BGAR(1) {#sec-bgar}
-------

We now discuss the first order autoregressive Beta-Gamma process of @lewis1989gamma, a stochastic process which is marginally Gamma distributed. The authors referred to the process as “BGAR(1)”. However, to the best of our knowledge, no extension to higher-order autoregressive processes exists in the time series literature. As such, from now on, we will simply refer to it as “BGAR”.

### Model

Consider $\alpha > 0$, $\beta > 0$, $\rho \in [0,1[$. The BGAR process is defined as: $$\begin{aligned}
h_1 & \sim \text{Gamma}(\alpha,\beta), \label{eq-gmc-bgar-def1} \\
h_n & = b_n h_{n-1} + \epsilon_n \qquad \text{for~} n \geq 2 \label{eq-gmc-bgar-def2},\end{aligned}$$ where $b_n$ and $\epsilon_n$ are i.i.d. random variables distributed as: $$\begin{aligned}
b_n & \sim \text{Beta}(\alpha \rho, \alpha(1-\rho)), \label{eq-gmc-bgar-def3} \\
\epsilon_n & \sim \text{Gamma}(\alpha(1-\rho), \beta) \label{eq-gmc-bgar-def4}.\end{aligned}$$ $(h_n)_{n \geq 0}$ is called the BGAR process. It is parametrized by $\alpha$, $\beta$ and $\rho$. We emphasize that the distribution $p(h_n|h_{n-1})$ is not known in closed form. Only $p(h_n|h_{n-1}, b_n)$ is known; it is a shifted Gamma distribution. The generative model may therefore be rewritten as $$\begin{aligned}
h_1 & \sim \text{Gamma}(\alpha,\beta), \label{eq-b1} \\
b_n & \sim \text{Beta}(\alpha \rho, \alpha(1-\rho)) \quad \text{for~}n \geq 2, \label{eq-b2} \\
h_n | b_n, h_{n-1} & \sim \text{Gamma}(\alpha(1-\rho), \beta, \text{loc} = b_n h_{n-1}) \label{eq-b3} \\
& \qquad \qquad \qquad \qquad \qquad~\text{for~}n \geq 2, \notag\end{aligned}$$ where the distribution in Eq.  is a shifted Gamma distribution with a location parameter “loc”.

We have $$\begin{aligned}
\mathbb{E}(h_n|h_{n-1}) & = \rho h_{n-1} + \frac{\alpha(1-\rho)}{\beta}, \label{eq-gmc-bgar-cond-moment1} \\
\text{var}(h_n|h_{n-1}) & = \frac{\rho(1-\rho)}{\alpha + 1}h^2_{n-1} + \frac{\alpha(1-\rho)}{\beta^2}. \label{eq-gmc-bgar-cond-moment2}\end{aligned}$$ As we can see, BGAR(1) already differs from the three previously presented models because the conditional expectation $\mathbb{E}(h_n|h_{n-1})$ is not proportional to $h_{n-1}$ (it is an affine transformation).

### Analysis

To study the marginal distribution of the process, we recall the following lemma.

If $X \sim \text{Beta}(a,b)$ and $Y \sim \text{Gamma}(a+b,c)$ are independent random variables, then $Z = XY$ is $\text{Gamma}(a,c)$ distributed. \[lemma-ch4\]

Let $(h_n)_{n \geq 1}$ be a BGAR process. Then $h_n$ is marginally $\text{Gamma}(\alpha,\beta)$ distributed.

Follows by induction. Consider $n$ such that $h_n$ is $\text{Gamma}(\alpha,\beta)$ distributed. Then, $\epsilon_{n+1}h_n$ is $\text{Gamma}(\alpha \rho,\beta)$ distributed (Lemma \[lemma-ch4\]). Finally, $h_{n+1} = \epsilon_{n+1}h_n + b_{n+1}$ is $\text{Gamma}(\alpha,\beta)$ distributed (sum of independent Gamma random variables), which concludes the proof.

Therefore the parameters $\alpha$ and $\beta$ control the marginal distribution. The parameter $\rho$ controls the correlation between successive values, as discussed in the following proposition.

Let $(h_n)_{n \geq 1}$ be a BGAR process. Let $n$ and $r$ be two integers such that $r > 1$. We have $\text{corr}(h_n, h_{n+r}) = \rho^{r}$. \[prp2\]

See Appendix \[app-b\] for $r=1$.

Proposition \[prp2\] implies that the BGAR(1) process admits a (second order) AR(1) representation. Two limit cases of BGAR can be exhibited:

-   When $\rho = 0$, the $h_n$ are i.i.d. random variables;

-   When $\rho \rightarrow 1$, the process is not random anymore, and $h_n = h_1$ for all $n$ (note that $\rho = 1$ is not an admissible value).

Finally, from Eq. , we have $$\bigg( \mathbb{E}(h_n|h_{n-1}) > h_{n-1} \bigg) \Leftrightarrow \bigg( h_{n-1} < \frac{\alpha}{\beta} \bigg).$$ If $h_{n-1}$ is below the mean of the marginal distribution ($\frac{\alpha}{\beta}$), then $h_n$ will be in expectation above $h_{n-1}$, and vice-versa.

Note that BGAR is not the only Markovian process with a marginal Gamma distribution considered in the literature. We mention the GAR(1) process (first-order autoregressive Gamma process) of @gaver1980first, which is also marginally Gamma distributed. However, this particular process is piecewise deterministic, and its parameters are “coupled”: the parameters of the marginal distribution also have an influence on other properties of the model. As such, it is less suited to our problem, and will not be considered here.

Figure \[fig-ch4-bgar\] displays three realizations of the BGAR process, with parameters fixed to $\alpha = 2$ and $\beta = 1$, and a different parameter $\rho$ in each subplot. The mean of the marginal distribution is displayed in red. When $\rho = 0.5$, the correlation is weak, and no particular structure is observed. However, as $\rho$ goes to 1, the correlation becomes stronger, and we typically observe piecewise constant trajectories.

MAP estimation in temporal NMF models {#sec-map}
=====================================

We now turn to the problem of maximum a posteriori (MAP) estimation in temporal NMF models. More precisely, we assume a Poisson likelihood, that is $$v_{fn} \sim \text{Poisson}([\mathbf{WH}]_{fn}),$$ and we also assume that $\mathbf{W}$ is a deterministic variable. We consider four different models corresponding to the four temporal structures on $\mathbf{H}$ presented in Section \[sec-study\]. As such, $\mathbf{V}$ and $\mathbf{H}$ define a hidden Markov model, as displayed on Figure \[figure-hmm\].

(0,0) circle (0.5); (1.5,0) circle (0.5); (3,0) circle (0.5); (0,-1.5) circle (0.5); (1.5,-1.5) circle (0.5); (3,-1.5) circle (0.5); (0.5,0) – (1,0); (2,0) – (2.5,0); (0, -1) – (0, -0.5); (1.5, -1) – (1.5, -0.5); (3, -1) – (3, -0.5); (3.5, 0) – (3.95, 0); (-0.5, 0) – (-0.95, 0); (0,0) node[[$\mathbf{h}_{n-1}$]{}]{}; (1.5,0) node[[$\mathbf{h}_{n}$]{}]{}; (3,0) node[[$\mathbf{h}_{n+1}$]{}]{}; (0,-1.5) node[[$\mathbf{v}_{n-1}$]{}]{}; (1.5,-1.5) node[[$\mathbf{v}_{n}$]{}]{}; (3,-1.5) node[[$\mathbf{v}_{n+1}$]{}]{};

(1.5,-2.5) node [$\bullet$]{}; (1.5,-2.8) node [[$\mathbf{W}$]{}]{}; (1.5,-2.5) – (1.5,-2); (1.5,-2.5) – (0,-2); (1.5,-2.5) – (3,-2); (1.5, -2.5) – (-0.75, -2.125); (1.5, -2.5) – (3.75, -2.125);

Generally speaking, joint MAP estimation in such models amounts to minimizing the following criterion $$\begin{aligned}
C(\mathbf{W},\mathbf{H}, \boldsymbol{\beta}) & = -\log p(\mathbf{V},\mathbf{H};\mathbf{W},\boldsymbol{\beta}) \label{eq-ccc} \\
& = -\log p(\mathbf{V}|\mathbf{H};\mathbf{W}) - \sum_k \left( \log p(h_{k1}) + \sum_{n \geq 2} \log p(h_{kn}|h_{k(n-1)};\beta_k) \right),\end{aligned}$$ that is to say that the factors $\mathbf{W}$ and $\mathbf{H}$, as well as the scale hyperparameters $\boldsymbol{\beta} = [\beta_1, \dotsc, \beta_K]^{\text{T}}$, are going to be estimated (shape hyperparameters $\alpha$ or $\rho$ will be treated as fixed).

The optimization of the function $C$ is carried out with a block coordinate descent scheme over the variables $\mathbf{W}$, $\mathbf{H}$, and $\boldsymbol{\beta}$.

For the first two steps, we resort to a majorization-minimization (MM) scheme, which consists in iteratively majorizing the function $C$ (by a so-called auxiliary function, tight for some $\tilde{\mathbf{W}}$ or $\tilde{\mathbf{H}}$), and minimizing this auxiliary function instead. We refer the reader to @hunter2004tutorial for a detailed tutorial. Under this scheme, the function $C$ is non-increasing. As it turns out, only the Poisson likelihood term $-\log p(\mathbf{V}|\mathbf{H};\mathbf{W})$ needs to be majorized. This is a well-studied issue in the NMF literature. As stated in @lee2000algorithms [@fevotte2011algorithms], the function $$G_1(\mathbf{H};\tilde{\mathbf{H}}) = - \sum_{k,n} {p}_{kn} \log (h_{kn}) + \sum_{k,n} q_{k} h_{kn},
\label{eq-g1}$$ with the notations $$p_{kn} = \tilde{h}_{kn} \sum_f w_{fk} \frac{v_{fn}}{[\mathbf{W\tilde{H}}]_{fn}}, \quad q_k = \sum_f w_{fk},$$ is a tight auxiliary function of $- \log p(\mathbf{V}|\mathbf{H};\mathbf{W})$ at $\mathbf{H} = \tilde{\mathbf{H}}$. Similarly the function $$G_2(\mathbf{W};\tilde{\mathbf{W}}) = - \sum_{f,k} p'_{fk} \log (w_{fk}) + \sum_{f,k} q'_{k} w_{fk},
\label{eq-g2}$$ with the notations $$p'_{fk} = \tilde{w}_{fk} \sum_n h_{kn} \frac{v_{fn}}{[\mathbf{\tilde{W}H}]_{fn}}, \quad q'_k = \sum_n h_{kn},$$ is a tight auxiliary function of $- \log p(\mathbf{V}|\mathbf{H};\mathbf{W})$ at $\mathbf{W} = \tilde{\mathbf{W}}$.

Finally, for all considered models the function $C$ can be minimized in closed form w.r.t. the variable $\beta_k$.

Minimization w.r.t. **W**
-------------------------

The optimization w.r.t. $\mathbf{W}$ is common to all algorithms, and amounts to minimizing $G_2(\mathbf{W};\tilde{\mathbf{W}})$ only. The scale of $\mathbf{W}$ must be however be fixed in order to prevent potential degenerate solutions such that $\mathbf{W} \rightarrow + \infty$ and $\mathbf{H} \rightarrow 0$. Indeed, consider $\mathbf{W}^{\star}$ and $\mathbf{H}^{\star}$ minimizers of Eq. , and let $\boldsymbol{\Lambda}$ be a diagonal matrix with non-negative entries. Then $$\begin{aligned}
C(\mathbf{W}^{\star} \boldsymbol{\Lambda}^{-1}, \boldsymbol{\Lambda} \mathbf{H}^{\star}) & = - \log p(\mathbf{V}|\boldsymbol{\Lambda} \mathbf{H}^{\star}; \mathbf{W}^{\star} \boldsymbol{\Lambda}^{-1}) - \log p(\boldsymbol{\Lambda} \mathbf{H}^{\star}) \\
& = - \log p(\mathbf{V}|\mathbf{H}^{\star}; \mathbf{W}^{\star}) - \log p(\boldsymbol{\Lambda} \mathbf{H}^{\star}),\end{aligned}$$ and depending on the choice of the prior distribution $p(\mathbf{H})$, we may obtain $C(\mathbf{W}^{\star} \boldsymbol{\Lambda}^{-1}, \boldsymbol{\Lambda} \mathbf{H}^{\star}) < C(\mathbf{W}^{\star}, \mathbf{H}^{\star})$, i.e., a contradiction. Therefore, in the following we impose that $||\mathbf{w}_k||_1 = 1$.

The constrained optimization is performed with the following update rule $$w_{fk} = \frac{p'_{fk}}{\sum_f p'_{fk}}, \label{eq-upw}$$ see Appendix \[app-c\] for the proof.

The following subsections detail the optimization w.r.t. $\mathbf{H}$ (and other variables when necessary), which amounts to the minimization of $G_1(\mathbf{H};\tilde{\mathbf{H}}) - \log p(\mathbf{H})$, as well as the minimization of $C$ w.r.t. $\beta_k$, for each considered model

Chaining on the rate parameter {#chaining-on-the-rate-parameter}
------------------------------

The transition distribution $p(h_{kn}|h_{k(n-1)})$ is given by Eq. . The optimization w.r.t. $h_{kn}$ amounts to solving an order-2 polynomial equation $$a_{2,{kn}} h_{kn}^2 + a_{1,{kn}} h_{kn} + a_{0,{kn}} = 0.
\label{eq-poly-1}$$ The coefficients of the polynomial equation are given in Table \[table-1\]. This bears resemblance with the methodology described in @fevotte2009nonnegative, where the authors aimed at retrieving MAP estimates with a EM-like algorithm (with an exponential likelihood).

[cccc]{} $n$ & $a_{2,{kn}}$ & $a_{1,{kn}}$ & $a_{0,{kn}}$\
$1$ & $q_{k}$ & $\alpha - p_{1k}$ & $-\beta_k h_{k2}$\
\
$2,\dotsc,N-1$ & $q_k + \frac{\beta_k}{h_{k(n-1)}}$ & $1 - p_{kn}$ & $- \beta_k h_{k(n+1)}$\
\
$N$ & 0 & $q_k$ + $\frac{\beta}{h_{k(N-1)}}$ & $1- \alpha - p_N$\

The update for the hyperparameter $\beta_k$ is given by $$\beta_k = \frac{\alpha (N-1)}{\sum_{n \geq 2} \frac{h_{kn}}{h_{k(n-1)}}}.$$

Hierarchical chaining with an auxiliary variable {#hierarchical-chaining-with-an-auxiliary-variable}
------------------------------------------------

In this case, since the transition distribution $p(h_{kn}|h_{k(n-1)})$ is not known in closed form, we resort to optimizing the slighlty more involved following criterion $$\begin{aligned}
& C(\mathbf{W},\mathbf{H},\mathbf{Z}, \boldsymbol{\beta}_h, \boldsymbol{\beta}_z) =  \notag \\
& -\log p(\mathbf{V}|\mathbf{H};\mathbf{W}) - \sum_k \bigg[ \log p(h_{k1}) + \sum_{n \geq 2} \left( \log p(z_{kn}|h_{k(n-1)};\beta_{z,k}) + \log p(h_{kn}|z_{kn};\beta_{h,k}) \right) \bigg],\end{aligned}$$ where $\boldsymbol{\beta}_h = [\beta_{h,1}, \dotsc, \beta_{h,K}]^{\text{T}}$ and $\boldsymbol{\beta}_z = [\beta_{z,1}, \dotsc, \beta_{z,K}]^{\text{T}}$. We recall that $p(z_{kn}|h_{k(n-1)})$ and $p(h_{kn}|z_{kn})$ are given by Eq.  and Eq. , respectively. Note that @cemgil2007conjugate proposed a Gibbs sampler and variational inference, and as such the development of the MAP algorithm is novel.

The update for $z_{kn}$ is given by $$z_{kn} = \frac{\alpha_{z} + \alpha_{h} - 1}{\beta_{z,k} h_{k(n-1)} + \beta_{h,k} h_{kn}}.$$ As for the updates for $h_{kn}$, they are given by $$\begin{aligned}
h_{k1} & = \frac{p_{k1} + \alpha_{z}}{q_k + \beta_{z,k} z_{k2}}, \\
h_{kn} & = \frac{p_{kn} + \alpha_{h} + \alpha_{z} - 1}{q_k + \beta_{h,k} z_{kn} + \beta_{z,k} z_{k(n+1)}} \qquad n \in \{ 2,\dotsc, N \}, \\
h_{kN} & = \frac{p_{kN} + \alpha_{h} - 1}{q_k + \beta_{h,k} z_{kN}}.\end{aligned}$$

As such, imposing $\alpha_{k} \geq 1$ is a sufficient condition to preserve the non-negativity of all the updates. Finally, the update for the parameters $\beta_{z,k}$ and $\beta_{h,k}$ are given by $$\beta_{z,k} = \frac{(N-1)\alpha_{z}}{\sum_{n \geq 2} h_{k(n-1)} z_{kn}}$$ and $$\beta_{h,k} = \frac{(N-1)\alpha_h}{\sum_{n \geq 2} h_{kn} z_{kn}}.$$

Chaining on the shape parameter {#chaining-on-the-shape-parameter}
-------------------------------

The transition distribution $p(h_{kn}|h_{k(n-1)})$ is given by Eq. . The optimization w.r.t. $h_{kn}$ amounts to solving the following equations $$-p_{k1} + (q_k - \alpha \log(\beta_k h_{k2}) + \alpha \Psi(\alpha h_{k1}) )h_{k1} = 0,$$ $$\begin{aligned}
(1-\alpha h_{k(n-1)}-p_{kn}) & + (q_k + \beta_k - \alpha \log(\beta_k h_{k(n+1)})) h_{kn} + \alpha \Psi(\alpha h_{kn}) h_{kn} = 0,\end{aligned}$$ where $\Psi$ denotes the digamma function. Solving such equations can be done numerically with Newton’s method. Finally the update for $h_{kN}$ is given by $$h_{kN} = \frac{p_{kn} + \alpha h_{k(N-1)} -1}{q_k + \beta_k}.$$ The update for $\beta_k$ is given by $$\beta_k = \alpha \frac{\sum_{n \geq 2} h_{k(n-1)}}{\sum_{n \geq 2} h_{kn}}.$$

Note that a Gibbs sampling procedure is proposed in @acharya15nonparametric [@schein2016poisson], and as such the development of the MAP algorithm is novel.

BGAR(1) {#alg-map-bgar}
-------

In this case, since the transition distribution $p(h_{kn}|h_{k(n-1)})$ is not known in closed form, we resort to optimizing the slightly more involved following criterion $$\begin{aligned}
& C(\mathbf{W},\mathbf{H},\mathbf{B}, \boldsymbol{\beta}) = \notag \\
& -\log p(\mathbf{V}|\mathbf{H};\mathbf{W}) - \sum_k \bigg( \log p(h_{k1})~+ \sum_{n \geq 2} \left( \log p(h_{kn}|h_{k(n-1)},b_{kn};\beta_k) + \log p(b_{kn}) \right) \bigg) \label{eq-map-bgar}.\end{aligned}$$

In the following, we will use the notations $\gamma_k = \alpha_k(1-\rho_k)$ and $\eta_k = \alpha_k \rho_k$.

### Constraints

By construction, the variables $h_{kn}$ and $b_{kn}$ must lie in a specific interval given the values of all the other variables. Indeed, as $h_{kn} = b_{kn} h_{k(n-1)} + \epsilon_{kn}$ (Eq. ), where $\epsilon_{kn}$ is a non-negative random variable, we obtain $h_{kn} \geq b_{kn} h_{k(n-1)}$, $b_{kn} \leq \frac{h_{kn}}{h_{k(n-1)}}$, and $h_{kn} \leq \frac{h_{k(n+1)}}{b_{k(n+1)}}$.

This leads to the following constraints $$\begin{aligned}
0 & \leq  h_{k1} \leq \frac{h_{k2}}{b_{k2}}, \\
b_{kn} h_{k(n-1)} & \leq h_{kn} \leq \frac{h_{k(n+1)}}{b_{k(n+1)}} \qquad 2 \leq n < N, \\
b_{kN} h_{k(N-1)} & \leq h_{kN},\end{aligned}$$ and $$0 \leq b_{kn} \leq \min \left( 1, \frac{h_{kn}}{h_{k(n-1)}} \right).$$

We therefore introduce the notations $$\begin{aligned}
c_{kn} = b_{kn} h_{k(n-1)}, \quad d_{kn} = \frac{h_{k(n+1)}}{b_{k(n+1)}}, \quad x_{kn} = \frac{h_{kn}}{h_{k(n-1)}},\end{aligned}$$ as these quantities arise naturally in our derivations.

### Minimization w.r.t. $h_{kn}$

The optimization of Eq.  w.r.t. $h_{kn}$ may give rise to intractable problems, due to the logarithmic terms in the objective function. To alleviate this issue, we propose to control the limit values of the auxiliary function, by restricting ourselves to certain values of the hyperparameters. In particular, choosing $(1-\gamma_k) < 0$ ensures the existence of at least one minimizer.

In all sub-cases, the optimization w.r.t. $h_{kn}$ amounts to solving an order-3 polynomial equation $$a_{3,{kn}} h_{kn}^3 + a_{2,{kn}} h_{kn}^2 + a_{1,{kn}} h_{kn} + a_{0,{kn}} = 0.
\label{eq-poly-2}$$ The coefficients of the polynomial equation are given in Table \[table-2\]. If several roots belong to the definition interval, we simply choose the root which gives the lowest objective value.

### Minimization w.r.t. $b_{kn}$

Similarly, logarithmic terms of the objective function may give rise to degenerate solutions. Using the same reasoning, we choose to impose $(1-\gamma_k) < 0$ and $(1-\eta_k) < 0$ to ensure the existence of at least one minimizer.

The minimization of the auxiliary function w.r.t. $b_{kn}$ amounts to solving the following order 3 polynomial over the interval $[0, \min(1, x_{kn})]$ $$a_{3,kn} b_{kn}^3 + a_{2,kn} b_{kn}^2 + a_{1,kn} b_{kn} + a_{0,kn} d_{kn},$$ where $$\begin{aligned}
a_{3,kn} & = -\beta_k h_{k(n-1)}, \\
a_{2,kn} & = 2(1-\gamma_k) + (1-\eta_k) + \beta_k h_{k(n-1)}(x_{kn}+1), \\
a_{1,kn} & = -(1-\gamma_k)(x_{kn}+1) - (1-\eta_k)(x_{kn} + 1) - \beta_k h_{k(n-1)} x_{kn}, \\
a_{0,kn} & = (1-\eta_k) x_{kn}.\end{aligned}$$

### Minization w.r.t. $\beta_k$

The minimization of $C$ w.r.t. $\beta_k$ can be done in closed form and results in the following update rule $$\beta_k = \frac{(N-1) \alpha (1-\rho)}{\sum_{n \geq 2} (h_{kn} - b_{kn}h_{k(n-1)})}.$$

### Admissible values of hyperparameters

To recap the discussion on admissible values of hyperparameters, to ensure the existence of minimizers of the auxiliary function, we have restricted ourselves to $$\left\{
\begin{array}{l}
\alpha_k(1-\rho_k) > 1 \\
\alpha_k \rho_k > 1
\end{array}
\right.$$ This set is graphically displayed on Figure \[fig-admv\]. As we can see, choosing the value of $\rho_k$ to be close to one (to ensure correlation) leads to high values of $\alpha_k$.

Experimental work {#sec-exp}
=================

We now compare the performance of all considered temporal NMF models on a prediction task on three real datasets. This task will consist in hiding random columns of the considered datasets and predicting those missing values. We will also include the performance of a naive baseline, which we detail in the following subsection. Adapting the MAP algorithms presented in Section \[sec-study\] in a setting with a mask of missing values only consist in a slight modification, presented in Appendix \[app-d\]. Python code will be made available upon acceptance.

Experimental protocol
---------------------

For each considered dataset, the experimental protocol is as follows.

First of all, a value of the factorization rank $K$ (which will be used for all considered methods) must be selected. To do so, we apply the standard KL-NMF algorithm [@lee2000algorithms; @fevotte2011algorithms] on 10 random training sets, which consist of $80 \%$ of the original data, with a pre-defined grid of values for $K$. We then select the value of $K$ which yields the lowest generalized Kullback-Leibler error (KLE) (see definition below) on the remaining $20 \%$ of the data.

For the prediction experiment itself, we create 5 random splits of the data matrix, where $80 \%$ corresponds to the training set, $10 \%$ to the validation set, and the remaining $10 \%$ to the test set. To do so, we randomly select non-adjacent columns of the data matrix (excluding the first one and the last one), half of which will make up the validation set and the other half the test set. We also consider 5 different random initializations.

Thus, for each split-initialization pair, all the algorithms are run from this initialization point on the training set until convergence (the algorithms are stopped when the relative decrease of the objection function falls under $10^{-5}$). For each method, a grid of shape hyperparameters is considered, and the selection of this parameter is based on the lowest KLE on the validation set. The predictive performance of each method is then computed on the test set by comparing the original value $v_{fn}$ and its associated estimate $\hat{v}_{fn} = [\mathbf{WH}]_{fn}$ with two different metrics. Denoting by $\mathcal{T}$ the test set, we consider

-   the generalized Kullback-Leibler error (KLE) $$\text{KLE} = \sum_{(f,n) \in \mathcal{T}} \left[ v_{fn} \log \left( \frac{v_{fn}}{\hat{v}_{fn}} \right) - v_{fn} + \hat{v}_{fn} \right];$$

-   the relative error, as in @schein2016poisson $$\text{RE} = \sum_{(f,n) \in \mathcal{T}} \frac{|v_{fn} - \hat{v}_{fn}|}{v_{fn} + 1}.$$

Finally, we compare the NMF-based approaches to the following naive baseline, based on a random guess. In this case, the values of the missing columns are simply estimated by drawing $\hat{v}_{fn}$ from the empirical distribution of the observed data coefficient in every row $f$.

Datasets
--------

The following datasets are considered

-   The `NIPS` dataset[^3], which contains word counts (with stop words removed) of all the articles published at the NIPS[^4] conference between 1987 and 2015. We grouped the articles per year, yielding an observation matrix of size $11463 \times 29$. We obtained $K = 3$.

-   The `ICEWS` dataset[^5], an international relations dataset, which contains the number of interactions between two countries for each day of the year 2003. The matrix is of size $6197 \times 365$. We obtained $K = 15$.

-   The `last.fm` dataset, based on the so-called “last.fm 1K” users[^6], which contains the listening history of users with timestamps information. We preprocessed this dataset to obtain the monthly evolution of the listening counts of artists with at least 20 different listeners. This yields a dataset of size $7017 \times 53$. We obtained $K = 6$.

Experimental results
--------------------

The averaged KLE and RE over the 25 split-initialization pairs are reported on Table \[table-nips\] for the `NIPS` dataset, on Table \[table-icews\] for the `ICEWS` dataset, and on Table \[table-last\] for the `last.fm` dataset.

First of all, on all the considered datasets, the naive baseline yields the worst performance results, for all metrics, as expected. Moreover, all temporal models achieve comparable predictive performance, and the slight advantage of some methods over the others being data-dependent.

          Model                                KLE                                                RelE
  --------------------- -------------------------------------------------- --------------------------------------------------
        Baseline             $14.8 \times 10^5 \pm 57.4 \times 10^4$            $7.45 \times 10^4 \pm 17.9 \times 10^3$
       Rate (II.A)       $\mathbf{1.07 \times 10^5 \pm 2.01 \times 10^4}$   $\mathbf{2.98 \times 10^4 \pm 2.53 \times 10^3}$
   Hierarchical (II.B)   $\mathbf{1.07 \times 10^5 \pm 2.02 \times 10^4}$   $\mathbf{2.98 \times 10^4 \pm 2.45 \times 10^3}$
      Shape (II.C)           $1.29 \times 10^5 \pm 2.75 \times 10^4$        $\mathbf{3.07 \times 10^4 \pm 4.07 \times 10^3}$
       BGAR (II.D)       $\mathbf{1.05 \times 10^5 \pm 1.74 \times 10^4}$   $\mathbf{3.02 \times 10^4 \pm 2.00 \times 10^3}$

          Model                                KLE                                                RelE
  --------------------- -------------------------------------------------- --------------------------------------------------
        Baseline             $99.1 \times 10^4 \pm 44.2 \times 10^3$            $3.33 \times 10^4 \pm 4.06 \times 10^2$
       Rate (II.A)       $\mathbf{8.68 \times 10^4 \pm 2.95 \times 10^3}$   $\mathbf{2.65 \times 10^4 \pm 5.93 \times 10^2}$
   Hierarchical (II.B)   $\mathbf{8.80 \times 10^4 \pm 3.39 \times 10^3}$   $\mathbf{2.62 \times 10^4 \pm 7.80 \times 10^2}$
      Shape (II.C)       $\mathbf{8.79 \times 10^4 \pm 3.31 \times 10^3}$   $\mathbf{2.61 \times 10^4 \pm 6.90 \times 10^2}$
       BGAR (II.D)           $9.09 \times 10^4 \pm 4.16 \times 10^3$        $\mathbf{2.54 \times 10^4 \pm 5.30 \times 10^2}$

          Model                                KLE                                                RelE
  --------------------- -------------------------------------------------- --------------------------------------------------
        Baseline             $66.5 \times 10^4 \pm 1790 \times 10^2$            $3.40 \times 10^4 \pm 89.9 \times 10^2$
       Rate (II.A)       $\mathbf{1.55 \times 10^4 \pm 7.13 \times 10^2}$   $\mathbf{1.18 \times 10^4 \pm 5.55 \times 10^2}$
   Hierarchical (II.B)   $\mathbf{1.57 \times 10^4 \pm 6.18 \times 10^2}$   $\mathbf{1.17 \times 10^4 \pm 6.06 \times 10^2}$
      Shape (II.C)           $2.01 \times 10^4 \pm 51.9 \times 10^2$        $\mathbf{1.19 \times 10^4 \pm 7.16 \times 10^2}$
       BGAR (II.D)       $\mathbf{1.59 \times 10^4 \pm 9.10 \times 10^2}$   $\mathbf{1.20 \times 10^4 \pm 6.86 \times 10^2}$

Conclusion {#sec-ccl}
==========

In this paper, we have reviewed existing temporal NMF models in a unified MAP framework and introduced a new one. These models differ by the choice of the Markov chain structure used on the activation coefficients to induce temporal correlation. We began by studying the previously proposed Gamma Markov chains of the NMF literature, only to find that they all share the same drawback, namely the absence of a well-defined stationary distribution. This leads to problematic behaviors from the generative perspective, because the realizations of the chains are degenerate (although this is not necessarily a problem in MAP estimation). We then introduced a Markovian process from the time series literature, called BGAR(1), which overcomes this limitation, and which, to the best of our knowledge, had never been exploited for learning tasks.

We then derived MAP estimation algorithms for all these models, in the context of a Poisson likelihood, which allowed for a comprehensive comparison on a prediction task on real datasets. As it turns out, we cannot claim that there is a single model which outperforms all the others. Their strengths and weaknesses appear to depend on the nature of the data at hand, which is reasonable.

Future work will focus on finding a way to perform inference with the BGAR prior for a less restrictive set of hyperparameters, which might increase the performance of this particular model. We will also work to derive similar algorithms in the context of different likelihoods, such as an exponential likelihood [@fevotte2009nonnegative], which can be easily done thanks to the MM framework.

Acknowledgments {#acknowledgments .unnumbered}
===============

This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grant agreement No 681839 (project FACTORY).

The Beta-Prime distribution {#app-a}
===========================

Distribution for a continuous random variable in $[0,+\infty[$, with parameters $\alpha > 0$, $\beta > 0$, $p > 0$ and $q > 0$. Its p.d.f. writes, for $x \geq 0$: $$f(x;\alpha,\beta,p,q) = \frac{p \left( \frac{x}{q} \right)^{\alpha p -1} \left(1 + \left(\frac{x}{q}\right)^p \right)^{-\alpha - \beta}}{q \text{B}(\alpha, \beta)}.$$

BGAR(1) linear correlation {#app-b}
==========================

We have between two successive values $h_n$ and $h_{n+1}$: $$\begin{aligned}
& \text{corr}(h_n, h_{n+1}) \\
& = \frac{\mathbb{E}(h_n h_{n+1}) - \mathbb{E}(h_n)\mathbb{E}(h_{n+1})}{\sigma(h_n)\sigma(h_{n+1})} \\
& = \frac{\mathbb{E}(h_n(b_{n+1}h_n + \epsilon_{n+1})) - \mathbb{E}(h_n)\mathbb{E}(h_{n+1})}{\sigma(h_n)\sigma(h_{n+1})} \\
& = \frac{\mathbb{E}(b_{n+1})\mathbb{E}(h_n^2) + \mathbb{E}(h_n)\mathbb{E}(\epsilon_{n+1}) - \mathbb{E}(h_n)\mathbb{E}(h_{n+1})}{\sigma(h_n)\sigma(h_{n+1})} \\
& = \frac{\frac{\alpha \rho}{\alpha \rho + \alpha(1-\rho)} \frac{\alpha(\alpha+1)}{\beta^2} + \frac{\alpha}{\beta} \frac{\alpha(1-\rho)}{\beta} - \frac{\alpha}{\beta} \frac{\alpha}{\beta} }{\frac{\alpha}{\beta^2}} \\
& = \rho.\end{aligned}$$

Constrained optimization {#app-c}
========================

We want to optimize $G_2(\mathbf{W};\tilde{\mathbf{W}})$ w.r.t. $\mathbf{W}$ s.t. $\sum_f{w_{fk}} = 1$. Rewriting this with Lagrange multipliers $\boldsymbol{\lambda} = [\lambda_1,\dotsc,\lambda_{K}]^{\text{T}}$, this is tantamount to $$\min_{\mathbf{W}, \boldsymbol{\lambda}} G_2(\mathbf{W};\tilde{\mathbf{W}}) + \sum_k \lambda_k (||\mathbf{w}_k||_1 - 1).$$ Deriving w.r.t $w_{fk}$ yields $$w_{fk} = \frac{p'_{fk}}{q'_k + \lambda_k}. \label{eq-lagr}$$ We retrieve the constraint by summing this expression over $f$. This gives the expression of the Lagrange multiplier: $\lambda_k = \sum_f p'_{fk} - q'_k$. Substituting this expression into Eq. , we obtain the following update rule $$w_{fk} = \frac{p'_{fk}}{\sum_f p'_{fk}}.$$

Algorithms with missing values {#app-d}
==============================

In the context of missing values, let us consider a mask matrix $\mathbf{M}$ of size $F \times N$ such that $m_{fn} = 1$ if the entry $v_{fn}$ is observed and 0 otherwise. The likelihood term can then be written as $$-\log p(\mathbf{V}|\mathbf{H};\mathbf{W}) = - \sum_{f,n} m_{fn} \log p(v_{fn}|[\mathbf{WH}]_{fn}).$$ The auxiliary function $G_1$ of Eq.  and $G_2$ of Eq.  can then be written is the same way, with $$p_{kn} = \tilde{h}_{kn} \sum_f w_{fk} \frac{m_{fn} v_{fn}}{[\mathbf{W\tilde{H}}]_{fn}}, \quad q_{kn} = \sum_f m_{fn} w_{fk},$$ for $G_1$, and $$p'_{fk} = \tilde{w}_{fk} \sum_n h_{kn} \frac{m_{fn} v_{fn}}{[\mathbf{\tilde{W}H}]_{fn}}, \quad q'_{kn} = \sum_n m_{fn} h_{kn},$$ for $G_2$.

[^1]: These models are sometimes generically referred to as “Poisson factorization” or “Poisson factor analysis”.

[^2]: Note that this particular mapping is actually slightly more complex, as the $K$-dimensional real vector must be mapped to the $(K-1)$ simplex due to further constraints in the model.

[^3]: <https://archive.ics.uci.edu/ml/datasets/NIPS+Conference+Papers+1987-2015>

[^4]: Now called NeurIPS.

[^5]: <https://github.com/aschein/pgds>

[^6]: <http://ocelma.net/MusicRecommendationDataset/>
