Example: Estimation of a Gaussian Vector

From Gampmatlab

Jump to: navigation, search

Main -> Users' Guide and Examples -> Estimation of a Gaussian Vector.

Next: Estimating a sparse vector.



The code for the example is in gampmatlab/trunk/code/examples/basic/simpleAWGN.m.

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.

Problem Description

In this first problem, we will show how to use the gampEst function to estimate a Gaussian random vector LaTeX: \mathbf{x} from random linear measurements of the form

\quad \mathbf{y}=\mathbf{z} + \mathbf{w}, \quad \mathbf{z}= \mathbf{A}\mathbf{x},

where LaTeX: \mathbf{A} is a Gaussian random matrix and LaTeX: \mathbf{w} is Gaussian noise. The solution to this problem is, of course, standard and be performed easily without GAMP: For example, if the components of LaTeX: \mathbf{x} and LaTeX: \mathbf{w} are i.i.d., Gaussian zero-mean with variances LaTeX: \sigma^2_x and LaTeX: \sigma^2_w, respectively, then the MMSE estimate is given by the least squares solution:

\quad \widehat{\mathbf{x}}_{\rm LS} = \sigma^2_x\left( \sigma^2_x\mathbf{A}'\mathbf{A}
<pre>   + \sigma^2_wI\right)^{-1}\mathbf{A}'\mathbf{y}.   \quad (1)

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
% 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 LaTeX: \mathbf{x}:

% Create a random Gaussian vector
xmean0 = 0;
xvar0 = 1;
x = normrnd(xmean0, sqrt(xvar0),nx,1);

The matrix LaTeX: \mathbf{A} and the noise vector LaTeX: \mathbf{w} are generated with Gaussian i.i.d. components, with the noise scaled according to the parameter snr:

% 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 LaTeX: \mathbf{x} from the measurement vector LaTeX: \mathbf{y} appears near the bottom of the file:

[xhat, xvar] = gampEst(inputEst, outputEst, A, opt);

The gampEst function takes as input arguments:

  • inputEst is a class of type EstimIn that describes the prior LaTeX: p(x_j) on the components of the unknown vector LaTeX: \mathbf{x};
  • A is a matrix for the transform for the relation LaTeX: \mathbf{z}=\mathbf{A}\mathbf{x}. This can also be specified by a LinTrans class as we will see in subsequent examples.
  • outputEst is class of type EstimOut describing the transition probabilities for the measurements LaTeX: p(y_i|z_i); and
  • opt is a set of options.

The outputs xhat and xvar are the estimates of the mean and variance of the vector LaTeX: \mathbf{x}. 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 inputEst and 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. However, the 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, the class inputEst is created using the pre-built AwgnEstimIn 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);

creates 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 gampEst function.

With the estimator classes defined, the program then calls gampEst which returns an estimate vector xhat.

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.

Personal tools