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