---
kind: concept
slug: non-negative-matrix-factorization
title: Non-negative matrix factorization
summary: >-
  Fitting non-negative matrix factorisation to count data, especially Poisson NMF, by alternating Poisson regressions with EM, coordinate-descent or extrapolated updates. The same fits give maximum-likelihood topic models after a change of variables. One card covers it; speed results are author-reported.
tags: [matrix-factorization, optimization, topic-models]
papers: ["2105.13440"]
related: [topic-modelling, expectation-maximization, coordinate-descent]
status: draft
updated_at: 2026-11-08
---

## Overview

Non-negative matrix factorisation (NMF) writes a non-negative matrix as the
product of two non-negative factor matrices with a small inner dimension K.
For count data, the Poisson version models each count x_ij as Poisson with a
rate equal to the sum over k of h_ik w_jk [arxiv:2105.13440 §2]. Fitting it
means finding non-negative H and W that maximise the likelihood.

The core idea of the cited paper is a link to topic models. A multinomial
topic model, whose rows and columns each sum to one, has a maximum-likelihood
estimate that matches the Poisson NMF estimate after a change of variables
[arxiv:2105.13440 §2]. The mapping is: normalise the columns of W and keep
their totals; scale the columns of H by those totals; normalise the rows and
keep their totals [arxiv:2105.13440 §2]. Because the topic-model constraints
disappear in the NMF form, NMF optimisers can fit topic models directly
[arxiv:2105.13440 §2] [arxiv:2105.13440 §3].

Most solvers fit one factor matrix at a time. Each row of H is the maximum-
likelihood estimate of an additive Poisson regression whose design matrix is
W; then each row of W is fitted with the new H fixed. The row problems are
independent [arxiv:2105.13440 §3]. Variants differ only in the inner solver
and in any acceleration added to the outer loop.

Raw H and W are unique only up to rescaling each column pair, so raw factors
from two fits cannot be compared column by column. The normalised L and F from
the mapping can be compared [arxiv:2105.13440 §2].

## Variants

### Multiplicative updates (EM form)

The classic multiplicative update for Poisson NMF is, after rewriting, the EM
algorithm for the Poisson regression inside each block. The E-step computes
expected counts for each component; the M-step divides their sum by the column
sum of the design matrix [arxiv:2105.13440 §3]. Applied through the mapping,
the same updates become EM for the topic model [arxiv:2105.13440 §3]
[arxiv:2105.13440 tab. 1]. The update has no step size to tune.

Its weakness, in the paper's example, is slow progress far from the optimum.
On MCF-7 (three topics), 200 EM updates from the same start left the estimate
more than 150,000 log-likelihood units from the maximum-likelihood estimate.
EM eventually reached it after many more updates [arxiv:2105.13440 §4].

Card: [Carbonetto et al., 2021](/papers/2105.13440).

### Coordinate descent (CCD and SCD)

Each component is updated with a Newton-type step, using the first and second
derivatives of the negative log-likelihood. The step size is fixed at one and
the result is clipped at zero. There is no line search, so a single step can
increase the objective [arxiv:2105.13440 §3]. The paper runs one EM update
before each coordinate-descent update to reduce such failures
[arxiv:2105.13440 app. B].

In the same MCF-7 example, 200 coordinate-descent updates were within 0.079
log-likelihood units of the maximum-likelihood estimate [arxiv:2105.13440 §4].
This is an author-reported result on one data set.

Card: [Carbonetto et al., 2021](/papers/2105.13440).

### Extrapolated updates

An extrapolation step moves beyond the new estimate by a fraction β of its
last change, then clips at zero. β lies in [0, 1], and β = 0 gives the plain
update. The scheme follows Ang and Gillis (2019) [arxiv:2105.13440 app. B].
The paper turns it on at iteration 50 of coordinate descent on MCF-7
[arxiv:2105.13440 §4]. The step size is adapted to performance, but the paper
does not give the schedule [arxiv:2105.13440 app. B]. The paper says the extra
cost is minimal [arxiv:2105.13440 §3].

In the authors' comparisons on four count data sets, with K from 2 to 12,
extrapolated coordinate descent converged at least as fast as the other
variants, often much faster, or reached the best fit in the time allowed. The
paper gives these results as plots only, not as a table of speed-ups
[arxiv:2105.13440 §4].

Card: [Carbonetto et al., 2021](/papers/2105.13440).

## When to use which

| Situation | Use | Why |
|---|---|---|
| Simple code, no step size to tune, and slow convergence is acceptable | Multiplicative or EM updates | No step size; the same update is EM for the topic model [arxiv:2105.13440 §3]. Slow far from the optimum [arxiv:2105.13440 §4]. |
| Topic-model fit needed quickly, and extra code is acceptable | Extrapolated coordinate descent | Usually as fast as or faster than the alternatives on the paper's four data sets. Author-reported [arxiv:2105.13440 §4]. |
| Plain coordinate descent with no line search | Coordinate descent, with an EM step before each sweep | The objective need not fall at every step [arxiv:2105.13440 §3]. The paper pairs it with one EM update [arxiv:2105.13440 app. B]. |
| Sparse, large count matrices | Any variant, with a loss that scales with non-zeros | The paper's loss costs about (N + n + m)K operations for N non-zeros, against nmK when dense. Its data sets had 0.2% to 9.3% non-zeros [arxiv:2105.13440 app. B] [arxiv:2105.13440 tab. 2]. |
| Posterior uncertainty or full Bayesian inference | Not covered by this card | The paper treats point estimates only: the maximum-likelihood estimate, and a MAP estimate with Dirichlet priors on F [arxiv:2105.13440 §5] [arxiv:2105.13440 app. A]. |
| Comparing topics from two fits | Compare the mapped L and F | Raw H and W are determined only up to column scaling [arxiv:2105.13440 §2]. |

