Index exceeds matrix dimensions.错误_百度知道From Wikimization
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.
for solving larger completion problems is presently under development.
It is written in C and Fortran and can be compiled in Matlab.
% Written by: Emmanuel Candes
% Email: emmanuel@acm.caltech.edu
% Created: October 2008
%% Setup a matrix
randn('state',2008);
rand('state',2008);
M = randn(n,r)*randn(r,n);
df = r*(2*n-r);
oversampling = 5;
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
[U,S,V,numiter] = SVT(n,Omega,data,delta,maxiter,tol);
%% 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)))
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
[U,S,V,numiter] = SVT(n,Omega,b,delta,maxiter,opts)
% 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);
fprintf('\nIteration:
for k = 1:maxiter
fprintf('\b\b\b%3d',k);
s = r + 1;
[U,Sigma,V] = lansvd(Y,s,'L');
OK = Sigma(s,s) &=
sigma = diag(Sigma);
r = sum(sigma & tau);
U = U(:,1:r);
V = V(:,1:r);
sigma = sigma(1:r) -
Sigma = diag(sigma);
A = U*diag(sigma)*V';
x = A(Omega);
if norm(x-b)/normb & tol
y = y + delta*(b-x);
updateSparse(Y,y,indx);
fprintf('\n');
function [U,S,V,bnd,j] = lansvd(varargin)
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
OPTIONS.tol
Convergence tolerance
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'&: iterated modified Gram-Schmidt
'1'&: iterated classical Gram-Schmidt
OPTIONS.elr
If equal to 1 then extended local
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.');
A = varargin{1};
if ~isstr(A)
if ~isreal(A)
error('A must be real')
[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';
sigma=varargin{3}; end
if length(varargin) & 4, options = [];
options=varargin{4}; end
if length(varargin)&4
error('Not enough input arguments.');
Atrans = varargin{2};
if ~isstr(Atrans)
error('Atransfunc must be the name of a function')
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
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.')
% 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);
U = eye(m,k); S = zeros(k,k);
V = eye(n,k);
bnd = zeros(k,1);
elseif min(n,m) == 1 & k&0
if isstr(A)
&% Extract the single column or row of A
A = feval(A,1);
A = feval(Atrans,1)';
if nargout==1
U = norm(A);
[U,S,V] = svd(full(A));
% A is the matrix of all zeros (not detectable if A is defined by an m-file)
if isnumeric(A)
if nargout&3
U = zeros(k,1);
U = eye(m,k); S = zeros(k,k);
V = eye(n,k);
bnd = zeros(k,1);
lanmax = min(m,n);
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
% 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')
lanmax = min(lanmax,min(m,n));
if k&lanmax
error('K must satisfy
K &= LANMAX &= MIN(M,N).');
%%%%%%%%%%%%%%%%%%%%% Here begins the computation &%%%%%%%%%%%%%%%%%%%%%%
if strcmp(sigma,'S')
if isstr(A)
error('Shift-and-invert works only when the matrix A is given explicitly.');
&% Prepare for shift-and-invert Lanczos.
if issparse(A)
pmmd = colmmd(A);
A.A = A(:,pmmd);
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)
[A.Q,A.R] = qr(A.A,0);
A.Rt = A.R';
p = A.Q'*p;&% project starting vector on span(Q1)
error('Sorry, shift-and-invert for m&n not implemented yet!')
A.R = qr(A.A',0);
A.Rt = A.R';
condR = condest(A.R);
if condR & 1/eps
error(['A is rank deficient or too ill-conditioned to do shift-and-' ...
' invert.'])
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);
[U,B,V,p,ierr,w] = lanbpro(A,Atrans,m,n,j,p,options,U,B,V,anorm);
work= work +
if ierr&0&% Invariant subspace of dimension -ierr found.
&%%%%%%%%%%%%%%%%%% 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&%%%%%%%%%%%%%%%%%%%%
while i&=min(j,k)
if (bnd(i) &= tol*abs(S(i)))
neig = neig + 1;
i = min(j,k)+1;
&%%%%%%%%%% Check whether to stop or to extend the Krylov basis?&%%%%%%%%%%
if ierr&0&% Invariant subspace found
warning(['Invariant subspace of dimension ',num2str(j-1),' found.'])
if j&=lanmax&% Maximal dimension of Krylov subspace reached. Bail out
if j&=min(m,n)
if neig&ksave
warning(['Maximum dimension of Krylov subspace exceeded prior',...
' to convergence.']);
&% Increase dimension of Krylov subspace
&% 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)));
&% 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);
j = ceil(min(j+1,lanmax));
nrestart = nrestart + 1;
%%%%%%%%%%%%%%%% 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);
&% 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,:);
U = U*P(1:j,:);
nq = norm(V(:,i));
if isfinite(nq) & nq~=0 & nq~=1
V(:,i) = V(:,i)/
nq = norm(U(:,i));
if isfinite(nq) & nq~=0 & nq~=1
U(:,i) = U(:,i)/
% Pick out desired part the spectrum
S = S(1:k);
bnd = bnd(1:k);
if strcmp(sigma,'S')
[S,p] = sort(-1./S);
bnd = bnd(p);
if nargout&2
if issparse(A.A)
U = A.A*(A.R\U(:,p));
V(pmmd,:) = V(:,p);
U = A.Q(:,1:min(m,n))*U(:,p);
V = V(:,p);
if nargout&3
S = B;&% Undocumented feature -
for checking B.
S = diag(S);
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.
&: Sub-diagonal elements.
% Output parameters:
sigma(1:k) &: Computed eigenvalues.
&: 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')
B = spdiags([alpha(:),beta(:)],[0,-1],k+1,k);
[U,S,V] = svd(full(B),0);
sigma = diag(S);
bnd = U(end,1:k)';
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
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
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.')
switch strategy
I0 = find(abs(mu(1:j))&=delta);
if length(I0)==0
[mm,I0] = max(abs(mu(1:j)));
int = zeros(j,1);
for i = 1:length(I0)
for r=I0(i):-1:1
if abs(mu(r))&eta | int(r)==1
int(r) = 1;
int(max(1,r-extra+1):r) = 1;
for s=I0(i)+1:j
if abs(mu(s))&eta | int(s)==1
int(s) = 1;
int(s:min(j,s+extra-1)) = 1;
int(1:LL) = 0;
int = find(int);
int=find(abs(mu(1:j))&eta);
int = max(LL+1,min(int)-extra):min(max(int)+extra,j);
int=find(abs(mu(1:j))&=eta);
int = int(:);
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
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
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
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
if nargin&1 | length(varargin)&2
error('Not enough input arguments.');
narg=length(varargin);
A = varargin{1};
if isnumeric(A) | isstruct(A)
if isnumeric(A)
if ~isreal(A)
error('A must be real')
[m n] = size(A);
elseif isstruct(A)
[m n] = size(A.R);
k=varargin{2};
if narg &= 3 & ~isempty(varargin{3});
p = varargin{3};
p = rand(m,1)-0.5;
if narg & 4, options = []; else options=varargin{4}; end
if narg & 4
error('All or none of U_old, B_old and V_old must be provided.')
U = varargin{5}; B_k = varargin{6}; V = varargin{7};
U = []; B_k = []; V = [];
if narg & 7, anorm=varargin{8}; else anorm = []; end
error('Not enough input arguments.');
Atrans = varargin{2};
if ~isstr(Atrans)
error('Afunc and Atransfunc must be names of m-files')
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.')
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
error('All or none of U_old, B_old and V_old must be provided.')
U = varargin{8}; B_k = varargin{9}; V = varargin{10};
U = []; B_k = []; V=[];
if narg & 10, anorm=varargin{11}; else anorm = [];
% Quick return for min(m,n) equal to 0 or 1.
if min(m,n) == 0
work = zeros(2,2);
min(m,n) == 1
if isnumeric(A)
p = 0; ierr = 0; work = zeros(2,2);
B_k = feval(A,1); V = 1; p = 0; ierr = 0; work = zeros(2,2);
if nargout&3
% Set options.
%m2 = 3/2*(sqrt(m)+1);
%n2 = 3/2*(sqrt(n)+1);
delta = sqrt(eps/k);&% Desired level of orthogonality.
eta = eps^(3/4)/sqrt(k);
&% Level of orth. after reorthogonalization.
&% Flag for switching between iterated MGS and CGS.
&% 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');
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
waitbarh = waitbar(0,'Lanczos bidiagonalization in progress...');
if isempty(anorm)
anorm = []; est_anorm=1;
est_anorm=0;
% 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;
% 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;
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,
int = [1:j]';
[p,beta(j+1),rr] = reorth(U,p,beta(j+1),int,gamma,cgs);
nreorthu = 1;
force_reorth= 1;
&% Compute Gerscgorin bound on ||B_k||_2
if est_anorm
anorm = FUDGE*sqrt(norm(B_k'*B_k,1));
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;
if isnumeric(A)
if delta==0
fro = 1;&% The user has requested full reorthogonalization.
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);
% Perform Lanczos bidiagonalization with partial reorthogonalization.
for j=j0:k
waitbar(j/k,waitbarh)
if beta(j) ~= 0
U(:,j) = p/beta(j);
&% Replace norm estimate with largest Ritz value.
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;
&%%%%%%%%%% Lanczos step to generate v_j.&%%%%%%%%%%%%%
if isnumeric(A)
r = At*U(:,1);
elseif isstruct(A)
r = A.R\U(:,1);
r = feval(Atrans,U(:,1));
alpha(1) = norm(r);
if est_anorm
anorm = FUDGE*alpha(1);
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);
r = feval(Atrans,U(:,j))
- beta(j)*V(:,j-1);
alpha(j) = norm(r);
&% Extended local reorthogonalization
if alpha(j)&gamma*beta(j) & elr & ~fro
normold = alpha(j);
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) +
if alpha(j)&=gamma*normold
normold = alpha(j);
if est_anorm
anorm = max(anorm,FUDGE*sqrt(alpha(1)^2+beta(2)^2+alpha(2)*beta(2)));
anorm = max(anorm,FUDGE*sqrt(alpha(j-1)^2+beta(j)^2+alpha(j-1)* ...
beta(j-1) + alpha(j)*beta(j)));
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)));
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);
nu(j-1) = n2*
&% 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);
&% 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* &% Reset nu for orthogonalized vectors.
&% If necessary force reorthogonalization of u_{j+1}
&% to avoid spillover
if force_reorth==0
force_reorth = 1;
force_reorth = 0;
nreorthv = nreorthv + 1;
&% 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:
&%disp('restarting, alpha = 0')
alpha(j) = 0;
bailout = 1;
for attempt=1:3
r = rand(m,1)-0.5;
if isnumeric(A)
elseif isstruct(A)
r = A.R\r;
r = feval(Atrans,r);
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*
if nrmnew & 0
% A vector numerically orthogonal to span(Q_k(:,1:j)) was found.
% Continue iteration.
bailout=0;
if bailout
ierr = -j;
r=r/&% Continue with new normalized r as starting vector.
force_reorth = 1;
if delta&0
&% Turn off full reorthogonalization.
j&k & ~fro & anorm*eps & delta*alpha(j)
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);
if alpha(j) ~= 0
V(:,j) = r/alpha(j);
&%%%%%%%%%% Lanczos step to generate u_{j+1}.&%%%%%%%%%%%%%
waitbar((2*j+1)/(2*k),waitbarh)
if isnumeric(A)
p = A*V(:,j) - alpha(j)*U(:,j);
elseif isstruct(A)
p = A.Rt\V(:,j) - alpha(j)*U(:,j);
p = feval(A,V(:,j)) - alpha(j)*U(:,j);
beta(j+1) = norm(p);
&% Extended local reorthogonalization
if beta(j+1)&gamma*alpha(j) & elr & ~fro
normold = beta(j+1);
while ~stop
t = U(:,j)'*p;
p = p - U(:,j)*t;
beta(j+1) = norm(p);
if alpha(j) ~= 0
alpha(j) = alpha(j) +
if beta(j+1) &= gamma*normold
normold = beta(j+1);
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!
anorm = max(anorm,FUDGE*pythag(alpha(1),beta(2)));
anorm = max(anorm,FUDGE*sqrt(alpha(j)^2+beta(j+1)^2 + alpha(j)*beta(j)));
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)));
if LANBPRO_TRUTH==1
MU(1:j,j) = mu(1:j);
MUTRUE(1:j,j) = U(:,1:j)'*p/beta(j+1);
mu(j) = m2*
&% 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);
int = [ max(int)+1];
&% 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)
&% 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*
if force_reorth==0
force_reorth = 1;&% Force reorthogonalization of v_{j+1}.
force_reorth = 0;
&% Check for convergence or failure to maintain semiorthogonality
if beta(j+1) & max(m,n)*anorm*eps
&% 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:
&%disp('restarting, beta = 0')
beta(j+1) = 0;
bailout = 1;
for attempt=1:3
p = rand(n,1)-0.5;
if isnumeric(A)
elseif isstruct(A)
p = A.Rt\p;
p = feval(A,p);
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*
if nrmnew & 0
% A vector numerically orthogonal to span(Q_k(:,1:j)) was found.
% Continue iteration.
bailout=0;
if bailout
ierr = -j;
p=p/&% Continue with new normalized p as starting vector.
force_reorth = 1;
if delta&0
&% Turn off full reorthogonalization.
j&k & ~fro & anorm*eps & delta*beta(j+1)
if LANBPRO_TRUTH==1
MU_AFTER(1:j,j) = mu(1:j);
MUTRUE_AFTER(1:j,j) = U(:,1:j)'*p/beta(j+1);
close(waitbarh)
B_k = spdiags([alpha(1:k) [beta(2:k);0]],[0 -1],k,k);
if nargout==1
elseif k~=size(U,2) | k~=size(V,2)
U = U(:,1:k);
V = V(:,1:k);
if nargout&5
work = [[nreorthu,npu];[nreorthv,npv]];
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);
eps1 = 100*eps/2;
T = eps1*(pythag(alpha(1),beta(2)) + pythag(alpha(1),beta(1)));
T = T + eps1*
mu(1) = T / beta(2);
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*
mu(1) = (mu(1) + sign(mu(1))*T) / beta(j+1);
&% Vectorized version of loop:
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*
mu(k) = binv*(mu(k) + sign(mu(k)).*T);
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*
mu(j) = beta(j)*nu(j-1);
mu(j) = (mu(j) + sign(mu(j))*T) / beta(j+1);
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.
ainv = 1/alpha(j);
eps1 = 100*eps/2;
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*
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);
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 = rmax.*sqrt((y./rmax).^2 + (z./rmax).^2);
x = reshape(x,m,n);
rmax = max(abs([y;z]));
if (rmax==0)
x = rmax*sqrt((y/rmax)^2 + (z/rmax)^2);
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);
% 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;
% 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) =
function [r,normr,nre,s] = reorth(Q,r,normr,index,alpha,method)
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.')
[n k1] = size(Q);
if nargin&3 | isempty(normr)
normr = norm(r);
normr = sqrt(r'*r);
if nargin&4 | isempty(index)
index = [1:k]';
simple = 1;
k = length(index);
if k==k1 & index(:)==[1:k]'
simple = 1;
simple = 0;
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.
if nargin&6 | isempty(method)
method = 0;
if k==0 | n==0
if nargout&3
s = zeros(k,1);
normr_old = 0;
while normr & alpha*normr_old | nre==0
if method==1
r = r - Q*t;
t = Q(:,index)'*r;
r = r - Q(:,index)*t;
for i=index,
t = Q(:,i)'*r;
r = r - Q(:,i)*t;
if nargout&3
normr_old =
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;
A precompiled .
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
#ifndef false
#define false 0
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(
mxArray *plhs[],
int nrhs, const mxArray *prhs[]
/* Declare variable */
int M, N, i, j, m,
double *b, *S, *
int SORTED =
/* 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 =
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&:
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];
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}