Example: Estimation of a Gaussian Vector
The code for the example is in
To run the code, simply
cd to that directory and run
simpleAWGN at the MATLAB prompty. You will see a plot of a sparse vector and its estimate. In addition, the mean squared error (MSE) for both the least-squares and GAMP methods will be printed.
In this first problem, we will show how to use the
gampEst function to estimate a Gaussian random vector from random linear measurements of the form
where is a Gaussian random matrix and is Gaussian noise. The solution to this problem is, of course, standard and be performed easily without GAMP: For example, if the components of and are i.i.d., Gaussian zero-mean with variances and , respectively, then the MMSE estimate is given by the least squares solution:
Nevertheless, we will show how to get an approximate solution to this problem using GAMP simply to illustrate the syntax and usage of the MATLAB command.
MATLAB Program Description
Generation of the Data
At the top of the program
simpleAWGN.m, you will see the lines:
% Set path to the main directory addpath('../../main/'); % Parameters nx = 1000; % Number of input components (dimension of x) nz = 2000; % Number of output components (dimension of y) snr = 20; % SNR in dB.
The first line adds the necessary directories to the MATLAB path. You can also set the MATLAB path permanently if you intend to use GAMP significantly. The remaining lines set the parameters of the problem. Feel free to change these to experiment.
The next lines generate a random Gaussian vector :
% Create a random Gaussian vector xmean0 = 0; xvar0 = 1; x = normrnd(xmean0, sqrt(xvar0),nx,1);
The matrix and the noise vector are generated with Gaussian
i.i.d. components, with the noise scaled according to the parameter
% Create a random measurement matrix A = (1/sqrt(nx))*randn(nz,nx); % Compute the noise level based on the specified SNR. Since the components % of A have power 1/nx, the E|y(i)|^2 = E|x(j)|^2 = xvar. wvar = 10^(-0.1*snr)*xvar0; % Generate the noise w = normrnd(0, sqrt(wvar), nz, 1); y = A*x + w;
Solution via GAMP
The main GAMP command to perform the estimation of from the measurement vector appears near the bottom of the file:
[xhat, xvar] = gampEst(inputEst, outputEst, A, opt);
gampEst function takes as input arguments:
inputEstis a class of type
EstimInthat describes the prior on the components of the unknown vector ;
Ais a matrix for the transform for the relation . This can also be specified by a
LinTransclass as we will see in subsequent examples.
outputEstis class of type
EstimOutdescribing the transition probabilities for the measurements ; and
optis a set of options.
xvar are the estimates of the mean and variance
of the vector . You can, in fact, get other outputs to estimate the full marginal posteriors
of each the variables -- we will illustrate this in a later example.
The idea of specifying the prior and measurement channel
via the MATLAB classes
outputEst may be confusing at first.
However, as we will see in this sequence of examples, the use of classes
allows the algorithm to be very general, modular and capable of handling many different
types of distributions. If you are not familiar with MATLAB classes or object-oriented programming, you can read more on the
MATLAB object-oriented webpage.
gampmatlab package contains pre-built
classes for many common input and output distributions.
For most problems, you can just use one of these
pre-built classes and not have to build one from scratch.
In this particular example, since the input is a Gaussian random vector,
inputEst is created using the pre-built
class with the line:
inputEst = AwgnEstimIn(xmean0, xvar0);
The output estimator class is also created via a pre-built class. Specifically, the line
outputEst = AwgnEstimOut(y, wvar);
outputEst class using the function
AwgnEstimOut for an AWGN measurement
channel. Note that the observation vector
y is passed to the estimator class -- not the
With the estimator classes defined, the program then calls
gampEst which returns
an estimate vector
The GAMP estimate is compared to the MMSE estimator. For this problem, the MMSE estimate is given by the solution in equation (1) above to the least-squares problem. This estimate is computed in the MATLAB code with the lines:
xhatLS = xvar0*(A'*A*xvar0 + wvar*eye(nx))\(A'*y);
The program concludes by plotting
x, the GAMP estimate
xhat and LS estimate
xhatLS. The MSE for the GAMP and LS estimates are also printed. You should see that GAMP estimate obtains an error close to the LS estimate.