## Analysis Operator

The analysis operator W is p -by- n iid Gaussian and thus a nearly uniform and nearly tight frame

```n=1e3; % signal dimension
overcompleteness = 2.0; % ratio of analysis dimension to signal dimension
p=round(n*overcompleteness); % analysis dimension
W = randn(p,n) / sqrt(n); % analysis operator with approximately unit-norm rows
```

## Measurement Matrix

The measurement operator M is m -by- n iid Gaussian

```delta = 0.5 + (overcompleteness-1)*0.25; % sampling ratio
m=round(delta*n); % number of measurements
M = randn(m,n)/sqrt(n);  % measurement operator with approximately unit-norm rows
```

## Cosparse Signal

The signal x is l-cosparse on W, i.e., W * x has l zeros, located uniformly at random

```rho = .3;  % rho=(n-l)/m: an uncertainty parameter between 0 (easy) and 1 (difficult)
l=n - round(rho*m);  % cosparsity: the number of zeros in W*x
cosup = randperm(p,l); % cosupport: a random selection of l rows
P = null(W(cosup,:)); % the rows of P form a basis for the nullspace of the cosupported rows of W
x = P*P'*randn(n,1); % a random cosparse signal created by projecting onto the nullspace
x = x/norm(x);
z = M*x; % true noiseless measurements
v = W*x; % true analysis outputs
```

## Noisy Compressive Linear Measurements

The measurements are constructed as y = M x + w with AWGN w

```SNRdB = 80; % measurement signal-to-noise ratio in dB
wvar = (norm(z)^2/m) * 10^(-SNRdB/10); % AWGN variance
y = z + randn(m,1)*sqrt(wvar); % noisy measurements
```

## Configuring and Running GrAMPA

GrAMPA is based on the generalized approximate message passing (GAMP) algorithm, a powerful Bayesian framework for large-scale generalized-linear inference exploiting a separable signal prior, a separable likelihood, and a sufficiently large, dense, linear transform.

To get started using GAMP, download the latest GAMPmatlab release. You can find some introductory examples described on the wikia site as well as intructions on how to get the latest development code.

The workhorse of GAMPmatlab is the gampEst function. To implement GrAMPA, we need to specify its four input arguments: the signal prior GrampaEstimIn, the measurement likelihood GrampaEstimOut, the linear transform GrampaLinTrans, and the algorithmic options GrampaOptions.

```% When using GAMP to implement GrAMPA, we recommend the following options.
GrampaOptions.xvar0 = (norm(y)^2-wvar*m)/norm(M,'fro')^2; % guess at signal variance
GrampaOptions.step = 0.1; % initial stepsize
GrampaOptions.stepMax = 0.5; % maximum stepsize
GrampaOptions.stepIncr = 1.01; % stepsize increase rate
GrampaOptions.tol = 1e-6; % stopping tolerance
GrampaOptions.nit = 1000; % maximum number of iterations
GrampaOptions.varNorm = true; % turn on internal normalization
GrampaOptions.zvarToPvarMax = inf; % analysis variance clipping ratio

% For this example, we assume nothing about the signal, and so we use the trivial
% signal prior.  If we knew the signal was bounded on an interval, we would
GrampaEstimIn = NullEstimIn(0,1);

% GrampaLinTrans is constructed by "stacking" the measurement operator
% on top of the analysis operator.  Since our operators are explicit matrices,
% we can simply stack the matrices, but in general we'd use *LinTransConcat()*.
% Note that, when stacking, the operators must have the same (average) row norms,
% which in this example is ensured by the way that we built the operators.
GrampaLinTrans = [M;W];
iMeas = (1:m)';	% measurement output indices
iAnal = (1:p)'+m; % analysis output indices

% The likelihood is then constructed by stacking the measurement likelihood
% on top of the analysis prior.  In this example, the measurement likelihood is AWGN
% with variance *wvar*.  For the analysis prior, we use ``SNIPEstim,'' which depends on
% a parameter *omega*.  Below, we guess the value *omega*=3.
MeasEstimOut = AwgnEstimOut(y,wvar);
omega = 3;
AnaEstimOut = SNIPEstim(omega);
GrampaEstimOut = EstimOutConcat({MeasEstimOut;AnaEstimOut},[m,p]);

% Finally, we call *gampEst* using the objects contructed above, and measure the
% resulting normalized output SNR.
estFin1 = gampEst(GrampaEstimIn,GrampaEstimOut,GrampaLinTrans,GrampaOptions);
snr_gamp1 = 20*log10( norm(x)/norm(x-estFin1.xhat) );
fprintf('With guessed omega=%.1f, GrAMPA gives NMSE=%.2fdB\n',AnaEstimOut.omega,-snr_gamp1);
```
```With guessed omega=3.0, GrAMPA gives NMSE=-81.87dB
```

