Sparse non-negative matrix factorization using 1) adaptive overrelaxed multiplicative updates and 2) non-negative least squares with block principal pivoting.

function [A,B] = fastnmf(X,N,varargin)
FASTNMF Fast least-squares non-negative matrix factorization 
Approximates X by A*B s.t. A,B>=0
Usage
[A,B] = fastnmf(X,N,[options])
Input
  X           Data matrix (I x J)
  N           Number of components 
  options
    maxiter   Number of iterations (default 100)
    iiter     Number of inner iterations (default 10)
    A         Initial value for A (I x N)
    B         Initial value for B (N x J)
Output
  A           Samples of A (I x N x M)
  B           Samples of B (N x J x M)

function [W,H,cost] = isp_snmf(V, d, varargin)
ISP_SNMF: SPARSE NON-NEGATIVE MATRIX FACTORIZATION 
Usage
  [W,H] = isp_snmf(V,d,[options])
Input
  V                 M x N data matrix
  d                 Number of factors
  options
    .costfcn        Cost function to optimize
                      'ls': Least squares (default)
                      'kl': Kullback Leibler
    .W              Initial W, array of size M x d
    .H              Initial H, array of size d x N 
    .lambda         Sparsity weight on H
    .updateW        Update W [<on> | off]
    .maxiter        Maximum number of iterations (default 100)
    .conv_criteria  Function exits when cost/delta_cost exceeds this
    .plotfcn        Function handle for plot function
    .plotiter       Plot only every i'th iteration
    .accel          Wild driver accelleration parameter (default 1)
    .displaylevel   Level of display: [off | final | <iter>]
Output
  W                 M x d
  H                 d x N
Example I, Standard NMF:
  d = 4;                                % Four components
  [W,H] = isp_snmf(V,d);
Example I, Sparse NMF:
  d = 2;                                % Two components
  opts.costfcn = 'kl';                  % Kullback Leibler cost function
  opts.lambda = 0.1;                    % Sparsity
  [W,H] = isp_snmf(V,d,opts);