Matrix Completion.m
From Wikimization
(Difference between revisions)
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 03: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)';