Regarding the problem of Massive MIMO uplink detection [HBD13], the
architecture serves tens of users by employing hundreds of antennas,

where the channel has its elements
sampled from , ,
is the received signal, AWGN noise components
are i.i.d with ;
regarding the transmitted , we only assume that it’s zero
mean and finite variance .

Before incorporating the AMP algorithm, we should be well aware of
two facts: 1. directly using maximum a priori (MAP)
or MMSE estimation
to work with the exact prior
degrade the necessity of employing AMP, because achieving a full
diversity requires an extremely large set of constellation points, in
which AMP works slowly while doing the moment matching process, not
to mention problems about its inability to converge to the lowest
fixed point. 2. In the CDMA multiuser detection theory [Verdu98,
etc.], their “MMSE” detector does not mean the one working with
exact prior , but rather the one assuming a Gaussian prior.

So we use a proxy prior for detecting , i.e., assuming that
, and it will result in a MAP function that is statistically extremely close to the true MAP in Massive MIMO. In this occurrence, we have the signal power
in QPSK, in 16QAM, etc. So
the target function becomes:

The AMP algorithm to solve the above problem only requires three
lines (see [BM11], or our next paper which delivers a universal
version of this algorithm):

(1)  

(2)  

(3)  

where the initialization is to let , ,
. In terms of complexity, it only costs
. Also, according to the second
equation of the algorithm, it is converging extremely fast. On the
contrary, MMSE has complexity . It is noteworthy that
known approximation methods to MMSE, such as Richardson’s method or
Newman series approximation, both fall behind the complexity-performance
trade-off of AMP according to our simulations.

Now we share simulation results as well as MATLAB source codes.

