Matrix Completion.m

From Wikimization

Revision as of 03:03, 16 February 2009

Matlab demonstration of Cai, Candès, & Shen

A Singular Value Thresholding Algorithm for Matrix Completion, 2008

test_SVT.m

% Written by: Emmanuel Candes
% Email: emmanuel@acm.caltech.edu
% Created: October 2008

%% Set path and global variables
    global SRB 
    SRB = true;
  
%% Setup a matrix
randn('state',2008);
rand('state',2008);

n = 1000; r = 10;
M = randn(n,r)*randn(r,n);

df = r*(2*n-r);
oversampling = 5;  m = 5*df; 

Omega = randsample(n^2,m); 
data = M(Omega);

%% Set parameters and solve

p  = m/n^2;  delta = 1.2/p;      
maxiter = 500;
tol = 1e-4;

%% Approximate minimum nuclear norm solution by SVT algorithm

    tic
    [U,S,V,numiter] = SVT(n,Omega,data,delta,maxiter,tol);
    toc 
    
%% Show results

 X = U*S*V';

disp(sprintf('The relative error on Omega is: %d ', norm(data-X(Omega))/norm(data)))
disp(sprintf('The relative recovery error is: %d ', norm(M-X,'fro')/norm(M,'fro')))
disp(sprintf('The relative recovery in the spectral norm is: %d ', norm(M-X)/norm(M)))

SVT()

function [U,Sigma,V,numiter]  = SVT(n,Omega,b,delta,maxiter,tol)
% 
% Finds the minimum of tau ||X||_* + .5 || X ||_F^2 
%
% subject to P_Omega(X) = P_Omega(M)
%
% using linear Bregman iterations
%
% Usage:  [U,S,V,numiter]  = SVT(n,Omega,b,delta,maxiter,opts)
%
% Inputs:
%
% n - size of the matrix X assumed n by n
%
% Omega - set of observed entries
%
% b - data vector of the form M(Omega)
%
% delta - step size 
%
% maxiter - maximum number of iterations
%
% Outputs: matrix X stored in SVD format X = U*diag(S)*V'
% 
% U - nxr left singular vectors 
% 
% S - rx1 singular values
%
% V - nxr right singular vectors 
%
% numiter - number of iterations to achieve convergence

% Description: 
% Reference:
%
%    Cai, Candes and Shen
%    A singular value thresholding algorithm for matrix completion
%    Submitted for publication, October 2008.
%
% Written by: Emmanuel Candes
% Email: emmanuel@acm.caltech.edu
% Created: October 2008

m = length(Omega); [temp,indx] = sort(Omega); 
tau = 5*n; incre = 5; 

[i, j] = ind2sub([n,n], Omega);
ProjM = sparse(i,j,b,n,n,m);

normProjM = normest(ProjM);
k0 = ceil(tau/(delta*normProjM));

normb = norm(b);

y = k0*delta*b;
Y = sparse(i,j,y,n,n,m);
r = 0;

fprintf('\nIteration:   ');
for k = 1:maxiter,
    fprintf('\b\b\b%3d',k); 
    s = r + 1;
    
    OK = 0;
    while ~OK 
        [U,Sigma,V] = lansvd(Y,s,'L');
        OK = (Sigma(s,s) <= tau);
        s = s + incre;
    end
   
    sigma = diag(Sigma); r = sum(sigma > tau);
    U = U(:,1:r); V = V(:,1:r); sigma = sigma(1:r) - tau; Sigma = diag(sigma);
    
    x = XonOmega(U*diag(sigma),V,Omega);
    
    if (norm(x-b)/normb < tol)
        break
    end
    
   y = y + delta*(b-x);
   updateSparse(Y,y,indx);   
end

fprintf('\n');
numiter = k;

bdsqr()

function [sigma,bnd] = bdsqr(alpha,beta)

% BDSQR: Compute the singular values and bottom element of
%        the left singular vectors of a (k+1) x k lower bidiagonal 
%        matrix with diagonal alpha(1:k) and lower bidiagonal beta(1:k),
%        where length(alpha) = length(beta) = k.
%
% [sigma,bnd] = bdsqr(alpha,beta)
%
% Input parameters:
%   alpha(1:k)   : Diagonal elements.
%   beta(1:k)    : Sub-diagonal elements.
% Output parameters:
%   sigma(1:k)  : Computed eigenvalues.
%   bnd(1:k)    : Bottom elements in left singular vectors.

% Below is a very slow replacement for the BDSQR MEX-file.

%warning('PROPACK:NotUsingMex','Using slow matlab code for bdsqr.')
k = length(alpha);
if min(size(alpha)') ~= 1  | min(size(beta)') ~= 1
  error('alpha and beta must be vectors')
elseif length(beta) ~= k
  error('alpha and beta must have the same lenght')
end    
B = spdiags([alpha(:),beta(:)],[0,-1],k+1,k);
[U,S,V] = svd(full(B),0);
sigma = diag(S);
bnd = U(end,1:k)';