% Steve Penny % AMSC 600 % 10/18/07 % HW 2_2c west0479 = load('west0479.mat'); C = west0479.west0479; n = length(C); k = 100; %%%%%%%%%%%%%% % Initialize % %%%%%%%%%%%%%% Z = zeros(n,k); W = zeros(n,k); B = zeros(n,k); Z(:,1) = ones(length(Z),1)/norm(ones(length(Z),1),2); % iterate algorithm from 2b % to generate sigular values of C for j=1:k % Part 1 if j==1 W(:,j) = C*Z(:,j); else W(:,j) = C*Z(:,j) - W(:,j-1)*B(j-1,j); end B(j,j) = norm(W(:,j),2); W(:,j) = W(:,j)/B(j,j); % Part 2 Z(:,j+1) = C'*W(:,j) - Z(:,j)*B(j,j); B(j,j+1) = norm(Z(:,j+1),2); Z(:,j+1) = Z(:,j+1)/B(j,j+1); end Sig_C_approx = svd(B); Sig_C = svd(full(C)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Display Results: %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Print some of the singular values Sig_C(1:10) Sig_C(10:20) Sig_C(20:30) Sig_C(end-10:end) Sig_C_approx(1:10) Sig_C_approx(10:20) Sig_C_approx(20:30) Sig_C_approx(end-10:end) % plot the singular values figure(1) plot(Sig_C,Sig_C,'r.') hold on plot(Sig_C_approx,Sig_C_approx,'bo','MarkerSize',10.0) hold off ylabel('singular values') xlabel('singular values') set(gcf,'color',[1 1 1]) set(gca,'box','off') % plot the upper singular values figure(2) upper_cut = 3*10^5; plot(Sig_C(find(Sig_C >= upper_cut)),Sig_C(find(Sig_C >= upper_cut)),'r.') hold on plot(Sig_C_approx(find(Sig_C_approx >= upper_cut)),Sig_C_approx(find(Sig_C_approx >= upper_cut)),'bo','MarkerSize',10.0) hold off ylabel('singular values') xlabel('singular values') set(gcf,'color',[1 1 1]) set(gca,'box','off') % plot the lower singular values figure(3) lower_cut = 0.5*10^5; plot(Sig_C(find(Sig_C <= lower_cut)),Sig_C(find(Sig_C <= lower_cut)),'r.') hold on plot(Sig_C_approx(find(Sig_C_approx <= lower_cut)),Sig_C_approx(find(Sig_C_approx <= lower_cut)),'bo','MarkerSize',10.0) hold off ylabel('singular values') xlabel('singular values') set(gcf,'color',[1 1 1]) set(gca,'box','off') return