%   AMP detector in Massive MIMO
%   written by Shanxiang Lyu (s.lyu14@imperial.ac.uk)
%   Last updated on 04/12/2015
function main()
clc;clear; close all; 
m=128;% # of received antennas
n=16;% # of users
SNRrange=[1:20];
count=0;
for s=SNRrange
SNRdb=s;
    for monte=1:1000
    x=(2*randi([0,1],n,1)-ones(n,1))+sqrt(-1)*(2*randi([0,1],n,1)-ones(n,1));
    sigmas2=2;%signal variance in QPSK
    H=1/sqrt(2*m)*randn(m,n)+sqrt(-1)/sqrt(2*m)*randn(m,n);
    sigma2=2*n/m*10^(-SNRdb/10); %noise variance in control by SNR in DB
    w=sqrt(2*sigma2)*randn(m,1)+sqrt(-1)*sqrt(2*sigma2)*randn(m,1);
    y=H*x+w; %channel model

    %iterAMP is # of iterations in AMP
    iterAMP1=2;
    xhat1=AMP(y,H,sigma2,sigmas2,iterAMP1,m,n);
     iterAMP2=4;
    xhat2=AMP(y,H,sigma2,sigmas2,iterAMP2,m,n);
    
     x_mmse=(sigma2/sigmas2*eye(n)+H'*H)^(-1)*H'*y;
     x_mmse=sign(real(x_mmse))+sqrt(-1)*sign(imag(x_mmse));
    errorAMP1(monte)=sum(x~=xhat1);
    errorAMP2(monte)=sum(x~=xhat2);
    errorMMSE(monte)=sum(x~=x_mmse);
    end
    count=count+1;
serAMP1(count)=mean(errorAMP1);
serAMP2(count)=mean(errorAMP2);
serMMSE(count)=mean(errorMMSE);
end
figure(1)% plot the SER
semilogy(SNRrange,serAMP1,'-+r', SNRrange,serAMP2,'-pk',SNRrange, serMMSE,'-ob'); 
grid on;
legend(['AMP iteration=' int2str(iterAMP1)], ['AMP iteration=' int2str(iterAMP2)], 'MMSE');
xlabel('SNR (dB)'); ylabel('SER');
title(['BER performance comparison in system m= ' int2str(m)  '  n=' int2str(n)]);
end
function xhat=AMP(y,H,sigma2,sigmas2,iterAMP,m,n)
%   AMP detector in Massive MIMO
%   written by Shanxiang Lyu (s.lyu14@imperial.ac.uk)
%   Last updated on 04/12/2015
    r=zeros(m,1);
    xhat=zeros(n,1);
    alpha=sigmas2;%initial estimation variance
        for t=1:iterAMP
        r=y-H*xhat+(n/m)*sigmas2/(sigmas2+alpha)*r;
        alpha=sigma2+(n/m)*sigmas2*alpha/(sigmas2+alpha);
        xhat=(sigmas2/(sigmas2+alpha))*(H'*r+xhat);
        end
    xhat=sign(real(xhat))+sqrt(-1)*sign(imag(xhat));
end

(pp. 2799, para 2 on the left column)”Since  only needs to be orthogonal to , we must have  and hence “.

In this post, we try to establish a proof about ““.

First of all, necessary notations are given, as those in Ref. [1]. Let  be the Euclidean distance from point  to the -th facet of the decision regions of ZF and SIC, respectively. Then the distance spectrum of them are given by:

where  denotes the angle between  and the space span by all the rest columns of matrix , and  is the angle between  and the space span by the columns before . We dubs the two notations as

In order to show that , it is indeed showing

where the “Proj” operator means taking the projection of vectors onto their orthogonal complement. Let the normalized Gram-Schmidt basis of  note as , then we have

where  and , among which the vectors are pair-wise orthogonal. Thus , and we arrive in the conclusion that

So the decoding radius of SIC is larger.

References:

[1] C. Ling, “On the proximity factors of lattice reduction aided decoding”, IEEE Transactions on signal processing, Vol.59, No.6, pp. 2795–2808, 2011.

[2] Y.H. Gan, C. Ling and W.H. Mow, “Complex lattice reduction algorithm for low-complexity full-diversity MIMO detection”, IEEE Transactions on signal processing, Vol.57, No.7, pp. 2701–2710, 2009.

In this post, we firstly try to classify the spectral graph model, and then given out the matlab code for partitioning a bipartite graph.

A general model for structural graph is:

: let  which meaning the mapping from vertices to classes, and  be its generating matrix of size , where . So the connecting probability of any edge  is .

So the weighted adjacency matrix of  looks like

,

where the symmetric property of this matrix descents from that it is an in-directed graph.

We further specialize the bisection model hereby

• Planted Bisection : the connecting probability of connecting to its own section is  and  respectively, while that of crossing is .

Suppose we are now given a graph of this form, and want to simulate the process of partitioning all these nodes into two sections, with least possible crossing edges. Before delving into those technical details, we briefly describe a spectral partitioning principle for our simulation in this post:

The eigenvector of the second largest eigenvalue of the adjacency matrix  is a slightly perturbed version of its generator matrix , which is of the form

.

Now we run a MATLAB simulation to simulate how that spectral partitioning method works.

1. Generate a graph with  vertexes, and  vertexes are connecting each other with probability , while the other  vertexes are self connecting with probability . And that of their crossing is .

clc;clear all;
n=800;
x=randperm(n);%put the sets among random positions
n_1=round(n/4);
ind_set1=x(1:n_1);
ind_set2=x(n_1+1:end);
p_1=0.8;
p_2=0.8;
q=0.001;
A(ind_set1,ind_set1)=rand(n_1,n_1)<p_1;
A(ind_set2,ind_set2)=rand(n-n_1,n-n_1)<p_2;
A(ind_set1,ind_set2)=rand(n_1,n-n_1)<q;
A(ind_set2,ind_set1)=rand(n-n_1,n_1)<q;
A=triu(A,1);%to make the adjacency matrix symmetric
A=A+A';
figure(1)
spy(A);

When this graph is generated, we can hardly tell any relation among its connection according to figure 1.

Again, we should bear in mind that the second eigenvector  of  is a perturbed version of  of the generator matrix . (means the values of  are just a little deviated from ) And why is  this important? because it is of the form , which reflects the relation of all these vertexes. So our second step of the code is

2. Get the second eigenvector of , and re-sort the positions of these  vertexes.

[V D] = eigs(A, 2);
[useless ind_recover] = sort(V(:,2));
A_hat=A(ind_recover,ind_recover);
figure(2)
spy(A_hat);

Then we can see from figure 2 below that the relation it reflects is really a bisection.

References:

1. Daniel A. Spielman, “Spectral partitioning in the planted partition model”, Lecture notes(Lecture 21), 2009.

2. Frank McSherry, “Spectral partitioning of random graphs”, STOC, 2001.

3. David Gleich, “Spectral Graph Partitioning and the Laplacian with Matlab”, 2006. (Online resources)