# Matrix Completion.m

(Difference between revisions)
 Revision as of 02:06, 16 February 2009 (edit)← Previous diff Revision as of 02:15, 16 February 2009 (edit) (undo)Next diff → Line 928: Line 928: end end end end + + + === 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.
+                                                                %
+
+                                                                % 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 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=lanmax % Maximal dimension of Krylov subspace reached. Bail out
+                                                                if j>=min(m,n)
+                                                                neig = ksave;
+                                                                break;
+                                                                end
+                                                                if neig0
+                                                                % 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
+
+ + === 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 === + The precompiled Intel Windows-32bit Matlab [http://convexoptimization.com/TOOLS/updateSparse.mexw32 mex file is here]. Otherwise, compile this C program in Matlab via command
+
+                                                                /*
+                                                                * Stephen Becker, 11/10/08
+                                                                * Updates a sparse vector very quickly
+                                                                * calling format:
+                                                                * 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("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 */
+                                                                }
+                                                                }
+
+
+                                                                }

## Matlab demonstration of Cai, Candès, & Shen

### test_SVT.m

% Written by: Emmanuel Candes
% Email: emmanuel@acm.caltech.edu
% Created: October 2008

%% Set path and global variables
global SRB
SRB = true;

%% Setup a matrix
randn('state',2008);
rand('state',2008);

n = 1000; r = 10;
M = randn(n,r)*randn(r,n);

df = r*(2*n-r);
oversampling = 5;  m = 5*df;

Omega = randsample(n^2,m);
data = M(Omega);

%% Set parameters and solve

p  = m/n^2;  delta = 1.2/p;
maxiter = 500;
tol = 1e-4;

%% Approximate minimum nuclear norm solution by SVT algorithm

tic
[U,S,V,numiter] = SVT(n,Omega,data,delta,maxiter,tol);
toc

%% Show results

X = U*S*V';

disp(sprintf('The relative error on Omega is: %d ', norm(data-X(Omega))/norm(data)))
disp(sprintf('The relative recovery error is: %d ', norm(M-X,'fro')/norm(M,'fro')))
disp(sprintf('The relative recovery in the spectral norm is: %d ', norm(M-X)/norm(M)))


### SVT()

function [U,Sigma,V,numiter]  = SVT(n,Omega,b,delta,maxiter,tol)
%
% Finds the minimum of tau ||X||_* + .5 || X ||_F^2
%
% subject to P_Omega(X) = P_Omega(M)
%
% using linear Bregman iterations
%
% Usage:  [U,S,V,numiter]  = SVT(n,Omega,b,delta,maxiter,opts)
%
% Inputs:
%
% n - size of the matrix X assumed n by n
%
% Omega - set of observed entries
%
% b - data vector of the form M(Omega)
%
% delta - step size
%
% maxiter - maximum number of iterations
%
% Outputs: matrix X stored in SVD format X = U*diag(S)*V'
%
% U - nxr left singular vectors
%
% S - rx1 singular values
%
% V - nxr right singular vectors
%
% numiter - number of iterations to achieve convergence

% Description:
% Reference:
%
%    Cai, Candes and Shen
%    A singular value thresholding algorithm for matrix completion
%    Submitted for publication, October 2008.
%
% Written by: Emmanuel Candes
% Email: emmanuel@acm.caltech.edu
% Created: October 2008

m = length(Omega); [temp,indx] = sort(Omega);
tau = 5*n; incre = 5;

[i, j] = ind2sub([n,n], Omega);
ProjM = sparse(i,j,b,n,n,m);

normProjM = normest(ProjM);
k0 = ceil(tau/(delta*normProjM));

normb = norm(b);

y = k0*delta*b;
Y = sparse(i,j,y,n,n,m);
r = 0;

fprintf('\nIteration:   ');
for k = 1:maxiter,
fprintf('\b\b\b%3d',k);
s = r + 1;

OK = 0;
while ~OK
[U,Sigma,V] = lansvd(Y,s,'L');
OK = (Sigma(s,s) <= tau);
s = s + incre;
end

sigma = diag(Sigma); r = sum(sigma > tau);
U = U(:,1:r); V = V(:,1:r); sigma = sigma(1:r) - tau; Sigma = diag(sigma);

x = XonOmega(U*diag(sigma),V,Omega);

if (norm(x-b)/normb < tol)
break
end

y = y + delta*(b-x);
end

fprintf('\n');
numiter = k;


### bdsqr()

function [sigma,bnd] = bdsqr(alpha,beta)

% BDSQR: Compute the singular values and bottom element of
%        the left singular vectors of a (k+1) x k lower bidiagonal
%        matrix with diagonal alpha(1:k) and lower bidiagonal beta(1:k),
%        where length(alpha) = length(beta) = k.
%
% [sigma,bnd] = bdsqr(alpha,beta)
%
% Input parameters:
%   alpha(1:k)   : Diagonal elements.
%   beta(1:k)    : Sub-diagonal elements.
% Output parameters:
%   sigma(1:k)  : Computed eigenvalues.
%   bnd(1:k)    : Bottom elements in left singular vectors.

% Below is a very slow replacement for the BDSQR MEX-file.

%warning('PROPACK:NotUsingMex','Using slow matlab code for bdsqr.')
k = length(alpha);
if min(size(alpha)') ~= 1  | min(size(beta)') ~= 1
error('alpha and beta must be vectors')
elseif length(beta) ~= k
error('alpha and beta must have the same lenght')
end
B = spdiags([alpha(:),beta(:)],[0,-1],k+1,k);
[U,S,V] = svd(full(B),0);
sigma = diag(S);
bnd = U(end,1:k)';


### 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


### 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.
%

% 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


### 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


The 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:
* 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(" 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("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 */
}
}

}