% Matlab code RankDeficient.m % For "Applied Numerical Linear Algebra", Figures 3.4 and 3.5 % Written by James Demmel, Feb 4, 1996 % Modified, Jun 2, 1997 % % Illustrate solution of Rank Deficient Least Squares Problem using % Truncated SVD % % Given any m-by-n input matrix A and use supplied tolerance tol, % we produce a matrix Al with a range of singular values from % small to large. Then we compute a rank-deficient Ard by setting % the singular values of Al less than tol to zero. % Then we generate a rank-deficient least % squares problem whose solution we know exactly as follows: % % 1) Compute the SVD Ard = U*S*V', where U and V have rnk columns, % where rnk is the rank of Ard. % 2) Pick a random xrd of dimension rnk, let x=V*xrd, and % let b = Ard*x (= U*S*xrd). Then x is the exactly, minimum norm % solution of the rank-deficient least squares problems % min(norm(A*y-b,2)). % 3) Add random noise dA of varying norm to A, and solve the perturbed % least squares problem using the truncated SVD, getting perturbed % solution y. % 4) Plot norm(y) and norm(y-x) versus norm(dA). % % Inputs: % dimensions m (number of rows) and n (number of columns), m >= n % type = 1 to computed matrix with floor(n/2) zero singular values; % for figure 3.4 % = 2 to computed matrix with singular values roughtly geometrically % distributed between 1e-15 and 1; for figure 3.5 % tol = tolerance, below which to set small singular value to 0 % % % Here are two ways to compute Al. One should be commented out. A = randn(m,n); if (type == 1) % Compute Al by zeroing out the last half of the columns of A: Al = randn(m,m) * ... [diag([ones(floor(n/2),1);zeros(ceil(n/2),1)]);zeros(m-n,n)] * ... randn(n,n); end if (type == 2) % Compute Al by multiplying columns by geometric sequence from 1e-16 to 1. % As a result Al will have singular values which form an approximate % geometric sequence in the same range. Al = A * diag(exp(log(10)*(0:n-1)*(-16/n))); end % [Uorig,Sorig,Vorig]=svd(Al,0); diag(Sorig); rnk = length(find(diag(Sorig) >= tol)); disp(['rank(A) = ',int2str(rnk)]) U=Uorig(:,1:rnk); S=Sorig(1:rnk,1:rnk); V=Vorig(:,1:rnk); Ard = U*S*V'; xrd = randn(rnk,1); x = V * xrd; b = Ard*x; xnorm = norm(xrd); nrm = sort([(10*ones(1,17)).^(-16:0),tol*(.1:.1:.9),tol*(.9:.01:1),... tol*(1:.1:1.9),tol*(2:10)]); rnkdsav=[]; ynorm=[]; xmynorm=[]; for nn=nrm; dA = randn(m,n); dA = dA/norm(dA); dA = nn*dA; ArdpdA = Ard + dA; [Ud,Sd,Vd] = svd(ArdpdA); rnkd = length(find(diag(Sd) >= tol)); rnkdsav=[rnkdsav;rnkd]; Udd=Ud(:,1:rnkd); Sdd=Sd(1:rnkd,1:rnkd); Vdd=Vd(:,1:rnkd); y = Vdd*(Sdd\(Udd'*b)); ynorm = [ynorm;norm(y)/xnorm]; xmynorm = [xmynorm;norm(x-y)/xnorm]; end figure(1), hold off, clg % subplot(2,1,1) handl=semilogx(nrm',rnkdsav,'b'); set(handl,'LineWidth',2) hold on dS = min(max([diag(S)],1e-16),1); handl=semilogx(dS,ones(size(dS)),'rx');set(handl,'MarkerSize',10) axis([1e-16 1 0 min(m,n)]), grid title(['Rank of perturbed A, original singular values are x''s, tol=',num2str(tol)]) xlabel('Norm of perturbation') % subplot(2,1,2) handl=loglog(nrm',xmynorm,'b');set(handl,'LineWidth',2) hold on handl=semilogx(dS,ones(size(dS))*1e-15,'rx');set(handl,'MarkerSize',10) axis([1e-16 1 1e-16 1]), grid title('Norm(solution-perturbed solution)/norm(solution)') xlabel('Norm of perturbation')