Data Science Part 24: Approximate Message Passing

Approximate message passing is an iterative signal recovery algorithm that decouples the linear inverse problem y = Ax + z into N parallel scallar problems. x is independent and identically distributed and we add some Gaussian noise in the mix, denoted by the symbol Z. X follows the probability and density function fX.

Sometimes, we don't know we have a particular signal; however, we do have noisy estimates of it. Extra Gaussian Noise can often be difficult to deal with. We are beginning with a linear inverse problem y =h + z =  Ax + z. Every single measurement will be performed through a scalar estimation scheme.  v = x + noise, and we approximate a precise belief propagation approach. The estimate is shown through this bipartite graph: 

Step 1 is to calculate the residual which is the (noisy measurements) y - (predicted measurements) AxT and adding a correction term. We define this as rt. 

Then we calculate the pseudo data, which is x + AT * residual. Finally in the denoising step, X is a denoising function of the pseudo data Vt. We use a slightly modified denoiser as we iterate through. We can utilize structure inside of the image to reduce the affect that the noise has.

The correction term, called the OnSager Correction term, calculates the derivative of the denoiser and has a few rules. 

1. White error means that the Error (x^t - x) will be uncorrelated with x

2. Error will be Gaussian.

Without the Onsanger term and these rules, the error will quickly be correlated with x and then the denoiser barely works. In a standard AMP because everything is i.i.d we will use a scalar denoiser which is element by element. 

Another aspect that we need to consider is damping. Sometimes imprecisions in an algorithm can make it proceed a bit too assertively and it diverges so we can try to slow down the denoising process with λ as the damping parameter and the damping procedure can be defined as 

We can result in a state evolution for AMP, and with separate denoisers we have a characterization of the variance in the pseudodata in the next iteration conditioned on the current iteration. M/N is the aspect ratio. The exapected value in the below equation is the mean squared error of the denoising prpoblem. So we do 1/aspect ratio multiplied by the mean squared error in the denoising function.  This would be the state evolution of Approximate message passing. We can also generalize thie conjecture to the broad class of denoisers. 

Here is an example for a Markovian input. we typically have multiple zero's and one's in a row. If we look at our neighbors in a window of size 3, Vt-1, Vt, Vt+1 we have a better prediction in contrast to if we have a sliding window size of 1. The red cross marks are the empirical results of simulation. 

An improved approach is the following in regards to the scalar channel si:

After several iteration, the scalar channel becomes clean. As a result, we can use a thresholding operator. X is sparse and many of the V's will be close to 0. We use the thresholding operator in the last iteration, which tunes the tradeoff between false positive and negatives and serves as a near-optimal heuristic. 

Let's do an example. A measurement matrix AϵℝM*N. The measurement rate is the aspect ratio of a matrix. The sparsity of a matrix can be quantified with a score, which is the number of zero values in the matrix divided by the total number of elements in the matrix. Again, N is the length of the signal and M is the number of measurements in the signal. 

To denoise the signal, I have to computer the conditional expectation of x given the noisy measurements V which is just basically E[X|V]. The true value X is some scalar random variable. In the measurement process, a measurement matrix is used to sample only the components that best represent the signal. A Wiener filter is one of the filters that can be used to measure matrices and is denoted as follows (also detnoted as the Minimum Mean Squared Error Estimator): 

The Weiner filter can also be estimated as the conditional denoiser. 

And here is some estimation of the corresponding X signal: 

We want to multiply this Weiner filter by a probability density function where the distribution of X 

Here is some code for AMP: 

AMP Algorithm:

#---------

# AMP algorithm

#---------

# initialization

mse = np.zeros((max_iter,1)) # store mean square error

xt = np.zeros((N,1))# estimate of signal

dt = np.zeros((N,1))# derivative of denoiser

rt = np.zeros((M,1))# residual

for iter in range(0,max_iter):

 # update residual

 rt = y - np.dot(A,xt) + 1/ratio*np.mean(dt)*rt

 # compute pseudo-data

 vt = xt + np.dot(Atrans,rt)

 # estimate scalar channel noise variance estimator is due to 

Montanari

 var_t = np.mean(rt**2)

 # denoising

 xt1,dt = denoise(vt,var_x,var_t,epsilon)

 # damping step

 xt=lamda*xt1+(1-lamda)*xt

 mse[iter] = np.mean((xt-x)**2)

Now for Addendum Say we want to discuss Tanaka's fixed point equation which in turn measures the inverse of the degradation of a scalar channel. It shows that the larger the signal to noise ratio there is less degradation. The greater the aspect ratio (number of measurements in the measurement matrix) the greater the degradation. The equation is as follows: 

where the η is degradation of the scalar channel, , δ = M/N > 0 is the aspect ratio of the measurement matrix, MMSE() is the minimum mean squared error obtained by a denoiser, and γ > 0 is the signal to noise ratio and ε, It basically calculates the degradation. We can approximate this equation by calculating the mean square value of the denoiser where V = X + W as 

MMSE = E[(X - E[X|V])^2]. 

An alternative is to simulate many signals of X and scalar channel, denoise it, and compute the Mean squared error. 

Let's consider a Rademacher signal where each entrie of the vector x takes values -1 and 1 with probability ∈, else 0 with probability 1 - 2∈. Estimating x from the pseudo data can result in v = x + N(0, σ^2) where conditional probability density function is 

.

We can compute the the Conditional Expectation E[Xn|Vn] = Pr(Xn = 1|Vn) - Pr (Xn = -1)|Vn. 

Let's consider a communication signal that modulates a binary vector. Again a lot of applications measure x with additive white Gaussian noise in the form of v = x + q where the random variable Q is identically and independently distributed. Let's develop a denoising function for x, or E[X|V = v]. We want to make sure that the probability density function for each entries satisfies the following: 

. 

If we want to get a denoising function for a binary vector lets say x ε (-2, +3) ^ N, then we have the following: 

.