Matrix Completion.m

From Wikimization

(Difference between revisions)
Jump to: navigation, search
Line 44: Line 44:
disp(sprintf('The relative recovery error is: %d ', norm(M-X,'fro')/norm(M,'fro')))
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)))
disp(sprintf('The relative recovery in the spectral norm is: %d ', norm(M-X)/norm(M)))
 +
</pre>
 +
 +
 +
=== SVT() ===
 +
<pre>
 +
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;
 +
</pre>
 +
 +
 +
=== bdsqr() ===
 +
<pre>
 +
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)';
</pre>
</pre>

Revision as of 02:03, 16 February 2009

Contents

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)';
Personal tools