www.pudn.com > SURE_SVT_code.zip > demo_sure_vs_montecarlo.m, change:2012-10-13,size:4830b


% 
%   DEMO #1 
%   UNBIASED RISK ESTIMATES FOR SINGULAR VALUE THRESHOLDING 
%   E. J. Candes, C. A. Sing-Long and J. D. Trzasko 
% 
%   DESCRIPTION:    This script generates Figure 1, which compares Stein's 
%                   Unbiased Risk Estimate (SURE) against an estimation of 
%                   the risk using Monte Carlo, when using Singular Value 
%                   Thresholding, and a Gaussian additive noise model. 
% 
%   Oct. 2012. 
 
close all 
clear all 
drawnow 
clc 
 
%   set the state of the random number generator 
%   (this is the preferred method for newer versions of Matlab) 
rng(1234) 
 
%   experiment parameters 
M           =   200;            %   number of rows 
N           =   500;            %   number of columns 
Ns          =   50;             %   number of Monte Carlo samples 
Nl          =   50;             %   number of values for lambda 
lambda_max  =   50;             %   maximum value for lambda 
SNR         =   [0.5; 1; 2; 4];     %   signal-to-noise ratio 
 
%   matrices (mean) 
%       rank(X1) = 200 (full-rank) 
X1          =   normrnd(0, 1, [M N]); 
%       rank(X1) = 100 (50%) 
X2          =   normrnd(0, 1, [M N]); 
[U S V]     =   svd(X2); 
idx         =   round(0.5*M); 
S((idx+1):M, (idx+1):M)     =   zeros(M - idx); 
X2          =   U*S*V'; 
%       rank(X1) = 5 (5%) 
X3          =   normrnd(0, 1, [M N]); 
[U S V]     =   svd(X3); 
idx         =   round(0.05*M); 
S((idx+1):M, (idx+1):M)     =   zeros(M - idx); 
X3          =   U*S*V'; 
%       \sigma_i = \sqrt{M}/(1 + exp{(i - M/2)/20}), i = 1, ..., 200 
X4          =   normrnd(0, 1, [M N]); 
[U S V]     =   svd(X4); 
S(1:M, 1:M) =   diag(sqrt(M)./(1 + exp(((1:M)' - 0.5*M)/20))); 
X4          =   U*S*V'; 
 
%       normalize 
X1          =   X1/norm(X1(:)); 
X2          =   X2/norm(X2(:)); 
X3          =   X3/norm(X3(:)); 
X4          =   X4/norm(X4(:)); 
 
%   variables 
X0      =   {X1 X2 X3 X4};          %   matrices to consider 
lambda  =   cell(4, 4);             %   thresholds 
tau_w   =   zeros(4, 4);            %   noise standard-deviation 
MCR     =   zeros(Nl, 4, 4);        %   Monte Carlo estimated risk 
MCS     =   zeros(Nl, 4, 4, Ns);    %   Monte Carlo samples 
SURE    =   zeros(Nl, 4, 4);        %   SURE estimate 
for Ik = 1:4,                   	%   loop over matrices 
    X           =   X0{Ik};         %   matrix 
    for In = 1:4,                   %   loop over SNR 
        tau_w(Ik, In)   =   1/(SNR(In)*sqrt(M*N));                      %   noise standard-deviation 
        lambda{Ik, In}  =   linspace(0, lambda_max*tau_w(Ik, In), Nl)'; %   thresholds 
        for Il = 1:Nl,                      %   loop over lambda 
            for Is = 1:Ns,                  %   Monte Carlo sample 
                Y       =   X + tau_w(Ik, In)*normrnd(0, 1, [M N]); 
                %   singular value thresholding 
                [U, Sy, V]  =   svd(Y); 
                Sy          =   max(0, Sy - lambda{Ik, In}(Il)); 
                SVT_Y       =   U*Sy*V'; 
                clear U Sy V 
                 
                %   Monte Carlo estimate 
                MCS(Il, In, Ik, Is) =   sum((SVT_Y(:) - X(:)).^2); 
                MCR(Il, In, Ik)     =   MCR(Il, In, Ik) + MCS(Il, In, Ik, Is); 
                clear SVT_Y 
            end 
            %   Monte Carlo Estimate             
            MCR(Il, In, Ik)     =   MCR(Il, In, Ik)/Ns; 
             
            %   SURE 
            Y                   =   X + tau_w(Ik, In)*normrnd(0, 1, [M N]); 
            SURE(Il, In, Ik)    =   sure_svt(lambda{Ik, In}(Il), tau_w(Ik, In), Y); 
        end 
    end 
end 
 
save('data_sure_vs_montecarlo.mat', 'SURE', 'MCR', 'MCS', 'SNR', 'lambda', 'tau_w'); 
%load('data_sure_vs_montecarlo.mat', 'SURE', 'MCR', 'MCS', 'SNR', 'lambda', 'tau_w'); 
 
cmap    =   jet(4); 
cmap(2, :)  =   [0 1 0.5]; 
cmap(3, :)  =   [1 0.5 0]; 
lfsz    =   20; 
afsz    =   24; 
 
for In = 1:4, 
    ymax        =   0; 
    xmax        =   0; 
    figure('Renderer', 'painters', 'Position', [100 100 800 800]) 
    hold on 
    for Ik = 1:4, 
        plot(tau_w(Ik, In)*lambda{Ik, In}, MCR(:, In, Ik), ... 
             'Color', cmap(Ik, :), 'LineWidth', 2) 
        ymax    =   max([MCR(:, In, Ik); ymax]); 
        xmax    =   max([tau_w(Ik, In)*lambda{Ik, In}; xmax]); 
    end 
    for Ik = 1:4, 
        plot(tau_w(Ik, In)*lambda{Ik, In}(1:2:Nl), SURE(1:2:Nl, In, Ik), ... 
             'LineStyle', 'none', 'Marker', '+', 'MarkerSize', 15, 'Color', cmap(Ik, :)) 
    end 
    hold off 
    axis square; box on; grid on;  
    xlim([0 xmax]); ylim([-0.1*ymax 1.05*ymax]); 
    legend({'X_0^{(1)}', 'X_0^{(2)}', 'X_0^{(3)}', 'X_0^{(4)}'}, 'Location', 'SouthEast', 'FontSize', lfsz); 
    xlabel('\lambda\tau', 'FontSize', afsz); ylabel('average squared error', 'FontSize', afsz); 
end