nonnegMatrixFactorization

Given an $m \times n$ real matrix $A \geq 0$ , and an integer $k \leq \min (m, n)$ , compute a factorization $A \cong FG$. The matrix factors $F_{m \times k} \geq 0, G_{k \times n} \geq 0$ are computed to minimize the Frobenius, or sum of squares, norm of the error matrix: $E = \left\{ e_{i,j} \right\} = A - FG$ .

Synopsis

nonnegMatrixFactorization (a, f, g)

Required Arguments

float a[[]] (Input)

An array of length m × n containing the A matrix.

float f[[]] (Input/Output)

An array of length m × k containing the F matrix. If initialFactors is used, the sweeps begin using the input values for $F_{m \times k} \geq 0$ .

float g[[]] (Output)

An array of length k × n containing the G matrix.

Return Value

A scalar containing the Frobenius norm of the error matrix

$$E:\mathit{error} = \left( \sum_{i,j} e_{i,j}^2 \right)^{1/2}$$

Optional Arguments

weight, float (Input)

An array of length m × n containing the matrix W ≥ 0 of weights that will be applied to the entries of A ≥ 0 during the solution sweeps. The factorization obtained is FG ≅W ∘ A, where the weights are applied element-wise.

Default: Weights are not applied, or equivalently, the weights all have value 1.

initialFactors, int (Input)

A flag that signifies if the matrix F is given an input estimate. If initialFactors = 0, start sweeps using

$$\begin{split}F = \begin{bmatrix} I_k \\ 0 \end{bmatrix}\end{split}$$

Otherwise, use initial values in f as the matrix F to start the sweeps.

Default: initialFactors = 0

itmax, int (Input)

The maximum number of sweeps allowed for alternately solving for G ≥ 0, then F ≥ 0.

Default: itmax = 2 * (m + n + 1)

residualError, float (Input)

A scalar that will stop the sweeps at the first one satisfying error ≤ residualError.

Default: residualError = 0

relativeError, float (Input)

A scalar that will stop the sweeps at the first one satisfying

$$\mathit{error}_{iter‑2} ‑ \mathit{error}_{iter‑1} \leq \mathit{relativeError} \times \mathit{error}_{iter}, iter > 2.$$

This test is made after three values of the error matrix norm have been computed. The sequence {$error_{iter}$} is decreasing with increasing values of the iteration counter, iter. If $error_{iter} \geq error_{iter‑1}$ occurs, the sweeps stop.

Default: relativeError = (machine$(3))^{0.4}$.

stoppingCriterion (Output)

This flag has the value 0,1,2 or 3 depending on which of the following conditions stopped the sweeps:

stoppingCriterion	Description
0	Errors in user input occurred
1	Reached maximum iterations
2	Residual norm is small
3	Relative error convergence

nstepsTaken (Output)

The last value of the iteration count that gives the number of sweeps.

Description

Function nonnegMatrixFactorization computes an approximation $A \cong FG$ , or with weights, W ∘ A ≅ FG; the factors are constrained: $F_{m \times k} \geq 0, G_{k \times n} \geq 0$ . The matrix factors $F_{m \times k} \geq 0, G_{k \times n} \geq 0$ are computed to minimize the Frobenius or sum of squares, norm of the error matrix: $E = \left\{ e_{i,j} \right\} = A - FG$ .

The algorithm is based on Alternating Least Squares, presented by P. Paatero and U. Tapper, “Positive Matrix Factorization, etc.” Environmetrics, (5), p. 111-126 (1994).

Each constrained least squares problem is solved using nonnegLeastSquares. This process alternates between computing the batch of $n$ columns of $G$ and then the batch of $m$ rows of $F$. This constitutes a “sweep.”

There is no restriction on the relative sizes of $m$ and $n$ . The restrictions on the integer $k$ are $0 < k \leq \min (m,n)$. When an initial matrix $G$ is to be used, instead of an initial $F$, repose the factorization in transposed form $A^T \cong G^T F^T$, or with weights, $A^T \circ W^T \cong G^T F^T$.

The matrix factors $F, G$ are not unique. In the function, the output rows of $G$ are scaled to have sum equal to the value 1. The scaled columns of $F$ are sorted so the column sums are non-increasing. This sort order is then applied to the rows of $G$ .

Example

Five customers, Beth, Dick, Fred, Joe and Kay make purchases at a convenience store.

	Flour	Balloons	Beer	Sugar	Chips
Beth		3	8		1
Dick		2	5	1
Fred	5		1	10
Joe		20	40	2	1
Kay	10		1	10	1

This matrix $A_{5 \times 5}$ of customers versus items purchased is approximated by a non-negative matrix factorization, using $k = 2$ : $A \cong FG$ . The example is taken from one due to H. Jin and M. Saunders, “Exploring Nonnegative Matrix Factorization,” A Workshop on Algorithms for Massive Data Sets, Stanford University, June 25-28, (2008).

from __future__ import print_function
from numpy import *
from pyimsl.math.nonnegMatrixFactorization import nonnegMatrixFactorization
from pyimsl.math.writeMatrix import writeMatrix

m = 5
n = 5
k = 2
a = array([[0, 3, 8, 0, 1],
           [0, 2, 5, 1, 0],
           [5, 0, 1, 10, 0],
           [0, 20, 40, 2, 1],
           [10, 0, 1, 10, 1]])
f = zeros((m, k))
g = zeros((k, n))
reason = []
nsteps = []

# Solve for factors, constraining values to be non-negative.
# Get reason for stopping and number of sweeps.
error = nonnegMatrixFactorization(a, f, g,
                                  stoppingCriterion=reason,
                                  nstepsTaken=nsteps)

writeMatrix("Matrix Factor F", f)
writeMatrix("Matrix Factor G", g)
print("\nFrobenius Norm of E=A-F*G is %e" % error)
print("Reason for stopping sweeps: %d" % reason[0])
print("Number of sweeps taken: %d" % nsteps[0])

Output

Frobenius Norm of E=A-F*G is 3.195012e+00
Reason for stopping sweeps: 3
Number of sweeps taken: 18
 
      Matrix Factor F
             1            2
1        11.98         0.00
2         7.52         0.93
3         0.34        16.60
4        62.97         0.04
5         0.00        21.35
 
                          Matrix Factor G
             1            2            3            4            5
1       0.0000       0.3146       0.6366       0.0309       0.0179
2       0.4049       0.0000       0.0472       0.5191       0.0288