Matrix Completion.m
From Wikimization
(Difference between revisions)
(adding link to Software Package) |
(→Matlab demonstration of [http://www.temasek-lab.nus.edu.sg/~tslcaij Cai], [http://www.acm.caltech.edu/~emmanuel Candès], & [http://www.math.nus.edu.sg/~matzuows Shen]) |
||
Line 9: | Line 9: | ||
% Created: October 2008 | % Created: October 2008 | ||
- | %% Set path and global variables | ||
- | global SRB | ||
- | SRB = true; | ||
- | |||
%% Setup a matrix | %% Setup a matrix | ||
randn('state',2008); | randn('state',2008); |
Revision as of 03:14, 24 March 2009
Contents |
Matlab demonstration of Cai, Candès, & Shen
A Singular Value Thresholding Algorithm for Matrix Completion, 2008
The latest code can be found at [1]
% Written by: Emmanuel Candes % Email: emmanuel@acm.caltech.edu % Created: October 2008 %% 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); A = U*diag(sigma)*V'; x = A(Omega); if norm(x-b)/normb < tol break end y = y + delta*(b-x); updateSparse(Y,y,indx); end fprintf('\n'); numiter = k;
subroutines
lansvd()
function [U,S,V,bnd,j] = lansvd(varargin) %LANSVD Compute a few singular values and singular vectors. % LANSVD computes singular triplets (u,v,sigma) such that % A*u = sigma*v and A'*v = sigma*u. Only a few singular values % and singular vectors are computed using the Lanczos % bidiagonalization algorithm with partial reorthogonalization (BPRO). % % S = LANSVD(A) % S = LANSVD('Afun','Atransfun',M,N) % % The first input argument is either a matrix or a % string containing the name of an M-file which applies a linear % operator to the columns of a given matrix. In the latter case, % the second input must be the name of an M-file which applies the % transpose of the same operator to the columns of a given matrix, % and the third and fourth arguments must be M and N, the dimensions % of the problem. % % [U,S,V] = LANSVD(A,K,'L',...) computes the K largest singular values. % % [U,S,V] = LANSVD(A,K,'S',...) computes the K smallest singular values. % % The full calling sequence is % % [U,S,V] = LANSVD(A,K,SIGMA,OPTIONS) % [U,S,V] = LANSVD('Afun','Atransfun',M,N,K,SIGMA,OPTIONS) % % where K is the number of singular values desired and % SIGMA is 'L' or 'S'. % % The OPTIONS structure specifies certain parameters in the algorithm. % Field name Parameter Default % % OPTIONS.tol Convergence tolerance 16*eps % OPTIONS.lanmax Dimension of the Lanczos basis. % OPTIONS.p0 Starting vector for the Lanczos rand(n,1)-0.5 % iteration. % OPTIONS.delta Level of orthogonality among the sqrt(eps/K) % Lanczos vectors. % OPTIONS.eta Level of orthogonality after 10*eps^(3/4) % reorthogonalization. % OPTIONS.cgs reorthogonalization method used 0 % '0' : iterated modified Gram-Schmidt % '1' : iterated classical Gram-Schmidt % OPTIONS.elr If equal to 1 then extended local 1 % reorthogonalization is enforced. % % See also LANBPRO, SVDS, SVD % References: % R.M. Larsen, Ph.D. Thesis, Aarhus University, 1998. % % B. N. Parlett, ``The Symmetric Eigenvalue Problem'', % Prentice-Hall, Englewood Cliffs, NJ, 1980. % % H. D. Simon, ``The Lanczos algorithm with partial reorthogonalization'', % Math. Comp. 42 (1984), no. 165, 115--142. % Rasmus Munk Larsen, DAIMI, 1998 %%%%%%%%%%%%%%%%%%%%% Parse and check input arguments. %%%%%%%%%%%%%%%%%%%%%% if nargin<1 | length(varargin)<1 error('Not enough input arguments.'); end A = varargin{1}; if ~isstr(A) if ~isreal(A) error('A must be real') end [m n] = size(A); if length(varargin) < 2, k=min(min(m,n),6); else k=varargin{2}; end if length(varargin) < 3, sigma = 'L'; else sigma=varargin{3}; end if length(varargin) < 4, options = []; else options=varargin{4}; end else if length(varargin)<4 error('Not enough input arguments.'); end Atrans = varargin{2}; if ~isstr(Atrans) error('Atransfunc must be the name of a function') end m = varargin{3}; n = varargin{4}; if length(varargin) < 5, k=min(min(m,n),6); else k=varargin{5}; end if length(varargin) < 6, sigma = 'L'; else sigma=varargin{6}; end if length(varargin) < 7, options = []; else options=varargin{7}; end end if ~isnumeric(n) | real(abs(fix(n))) ~= n | ~isnumeric(m) | ... real(abs(fix(m))) ~= m | ~isnumeric(k) | real(abs(fix(k))) ~= k error('M, N and K must be positive integers.') end % Quick return for min(m,n) equal to 0 or 1 or for zero A. if min(n,m) < 1 | k<1 if nargout<3 U = zeros(k,1); else U = eye(m,k); S = zeros(k,k); V = eye(n,k); bnd = zeros(k,1); end return elseif min(n,m) == 1 & k>0 if isstr(A) % Extract the single column or row of A if n==1 A = feval(A,1); else A = feval(Atrans,1)'; end end if nargout==1 U = norm(A); else [U,S,V] = svd(full(A)); bnd = 0; end return end % A is the matrix of all zeros (not detectable if A is defined by an m-file) if isnumeric(A) if nnz(A)==0 if nargout<3 U = zeros(k,1); else U = eye(m,k); S = zeros(k,k); V = eye(n,k); bnd = zeros(k,1); end return end end lanmax = min(m,n); tol = 16*eps; p = rand(m,1)-0.5; % Parse options struct if isstruct(options) c = fieldnames(options); for i=1:length(c) if any(strcmp(c(i),'p0')), p = getfield(options,'p0'); p=p(:); end if any(strcmp(c(i),'tol')), tol = getfield(options,'tol'); end if any(strcmp(c(i),'lanmax')), lanmax = getfield(options,'lanmax'); end end end % Protect against absurd options. tol = max(tol,eps); lanmax = min(lanmax,min(m,n)); if size(p,1)~=m error('p0 must be a vector of length m') end lanmax = min(lanmax,min(m,n)); if k>lanmax error('K must satisfy K <= LANMAX <= MIN(M,N).'); end %%%%%%%%%%%%%%%%%%%%% Here begins the computation %%%%%%%%%%%%%%%%%%%%%% if strcmp(sigma,'S') if isstr(A) error('Shift-and-invert works only when the matrix A is given explicitly.'); else % Prepare for shift-and-invert Lanczos. if issparse(A) pmmd = colmmd(A); A.A = A(:,pmmd); else A.A = A; end if m>=n if issparse(A.A) A.R = qr(A.A,0); A.Rt = A.R'; p = A.Rt\(A.A'*p); % project starting vector on span(Q1) else [A.Q,A.R] = qr(A.A,0); A.Rt = A.R'; p = A.Q'*p; % project starting vector on span(Q1) end else error('Sorry, shift-and-invert for m<n not implemented yet!') A.R = qr(A.A',0); A.Rt = A.R'; end condR = condest(A.R); if condR > 1/eps error(['A is rank deficient or too ill-conditioned to do shift-and-' ... ' invert.']) end end end ksave = k; neig = 0; nrestart=-1; j = min(k+max(8,k)+1,lanmax); U = []; V = []; B = []; anorm = []; work = zeros(2,2); while neig < k %%%%%%%%%%%%%%%%%%%%% Compute Lanczos bidiagonalization %%%%%%%%%%%%%%%%% if ~isstr(A) [U,B,V,p,ierr,w] = lanbpro(A,j,p,options,U,B,V,anorm); else [U,B,V,p,ierr,w] = lanbpro(A,Atrans,m,n,j,p,options,U,B,V,anorm); end work= work + w; if ierr<0 % Invariant subspace of dimension -ierr found. j = -ierr; end %%%%%%%%%%%%%%%%%% Compute singular values and error bounds %%%%%%%%%%%%%%%% % Analyze B resnrm = norm(p); % We might as well use the extra info. in p. % S = svd(full([B;[zeros(1,j-1),resnrm]]),0); % [P,S,Q] = svd(full([B;[zeros(1,j-1),resnrm]]),0); % S = diag(S); % bot = min(abs([P(end,1:j);Q(end,1:j)]))'; [S,bot] = bdsqr(diag(B),[diag(B,-1); resnrm]); % Use Largest Ritz value to estimate ||A||_2. This might save some % reorth. in case of restart. anorm=S(1); % Set simple error bounds bnd = resnrm*abs(bot); % Examine gap structure and refine error bounds bnd = refinebounds(S.^2,bnd,n*eps*anorm); %%%%%%%%%%%%%%%%%%% Check convergence criterion %%%%%%%%%%%%%%%%%%%% i=1; neig = 0; while i<=min(j,k) if (bnd(i) <= tol*abs(S(i))) neig = neig + 1; i = i+1; else i = min(j,k)+1; end end %%%%%%%%%% Check whether to stop or to extend the Krylov basis? %%%%%%%%%% if ierr<0 % Invariant subspace found if j<k warning(['Invariant subspace of dimension ',num2str(j-1),' found.']) end j = j-1; break; end if j>=lanmax % Maximal dimension of Krylov subspace reached. Bail out if j>=min(m,n) neig = ksave; break; end if neig<ksave warning(['Maximum dimension of Krylov subspace exceeded prior',... ' to convergence.']); end break; end % Increase dimension of Krylov subspace if neig>0 % increase j by approx. half the average number of steps pr. converged % singular value (j/neig) times the number of remaining ones (k-neig). j = j + min(100,max(2,0.5*(k-neig)*j/(neig+1))); else % As long a very few singular values have converged, increase j rapidly. % j = j + ceil(min(100,max(8,2^nrestart*k))); j = max(1.5*j,j+10); end j = ceil(min(j+1,lanmax)); nrestart = nrestart + 1; end %%%%%%%%%%%%%%%% Lanczos converged (or failed). Prepare output %%%%%%%%%%%%%%% k = min(ksave,j); if nargout>2 j = size(B,2); % Compute singular vectors [P,S,Q] = svd(full([B;[zeros(1,j-1),resnrm]]),0); S = diag(S); if size(Q,2)~=k Q = Q(:,1:k); P = P(:,1:k); end % Compute and normalize Ritz vectors (overwrites U and V to save memory). if resnrm~=0 U = U*P(1:j,:) + (p/resnrm)*P(j+1,:); else U = U*P(1:j,:); end V = V*Q; for i=1:k nq = norm(V(:,i)); if isfinite(nq) & nq~=0 & nq~=1 V(:,i) = V(:,i)/nq; end nq = norm(U(:,i)); if isfinite(nq) & nq~=0 & nq~=1 U(:,i) = U(:,i)/nq; end end end % Pick out desired part the spectrum S = S(1:k); bnd = bnd(1:k); if strcmp(sigma,'S') [S,p] = sort(-1./S); S = -S; bnd = bnd(p); if nargout>2 if issparse(A.A) U = A.A*(A.R\U(:,p)); V(pmmd,:) = V(:,p); else U = A.Q(:,1:min(m,n))*U(:,p); V = V(:,p); end end end if nargout<3 U = S; S = B; % Undocumented feature - for checking B. else S = diag(S); end
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)';
compute_int()
function int = compute_int(mu,j,delta,eta,LL,strategy,extra) %COMPUTE_INT: Determine which Lanczos vectors to reorthogonalize against. % % int = compute_int(mu,eta,LL,strategy,extra)) % % Strategy 0: Orthogonalize vectors v_{i-r-extra},...,v_{i},...v_{i+s+extra} % with nu>eta, where v_{i} are the vectors with mu>delta. % Strategy 1: Orthogonalize all vectors v_{r-extra},...,v_{s+extra} where % v_{r} is the first and v_{s} the last Lanczos vector with % mu > eta. % Strategy 2: Orthogonalize all vectors with mu > eta. % % Notice: The first LL vectors are excluded since the new Lanczos % vector is already orthogonalized against them in the main iteration. % Rasmus Munk Larsen, DAIMI, 1998. if (delta<eta) error('DELTA should satisfy DELTA >= ETA.') end switch strategy case 0 I0 = find(abs(mu(1:j))>=delta); if length(I0)==0 [mm,I0] = max(abs(mu(1:j))); end int = zeros(j,1); for i = 1:length(I0) for r=I0(i):-1:1 if abs(mu(r))<eta | int(r)==1 break; else int(r) = 1; end end int(max(1,r-extra+1):r) = 1; for s=I0(i)+1:j if abs(mu(s))<eta | int(s)==1 break; else int(s) = 1; end end int(s:min(j,s+extra-1)) = 1; end if LL>0 int(1:LL) = 0; end int = find(int); case 1 int=find(abs(mu(1:j))>eta); int = max(LL+1,min(int)-extra):min(max(int)+extra,j); case 2 int=find(abs(mu(1:j))>=eta); end int = int(:);
lanbpro()
function [U,B_k,V,p,ierr,work] = lanbpro(varargin) %LANBPRO Lanczos bidiagonalization with partial reorthogonalization. % LANBPRO computes the Lanczos bidiagonalization of a real % matrix using the with partial reorthogonalization. % % [U_k,B_k,V_k,R,ierr,work] = LANBPRO(A,K,R0,OPTIONS,U_old,B_old,V_old) % [U_k,B_k,V_k,R,ierr,work] = LANBPRO('Afun','Atransfun',M,N,K,R0, ... % OPTIONS,U_old,B_old,V_old) % % Computes K steps of the Lanczos bidiagonalization algorithm with partial % reorthogonalization (BPRO) with M-by-1 starting vector R0, producing a % lower bidiagonal K-by-K matrix B_k, an N-by-K matrix V_k, an M-by-K % matrix U_k and an M-by-1 vector R such that % A*V_k = U_k*B_k + R % Partial reorthogonalization is used to keep the columns of V_K and U_k % semiorthogonal: % MAX(DIAG((EYE(K) - V_K'*V_K))) <= OPTIONS.delta % and % MAX(DIAG((EYE(K) - U_K'*U_K))) <= OPTIONS.delta. % % B_k = LANBPRO(...) returns the bidiagonal matrix only. % % The first input argument is either a real matrix, or a string % containing the name of an M-file which applies a linear operator % to the columns of a given matrix. In the latter case, the second % input must be the name of an M-file which applies the transpose of % the same linear operator to the columns of a given matrix, % and the third and fourth arguments must be M and N, the dimensions % of then problem. % % The OPTIONS structure is used to control the reorthogonalization: % OPTIONS.delta: Desired level of orthogonality % (default = sqrt(eps/K)). % OPTIONS.eta : Level of orthogonality after reorthogonalization % (default = eps^(3/4)/sqrt(K)). % OPTIONS.cgs : Flag for switching between different reorthogonalization % algorithms: % 0 = iterated modified Gram-Schmidt (default) % 1 = iterated classical Gram-Schmidt % OPTIONS.elr : If OPTIONS.elr = 1 (default) then extended local % reorthogonalization is enforced. % OPTIONS.onesided % : If OPTIONS.onesided = 0 (default) then both the left % (U) and right (V) Lanczos vectors are kept % semiorthogonal. % OPTIONS.onesided = 1 then only the columns of U are % are reorthogonalized. % OPTIONS.onesided = -1 then only the columns of V are % are reorthogonalized. % OPTIONS.waitbar % : The progress of the algorithm is display graphically. % % If both R0, U_old, B_old, and V_old are provided, they must % contain a partial Lanczos bidiagonalization of A on the form % % A V_old = U_old B_old + R0 . % % In this case the factorization is extended to dimension K x K by % continuing the Lanczos bidiagonalization algorithm with R0 as a % starting vector. % % The output array work contains information about the work used in % reorthogonalizing the u- and v-vectors. % work = [ RU PU ] % [ RV PV ] % where % RU = Number of reorthogonalizations of U. % PU = Number of inner products used in reorthogonalizing U. % RV = Number of reorthogonalizations of V. % PV = Number of inner products used in reorthogonalizing V. % References: % R.M. Larsen, Ph.D. Thesis, Aarhus University, 1998. % % G. H. Golub & C. F. Van Loan, "Matrix Computations", % 3. Ed., Johns Hopkins, 1996. Section 9.3.4. % % B. N. Parlett, ``The Symmetric Eigenvalue Problem'', % Prentice-Hall, Englewood Cliffs, NJ, 1980. % % H. D. Simon, ``The Lanczos algorithm with partial reorthogonalization'', % Math. Comp. 42 (1984), no. 165, 115--142. % % Rasmus Munk Larsen, DAIMI, 1998. % Check input arguments. global LANBPRO_TRUTH LANBPRO_TRUTH=0; if LANBPRO_TRUTH==1 global MU NU MUTRUE NUTRUE global MU_AFTER NU_AFTER MUTRUE_AFTER NUTRUE_AFTER end if nargin<1 | length(varargin)<2 error('Not enough input arguments.'); end narg=length(varargin); A = varargin{1}; if isnumeric(A) | isstruct(A) if isnumeric(A) if ~isreal(A) error('A must be real') end [m n] = size(A); elseif isstruct(A) [m n] = size(A.R); end k=varargin{2}; if narg >= 3 & ~isempty(varargin{3}); p = varargin{3}; else p = rand(m,1)-0.5; end if narg < 4, options = []; else options=varargin{4}; end if narg > 4 if narg<7 error('All or none of U_old, B_old and V_old must be provided.') else U = varargin{5}; B_k = varargin{6}; V = varargin{7}; end else U = []; B_k = []; V = []; end if narg > 7, anorm=varargin{8}; else anorm = []; end else if narg<5 error('Not enough input arguments.'); end Atrans = varargin{2}; if ~isstr(Atrans) error('Afunc and Atransfunc must be names of m-files') end m = varargin{3}; n = varargin{4}; if ~isreal(n) | abs(fix(n)) ~= n | ~isreal(m) | abs(fix(m)) ~= m error('M and N must be positive integers.') end k=varargin{5}; if narg < 6, p = rand(m,1)-0.5; else p=varargin{6}; end if narg < 7, options = []; else options=varargin{7}; end if narg > 7 if narg < 10 error('All or none of U_old, B_old and V_old must be provided.') else U = varargin{8}; B_k = varargin{9}; V = varargin{10}; end else U = []; B_k = []; V=[]; end if narg > 10, anorm=varargin{11}; else anorm = []; end end % Quick return for min(m,n) equal to 0 or 1. if min(m,n) == 0 U = []; B_k = []; V = []; p = []; ierr = 0; work = zeros(2,2); return elseif min(m,n) == 1 if isnumeric(A) U = 1; B_k = A; V = 1; p = 0; ierr = 0; work = zeros(2,2); else U = 1; B_k = feval(A,1); V = 1; p = 0; ierr = 0; work = zeros(2,2); end if nargout<3 U = B_k; end return end % Set options. %m2 = 3/2*(sqrt(m)+1); %n2 = 3/2*(sqrt(n)+1); m2 = 3/2; n2 = 3/2; delta = sqrt(eps/k); % Desired level of orthogonality. eta = eps^(3/4)/sqrt(k); % Level of orth. after reorthogonalization. cgs = 0; % Flag for switching between iterated MGS and CGS. elr = 2; % Flag for switching extended local % reorthogonalization on and off. gamma = 1/sqrt(2); % Tolerance for iterated Gram-Schmidt. onesided = 0; t = 0; waitb = 0; % Parse options struct if ~isempty(options) & isstruct(options) c = fieldnames(options); for i=1:length(c) if strmatch(c(i),'delta'), delta = getfield(options,'delta'); end if strmatch(c(i),'eta'), eta = getfield(options,'eta'); end if strmatch(c(i),'cgs'), cgs = getfield(options,'cgs'); end if strmatch(c(i),'elr'), elr = getfield(options,'elr'); end if strmatch(c(i),'gamma'), gamma = getfield(options,'gamma'); end if strmatch(c(i),'onesided'), onesided = getfield(options,'onesided'); end if strmatch(c(i),'waitbar'), waitb=1; end end end if waitb waitbarh = waitbar(0,'Lanczos bidiagonalization in progress...'); end if isempty(anorm) anorm = []; est_anorm=1; else est_anorm=0; end % Conservative statistical estimate on the size of round-off terms. % Notice that {\bf u} == eps/2. FUDGE = 1.01; % Fudge factor for ||A||_2 estimate. npu = 0; npv = 0; ierr = 0; p = p(:); % Prepare for Lanczos iteration. if isempty(U) V = zeros(n,k); U = zeros(m,k); beta = zeros(k+1,1); alpha = zeros(k,1); beta(1) = norm(p); % Initialize MU/NU-recurrences for monitoring loss of orthogonality. nu = zeros(k,1); mu = zeros(k+1,1); mu(1)=1; nu(1)=1; numax = zeros(k,1); mumax = zeros(k,1); force_reorth = 0; nreorthu = 0; nreorthv = 0; j0 = 1; else j = size(U,2); % Size of existing factorization % Allocate space for Lanczos vectors U = [U, zeros(m,k-j)]; V = [V, zeros(n,k-j)]; alpha = zeros(k+1,1); beta = zeros(k+1,1); alpha(1:j) = diag(B_k); if j>1 beta(2:j) = diag(B_k,-1); end beta(j+1) = norm(p); % Reorthogonalize p. if j<k & beta(j+1)*delta < anorm*eps, fro = 1; ierr = j; end int = [1:j]'; [p,beta(j+1),rr] = reorth(U,p,beta(j+1),int,gamma,cgs); npu = rr*j; nreorthu = 1; force_reorth= 1; % Compute Gerscgorin bound on ||B_k||_2 if est_anorm anorm = FUDGE*sqrt(norm(B_k'*B_k,1)); end mu = m2*eps*ones(k+1,1); nu = zeros(k,1); numax = zeros(k,1); mumax = zeros(k,1); force_reorth = 1; nreorthu = 0; nreorthv = 0; j0 = j+1; end if isnumeric(A) At = A'; end if delta==0 fro = 1; % The user has requested full reorthogonalization. else fro = 0; end if LANBPRO_TRUTH==1 MUTRUE = zeros(k,k); NUTRUE = zeros(k-1,k-1); MU = zeros(k,k); NU = zeros(k-1,k-1); MUTRUE_AFTER = zeros(k,k); NUTRUE_AFTER = zeros(k-1,k-1); MU_AFTER = zeros(k,k); NU_AFTER = zeros(k-1,k-1); end % Perform Lanczos bidiagonalization with partial reorthogonalization. for j=j0:k if waitb waitbar(j/k,waitbarh) end if beta(j) ~= 0 U(:,j) = p/beta(j); else U(:,j) = p; end % Replace norm estimate with largest Ritz value. if j==6 B = [[diag(alpha(1:j-1))+diag(beta(2:j-1),-1)]; ... [zeros(1,j-2),beta(j)]]; anorm = FUDGE*norm(B); est_anorm = 0; end %%%%%%%%%% Lanczos step to generate v_j. %%%%%%%%%%%%% if j==1 if isnumeric(A) r = At*U(:,1); elseif isstruct(A) r = A.R\U(:,1); else r = feval(Atrans,U(:,1)); end alpha(1) = norm(r); if est_anorm anorm = FUDGE*alpha(1); end else if isnumeric(A) r = At*U(:,j) - beta(j)*V(:,j-1); elseif isstruct(A) r = A.R\U(:,j) - beta(j)*V(:,j-1); else r = feval(Atrans,U(:,j)) - beta(j)*V(:,j-1); end alpha(j) = norm(r); % Extended local reorthogonalization if alpha(j)<gamma*beta(j) & elr & ~fro normold = alpha(j); stop = 0; while ~stop t = V(:,j-1)'*r; r = r - V(:,j-1)*t; alpha(j) = norm(r); if beta(j) ~= 0 beta(j) = beta(j) + t; end if alpha(j)>=gamma*normold stop = 1; else normold = alpha(j); end end end if est_anorm if j==2 anorm = max(anorm,FUDGE*sqrt(alpha(1)^2+beta(2)^2+alpha(2)*beta(2))); else anorm = max(anorm,FUDGE*sqrt(alpha(j-1)^2+beta(j)^2+alpha(j-1)* ... beta(j-1) + alpha(j)*beta(j))); end end if ~fro & alpha(j) ~= 0 % Update estimates of the level of orthogonality for the % columns 1 through j-1 in V. nu = update_nu(nu,mu,j,alpha,beta,anorm); numax(j) = max(abs(nu(1:j-1))); end if j>1 & LANBPRO_TRUTH NU(1:j-1,j-1) = nu(1:j-1); NUTRUE(1:j-1,j-1) = V(:,1:j-1)'*r/alpha(j); end if elr>0 nu(j-1) = n2*eps; end % IF level of orthogonality is worse than delta THEN % Reorthogonalize v_j against some previous v_i's, 0<=i<j. if onesided~=-1 & ( fro | numax(j) > delta | force_reorth ) & alpha(j)~=0 % Decide which vectors to orthogonalize against: if fro | eta==0 int = [1:j-1]'; elseif force_reorth==0 int = compute_int(nu,j-1,delta,eta,0,0,0); end % Else use int from last reorth. to avoid spillover from mu_{j-1} % to nu_j. % Reorthogonalize v_j [r,alpha(j),rr] = reorth(V,r,alpha(j),int,gamma,cgs); npv = npv + rr*length(int); % number of inner products. nu(int) = n2*eps; % Reset nu for orthogonalized vectors. % If necessary force reorthogonalization of u_{j+1} % to avoid spillover if force_reorth==0 force_reorth = 1; else force_reorth = 0; end nreorthv = nreorthv + 1; end end % Check for convergence or failure to maintain semiorthogonality if alpha(j) < max(n,m)*anorm*eps & j<k, % If alpha is "small" we deflate by setting it % to 0 and attempt to restart with a basis for a new % invariant subspace by replacing r with a random starting vector: %j %disp('restarting, alpha = 0') alpha(j) = 0; bailout = 1; for attempt=1:3 r = rand(m,1)-0.5; if isnumeric(A) r = At*r; elseif isstruct(A) r = A.R\r; else r = feval(Atrans,r); end nrm=sqrt(r'*r); % not necessary to compute the norm accurately here. int = [1:j-1]'; [r,nrmnew,rr] = reorth(V,r,nrm,int,gamma,cgs); npv = npv + rr*length(int(:)); nreorthv = nreorthv + 1; nu(int) = n2*eps; if nrmnew > 0 % A vector numerically orthogonal to span(Q_k(:,1:j)) was found. % Continue iteration. bailout=0; break; end end if bailout j = j-1; ierr = -j; break; else r=r/nrmnew; % Continue with new normalized r as starting vector. force_reorth = 1; if delta>0 fro = 0; % Turn off full reorthogonalization. end end elseif j<k & ~fro & anorm*eps > delta*alpha(j) % fro = 1; ierr = j; end if j>1 & LANBPRO_TRUTH NU_AFTER(1:j-1,j-1) = nu(1:j-1); NUTRUE_AFTER(1:j-1,j-1) = V(:,1:j-1)'*r/alpha(j); end if alpha(j) ~= 0 V(:,j) = r/alpha(j); else V(:,j) = r; end %%%%%%%%%% Lanczos step to generate u_{j+1}. %%%%%%%%%%%%% if waitb waitbar((2*j+1)/(2*k),waitbarh) end if isnumeric(A) p = A*V(:,j) - alpha(j)*U(:,j); elseif isstruct(A) p = A.Rt\V(:,j) - alpha(j)*U(:,j); else p = feval(A,V(:,j)) - alpha(j)*U(:,j); end beta(j+1) = norm(p); % Extended local reorthogonalization if beta(j+1)<gamma*alpha(j) & elr & ~fro normold = beta(j+1); stop = 0; while ~stop t = U(:,j)'*p; p = p - U(:,j)*t; beta(j+1) = norm(p); if alpha(j) ~= 0 alpha(j) = alpha(j) + t; end if beta(j+1) >= gamma*normold stop = 1; else normold = beta(j+1); end end end if est_anorm % We should update estimate of ||A|| before updating mu - especially % important in the first step for problems with large norm since alpha(1) % may be a severe underestimate! if j==1 anorm = max(anorm,FUDGE*pythag(alpha(1),beta(2))); else anorm = max(anorm,FUDGE*sqrt(alpha(j)^2+beta(j+1)^2 + alpha(j)*beta(j))); end end if ~fro & beta(j+1) ~= 0 % Update estimates of the level of orthogonality for the columns of V. mu = update_mu(mu,nu,j,alpha,beta,anorm); mumax(j) = max(abs(mu(1:j))); end if LANBPRO_TRUTH==1 MU(1:j,j) = mu(1:j); MUTRUE(1:j,j) = U(:,1:j)'*p/beta(j+1); end if elr>0 mu(j) = m2*eps; end % IF level of orthogonality is worse than delta THEN % Reorthogonalize u_{j+1} against some previous u_i's, 0<=i<=j. if onesided~=1 & (fro | mumax(j) > delta | force_reorth) & beta(j+1)~=0 % Decide which vectors to orthogonalize against. if fro | eta==0 int = [1:j]'; elseif force_reorth==0 int = compute_int(mu,j,delta,eta,0,0,0); else int = [int; max(int)+1]; end % Else use int from last reorth. to avoid spillover from nu to mu. % if onesided~=0 % fprintf('i = %i, nr = %i, fro = %i\n',j,size(int(:),1),fro) % end % Reorthogonalize u_{j+1} [p,beta(j+1),rr] = reorth(U,p,beta(j+1),int,gamma,cgs); npu = npu + rr*length(int); nreorthu = nreorthu + 1; % Reset mu to epsilon. mu(int) = m2*eps; if force_reorth==0 force_reorth = 1; % Force reorthogonalization of v_{j+1}. else force_reorth = 0; end end % Check for convergence or failure to maintain semiorthogonality if beta(j+1) < max(m,n)*anorm*eps & j<k, % If beta is "small" we deflate by setting it % to 0 and attempt to restart with a basis for a new % invariant subspace by replacing p with a random starting vector: %j %disp('restarting, beta = 0') beta(j+1) = 0; bailout = 1; for attempt=1:3 p = rand(n,1)-0.5; if isnumeric(A) p = A*p; elseif isstruct(A) p = A.Rt\p; else p = feval(A,p); end nrm=sqrt(p'*p); % not necessary to compute the norm accurately here. int = [1:j]'; [p,nrmnew,rr] = reorth(U,p,nrm,int,gamma,cgs); npu = npu + rr*length(int(:)); nreorthu = nreorthu + 1; mu(int) = m2*eps; if nrmnew > 0 % A vector numerically orthogonal to span(Q_k(:,1:j)) was found. % Continue iteration. bailout=0; break; end end if bailout ierr = -j; break; else p=p/nrmnew; % Continue with new normalized p as starting vector. force_reorth = 1; if delta>0 fro = 0; % Turn off full reorthogonalization. end end elseif j<k & ~fro & anorm*eps > delta*beta(j+1) % fro = 1; ierr = j; end if LANBPRO_TRUTH==1 MU_AFTER(1:j,j) = mu(1:j); MUTRUE_AFTER(1:j,j) = U(:,1:j)'*p/beta(j+1); end end if waitb close(waitbarh) end if j<k k = j; end B_k = spdiags([alpha(1:k) [beta(2:k);0]],[0 -1],k,k); if nargout==1 U = B_k; elseif k~=size(U,2) | k~=size(V,2) U = U(:,1:k); V = V(:,1:k); end if nargout>5 work = [[nreorthu,npu];[nreorthv,npv]]; end function mu = update_mu(muold,nu,j,alpha,beta,anorm) % UPDATE_MU: Update the mu-recurrence for the u-vectors. % % mu_new = update_mu(mu,nu,j,alpha,beta,anorm) % Rasmus Munk Larsen, DAIMI, 1998. binv = 1/beta(j+1); mu = muold; eps1 = 100*eps/2; if j==1 T = eps1*(pythag(alpha(1),beta(2)) + pythag(alpha(1),beta(1))); T = T + eps1*anorm; mu(1) = T / beta(2); else mu(1) = alpha(1)*nu(1) - alpha(j)*mu(1); % T = eps1*(pythag(alpha(j),beta(j+1)) + pythag(alpha(1),beta(1))); T = eps1*(sqrt(alpha(j).^2+beta(j+1).^2) + sqrt(alpha(1).^2+beta(1).^2)); T = T + eps1*anorm; mu(1) = (mu(1) + sign(mu(1))*T) / beta(j+1); % Vectorized version of loop: if j>2 k=2:j-1; mu(k) = alpha(k).*nu(k) + beta(k).*nu(k-1) - alpha(j)*mu(k); %T = eps1*(pythag(alpha(j),beta(j+1)) + pythag(alpha(k),beta(k))); T = eps1*(sqrt(alpha(j).^2+beta(j+1).^2) + sqrt(alpha(k).^2+beta(k).^2)); T = T + eps1*anorm; mu(k) = binv*(mu(k) + sign(mu(k)).*T); end % T = eps1*(pythag(alpha(j),beta(j+1)) + pythag(alpha(j),beta(j))); T = eps1*(sqrt(alpha(j).^2+beta(j+1).^2) + sqrt(alpha(j).^2+beta(j).^2)); T = T + eps1*anorm; mu(j) = beta(j)*nu(j-1); mu(j) = (mu(j) + sign(mu(j))*T) / beta(j+1); end mu(j+1) = 1; function nu = update_nu(nuold,mu,j,alpha,beta,anorm) % UPDATE_MU: Update the nu-recurrence for the v-vectors. % % nu_new = update_nu(nu,mu,j,alpha,beta,anorm) % Rasmus Munk Larsen, DAIMI, 1998. nu = nuold; ainv = 1/alpha(j); eps1 = 100*eps/2; if j>1 k = 1:(j-1); % T = eps1*(pythag(alpha(k),beta(k+1)) + pythag(alpha(j),beta(j))); T = eps1*(sqrt(alpha(k).^2+beta(k+1).^2) + sqrt(alpha(j).^2+beta(j).^2)); T = T + eps1*anorm; nu(k) = beta(k+1).*mu(k+1) + alpha(k).*mu(k) - beta(j)*nu(k); nu(k) = ainv*(nu(k) + sign(nu(k)).*T); end nu(j) = 1; function x = pythag(y,z) %PYTHAG Computes sqrt( y^2 + z^2 ). % % x = pythag(y,z) % % Returns sqrt(y^2 + z^2) but is careful to scale to avoid overflow. % Christian H. Bischof, Argonne National Laboratory, 03/31/89. [m n] = size(y); if m>1 | n>1 y = y(:); z=z(:); rmax = max(abs([y z]'))'; id=find(rmax==0); if length(id)>0 rmax(id) = 1; x = rmax.*sqrt((y./rmax).^2 + (z./rmax).^2); x(id)=0; else x = rmax.*sqrt((y./rmax).^2 + (z./rmax).^2); end x = reshape(x,m,n); else rmax = max(abs([y;z])); if (rmax==0) x = 0; else x = rmax*sqrt((y/rmax)^2 + (z/rmax)^2); end end
refinebounds()
function [bnd,gap] = refinebounds(D,bnd,tol1) %REFINEBONDS Refines error bounds for Ritz values based on gap-structure % % bnd = refinebounds(lambda,bnd,tol1) % % Treat eigenvalues closer than tol1 as a cluster. % Rasmus Munk Larsen, DAIMI, 1998 j = length(D); if j<=1 return end % Sort eigenvalues to use interlacing theorem correctly [D,PERM] = sort(D); bnd = bnd(PERM); % Massage error bounds for very close Ritz values eps34 = sqrt(eps*sqrt(eps)); [y,mid] = max(bnd); for l=[-1,1] for i=((j+1)-l*(j-1))/2:l:mid-l if abs(D(i+l)-D(i)) < eps34*abs(D(i)) if bnd(i)>tol1 & bnd(i+l)>tol1 bnd(i+l) = pythag(bnd(i),bnd(i+l)); bnd(i) = 0; end end end end % Refine error bounds gap = inf*ones(1,j); gap(1:j-1) = min([gap(1:j-1);[D(2:j)-bnd(2:j)-D(1:j-1)]']); gap(2:j) = min([gap(2:j);[D(2:j)-D(1:j-1)-bnd(1:j-1)]']); gap = gap(:); I = find(gap>bnd); bnd(I) = bnd(I).*(bnd(I)./gap(I)); bnd(PERM) = bnd;
reorth()
function [r,normr,nre,s] = reorth(Q,r,normr,index,alpha,method) %REORTH Reorthogonalize a vector using iterated Gram-Schmidt % % [R_NEW,NORMR_NEW,NRE] = reorth(Q,R,NORMR,INDEX,ALPHA,METHOD) % reorthogonalizes R against the subset of columns of Q given by INDEX. % If INDEX==[] then R is reorthogonalized all columns of Q. % If the result R_NEW has a small norm, i.e. if norm(R_NEW) < ALPHA*NORMR, % then a second reorthogonalization is performed. If the norm of R_NEW % is once more decreased by more than a factor of ALPHA then R is % numerically in span(Q(:,INDEX)) and a zero-vector is returned for R_NEW. % % If method==0 then iterated modified Gram-Schmidt is used. % If method==1 then iterated classical Gram-Schmidt is used. % % The default value for ALPHA is 0.5. % NRE is the number of reorthogonalizations performed (1 or 2). % References: % Aake Bjorck, "Numerical Methods for Least Squares Problems", % SIAM, Philadelphia, 1996, pp. 68-69. % % J.~W. Daniel, W.~B. Gragg, L. Kaufman and G.~W. Stewart, % ``Reorthogonalization and Stable Algorithms Updating the % Gram-Schmidt QR Factorization'', Math. Comp., 30 (1976), no. % 136, pp. 772-795. % % B. N. Parlett, ``The Symmetric Eigenvalue Problem'', % Prentice-Hall, Englewood Cliffs, NJ, 1980. pp. 105-109 % Rasmus Munk Larsen, DAIMI, 1998. % Check input arguments. %warning('PROPACK:NotUsingMex','Using slow matlab code for reorth.') if nargin<2 error('Not enough input arguments.') end [n k1] = size(Q); if nargin<3 | isempty(normr) % normr = norm(r); normr = sqrt(r'*r); end if nargin<4 | isempty(index) k=k1; index = [1:k]'; simple = 1; else k = length(index); if k==k1 & index(:)==[1:k]' simple = 1; else simple = 0; end end if nargin<5 | isempty(alpha) alpha=0.5; % This choice garanties that % || Q^T*r_new - e_{k+1} ||_2 <= 2*eps*||r_new||_2, % cf. Kahans ``twice is enough'' statement proved in % Parletts book. end if nargin<6 | isempty(method) method = 0; end if k==0 | n==0 return end if nargout>3 s = zeros(k,1); end normr_old = 0; nre = 0; while normr < alpha*normr_old | nre==0 if method==1 if simple t = Q'*r; r = r - Q*t; else t = Q(:,index)'*r; r = r - Q(:,index)*t; end else for i=index, t = Q(:,i)'*r; r = r - Q(:,i)*t; end end if nargout>3 s = s + t; end normr_old = normr; % normr = norm(r); normr = sqrt(r'*r); nre = nre + 1; if nre > 4 % r is in span(Q) to full accuracy => accept r = 0 as the new vector. r = zeros(n,1); normr = 0; return end end
updateSparse.c
A precompiled Intel Windows-32bit Matlab mex file is here.
Otherwise, compile this C program in Matlab via commandmex updateSparse.c
/* * Stephen Becker, 11/10/08 * Updates a sparse vector very quickly * calling format: * updateSparse(Y,b) * which updates the values of Y to be b * * Modified 11/12/08 to allow unsorted omega * (omega is the implicit index: in Matlab, what * we are doing is Y(omega) = b. So, if omega * is unsorted, then b must be re-ordered appropriately * */ #include "mex.h" #ifndef true #define true 1 #endif #ifndef false #define false 0 #endif void printUsage() { mexPrintf("usage:\tupdateSparse(Y,b)\nchanges the sparse Y matrix"); mexPrintf(" to have values b\non its nonzero elements. Be careful:\n\t"); mexPrintf("this assumes b is sorted in the appropriate order!\n"); mexPrintf("If b (i.e. the index omega, where we want to perform Y(omega)=b)\n"); mexPrintf(" is unsorted, then call the command as follows:\n"); mexPrintf("\tupdateSparse(Y,b,omegaIndx)\n"); mexPrintf("where [temp,omegaIndx] = sort(omega)\n"); } void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) { /* Declare variable */ int M, N, i, j, m, n; double *b, *S, *omega; int SORTED = true; /* Check for proper number of input and output arguments */ if ( (nrhs < 2) || (nrhs > 3) ) { printUsage(); mexErrMsgTxt("Needs 2 or 3 input arguments"); } if ( nrhs == 3 ) SORTED = false; if(nlhs > 0){ printUsage(); mexErrMsgTxt("No output arguments!"); } /* Check data type of input argument */ if (!(mxIsSparse(prhs[0])) || !((mxIsDouble(prhs[1]))) ){ printUsage(); mexErrMsgTxt("Input arguments wrong data-type (must be sparse, double)."); } /* Get the size and pointers to input data */ /* Check second input */ N = mxGetN( prhs[1] ); M = mxGetM( prhs[1] ); if ( (M>1) && (N>1) ) { printUsage(); mexErrMsgTxt("Second argument must be a vector"); } N = (N>M) ? N : M; /* Check first input */ M = mxGetNzmax( prhs[0] ); if ( M != N ) { printUsage(); mexErrMsgTxt("Length of second argument must match nnz of first argument"); } /* if 3rd argument provided, check that it agrees with 2nd argument */ if (!SORTED) { m = mxGetM( prhs[2] ); n = mxGetN( prhs[2] ); if ( (m>1) && (n>1) ) { printUsage(); mexErrMsgTxt("Third argument must be a vector"); } n = (n>m) ? n : m; if ( n != N ) { printUsage(); mexErrMsgTxt("Third argument must be same length as second argument"); } omega = mxGetPr( prhs[2] ); } b = mxGetPr( prhs[1] ); S = mxGetPr( prhs[0] ); if (SORTED) { /* And here's the really fancy part: */ for ( i=0 ; i < N ; i++ ) S[i] = b[i]; } else { for ( i=0 ; i < N ; i++ ) { /* this is a little slow, but I should really check * to make sure the index is not out-of-bounds, otherwise * Matlab could crash */ j = (int)omega[i]-1; /* the -1 because Matlab is 1-based */ if (j >= N){ printUsage(); mexErrMsgTxt("Third argument must have values < length of 2nd argument"); } /* S[ j ] = b[i]; */ /* this is incorrect */ S[ i ] = b[j]; /* this is the correct form */ } } }