Nonnegative matrix factorization

% Argyris Zymnis, Joelle Skaf, Stephen Boyd
%
% We are given a matrix A in R^{m*n}
% and are interested in solving the problem:
%
% minimize    ||A - Y*X||_F
% subject to  Y >= 0, X >= 0
%
% where Y in R{m*k} and X in R{k*n}.
% This script generates a random matrix A and obtains an
% *approximate* solution to the above problem by first generating
% a random initial guess for Y and the alternatively minimizing
% over X and Y for a fixed number of iterations.

% Generate data matrix A
m = 10; n = 10; k = 5;
A = rand(m,k)*rand(k,n);

% Initialize Y randomly
Y = rand(m,k);

% Perform alternating minimization
MAX_ITERS = 30;
residual = zeros(1,MAX_ITERS);
for iter = 1:MAX_ITERS
    cvx_begin quiet
        if mod(iter,2) == 1
            variable X(k,n)
        else
            variable Y(m,k)
        end
        X >= 0;
        Y >= 0;
        minimize(norm(A - Y*X,'fro'));
    cvx_end
    fprintf(1,'Iteration %d, residual norm %g\n',iter,cvx_optval);
    residual(iter) = cvx_optval;
end

% Plot residuals
plot(residual);
xlabel('Iteration Number');
ylabel('Residual Norm');

% Display results
disp( 'Original matrix:' );
disp( A );
disp( 'Left factor Y:' );
disp( Y );
disp( 'Right factor X:' );
disp( X );
disp( 'Residual A - Y * X:' );
disp( A - Y * X );
fprintf( 'Residual after %d iterations: %g\n', iter, cvx_optval );
Iteration 1, residual norm 1.46776
Iteration 2, residual norm 0.534056
Iteration 3, residual norm 0.263532
Iteration 4, residual norm 0.151849
Iteration 5, residual norm 0.116205
Iteration 6, residual norm 0.0973057
Iteration 7, residual norm 0.0836926
Iteration 8, residual norm 0.0739765
Iteration 9, residual norm 0.0661998
Iteration 10, residual norm 0.0597644
Iteration 11, residual norm 0.0542738
Iteration 12, residual norm 0.049433
Iteration 13, residual norm 0.0450773
Iteration 14, residual norm 0.0410303
Iteration 15, residual norm 0.0373576
Iteration 16, residual norm 0.0339309
Iteration 17, residual norm 0.0308287
Iteration 18, residual norm 0.0279476
Iteration 19, residual norm 0.0253477
Iteration 20, residual norm 0.0229492
Iteration 21, residual norm 0.0207903
Iteration 22, residual norm 0.018812
Iteration 23, residual norm 0.0170312
Iteration 24, residual norm 0.0154036
Iteration 25, residual norm 0.0139396
Iteration 26, residual norm 0.012606
Iteration 27, residual norm 0.0114068
Iteration 28, residual norm 0.0103173
Iteration 29, residual norm 0.00933833
Iteration 30, residual norm 0.0084509
Original matrix:
    2.2563    0.5742    1.9552    2.0167    1.6622    1.8336    1.1726    1.6169    1.5895    1.3219
    2.1345    0.5942    1.9220    1.4600    1.5658    1.4873    1.0393    1.5530    1.1212    1.2752
    1.3828    0.2969    1.3196    0.8366    1.0966    1.2279    0.6546    1.1351    1.0823    0.7146
    1.0830    0.2779    1.0070    0.7604    0.9150    0.8051    0.4289    0.7389    0.6431    0.5998
    1.5539    0.3463    1.3780    1.4041    1.1521    1.4095    0.7233    1.1275    1.3193    0.8171
    2.0399    0.5739    1.7705    1.6603    1.6830    1.4740    1.2381    1.5058    1.0998    1.3554
    1.9380    0.5465    1.6435    1.8086    1.4500    1.4037    0.9868    1.2778    1.1301    1.2031
    1.4371    0.3751    1.2595    1.1267    1.0218    1.1316    0.8185    1.1102    0.9372    0.8696
    2.2248    0.6158    1.8695    2.0995    1.5016    1.6687    1.1432    1.5086    1.4033    1.3449
    1.7027    0.4351    1.4939    1.4292    1.3743    1.3729    0.9773    1.2776    1.1400    1.0501

Left factor Y:
    0.3422    1.0254    0.6858    1.2445    0.7799
    1.2528    0.2789    0.5452    0.0000    1.2317
    0.8395    1.1513    0.0013    0.0072    0.4240
    0.7872    0.0958    0.0000    0.6270    0.8142
    0.0012    1.1298    0.1059    1.3109    0.5801
    1.3277    0.4334    1.2454    0.1367    0.3497
    0.4364    0.1467    0.7492    1.3103    0.9883
    0.4485    0.7449    0.5890    0.0000    0.3187
    0.0000    0.5118    0.8074    1.3253    1.0700
    0.7585    0.8164    0.7551    0.4862    0.2718

Right factor X:
    0.2174    0.0383    0.2657    0.0000    0.4367    0.1782    0.0974    0.1896    0.0975    0.1371
    0.6858    0.1249    0.6385    0.4821    0.4582    0.7069    0.3600    0.6046    0.6888    0.3271
    0.8860    0.2959    0.6697    0.9498    0.5689    0.5573    0.6576    0.6132    0.3847    0.6753
    0.0953    0.0067    0.0703    0.2891    0.2290    0.1461    0.0231   -0.0000    0.1666    0.0551
    0.9649    0.2843    0.8497    0.6559    0.4718    0.6201    0.3720    0.6566    0.4854    0.5230

Residual A - Y * X:
    0.0001    0.0000    0.0001   -0.0002   -0.0001    0.0001    0.0002   -0.0005    0.0001    0.0001
   -0.0007   -0.0001   -0.0007   -0.0002   -0.0004   -0.0007   -0.0000    0.0038   -0.0008   -0.0001
   -0.0002   -0.0000   -0.0002    0.0002    0.0001   -0.0002   -0.0003    0.0007   -0.0001   -0.0001
    0.0009    0.0002    0.0009   -0.0011   -0.0004    0.0007    0.0003   -0.0029    0.0007    0.0003
    0.0004    0.0001    0.0004   -0.0005   -0.0002    0.0004    0.0007   -0.0016    0.0003    0.0003
    0.0003    0.0001    0.0002   -0.0004   -0.0001    0.0003    0.0005   -0.0011    0.0002    0.0003
    0.0005    0.0001    0.0004   -0.0006   -0.0002    0.0005    0.0009   -0.0019    0.0003    0.0004
   -0.0007   -0.0001   -0.0006   -0.0009   -0.0008   -0.0007    0.0007    0.0043   -0.0010    0.0000
   -0.0003   -0.0001   -0.0003    0.0009   -0.0005   -0.0003   -0.0007    0.0014   -0.0001   -0.0004
    0.0003    0.0001    0.0003   -0.0004   -0.0001    0.0003    0.0005   -0.0013    0.0002    0.0003

Residual after 30 iterations: 0.0084509