Contents

function comparison_matrix_completion_3
% Comparison of different linearized Bregman approaches on recovering a low-rank matrix from its random samples

Generate problem data

rand('seed', 0); randn('seed', 0);

% --- matrix size, rank, and sample ratio ---
m = 60; n = 100; % matrix size
r = 5; % matrix rank
sr = 0.4; % sample ratio
p = round(m*n*sr); % samples

% --- random low-rank matrix ---
X_ref = rand(m,r)*randn(r,n);

% --- downsample operator and its adjoint ---
S = randsample(m*n,p); % sample positions
A = @(X) reshape(X(S),[],1);  % downsample X at S, takes a matrix and returns a vector
AT = @(y) downsample_adjoint(y,S,[m n]);  % the adjoint of A, takes a vector and returns a matrix

% --- matrix sampling ---
b = A(X_ref);

set parameters

alpha = 10*normest(X_ref); % don't need to be exact, roughly 1 - 10 times norm(X_ref,2) is fine
opts.tol = 1e-8;
opts.maxit = 2000;
opts.X_ref = X_ref;

LBreg: fixed stepsize ---

opts.stepsize = 2/alpha; % roughly 2/alpha/norm(A)^2, where in this example norm(A)^2 = norm(A*A') = norm(I)
t0 = tic;
[X,out] = lbreg_mtx_fixedstep(A,AT,b,alpha,opts);
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out.iter, time);
fprintf('solution relative error = %4.2e\n', norm(X - X_ref,'fro')/norm(X_ref,'fro'));
opts = rmfield(opts, 'stepsize');
iter = 2000, time = 3.31e+00, solution relative error = 2.44e-03

LBreg: Barzilai-Borwein and non-montone line search ---

opts.stepsize = 2/alpha; % roughly 2/alpha/norm(A)^2, where in this example norm(A)^2 = norm(A*A') = norm(I)
t0 = tic;
[X1,out1] = lbreg_mtx_bbls(A,AT,b,alpha,opts);
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out1.iter, time);
fprintf('solution relative error = %4.2e\n', norm(X1 - X_ref,'fro')/norm(X_ref,'fro'));
opts = rmfield(opts, 'stepsize');
iter = 1219, time = 2.90e+00, solution relative error = 3.35e-08

LBreg: accelerated, no restart ---

t0 = tic;
opts.reset = 0; % acceleration with no reset
[X2,out2] = lbreg_mtx_accel_w_reset(A,AT,b,alpha,opts);
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out2.iter, time);
fprintf('solution relative error = %4.2e\n', norm(X2 - X_ref,'fro')/norm(X_ref,'fro'));
opts = rmfield(opts, 'reset');
iter = 2000, time = 3.32e+00, solution relative error = 8.88e-06

LBreg: skip: monotonic gradient scheme ---

opts.reset = 6;
t0 = tic;
[X3,out3] = lbreg_mtx_accel_w_reset(A,AT,b,alpha,opts); % no restart
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out3.iter, time);
fprintf('solution relative error = %4.2e\n', norm(X3 - X_ref,'fro')/norm(X_ref,'fro'));

clear opts
iter = 530, time = 8.78e-01, solution relative error = 1.01e-08

Reporting

figure;
plot(1:out.iter,  out.hist_obj,  'k-', ...
     1:out1.iter, out1.hist_obj, 'b--', ...
     1:out2.iter, out2.hist_obj, 'g--', ...
     1:out3.iter, out3.hist_obj, 'r--', ...
     'LineWidth', 2);
legend('fixed stepsize', 'BB+line search', 'Nesterov accel', 'Nesterov+skip','Location','SouthEast');
title('Dual objective')
xlabel('iteration'); ylabel('dual objective');

figure;
semilogy(1:out.iter,  out.hist_err/norm(X_ref,'fro'),  'k-', ...
         1:out1.iter, out1.hist_err/norm(X_ref,'fro'), 'b--', ...
         1:out2.iter, out2.hist_err/norm(X_ref,'fro'), 'g--', ...
         1:out3.iter, out3.hist_err/norm(X_ref,'fro'), 'r--', ...
         'LineWidth', 2);
legend('fixed stepsize', 'BB+line search', 'Nesterov accel', 'Nesterov+skip','Location','NorthEast');
title('Primal solution relative error')
xlabel('iteration'); ylabel('||X - X_{ref}||_F/||X_{ref}||_F');
end

function res = downsample_adjoint(y,S,sizeX)
res = zeros(sizeX);
res(S) = y;
end

Download this m-file