# Avoid evaluation of a very large matrix in non-negative matrix factorization

This is somewhere in between a math and a programming question, so please send me back to SO if you think it's off-topic.

I'm implementing non-negative sparse coding, a regularized variant of non-negative matrix factorization. This entails finding $W \in \mathbb{R}^{n\times k}$ and $H \in \mathbb{R}^{k\times m}$ that minimize

$\frac{1}{2} ||X - WH||_2^2 + \lambda ||W||_1$.

$k$ is a free parameter, but $n,m$ depend on the data and can be very large, so that evaluating $WH$ is prohibitively expensive. I can avoid it by evaluating the first term, the squared Frobenius norm, as

$||X||_2^2 + \operatorname{tr}(((W^T W)H)H^T) - 2\operatorname{tr}((X H^T)W^T)$,

which can be evaluated cheaply by noticing that

$\operatorname{tr}(A B^T) = \operatorname{vec}(A)^T\operatorname{vec}(B)$.

Now comes the problem: as part of the algorithm, I have to evaluate a gradient update

$W - \mu (WH - X)H^T$,

which again involves $WH$. $\mu$ is a learning rate, and the subtraction is element-wise. However, I don't see how I can avoid evaluating $WH$ this time; is there a way to rewrite this so I don't have to?

(I could evaluate it piece by piece, but I'm really looking for an expression that lends itself to computation using a Matlab-like linear algebra library.)

• Salvation may lie in how you plan to use $W - \mu(WH-X)H^T$, in that it may be possible to use it, e.g. for matrix-vector multiplication, in its "expanded" form. – hardmath Nov 5 '13 at 10:39
• Is $k$ small, relative to $n,m$? – hardmath Nov 5 '13 at 10:50
• @hardmath: $n,m$ may be on the order of millions; $k$ is user-specified and typically not more than a few hundred. – Fred Foo Nov 5 '13 at 12:07

Forming each entry of $WH$ requires $k$ multiplications (and $k-1$ adds) using a conventional matrix multiplication, for a total of $kmn$ multiplications.

Alternatively one can form $HH^T$ with $k^2 m$ multiplications, and then $W(HH^T)$ requires a further $k^2 n$ multiplications.

So if you have to evaluate it, I'd try:

$$(WH - X)H^T = W(HH^T) - XH^T$$

since $m+n$ is much less than $mn$.

Still $XH^T$ looks expensive to evaluate. Perhaps this piece should be left unevaluated/accessed implicitly (by storing $X$ and $H^T$ stepwise).

Comparing the notation to the article linked in the Question, presumably $H$ (and therefore $H^T$) are sparse but in a random fashion. It's possible that this sparsity can be exploit to perform $XH^T$ quickly, but this is not likely the case with a high-level matrix library/interface similar to MATLAB.

• $H$ is initially very dense (it's initialized as a random, strictly positive matrix), but by keeping $X$ in a sparse format $XH^T$ is feasible. Thanks! – Fred Foo Nov 5 '13 at 13:37