Non-negative matrix factorization

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 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 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 2105.13440 §2. Because the topic-model constraints disappear in the NMF form, NMF optimisers can fit topic models directly 2105.13440 §2 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 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 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 2105.13440 §3. Applied through the mapping, the same updates become EM for the topic model 2105.13440 §3 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 2105.13440 §4.

Card: Carbonetto et al., 2021.

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 2105.13440 §3. The paper runs one EM update before each coordinate-descent update to reduce such failures 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 2105.13440 §4. This is an author-reported result on one data set.

Card: Carbonetto et al., 2021.

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) 2105.13440 app. B. The paper turns it on at iteration 50 of coordinate descent on MCF-7 2105.13440 §4. The step size is adapted to performance, but the paper does not give the schedule 2105.13440 app. B. The paper says the extra cost is minimal 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 2105.13440 §4.

Card: Carbonetto et al., 2021.

When to use which

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.

# 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 2105.13440 §3:

# 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 2105.13440 §3:

# 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 2105.13440 app. B:

# 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 2105.13440 §2:

# 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 2105.13440 §2. Near-zero fitted means make the updates unstable; the paper floors estimates at 1e-15 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 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 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

• Non-negative matrix factorization algorithms generally improve topic model fits (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.