## Implementation sketch

The common core is the alternating loop. The inner solver is the only part
that changes between variants. Python-flavoured, with citations above each
step.

```pseudocode
# Alternate block updates for Poisson NMF  [arxiv:2105.13440 §3]
# Input:  X (n, m) counts; H (n, K), W (m, K) positive start
# Output: H (n, K), W (m, K), both non-negative

def fit_pnmf(X, H, W, outer_iters, solver):
    for t in range(outer_iters):
        # Refit each row of H with W fixed  [arxiv:2105.13440 §3]
        for i in range(n):
            H[i] = solver(X[i, :], W, H[i])     # (K,)
        # Refit each row of W with the new H fixed  [arxiv:2105.13440 §3]
        for j in range(m):
            W[j] = solver(X[:, j], H, W[j])     # (K,)
    return H, W
```

Inner solver A, one EM step for an additive Poisson regression
[arxiv:2105.13440 §3]:

```pseudocode
# One EM step for y ~ Pois(A @ b), b >= 0  [arxiv:2105.13440 §3]
def em_step(y, A, b):
    # Fitted means  [arxiv:2105.13440 §3]
    mu = A @ b                                   # (n,)
    # Expected counts per component  [arxiv:2105.13440 §3]
    z = y[:, None] * A * b[None, :] / mu[:, None]  # (n, K)
    # New b: expected counts over the column sums of A  [arxiv:2105.13440 §3]
    return z.sum(axis=0) / A.sum(axis=0)         # (K,)
```

Inner solver B, one coordinate-descent sweep with a unit step and no line
search [arxiv:2105.13440 §3]:

```pseudocode
# One coordinate-descent sweep for y ~ Pois(A @ b)  [arxiv:2105.13440 §3]
def cd_sweep(y, A, b):
    for k in range(K):
        # Fitted means at the current b  [arxiv:2105.13440 §3]
        mu = A @ b                               # (n,)
        # First and second derivatives of the negative log-likelihood
        g = (A[:, k] * (1 - y / mu)).sum()       # scalar
        q = (y * A[:, k] ** 2 / mu ** 2).sum()   # scalar
        # Unit Newton step, clipped at zero  [arxiv:2105.13440 §3]
        b[k] = max(0.0, b[k] - g / q)
    return b
```

Extrapolation of a non-negative estimate [arxiv:2105.13440 app. B]:

```pseudocode
# Extrapolate past the new estimate, then project  [arxiv:2105.13440 app. B]
def extrapolate(new, previous, beta):
    # beta in [0, 1]; beta = 0 gives the plain update  [arxiv:2105.13440 app. B]
    candidate = new + beta * (new - previous)    # (n, K)
    return np.maximum(candidate, 0.0)            # (n, K)
```

Mapping to a topic model [arxiv:2105.13440 §2]:

```pseudocode
# Poisson NMF factors to topic-model parameters  [arxiv:2105.13440 §2]
def pnmf_to_topic_model(H, W):
    # Column-normalise W, keeping the column totals  [arxiv:2105.13440 §2]
    u = W.sum(axis=0)                            # (K,)
    F = W / u                                    # (m, K)
    # Row-normalise H * u, keeping the row totals  [arxiv:2105.13440 §2]
    scaled = H * u                               # (n, K)
    s = scaled.sum(axis=1)                       # (n,)
    L = scaled / s[:, None]                      # (n, K)
    return L, F, s
```

Implementation notes. Drop documents with zero total before the mapping,
because the inverse is unique only when every total is positive. Do not let
column sums of W reach zero, since the mapping divides by them
[arxiv:2105.13440 §2]. Near-zero fitted means make the updates unstable; the
paper floors estimates at 1e-15 [arxiv:2105.13440 app. B].

## Contradictions

None between papers: only one card lists this concept, so no two sources
disagree. Within the paper, the speed results are author-reported and shown
mainly as plots [arxiv:2105.13440 §4]. The paper says other acceleration
methods (DAAREM and quasi-Newton EM) bring little or no benefit, but it does
not show those results [arxiv:2105.13440 §5]. The released code differs from
the paper in its extrapolation schedule and default setting; the differences
are listed on the card. [Carbonetto et al., 2021](/papers/2105.13440)

## Papers

- [Non-negative matrix factorization algorithms generally improve topic model fits](/papers/2105.13440) (Carbonetto et al., 2021): shows that Poisson NMF and the multinomial topic model share maximum-likelihood fits under a change of variables, and compares EM with coordinate descent and extrapolation on four count data sets.