## An Oracle Bound

The "cosupport genie" magically knows the true cosupport and then finds the signal vector that best matches the measurements (in a least-squares sense) while yielding that cosupport. To implement this, we can use the previously constructed matrix P, whose rows must now be orthogonal to the genie estimate of x:

```xhat_genie = P * ( [M*P] \ y ); % solve LLS over that subspace
gMSEx = norm(x-xhat_genie)^2;
snr_genie = 10*log10( norm(x)^2/gMSEx );
fprintf('Meanwhile, the cosupport genie gives NMSE=%.2fdB\n',-snr_genie)
```
```Meanwhile, the cosupport genie gives NMSE=-83.26dB
```

## SURE Tuning of SNIPEstim:

Now let's see if we can improve performance by (genie) tuning the SNIPEstim parameter omega using Stein's Unbiased Risk Estimate (SURE) computed from the GAMP analysis outputs, which should be AWGN corrupted versions of true analysis outputs with known noise variance. We will also compute the true signal MSE to verify that SURE is working as expected.

```% First, let's compute the SURE estimates over a grid of omegas
fprintf('Tuning *omega* by computing SURE risk over a grid...\n')
omegaGrid = linspace(0,8,9)'; % typical range for SNIPEstim *omega*
SURE = zeros(size(omegaGrid));
MSEx = zeros(size(omegaGrid));
for i=1:length(omegaGrid)
% run GrAMPA at an omega on the grid
AnaEstimOut = SNIPEstim(omegaGrid(i));
GrampaEstimOut = EstimOutConcat({MeasEstimOut;AnaEstimOut},[m,p]);
estFin2 = gampEst(GrampaEstimIn, GrampaEstimOut, GrampaLinTrans, GrampaOptions);
phat = estFin2.phat;
pvar = estFin2.pvar;
zhat = estFin2.zhat;
zvar = estFin2.zvar;

% compute SURE approximation to analysis-coefficient MSE at hypothesized omega
SURE(i) = sum(abs(-pvar + abs(zhat - phat).^2 + 2*zvar));

% compute true signal-MSE at hypothesized omega (requiring knowledge of signal)
xhat = estFin2.xhat;
MSEx(i) = norm(xhat - x)^2;
MSEz(i) = norm(zhat - [z;v])^2;
end
[bestSURE,iBest] = min(SURE);
bestOmega = omegaGrid(iBest);
snr_gamp2 = 10*log10( norm(x)^2/MSEx(iBest) );
disp('')
fprintf('With SURE-tuned omega=%.1f, GrAMPA gives NMSE=%.2fdB\n',bestOmega,-snr_gamp2);

% Finally, let's plot GrAMPA's MSE versus *omega* along with the SURE-estimated risk
figure(1); clf;
plot(bestOmega,10*log10(bestSURE),'bo',...
omegaGrid,10*log10(SURE),'b+-',...
omegaGrid,10*log10(MSEz),'b.-',...
omegaGrid,10*log10(MSEx),'r.-',...
omegaGrid,10*log10(gMSEx)*ones(size(omegaGrid)),'r--');
axe1=axis; axis('tight'); axe2=axis; axis([axe1(1:2),axe2(3)-5,min(axe2(4),axe2(3)+50)])
grid on;
legend('SURE tuning','SURE on z','MSE on z','MSE on x','genie MSE on x','Location','Best')
ylabel('dB')
xlabel('omega')
```
```Tuning *omega* by computing SURE risk over a grid...
With SURE-tuned omega=8.0, GrAMPA gives NMSE=-83.30dB
```