A MATLAB comparison of ssAMP to recent algorithms for 1-D piecewise-constant recovery

path(path, './solver/flsa');

path(path, './solver/ssAMP_denoiser');

path(path, './solver/TVAL3_beta2.4');

path(path, './solver/GAMP_codes');

path(path, './solver');

path(path, './etc_tool');

path(path,genpath(pwd));

alpha=0.25;% undersampling ratio M/N

N=25*25; % the signal vector dimension

M=round (N*alpha); % the measurement vector dimension

% Additive noise variance

Delta=0.01;

1) Gaussian , 2) Bernoulli , and 3) Ternary matrices . H=randn(M,N)/sqrt(M);% 1) Gaussian matrices

%H=(2*(rand(M,N)<0.5)-1)/sqrt(M);% 2) Bernoulli matrices h \in {1,-1}

%H=randint2(M,N,[-1,1])/sqrt(dilution*M);% 3) Ternary matrices h \in {0,1,-1

% These parameters are applied to the spike-and-slab potential function in the ssAMP algorithm

q =0.05;% signal piecewise-constancy

sigma1 = 0.01;% signal level

sigma0 = 1;% variance for difference between neigboring values

% Piecewise-constant signal generation

X=piecewiseGaussian_Gibbsgen(N,q,sigma0);% Gaussian signal: X_i - X_i-1 \in Gaussian

%X=piecewiseBinary_Gibbsgen(N,q,sigma0);% Binary signal: X_i - X_i-1 \in {0,1}

%X=piecewiseTernary_Gibbsgen(N,q,sigma0);% Ternary signal: X_i - X_i-1 \in {0,1,-1}

disp('%------------------------------------------------------------------------------------------%');

disp('<Experiment codition>')

disp(sprintf('Piecewise-constancy, K/M = %8.3f',q/alpha));

disp(sprintf('Undersampling ratio, M/N = %8.3f',M/N));

disp(sprintf('Noise variance, delta = %8.3f',Delta));

disp('%------------------------------------------------------------------------------------------%');

disp('<The number of iteration expended>')

is generated.

noise=sqrt(Delta)*randn(M,1); % Additive measurement noise generation

y=H*X+noise;

% Common parameters

maxiter =200; iteration_stopping_tol=1E-6;

%--------------- option setting for EFLA solver ------------------------%

rho=0; % the regularization parameter;it is a ratio between (0,1), if .rFlag=1

opts_EFLA=[];

opts_EFLA.init=2; % starting from a zero point

opts_EFLA.tFlag=5; % run .maxIter iterations

opts_EFLA.maxIter=maxiter; % maximum number of iterations

opts_EFLA.nFlag=0; % without normalization

opts_EFLA.rFlag=1; % the input parameter 'rho' is a ratio in (0, 1)

opts_EFLA.fusedPenalty=0.1;

opts_EFLA.lFlag=0; % line search off

%-------------------- option setting for TVAL3 ---------------------------%

% The following options are recommended by the website "http://www.caam.rice.edu/~optimization/L1/TVAL3/"

opts_TVAL3.mu = 2^7;

opts_TVAL3.beta = 2^4;

opts_TVAL3.mu0 = opts_TVAL3.mu*0.1;

opts_TVAL3.beta0 = opts_TVAL3.beta*0.1;

opts_TVAL3.maxit=maxiter;

opts_TVAL3.maxcnt = 10;

opts_TVAL3.tol_inn = 1e-3;

opts_TVAL3.tol = iteration_stopping_tol;

opts_TVAL3.init=0;

opts_TVAL3.TVnorm = 1; % Anisortopic setting

%---------------------- Configuring GrAMPA ----------------------------%

D=zeros(N,N);% co-sparse matrix for TV penalty for GrAMPA -

for i=2:N-1

D(i,:)=[zeros(1,i-1) -1 1 zeros(1, N-i-1)];

end

D(1,:)=[-1 1 zeros(1, N-2)];

% The following options are recommended by the website "http://www2.ece.ohio-state.edu/~schniter/GrAMPA/html/grampa_example1.html"

GrampaOptions = GampOpt('legacyOut',false,'uniformVariance',true,'adaptStep',false);

GrampaOptions.xvar0 = (norm(y)^2-Delta*M)/norm(H,'fro')^2; % guess at signal variance

GrampaOptions.step = 0.1; % initial stepsize

GrampaOptions.stepMax = 0.5; % maximum stepsize

GrampaOptions.stepIncr = 5; % stepsize increase rate

GrampaOptions.tol = iteration_stopping_tol; % stopping tolerance

GrampaOptions.nit = maxiter; % maximum number of iterations

GrampaOptions.varNorm = true; % turn on internal normalization

GrampaOptions.zvarToPvarMax = inf; % analysis variance clipping ratio

GrampaEstimIn = NullEstimIn(0,1);

GrampaLinTrans = [H;D];

iMeas = (1:M)'; % measurement output indices

iAnal = (1:N)'+M; % analysis output indices

MeasEstimOut = AwgnEstimOut(y,Delta);

omega = 3;

AnaEstimOut = SNIPEstim(omega);

GrampaEstimOut = EstimOutConcat({MeasEstimOut;AnaEstimOut},[M,N]);

<Signal recovery>

% ssAMP solving [1] (Proposed)

tstart=tic;

[estX_fusedAMP,varX,end_iter,theta]= solve_FusedAMP_mat_ver5(H,y,sigma0,q,Delta,maxiter);

telapsed_fusedMP=toc(tstart);

% EFLA solving [2]

tstart=tic;

[estX_fusedlasso, funVal1, ValueL1]= fusedLeastR(H, y, rho, opts_EFLA);

telapsed_fusedlasso=toc(tstart);

% TV-AMP solving [3]

tstart=tic;

[estX_TVAMP,tau]=solve_TVAMP_mat_ver4(H,y,maxiter,X);

telapsed_TVAMP=toc(tstart);

% TVAL3 solving [4]

tstart=tic;

estX_TVAL3=TVAL3(H,y,sqrt(N),sqrt(N),opts_TVAL3);

telapsed_TVAL3=toc(tstart);

estX_TVAL3=estX_TVAL3(:);

% chambolle-pock TV [5]

tstart=tic;

estX_CP= solve_chambolle_pock_TV_ver2(H,y,H'*y, opts_EFLA.fusedPenalty,iteration_stopping_tol,maxiter);

telapsed_CP=toc(tstart);

% GrAMPA estimate [6]

tstart=tic;

estFin1 = gampEst(GrampaEstimIn,GrampaEstimOut,GrampaLinTrans,GrampaOptions);

telapsed_GrAMPA=toc(tstart);

estX_GrAMPA=estFin1.xhat;

disp(sprintf('Number of GrAMPA iterations = %d',estFin1.nit));

< Normalized MSE calculation and display>

%----------------- Normalized MSE calculation--------------------%

normX_sqr=norm(X)^2;-----

MSE_fusedMP= norm(estX_fusedAMP-X)^2/normX_sqr;

MSE_fusedlasso= norm(estX_fusedlasso-X)^2/normX_sqr;

MSE_TVAMP= norm(estX_TVAMP-X)^2/normX_sqr;

MSE_TVAL3=norm(estX_TVAL3-X)^2/normX_sqr;

MSE_CP=norm(estX_CP-X)^2/normX_sqr;

MSE_GrAMPA=norm(estX_GrAMPA-X)^2/normX_sqr;

%--------------------------- Display -------------------------------%

disp('%------------------------------------------------------------------------------------------%');

disp('<Recovery result>')

disp(sprintf('ssAMP: Nomalized MSE = %8.7f',MSE_fusedMP));

disp(sprintf('EFLA: Nomalized MSE = %8.7f',MSE_fusedlasso));

disp(sprintf('TVAMP: Nomalized MSE = %8.7f',MSE_TVAMP));

disp(sprintf('TVAL3: Nomalized MSE = %8.7f',MSE_TVAL3));

disp(sprintf('Chambolle-Pock TV:Nomalized MSE = %8.7f',MSE_CP));

disp(sprintf('GrAMPA: Nomalized MSE = %8.7f',MSE_GrAMPA));

disp('%------------------------------------------------------------------------------------------%');

disp(sprintf('Running Time of ssAMP = %8.7f sec',telapsed_fusedMP));

disp(sprintf('Running Time of EFLA = %8.7f sec',telapsed_fusedlasso));

disp(sprintf('Running Time of TVAMP= %8.7f sec',telapsed_TVAMP));

disp(sprintf('Running Time of TVAL3= %8.7f sec',telapsed_TVAL3));

disp(sprintf('Running Time of Chambolle-Pock TV= %8.7f sec',telapsed_CP));

disp(sprintf('Running Time of GrAMPA= %8.7f sec',telapsed_GrAMPA));

disp('%------------------------------------------------------------------------------------------%');

<Exemplary results>

- Case 1 (Noisy and insufficient measurements): , Gaussian signalsssAMP: Nomalized MSE = 0.2809199

EFLA: Nomalized MSE = 0.2888708

TVAMP: Nomalized MSE = 0.2803345

TVAL3: Nomalized MSE = 0.9737418

Chambolle-Pock TV:Nomalized MSE = 0.3674162

GrAMPA: Nomalized MSE = Inf

- Case 2 (Noisy and sufficient measurements): , Gaussian signalsssAMP: Nomalized MSE = 0.0037972

EFLA: Nomalized MSE = 0.0506075

TVAMP: Nomalized MSE = 0.0288965

TVAL3: Nomalized MSE = 0.6756543

Chambolle-Pock TV:Nomalized MSE = 0.0669542

GrAMPA: Nomalized MSE = 0.0043243

- Case 3 (Clean and highly sufficient measurements):

, Gaussian signals

ssAMP: Nomalized MSE = 0.0000001

EFLA: Nomalized MSE = 0.0030774

TVAMP: Nomalized MSE = 0.0000035

TVAL3: Nomalized MSE = 0.1173732

Chambolle-Pock TV:Nomalized MSE = 0.0004288

GrAMPA: Nomalized MSE = 0.0000000