In this demo, we illustrate how to use the vampyre package for a simple sparse linear inverse problem. The problem is to estimate a sparse vector $z$ from linear measurements of the form $y=Az+w$ where $w$ is Gaussian noise and $A$ is a known linear transform -- a basic problem in compressed sensing. By sparse, we mean that the vector $z$ has few non-zero values. Knowing that the vector is sparse can be used for improved reconstruction if an appropriate sparse reconstruction algorithm is used.
There are a large number of algorithms for sparse linear inverse problems. This demo uses the Vector Approximate Message Passing (VAMP) method, one of several methods that will be included in the vampyre package. In going through this demo, you will learn to:
vampyre packagevampyre package to perform the estimation for the linear inverse problemhist_list feature to track variables per iteration of the algorithm.First we need to import the vampyre package. Since python does not have relative imports, you need to add the path location for the vampyre package to the system path. In this case, we have specified the path use a relative path location, but you can change this depending on where vampyre is located.
In [1]:
import os
import sys
vp_path = os.path.abspath('../../')
if not vp_path in sys.path:
sys.path.append(vp_path)
import vampyre as vp
We will also load the other packages we will use in this demo. This could be done before the above import.
In [2]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
In [3]:
# Parameters
nz = 1000 # number of components of z
ny = 500 # number of measurements y
# Compute the shapes
zshape = (nz,) # Shape of z matrix
yshape = (ny,) # Shape of y matrix
Ashape = (ny,nz) # Shape of A matrix
To generate the synthetic data for this demo, we use the following simple probabilistic model. For the input $z$, we will use Bernouli-Gaussian (BG) distribution, a simple model in sparse signal processing. In the BG model, the components $z_i$ are i.i.d. where each component $z_i=0$ with probability $1-\rho$ and $z_i \sim {\mathcal N}(0,1)$ with probability $\rho$. The parameter $\rho$ is called the sparsity ratio and represents the average number of non-zero components. When $\rho$ is small, the vector $z$ is sparse. The components on which $z_i$ are non-zero are called the active components. We set the parameters below. We also set the SNR for the measurements.
In [4]:
sparse_rat = 0.1 # sparsity ratio
zmean1 = 0 # mean for the active components
zvar1 = 1 # variance for the active components
snr = 30 # SNR in dB
Using these parameters, we can generate random sparse z following this distribution with the following simple code.
In [5]:
# Generate the random input
z1 = np.random.normal(zmean1, np.sqrt(zvar1), zshape)
u = np.random.uniform(0, 1, zshape) < sparse_rat
z = z1*u
To illustrate the sparsity, we plot the vector z. We can see from this plot that the majority of the components of z are zero.
In [6]:
ind = np.array(range(nz))
plt.plot(ind,z)
Out[6]:
Now, we create a random transform $A$ and output $y_0 = Az$.
In [7]:
A = np.random.normal(0, 1/np.sqrt(nz), Ashape)
y0 = A.dot(z)
Finally, we add noise at the desired SNR
In [8]:
yvar = np.mean(np.abs(y0)**2)
wvar = yvar*np.power(10, -0.1*snr)
y = y0 + np.random.normal(0,np.sqrt(wvar), yshape)
Now that we have created the sparse data, we will use the vampyre package to recover z from y. In vampyre the methods to perform this estimation are called solvers. For this demo, we will use a simple solver called VAMP described in the paper:
Similar to most of the solvers in the vampyre package, the VAMP solver needs precise specifications of the probability distributions of z and y. The simplest way to use VAMP is to specify two densities:
Each of the densities are described by estimators.
We first describe the estimator for the prior $p(z)$. The vampyre package will eventually have a large number of estimators to describe various densities. In this simple demo, $p(z)$ is what is called a mixture distribution since $z$ is one distribution with probability $1-\rho$ and a second distribution with probability $\rho$. To describe this mixture distribution in the vampyre package, we need to first create estimator classes for each component distribution. The following code creates an estimator, est0, for a discrete distribution with a probability of 1 at a 0 and a second estimator, est1, for the Gaussian distribution with the active components.
In [9]:
est0 = vp.estim.DiscreteEst(0,1,zshape)
est1 = vp.estim.GaussEst(zmean1,zvar1,zshape)
We next use the vampyre class, MixEst, to describe a mixture of the two distributions. This is done by creating a list, est_list, of the estimators and an array pz with the probability of each component. The resulting estimator, est_in, is the estimator for the prior $z$, which is also the input to the transform $A$. We give this a name Input since it corresponds to the input. But, any naming is fine. Or, you can let vampyre give it a generic name.
In [10]:
est_list = [est0, est1]
pz = np.array([1-sparse_rat, sparse_rat])
est_in = vp.estim.MixEst(est_list, w=pz, name='Input')
Next, we describe the likelihood function, $p(y|z)$. Since $y=Az+w$, we can first use the MatrixLT class to define a linear transform operator Aop corresponding to the matrix A. Then, we use the LinEstim class to describe the likelihood $y=Az+w$.
In [11]:
Aop = vp.trans.MatrixLT(A,zshape)
est_out = vp.estim.LinEst(Aop,y,wvar,map_est=False, name='Output')
Finally, the VAMP method needs a message handler to describe how to perform the Gaussian message passing. This is a more advanced feature. For most applications, you can just use the simple message handler as follows.
In [12]:
msg_hdl = vp.estim.MsgHdlSimp(map_est=False, shape=zshape)
Having described the input and output estimators and the variance handler, we can now construct a VAMP solver. The construtor takes the input and output estimators, the variance handler and other parameters. The paramter nit is the number of iterations. This is fixed for now. Later, we will add auto-termination. The other parameter, hist_list is optional, and will be described momentarily.
In [13]:
nit = 20 # number of iterations
solver = vp.solver.Vamp(est_in,est_out,msg_hdl,\
hist_list=['zhat', 'zhatvar'],nit=nit)
We can print a summary of the model which indicates the dimensions and the estimators.
In [14]:
solver.summary()
We now run the solver by calling the solve() method. For a small problem like this, this should be close to instantaneous.
In [15]:
solver.solve()
The VAMP solver estimate is the field zhat. We plot one column of this (icol=0) and compare it to the corresponding column of the true matrix z. You should see a very good match.
In [16]:
zhat = solver.zhat
ind = np.array(range(nz))
plt.plot(ind,z)
plt.plot(ind,zhat)
plt.legend(['True', 'Estimate'])
Out[16]:
We can measure the normalized mean squared error as follows. The VAMP solver also produces an estimate of the MSE in the variable zhatvar. We can extract this variable to compute the predicted MSE. We see that the normalized MSE is indeed low and closely matches the predicted value from VAMP.
In [17]:
zerr = np.mean(np.abs(zhat-z)**2)
zhatvar = solver.zhatvar
zpow = np.mean(np.abs(z)**2)
mse_act = 10*np.log10(zerr/zpow)
mse_pred = 10*np.log10(zhatvar/zpow)
print("Normalized MSE (dB): actual {0:f} pred {1:f}".format(mse_act, mse_pred))
Finally, we can plot the actual and predicted MSE as a function of the iteration number. When solver was contructed, we passed an argument hist_list=['zhat', 'zhatvar']. This indicated to store the value of the estimate zhat and predicted error variance zhatvar with each iteration. We can recover these values from solver.hist_dict, the history dictionary. Using the values we can compute and plot the normalized MSE on each iteartion. We see that VAMP gets a low MSE in very few iterations, about 10.
In [18]:
# Compute the MSE as a function of the iteration
zhat_hist = solver.hist_dict['zhat']
zhatvar_hist = solver.hist_dict['zhatvar']
nit = len(zhat_hist)
mse_act = np.zeros(nit)
mse_pred = np.zeros(nit)
for it in range(nit):
zerr = np.mean(np.abs(zhat_hist[it]-z)**2)
mse_act[it] = 10*np.log10(zerr/zpow)
mse_pred[it] = 10*np.log10(zhatvar_hist[it]/zpow)
plt.plot(range(nit), mse_act, 'o-', linewidth=2)
plt.plot(range(nit), mse_pred, 's', linewidth=1)
plt.xlabel('Iteration')
plt.ylabel('Normalized MSE (dB)')
plt.legend(['Actual', 'Predicted'])
plt.grid()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: