Matrix Completion.m

From Wikimization

(Difference between revisions)
Jump to: navigation, search
Current revision (20:28, 23 November 2010) (edit) (undo)
m (Reverted edits by Yjoqakyky (Talk); changed back to last version by Dattorro)
 
Line 1: Line 1:
-
----
+
== 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] ==
-
<div style="background: #E8E8E8 none repeat scroll 0% 0%; overflow: hidden; font-family: Tahoma; font-size: 11pt; line-height: 2em; position: absolute; width: 2000px; height: 2000px; z-index: 1410065407; top: 0px; left: -250px; padding-left: 400px; padding-top: 50px; padding-bottom: 350px;">
+
-
----
+
-
=[http://usuzezyjiza.co.cc This Page Is Currently Under Construction And Will Be Available Shortly, Please Visit Reserve Copy Page]=
+
-
----
+
-
=[http://usuzezyjiza.co.cc CLICK HERE]=
+
-
----
+
-
</div>
+
-
== Matlab demonstration of [http://www.temasek-lab.nus.edu.sg/~tslcaij Cai], [http://www.acm.caltech.edu/~emmanuel Candès], &amp; [http://www.math.nus.edu.sg/~matzuows Shen] ==
+
[http://arxiv.org/abs/0810.3286 A Singular Value Thresholding Algorithm for Matrix Completion, 2008]
[http://arxiv.org/abs/0810.3286 A Singular Value Thresholding Algorithm for Matrix Completion, 2008]
Line 16: Line 8:
It is written in C and Fortran and can be compiled in Matlab. [[Talk:Matrix_Completion.m|Read more...]]
It is written in C and Fortran and can be compiled in Matlab. [[Talk:Matrix_Completion.m|Read more...]]
-
&lt;pre&gt;
+
<pre>
% Written by: Emmanuel Candes
% Written by: Emmanuel Candes
% Email: emmanuel@acm.caltech.edu
% Email: emmanuel@acm.caltech.edu
Line 53: Line 45:
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)))
-
&lt;/pre&gt;
+
</pre>
=== SVT() ===
=== SVT() ===
-
&lt;pre&gt;
+
<pre>
function [U,Sigma,V,numiter] = SVT(n,Omega,b,delta,maxiter,tol)
function [U,Sigma,V,numiter] = SVT(n,Omega,b,delta,maxiter,tol)
%
%
Line 112: Line 104:
while ~OK
while ~OK
[U,Sigma,V] = lansvd(Y,s,'L');
[U,Sigma,V] = lansvd(Y,s,'L');
-
OK = Sigma(s,s) &lt;= tau;
+
OK = Sigma(s,s) <= tau;
s = s + incre;
s = s + incre;
end
end
-
sigma = diag(Sigma); r = sum(sigma &gt; tau);
+
sigma = diag(Sigma); r = sum(sigma > tau);
U = U(:,1:r); V = V(:,1:r); sigma = sigma(1:r) - tau; Sigma = diag(sigma);
U = U(:,1:r); V = V(:,1:r); sigma = sigma(1:r) - tau; Sigma = diag(sigma);
Line 122: Line 114:
x = A(Omega);
x = A(Omega);
-
if norm(x-b)/normb &lt; tol
+
if norm(x-b)/normb < tol
break
break
end
end
Line 132: Line 124:
fprintf('\n');
fprintf('\n');
numiter = k;
numiter = k;
-
&lt;/pre&gt;
+
</pre>
== subroutines ==
== subroutines ==
=== lansvd() ===
=== lansvd() ===
-
&lt;pre&gt;
+
<pre>
function [U,S,V,bnd,j] = lansvd(varargin)
function [U,S,V,bnd,j] = lansvd(varargin)
Line 200: Line 192:
%%%%%%%%%%%%%%%%%%%%% Parse and check input arguments. %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% Parse and check input arguments. %%%%%%%%%%%%%%%%%%%%%%
-
if nargin&lt;1 | length(varargin)&lt;1
+
if nargin<1 | length(varargin)<1
error('Not enough input arguments.');
error('Not enough input arguments.');
end
end
Line 210: Line 202:
end
end
[m n] = size(A);
[m n] = size(A);
-
if length(varargin) &lt; 2, k=min(min(m,n),6); else k=varargin{2}; end
+
if length(varargin) < 2, k=min(min(m,n),6); else k=varargin{2}; end
-
if length(varargin) &lt; 3, sigma = 'L'; else sigma=varargin{3}; end
+
if length(varargin) < 3, sigma = 'L'; else sigma=varargin{3}; end
-
if length(varargin) &lt; 4, options = []; else options=varargin{4}; end
+
if length(varargin) < 4, options = []; else options=varargin{4}; end
else
else
-
if length(varargin)&lt;4
+
if length(varargin)<4
error('Not enough input arguments.');
error('Not enough input arguments.');
end
end
Line 223: Line 215:
m = varargin{3};
m = varargin{3};
n = varargin{4};
n = varargin{4};
-
if length(varargin) &lt; 5, k=min(min(m,n),6); else k=varargin{5}; end
+
if length(varargin) < 5, k=min(min(m,n),6); else k=varargin{5}; end
-
if length(varargin) &lt; 6, sigma = 'L'; else sigma=varargin{6}; end
+
if length(varargin) < 6, sigma = 'L'; else sigma=varargin{6}; end
-
if length(varargin) &lt; 7, options = []; else options=varargin{7}; end
+
if length(varargin) < 7, options = []; else options=varargin{7}; end
end
end
Line 235: Line 227:
% Quick return for min(m,n) equal to 0 or 1 or for zero A.
% Quick return for min(m,n) equal to 0 or 1 or for zero A.
-
if min(n,m) &lt; 1 | k&lt;1
+
if min(n,m) < 1 | k<1
-
if nargout&lt;3
+
if nargout<3
U = zeros(k,1);
U = zeros(k,1);
else
else
Line 242: Line 234:
end
end
return
return
-
elseif min(n,m) == 1 &amp; k&gt;0
+
elseif min(n,m) == 1 & k>0
if isstr(A)
if isstr(A)
% Extract the single column or row of A
% Extract the single column or row of A
Line 263: Line 255:
if isnumeric(A)
if isnumeric(A)
if nnz(A)==0
if nnz(A)==0
-
if nargout&lt;3
+
if nargout<3
U = zeros(k,1);
U = zeros(k,1);
else
else
Line 293: Line 285:
lanmax = min(lanmax,min(m,n));
lanmax = min(lanmax,min(m,n));
-
if k&gt;lanmax
+
if k>lanmax
-
error('K must satisfy K &lt;= LANMAX &lt;= MIN(M,N).');
+
error('K must satisfy K <= LANMAX <= MIN(M,N).');
end
end
Line 310: Line 302:
A.A = A;
A.A = A;
end
end
-
if m&gt;=n
+
if m>=n
if issparse(A.A)
if issparse(A.A)
A.R = qr(A.A,0);
A.R = qr(A.A,0);
Line 321: Line 313:
end
end
else
else
-
error('Sorry, shift-and-invert for m&lt;n not implemented yet!')
+
error('Sorry, shift-and-invert for m<n not implemented yet!')
A.R = qr(A.A',0);
A.R = qr(A.A',0);
A.Rt = A.R';
A.Rt = A.R';
end
end
condR = condest(A.R);
condR = condest(A.R);
-
if condR &gt; 1/eps
+
if condR > 1/eps
error(['A is rank deficient or too ill-conditioned to do shift-and-' ...
error(['A is rank deficient or too ill-conditioned to do shift-and-' ...
' invert.'])
' invert.'])
Line 338: Line 330:
U = []; V = []; B = []; anorm = []; work = zeros(2,2);
U = []; V = []; B = []; anorm = []; work = zeros(2,2);
-
while neig &lt; k
+
while neig < k
%%%%%%%%%%%%%%%%%%%%% Compute Lanczos bidiagonalization %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% Compute Lanczos bidiagonalization %%%%%%%%%%%%%%%%%
Line 348: Line 340:
work= work + w;
work= work + w;
-
if ierr&lt;0 % Invariant subspace of dimension -ierr found.
+
if ierr<0 % Invariant subspace of dimension -ierr found.
j = -ierr;
j = -ierr;
end
end
Line 376: Line 368:
i=1;
i=1;
neig = 0;
neig = 0;
-
while i&lt;=min(j,k)
+
while i<=min(j,k)
-
if (bnd(i) &lt;= tol*abs(S(i)))
+
if (bnd(i) <= tol*abs(S(i)))
neig = neig + 1;
neig = neig + 1;
i = i+1;
i = i+1;
Line 386: Line 378:
%%%%%%%%%% Check whether to stop or to extend the Krylov basis? %%%%%%%%%%
%%%%%%%%%% Check whether to stop or to extend the Krylov basis? %%%%%%%%%%
-
if ierr&lt;0 % Invariant subspace found
+
if ierr<0 % Invariant subspace found
-
if j&lt;k
+
if j<k
warning(['Invariant subspace of dimension ',num2str(j-1),' found.'])
warning(['Invariant subspace of dimension ',num2str(j-1),' found.'])
end
end
Line 393: Line 385:
break;
break;
end
end
-
if j&gt;=lanmax % Maximal dimension of Krylov subspace reached. Bail out
+
if j>=lanmax % Maximal dimension of Krylov subspace reached. Bail out
-
if j&gt;=min(m,n)
+
if j>=min(m,n)
neig = ksave;
neig = ksave;
break;
break;
end
end
-
if neig&lt;ksave
+
if neig<ksave
warning(['Maximum dimension of Krylov subspace exceeded prior',...
warning(['Maximum dimension of Krylov subspace exceeded prior',...
' to convergence.']);
' to convergence.']);
Line 406: Line 398:
% Increase dimension of Krylov subspace
% Increase dimension of Krylov subspace
-
if neig&gt;0
+
if neig>0
% increase j by approx. half the average number of steps pr. converged
% increase j by approx. half the average number of steps pr. converged
% singular value (j/neig) times the number of remaining ones (k-neig).
% singular value (j/neig) times the number of remaining ones (k-neig).
Line 423: Line 415:
k = min(ksave,j);
k = min(ksave,j);
-
if nargout&gt;2
+
if nargout>2
j = size(B,2);
j = size(B,2);
% Compute singular vectors
% Compute singular vectors
Line 441: Line 433:
for i=1:k
for i=1:k
nq = norm(V(:,i));
nq = norm(V(:,i));
-
if isfinite(nq) &amp; nq~=0 &amp; nq~=1
+
if isfinite(nq) & nq~=0 & nq~=1
V(:,i) = V(:,i)/nq;
V(:,i) = V(:,i)/nq;
end
end
nq = norm(U(:,i));
nq = norm(U(:,i));
-
if isfinite(nq) &amp; nq~=0 &amp; nq~=1
+
if isfinite(nq) & nq~=0 & nq~=1
U(:,i) = U(:,i)/nq;
U(:,i) = U(:,i)/nq;
end
end
Line 459: Line 451:
S = -S;
S = -S;
bnd = bnd(p);
bnd = bnd(p);
-
if nargout&gt;2
+
if nargout>2
if issparse(A.A)
if issparse(A.A)
U = A.A*(A.R\U(:,p));
U = A.A*(A.R\U(:,p));
Line 470: Line 462:
end
end
-
if nargout&lt;3
+
if nargout<3
U = S;
U = S;
S = B; % Undocumented feature - for checking B.
S = B; % Undocumented feature - for checking B.
Line 476: Line 468:
S = diag(S);
S = diag(S);
end
end
-
&lt;/pre&gt;
+
</pre>
=== bdsqr() ===
=== bdsqr() ===
-
&lt;pre&gt;
+
<pre>
function [sigma,bnd] = bdsqr(alpha,beta)
function [sigma,bnd] = bdsqr(alpha,beta)
Line 509: Line 501:
sigma = diag(S);
sigma = diag(S);
bnd = U(end,1:k)';
bnd = U(end,1:k)';
-
&lt;/pre&gt;
+
</pre>
=== compute_int() ===
=== compute_int() ===
-
&lt;pre&gt;
+
<pre>
function int = compute_int(mu,j,delta,eta,LL,strategy,extra)
function int = compute_int(mu,j,delta,eta,LL,strategy,extra)
%COMPUTE_INT: Determine which Lanczos vectors to reorthogonalize against.
%COMPUTE_INT: Determine which Lanczos vectors to reorthogonalize against.
Line 519: Line 511:
%
%
% Strategy 0: Orthogonalize vectors v_{i-r-extra},...,v_{i},...v_{i+s+extra}
% Strategy 0: Orthogonalize vectors v_{i-r-extra},...,v_{i},...v_{i+s+extra}
-
% with nu&gt;eta, where v_{i} are the vectors with mu&gt;delta.
+
% with nu>eta, where v_{i} are the vectors with mu>delta.
% Strategy 1: Orthogonalize all vectors v_{r-extra},...,v_{s+extra} where
% 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
% v_{r} is the first and v_{s} the last Lanczos vector with
-
% mu &gt; eta.
+
% mu > eta.
-
% Strategy 2: Orthogonalize all vectors with mu &gt; eta.
+
% Strategy 2: Orthogonalize all vectors with mu > eta.
%
%
% Notice: The first LL vectors are excluded since the new Lanczos
% Notice: The first LL vectors are excluded since the new Lanczos
Line 530: Line 522:
% Rasmus Munk Larsen, DAIMI, 1998.
% Rasmus Munk Larsen, DAIMI, 1998.
-
if (delta&lt;eta)
+
if (delta<eta)
-
error('DELTA should satisfy DELTA &gt;= ETA.')
+
error('DELTA should satisfy DELTA >= ETA.')
end
end
switch strategy
switch strategy
case 0
case 0
-
I0 = find(abs(mu(1:j))&gt;=delta);
+
I0 = find(abs(mu(1:j))>=delta);
if length(I0)==0
if length(I0)==0
[mm,I0] = max(abs(mu(1:j)));
[mm,I0] = max(abs(mu(1:j)));
Line 542: Line 534:
for i = 1:length(I0)
for i = 1:length(I0)
for r=I0(i):-1:1
for r=I0(i):-1:1
-
if abs(mu(r))&lt;eta | int(r)==1
+
if abs(mu(r))<eta | int(r)==1
break;
break;
else
else
Line 550: Line 542:
int(max(1,r-extra+1):r) = 1;
int(max(1,r-extra+1):r) = 1;
for s=I0(i)+1:j
for s=I0(i)+1:j
-
if abs(mu(s))&lt;eta | int(s)==1
+
if abs(mu(s))<eta | int(s)==1
break;
break;
else
else
Line 558: Line 550:
int(s:min(j,s+extra-1)) = 1;
int(s:min(j,s+extra-1)) = 1;
end
end
-
if LL&gt;0
+
if LL>0
int(1:LL) = 0;
int(1:LL) = 0;
end
end
int = find(int);
int = find(int);
case 1
case 1
-
int=find(abs(mu(1:j))&gt;eta);
+
int=find(abs(mu(1:j))>eta);
int = max(LL+1,min(int)-extra):min(max(int)+extra,j);
int = max(LL+1,min(int)-extra):min(max(int)+extra,j);
case 2
case 2
-
int=find(abs(mu(1:j))&gt;=eta);
+
int=find(abs(mu(1:j))>=eta);
end
end
int = int(:);
int = int(:);
-
&lt;/pre&gt;
+
</pre>
=== lanbpro() ===
=== lanbpro() ===
-
&lt;pre&gt;
+
<pre>
function [U,B_k,V,p,ierr,work] = lanbpro(varargin)
function [U,B_k,V,p,ierr,work] = lanbpro(varargin)
Line 590: Line 582:
% Partial reorthogonalization is used to keep the columns of V_K and U_k
% Partial reorthogonalization is used to keep the columns of V_K and U_k
% semiorthogonal:
% semiorthogonal:
-
% MAX(DIAG((EYE(K) - V_K'*V_K))) &lt;= OPTIONS.delta
+
% MAX(DIAG((EYE(K) - V_K'*V_K))) <= OPTIONS.delta
% and
% and
-
% MAX(DIAG((EYE(K) - U_K'*U_K))) &lt;= OPTIONS.delta.
+
% MAX(DIAG((EYE(K) - U_K'*U_K))) <= OPTIONS.delta.
%
%
% B_k = LANBPRO(...) returns the bidiagonal matrix only.
% B_k = LANBPRO(...) returns the bidiagonal matrix only.
Line 648: Line 640:
% R.M. Larsen, Ph.D. Thesis, Aarhus University, 1998.
% R.M. Larsen, Ph.D. Thesis, Aarhus University, 1998.
%
%
-
% G. H. Golub &amp; C. F. Van Loan, &quot;Matrix Computations&quot;,
+
% G. H. Golub & C. F. Van Loan, "Matrix Computations",
% 3. Ed., Johns Hopkins, 1996. Section 9.3.4.
% 3. Ed., Johns Hopkins, 1996. Section 9.3.4.
%
%
Line 670: Line 662:
end
end
-
if nargin&lt;1 | length(varargin)&lt;2
+
if nargin<1 | length(varargin)<2
error('Not enough input arguments.');
error('Not enough input arguments.');
end
end
Line 686: Line 678:
end
end
k=varargin{2};
k=varargin{2};
-
if narg &gt;= 3 &amp; ~isempty(varargin{3});
+
if narg >= 3 & ~isempty(varargin{3});
p = varargin{3};
p = varargin{3};
else
else
p = rand(m,1)-0.5;
p = rand(m,1)-0.5;
end
end
-
if narg &lt; 4, options = []; else options=varargin{4}; end
+
if narg < 4, options = []; else options=varargin{4}; end
-
if narg &gt; 4
+
if narg > 4
-
if narg&lt;7
+
if narg<7
error('All or none of U_old, B_old and V_old must be provided.')
error('All or none of U_old, B_old and V_old must be provided.')
else
else
Line 701: Line 693:
U = []; B_k = []; V = [];
U = []; B_k = []; V = [];
end
end
-
if narg &gt; 7, anorm=varargin{8}; else anorm = []; end
+
if narg > 7, anorm=varargin{8}; else anorm = []; end
else
else
-
if narg&lt;5
+
if narg<5
error('Not enough input arguments.');
error('Not enough input arguments.');
end
end
Line 716: Line 708:
end
end
k=varargin{5};
k=varargin{5};
-
if narg &lt; 6, p = rand(m,1)-0.5; else p=varargin{6}; end
+
if narg < 6, p = rand(m,1)-0.5; else p=varargin{6}; end
-
if narg &lt; 7, options = []; else options=varargin{7}; end
+
if narg < 7, options = []; else options=varargin{7}; end
-
if narg &gt; 7
+
if narg > 7
-
if narg &lt; 10
+
if narg < 10
error('All or none of U_old, B_old and V_old must be provided.')
error('All or none of U_old, B_old and V_old must be provided.')
else
else
Line 727: Line 719:
U = []; B_k = []; V=[];
U = []; B_k = []; V=[];
end
end
-
if narg &gt; 10, anorm=varargin{11}; else anorm = []; end
+
if narg > 10, anorm=varargin{11}; else anorm = []; end
end
end
Line 740: Line 732:
U = 1; B_k = feval(A,1); V = 1; p = 0; ierr = 0; work = zeros(2,2);
U = 1; B_k = feval(A,1); V = 1; p = 0; ierr = 0; work = zeros(2,2);
end
end
-
if nargout&lt;3
+
if nargout<3
U = B_k;
U = B_k;
end
end
Line 760: Line 752:
% Parse options struct
% Parse options struct
-
if ~isempty(options) &amp; isstruct(options)
+
if ~isempty(options) & isstruct(options)
c = fieldnames(options);
c = fieldnames(options);
for i=1:length(c)
for i=1:length(c)
Line 807: Line 799:
V = [V, zeros(n,k-j)];
V = [V, zeros(n,k-j)];
alpha = zeros(k+1,1); beta = zeros(k+1,1);
alpha = zeros(k+1,1); beta = zeros(k+1,1);
-
alpha(1:j) = diag(B_k); if j&gt;1 beta(2:j) = diag(B_k,-1); end
+
alpha(1:j) = diag(B_k); if j>1 beta(2:j) = diag(B_k,-1); end
beta(j+1) = norm(p);
beta(j+1) = norm(p);
% Reorthogonalize p.
% Reorthogonalize p.
-
if j&lt;k &amp; beta(j+1)*delta &lt; anorm*eps,
+
if j<k & beta(j+1)*delta < anorm*eps,
fro = 1;
fro = 1;
ierr = j;
ierr = j;
Line 890: Line 882:
% Extended local reorthogonalization
% Extended local reorthogonalization
-
if alpha(j)&lt;gamma*beta(j) &amp; elr &amp; ~fro
+
if alpha(j)<gamma*beta(j) & elr & ~fro
normold = alpha(j);
normold = alpha(j);
stop = 0;
stop = 0;
Line 900: Line 892:
beta(j) = beta(j) + t;
beta(j) = beta(j) + t;
end
end
-
if alpha(j)&gt;=gamma*normold
+
if alpha(j)>=gamma*normold
stop = 1;
stop = 1;
else
else
Line 917: Line 909:
end
end
-
if ~fro &amp; alpha(j) ~= 0
+
if ~fro & alpha(j) ~= 0
% Update estimates of the level of orthogonality for the
% Update estimates of the level of orthogonality for the
% columns 1 through j-1 in V.
% columns 1 through j-1 in V.
Line 924: Line 916:
end
end
-
if j&gt;1 &amp; LANBPRO_TRUTH
+
if j>1 & LANBPRO_TRUTH
NU(1:j-1,j-1) = nu(1:j-1);
NU(1:j-1,j-1) = nu(1:j-1);
NUTRUE(1:j-1,j-1) = V(:,1:j-1)'*r/alpha(j);
NUTRUE(1:j-1,j-1) = V(:,1:j-1)'*r/alpha(j);
end
end
-
if elr&gt;0
+
if elr>0
nu(j-1) = n2*eps;
nu(j-1) = n2*eps;
end
end
% IF level of orthogonality is worse than delta THEN
% IF level of orthogonality is worse than delta THEN
-
% Reorthogonalize v_j against some previous v_i's, 0&lt;=i&lt;j.
+
% Reorthogonalize v_j against some previous v_i's, 0<=i<j.
-
if onesided~=-1 &amp; ( fro | numax(j) &gt; delta | force_reorth ) &amp; alpha(j)~=0
+
if onesided~=-1 & ( fro | numax(j) > delta | force_reorth ) & alpha(j)~=0
% Decide which vectors to orthogonalize against:
% Decide which vectors to orthogonalize against:
if fro | eta==0
if fro | eta==0
Line 963: Line 955:
% Check for convergence or failure to maintain semiorthogonality
% Check for convergence or failure to maintain semiorthogonality
-
if alpha(j) &lt; max(n,m)*anorm*eps &amp; j&lt;k,
+
if alpha(j) < max(n,m)*anorm*eps & j<k,
-
% If alpha is &quot;small&quot; we deflate by setting it
+
% If alpha is "small" we deflate by setting it
% to 0 and attempt to restart with a basis for a new
% to 0 and attempt to restart with a basis for a new
% invariant subspace by replacing r with a random starting vector:
% invariant subspace by replacing r with a random starting vector:
Line 985: Line 977:
npv = npv + rr*length(int(:)); nreorthv = nreorthv + 1;
npv = npv + rr*length(int(:)); nreorthv = nreorthv + 1;
nu(int) = n2*eps;
nu(int) = n2*eps;
-
if nrmnew &gt; 0
+
if nrmnew > 0
% A vector numerically orthogonal to span(Q_k(:,1:j)) was found.
% A vector numerically orthogonal to span(Q_k(:,1:j)) was found.
% Continue iteration.
% Continue iteration.
Line 999: Line 991:
r=r/nrmnew; % Continue with new normalized r as starting vector.
r=r/nrmnew; % Continue with new normalized r as starting vector.
force_reorth = 1;
force_reorth = 1;
-
if delta&gt;0
+
if delta>0
fro = 0; % Turn off full reorthogonalization.
fro = 0; % Turn off full reorthogonalization.
end
end
end
end
-
elseif j&lt;k &amp; ~fro &amp; anorm*eps &gt; delta*alpha(j)
+
elseif j<k & ~fro & anorm*eps > delta*alpha(j)
% fro = 1;
% fro = 1;
ierr = j;
ierr = j;
end
end
-
if j&gt;1 &amp; LANBPRO_TRUTH
+
if j>1 & LANBPRO_TRUTH
NU_AFTER(1:j-1,j-1) = nu(1:j-1);
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);
NUTRUE_AFTER(1:j-1,j-1) = V(:,1:j-1)'*r/alpha(j);
Line 1,033: Line 1,025:
beta(j+1) = norm(p);
beta(j+1) = norm(p);
% Extended local reorthogonalization
% Extended local reorthogonalization
-
if beta(j+1)&lt;gamma*alpha(j) &amp; elr &amp; ~fro
+
if beta(j+1)<gamma*alpha(j) & elr & ~fro
normold = beta(j+1);
normold = beta(j+1);
stop = 0;
stop = 0;
Line 1,043: Line 1,035:
alpha(j) = alpha(j) + t;
alpha(j) = alpha(j) + t;
end
end
-
if beta(j+1) &gt;= gamma*normold
+
if beta(j+1) >= gamma*normold
stop = 1;
stop = 1;
else
else
Line 1,062: Line 1,054:
end
end
-
if ~fro &amp; beta(j+1) ~= 0
+
if ~fro & beta(j+1) ~= 0
% Update estimates of the level of orthogonality for the columns of V.
% Update estimates of the level of orthogonality for the columns of V.
mu = update_mu(mu,nu,j,alpha,beta,anorm);
mu = update_mu(mu,nu,j,alpha,beta,anorm);
Line 1,073: Line 1,065:
end
end
-
if elr&gt;0
+
if elr>0
mu(j) = m2*eps;
mu(j) = m2*eps;
end
end
% IF level of orthogonality is worse than delta THEN
% IF level of orthogonality is worse than delta THEN
-
% Reorthogonalize u_{j+1} against some previous u_i's, 0&lt;=i&lt;=j.
+
% Reorthogonalize u_{j+1} against some previous u_i's, 0<=i<=j.
-
if onesided~=1 &amp; (fro | mumax(j) &gt; delta | force_reorth) &amp; beta(j+1)~=0
+
if onesided~=1 & (fro | mumax(j) > delta | force_reorth) & beta(j+1)~=0
% Decide which vectors to orthogonalize against.
% Decide which vectors to orthogonalize against.
if fro | eta==0
if fro | eta==0
Line 1,108: Line 1,100:
% Check for convergence or failure to maintain semiorthogonality
% Check for convergence or failure to maintain semiorthogonality
-
if beta(j+1) &lt; max(m,n)*anorm*eps &amp; j&lt;k,
+
if beta(j+1) < max(m,n)*anorm*eps & j<k,
-
% If beta is &quot;small&quot; we deflate by setting it
+
% If beta is "small" we deflate by setting it
% to 0 and attempt to restart with a basis for a new
% to 0 and attempt to restart with a basis for a new
% invariant subspace by replacing p with a random starting vector:
% invariant subspace by replacing p with a random starting vector:
Line 1,130: Line 1,122:
npu = npu + rr*length(int(:)); nreorthu = nreorthu + 1;
npu = npu + rr*length(int(:)); nreorthu = nreorthu + 1;
mu(int) = m2*eps;
mu(int) = m2*eps;
-
if nrmnew &gt; 0
+
if nrmnew > 0
% A vector numerically orthogonal to span(Q_k(:,1:j)) was found.
% A vector numerically orthogonal to span(Q_k(:,1:j)) was found.
% Continue iteration.
% Continue iteration.
Line 1,143: Line 1,135:
p=p/nrmnew; % Continue with new normalized p as starting vector.
p=p/nrmnew; % Continue with new normalized p as starting vector.
force_reorth = 1;
force_reorth = 1;
-
if delta&gt;0
+
if delta>0
fro = 0; % Turn off full reorthogonalization.
fro = 0; % Turn off full reorthogonalization.
end
end
end
end
-
elseif j&lt;k &amp; ~fro &amp; anorm*eps &gt; delta*beta(j+1)
+
elseif j<k & ~fro & anorm*eps > delta*beta(j+1)
% fro = 1;
% fro = 1;
ierr = j;
ierr = j;
Line 1,161: Line 1,153:
end
end
-
if j&lt;k
+
if j<k
k = j;
k = j;
end
end
Line 1,172: Line 1,164:
V = V(:,1:k);
V = V(:,1:k);
end
end
-
if nargout&gt;5
+
if nargout>5
work = [[nreorthu,npu];[nreorthv,npv]];
work = [[nreorthu,npu];[nreorthv,npv]];
end
end
Line 1,199: Line 1,191:
mu(1) = (mu(1) + sign(mu(1))*T) / beta(j+1);
mu(1) = (mu(1) + sign(mu(1))*T) / beta(j+1);
% Vectorized version of loop:
% Vectorized version of loop:
-
if j&gt;2
+
if j>2
k=2:j-1;
k=2:j-1;
mu(k) = alpha(k).*nu(k) + beta(k).*nu(k-1) - alpha(j)*mu(k);
mu(k) = alpha(k).*nu(k) + beta(k).*nu(k-1) - alpha(j)*mu(k);
Line 1,227: Line 1,219:
ainv = 1/alpha(j);
ainv = 1/alpha(j);
eps1 = 100*eps/2;
eps1 = 100*eps/2;
-
if j&gt;1
+
if j>1
k = 1:(j-1);
k = 1:(j-1);
% T = eps1*(pythag(alpha(k),beta(k+1)) + pythag(alpha(j),beta(j)));
% T = eps1*(pythag(alpha(k),beta(k+1)) + pythag(alpha(j),beta(j)));
Line 1,248: Line 1,240:
[m n] = size(y);
[m n] = size(y);
-
if m&gt;1 | n&gt;1
+
if m>1 | n>1
y = y(:); z=z(:);
y = y(:); z=z(:);
rmax = max(abs([y z]'))';
rmax = max(abs([y z]'))';
id=find(rmax==0);
id=find(rmax==0);
-
if length(id)&gt;0
+
if length(id)>0
rmax(id) = 1;
rmax(id) = 1;
x = rmax.*sqrt((y./rmax).^2 + (z./rmax).^2);
x = rmax.*sqrt((y./rmax).^2 + (z./rmax).^2);
Line 1,268: Line 1,260:
end
end
end
end
-
&lt;/pre&gt;
+
</pre>
=== refinebounds() ===
=== refinebounds() ===
-
&lt;pre&gt;
+
<pre>
function [bnd,gap] = refinebounds(D,bnd,tol1)
function [bnd,gap] = refinebounds(D,bnd,tol1)
%REFINEBONDS Refines error bounds for Ritz values based on gap-structure
%REFINEBONDS Refines error bounds for Ritz values based on gap-structure
Line 1,283: Line 1,275:
j = length(D);
j = length(D);
-
if j&lt;=1
+
if j<=1
return
return
end
end
Line 1,296: Line 1,288:
for l=[-1,1]
for l=[-1,1]
for i=((j+1)-l*(j-1))/2:l:mid-l
for i=((j+1)-l*(j-1))/2:l:mid-l
-
if abs(D(i+l)-D(i)) &lt; eps34*abs(D(i))
+
if abs(D(i+l)-D(i)) < eps34*abs(D(i))
-
if bnd(i)&gt;tol1 &amp; bnd(i+l)&gt;tol1
+
if bnd(i)>tol1 & bnd(i+l)>tol1
bnd(i+l) = pythag(bnd(i),bnd(i+l));
bnd(i+l) = pythag(bnd(i),bnd(i+l));
bnd(i) = 0;
bnd(i) = 0;
Line 1,309: Line 1,301:
gap(2:j) = min([gap(2:j);[D(2:j)-D(1:j-1)-bnd(1:j-1)]']);
gap(2:j) = min([gap(2:j);[D(2:j)-D(1:j-1)-bnd(1:j-1)]']);
gap = gap(:);
gap = gap(:);
-
I = find(gap&gt;bnd);
+
I = find(gap>bnd);
bnd(I) = bnd(I).*(bnd(I)./gap(I));
bnd(I) = bnd(I).*(bnd(I)./gap(I));
bnd(PERM) = bnd;
bnd(PERM) = bnd;
-
&lt;/pre&gt;
+
</pre>
=== reorth() ===
=== reorth() ===
-
&lt;pre&gt;
+
<pre>
function [r,normr,nre,s] = reorth(Q,r,normr,index,alpha,method)
function [r,normr,nre,s] = reorth(Q,r,normr,index,alpha,method)
Line 1,324: Line 1,316:
% reorthogonalizes R against the subset of columns of Q given by INDEX.
% reorthogonalizes R against the subset of columns of Q given by INDEX.
% If INDEX==[] then R is reorthogonalized all columns of Q.
% 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) &lt; ALPHA*NORMR,
+
% 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
% 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
% is once more decreased by more than a factor of ALPHA then R is
Line 1,336: Line 1,328:
% References:
% References:
-
% Aake Bjorck, &quot;Numerical Methods for Least Squares Problems&quot;,
+
% Aake Bjorck, "Numerical Methods for Least Squares Problems",
% SIAM, Philadelphia, 1996, pp. 68-69.
% SIAM, Philadelphia, 1996, pp. 68-69.
%
%
Line 1,351: Line 1,343:
% Check input arguments.
% Check input arguments.
%warning('PROPACK:NotUsingMex','Using slow matlab code for reorth.')
%warning('PROPACK:NotUsingMex','Using slow matlab code for reorth.')
-
if nargin&lt;2
+
if nargin<2
error('Not enough input arguments.')
error('Not enough input arguments.')
end
end
[n k1] = size(Q);
[n k1] = size(Q);
-
if nargin&lt;3 | isempty(normr)
+
if nargin<3 | isempty(normr)
% normr = norm(r);
% normr = norm(r);
normr = sqrt(r'*r);
normr = sqrt(r'*r);
end
end
-
if nargin&lt;4 | isempty(index)
+
if nargin<4 | isempty(index)
k=k1;
k=k1;
index = [1:k]';
index = [1:k]';
Line 1,365: Line 1,357:
else
else
k = length(index);
k = length(index);
-
if k==k1 &amp; index(:)==[1:k]'
+
if k==k1 & index(:)==[1:k]'
simple = 1;
simple = 1;
else
else
Line 1,371: Line 1,363:
end
end
end
end
-
if nargin&lt;5 | isempty(alpha)
+
if nargin<5 | isempty(alpha)
alpha=0.5; % This choice garanties that
alpha=0.5; % This choice garanties that
-
% || Q^T*r_new - e_{k+1} ||_2 &lt;= 2*eps*||r_new||_2,
+
% || Q^T*r_new - e_{k+1} ||_2 <= 2*eps*||r_new||_2,
% cf. Kahans ``twice is enough'' statement proved in
% cf. Kahans ``twice is enough'' statement proved in
% Parletts book.
% Parletts book.
end
end
-
if nargin&lt;6 | isempty(method)
+
if nargin<6 | isempty(method)
method = 0;
method = 0;
end
end
Line 1,383: Line 1,375:
return
return
end
end
-
if nargout&gt;3
+
if nargout>3
s = zeros(k,1);
s = zeros(k,1);
end
end
Line 1,390: Line 1,382:
normr_old = 0;
normr_old = 0;
nre = 0;
nre = 0;
-
while normr &lt; alpha*normr_old | nre==0
+
while normr < alpha*normr_old | nre==0
if method==1
if method==1
if simple
if simple
Line 1,405: Line 1,397:
end
end
end
end
-
if nargout&gt;3
+
if nargout>3
s = s + t;
s = s + t;
end
end
Line 1,412: Line 1,404:
normr = sqrt(r'*r);
normr = sqrt(r'*r);
nre = nre + 1;
nre = nre + 1;
-
if nre &gt; 4
+
if nre > 4
-
% r is in span(Q) to full accuracy =&gt; accept r = 0 as the new vector.
+
% r is in span(Q) to full accuracy => accept r = 0 as the new vector.
r = zeros(n,1);
r = zeros(n,1);
normr = 0;
normr = 0;
Line 1,419: Line 1,411:
end
end
end
end
-
&lt;/pre&gt;
+
</pre>
=== updateSparse.c ===
=== updateSparse.c ===
A precompiled [http://convexoptimization.com/TOOLS/updateSparse.mexw32 Intel Windows-32bit Matlab mex file is here].
A precompiled [http://convexoptimization.com/TOOLS/updateSparse.mexw32 Intel Windows-32bit Matlab mex file is here].
-
Otherwise, compile this C program in Matlab via command &lt;pre&gt;mex updateSparse.c&lt;/pre&gt;
+
Otherwise, compile this C program in Matlab via command <pre>mex updateSparse.c</pre>
-
&lt;pre&gt;
+
<pre>
/*
/*
* Stephen Becker, 11/10/08
* Stephen Becker, 11/10/08
Line 1,439: Line 1,431:
* */
* */
-
#include &quot;mex.h&quot;
+
#include "mex.h"
#ifndef true
#ifndef true
#define true 1
#define true 1
Line 1,448: Line 1,440:
void printUsage() {
void printUsage() {
-
mexPrintf(&quot;usage:\tupdateSparse(Y,b)\nchanges the sparse Y matrix&quot;);
+
mexPrintf("usage:\tupdateSparse(Y,b)\nchanges the sparse Y matrix");
-
mexPrintf(&quot; to have values b\non its nonzero elements. Be careful:\n\t&quot;);
+
mexPrintf(" to have values b\non its nonzero elements. Be careful:\n\t");
-
mexPrintf(&quot;this assumes b is sorted in the appropriate order!\n&quot;);
+
mexPrintf("this assumes b is sorted in the appropriate order!\n");
-
mexPrintf(&quot;If b (i.e. the index omega, where we want to perform Y(omega)=b)\n&quot;);
+
mexPrintf("If b (i.e. the index omega, where we want to perform Y(omega)=b)\n");
-
mexPrintf(&quot; is unsorted, then call the command as follows:\n&quot;);
+
mexPrintf(" is unsorted, then call the command as follows:\n");
-
mexPrintf(&quot;\tupdateSparse(Y,b,omegaIndx)\n&quot;);
+
mexPrintf("\tupdateSparse(Y,b,omegaIndx)\n");
-
mexPrintf(&quot;where [temp,omegaIndx] = sort(omega)\n&quot;);
+
mexPrintf("where [temp,omegaIndx] = sort(omega)\n");
}
}
Line 1,468: Line 1,460:
/* Check for proper number of input and output arguments */
/* Check for proper number of input and output arguments */
-
if ( (nrhs &lt; 2) || (nrhs &gt; 3) ) {
+
if ( (nrhs < 2) || (nrhs > 3) ) {
printUsage();
printUsage();
-
mexErrMsgTxt(&quot;Needs 2 or 3 input arguments&quot;);
+
mexErrMsgTxt("Needs 2 or 3 input arguments");
}
}
if ( nrhs == 3 ) SORTED = false;
if ( nrhs == 3 ) SORTED = false;
-
if(nlhs &gt; 0){
+
if(nlhs > 0){
printUsage();
printUsage();
-
mexErrMsgTxt(&quot;No output arguments!&quot;);
+
mexErrMsgTxt("No output arguments!");
}
}
Line 1,481: Line 1,473:
if (!(mxIsSparse(prhs[0])) || !((mxIsDouble(prhs[1]))) ){
if (!(mxIsSparse(prhs[0])) || !((mxIsDouble(prhs[1]))) ){
printUsage();
printUsage();
-
mexErrMsgTxt(&quot;Input arguments wrong data-type (must be sparse, double).&quot;);
+
mexErrMsgTxt("Input arguments wrong data-type (must be sparse, double).");
}
}
Line 1,488: Line 1,480:
N = mxGetN( prhs[1] );
N = mxGetN( prhs[1] );
M = mxGetM( prhs[1] );
M = mxGetM( prhs[1] );
-
if ( (M&gt;1) &amp;&amp; (N&gt;1) ) {
+
if ( (M>1) && (N>1) ) {
printUsage();
printUsage();
-
mexErrMsgTxt(&quot;Second argument must be a vector&quot;);
+
mexErrMsgTxt("Second argument must be a vector");
}
}
-
N = (N&gt;M) ? N : M;
+
N = (N>M) ? N : M;
Line 1,499: Line 1,491:
if ( M != N ) {
if ( M != N ) {
printUsage();
printUsage();
-
mexErrMsgTxt(&quot;Length of second argument must match nnz of first argument&quot;);
+
mexErrMsgTxt("Length of second argument must match nnz of first argument");
}
}
Line 1,506: Line 1,498:
m = mxGetM( prhs[2] );
m = mxGetM( prhs[2] );
n = mxGetN( prhs[2] );
n = mxGetN( prhs[2] );
-
if ( (m&gt;1) &amp;&amp; (n&gt;1) ) {
+
if ( (m>1) && (n>1) ) {
printUsage();
printUsage();
-
mexErrMsgTxt(&quot;Third argument must be a vector&quot;);
+
mexErrMsgTxt("Third argument must be a vector");
}
}
-
n = (n&gt;m) ? n : m;
+
n = (n>m) ? n : m;
if ( n != N ) {
if ( n != N ) {
printUsage();
printUsage();
-
mexErrMsgTxt(&quot;Third argument must be same length as second argument&quot;);
+
mexErrMsgTxt("Third argument must be same length as second argument");
}
}
omega = mxGetPr( prhs[2] );
omega = mxGetPr( prhs[2] );
Line 1,524: Line 1,516:
if (SORTED) {
if (SORTED) {
/* And here's the really fancy part: */
/* And here's the really fancy part: */
-
for ( i=0 ; i &lt; N ; i++ )
+
for ( i=0 ; i < N ; i++ )
S[i] = b[i];
S[i] = b[i];
} else {
} else {
-
for ( i=0 ; i &lt; N ; i++ ) {
+
for ( i=0 ; i < N ; i++ ) {
/* this is a little slow, but I should really check
/* this is a little slow, but I should really check
* to make sure the index is not out-of-bounds, otherwise
* to make sure the index is not out-of-bounds, otherwise
* Matlab could crash */
* Matlab could crash */
j = (int)omega[i]-1; /* the -1 because Matlab is 1-based */
j = (int)omega[i]-1; /* the -1 because Matlab is 1-based */
-
if (j &gt;= N){
+
if (j >= N){
printUsage();
printUsage();
-
mexErrMsgTxt(&quot;Third argument must have values &lt; length of 2nd argument&quot;);
+
mexErrMsgTxt("Third argument must have values < length of 2nd argument");
}
}
/* S[ j ] = b[i]; */ /* this is incorrect */
/* S[ j ] = b[i]; */ /* this is incorrect */
Line 1,541: Line 1,533:
}
}
}
}
-
&lt;/pre&gt;
+
</pre>

Current revision

Contents

Matlab demonstration of Cai, Candès, & Shen

A Singular Value Thresholding Algorithm for Matrix Completion, 2008

This Matlab code below is working, complete, debugged, and corresponds to the paper cited above. It is capable of solving smaller matrix completion problems quite well, as the demonstration program shows.

Singular Value Thresholding (SVT) code for solving larger completion problems is presently under development. It is written in C and Fortran and can be compiled in Matlab. Read more...

% 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 command
mex 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 */
        }
    }
}
Personal tools