Code covered by the BSD License

# Grassmannian Design Package

by

### Ahmed (view profile)

This software can be used to generate Grassmannian line/subspace packings

[x fval]=OrthogonalAlgGenMatX3fast()
```function [x fval]=OrthogonalAlgGenMatX3fast()
clc
close all
clear all
global M
global Mt
global Qpre
global skip
global power
global Index
disp('======================================================================================')
disp('This code generates Complex Grassmannian Subspace Packings with Fubini-Study distance ')
disp('======================================================================================')
skip=0;
Qpre=input('Please enter the number of Precoding matrices (Packings)');
QLeft=Qpre;
Mt=input('Please enter the number of rows for each Precoder matrix ');
M=input('Please enter the number of columns for each Precoder matrix ');
disp(['You selected to design ',num2str(Qpre),' ',num2str(Mt),'x',num2str(M),' Precoders'])
ss=input('To start Press 1 or any other key to abort ');
if ss~=1
disp('Stopped by User')
x=[];
return
end
if isempty(ss)==1
disp('Stopped by User')
x=[];
return
end
x=[genorth(Mt,M,QLeft)];
tol=10^-4;
mag=100;
ma=10;
power=2;
Index=genindex(Qpre);
disp('=================')
disp('Optimizer Started')
disp('=================')
for jj=1:50
for ii=1:10
[x fval]=orthcon(x,@objectivefn,mag);
[x fval]=orthcon(x,@objectivefn2,ma);
mag=mag*10000;
power=2+power;
if power==20
power=2;
mag=10;
end

end
if mod(jj,1)==0
ss=input('Do you wish to continue? Press 1 to stop or Enter to continue');
if ss==1
break
end
close all
end
end
disp('End of Optimzation ')
disp('===================')
DispRes(x)
%==========================================================================
function [x fval]=orthcon(x,objectivefn,mag)
%====init=====
global  M Mt
global Qpre
global skip
tol=10^-2;
maxiter=10;
xbest=x;
fbest=objectivefn(x);
direction=[zeros(1,skip) ones(1,(Qpre-skip))];
rate=0;

%=========================================================================
for kkk=1:maxiter
gamma=ones(1,Qpre);
Z=[];
ttrace=[];
cont=[];
for jj=1:Qpre

f=x(:,(jj-1)*M+1:jj*M);
h=f*f';
I1=eye(size(h));
con1=(h-I1);

cont=[cont con1];
z1=direction(jj)*con1*Dx(:,(jj-1)*M+1:jj*M);

Z=[Z z1];
ttrace=[ttrace trace(z1'*z1)];
end

XX=x+Z;
xnn=x;
xn=xnn;
dd=feval(objectivefn,x);
for jj=(skip+1):Qpre
xn=xnn;
XX1=XX(:,(jj-1)*M+1:jj*M);

[Q1 R1]=qr(XX1);
xn(:,(jj-1)*M+1:jj*M)=Q1(:,1:M);
QRfac=feval(objectivefn,xn);
ii=0;
%     if (dd-QRfac)<0
%         [Q1 R1]=qr(YY1);
%         xn(:,(jj-1)*M+1:jj*M)=Q1(:,1:M);
%         QRfac=feval(objectivefn,xn);
% %         continue
%     end
QRfacinit=QRfac;
fact=1;
while (dd-QRfac)>= 0 && ii<5
gamma(jj)=2*gamma(jj);
temp=x(:,(jj-1)*M+1:jj*M)+gamma(jj)*Z(:,(jj-1)*M+1:jj*M);
[Q1 R1]=qr(temp);
xn(:,(jj-1)*M+1:jj*M)=Q1(:,1:M);
QRfac=feval(objectivefn,xn);
if QRfac>QRfacinit
break
end
QRfacinit=QRfac;
ii=ii+1;
fact=2;
end
gamma(jj)=gamma(jj)/fact;
temp=x(:,(jj-1)*M+1:jj*M)+gamma(jj)*Z(:,(jj-1)*M+1:jj*M);
[Q1 R1]=qr(temp);
xn(:,(jj-1)*M+1:jj*M)=Q1(:,1:M);
QRfac=feval(objectivefn,xn);
ii=0;
fact=1;
while (dd-QRfac)<0 && ii<5
gamma(jj)=gamma(jj)*0.25;
temp=x(:,(jj-1)*M+1:jj*M)+gamma(jj)*Z(:,(jj-1)*M+1:jj*M);
[Q1 R1]=qr(temp);
xn(:,(jj-1)*M+1:jj*M)=Q1(:,1:M);
QRfac=feval(objectivefn,xn);
ii=ii+1;

end
gamma(jj)=gamma(jj)/fact;
temp=x(:,(jj-1)*M+1:jj*M)+gamma(jj)*Z(:,(jj-1)*M+1:jj*M);
[Q1 R1]=qr(temp);
xnn(:,(jj-1)*M+1:jj*M)=Q1(:,1:M);
dd=feval(objectivefn,xnn);
end
xn=x;
for jj=(skip+1):Qpre
XX1=xnn(:,(jj-1)*M+1:jj*M);
[Q1 R1]=qr(XX1);
xn(:,(jj-1)*M+1:jj*M)=Q1(:,1:M);
end
x=xn;
fval=feval(objectivefn,x);
disp(['iteration number ',num2str(kkk),'  Function value ',num2str(fval)]);
%=========================================================================

mm=fval;
if mm<fbest
%     disp('fbest changed')
fbest=mm;
xbest=x;
rate=0;

elseif mm>=fbest

rate=rate+1;
if rate==2
x=xbest;
mag=mag/10;
end
if rate==4
x=xbest;
mag=mag*100;
end
if rate==6
x=xbest;
mag=mag/1000;
end

if rate==10
x=xbest;
break;
end
end
%==========================================================================

nn=isnan(x);
if nn>0
disp('Nan')

return
end
end

disp('Maximum iteration reached!')
DispRes(xbest)
%=========================================================================
beep
return
%=========================================================================
f0 = feval(fcn,x);
t=zeros(size(x));
delta = 1e-4;
for ii=1:1:size(x,1)
for kk=1:1:size(x,2)
z=t;z(ii,kk)=1*delta;
xx=x+z;
f1=feval(fcn,xx);
end
end
for ii=1:1:size(x,1)
for kk=1:1:size(x,2)
z=t;z(ii,kk)=1i*delta;
xx=x+z;
f1=feval(fcn,xx);
end
end
%==========================================================================

function f=objectivefn(F)
global M power
global Index
Dn=F'*F;
lenn=size(Dn,2)/M;
s=Dn(:,1:M);
for ii=2:lenn
s=[s;Dn(:,(ii-1)*M+1:ii*M)];
end
s=s.';
lenn=length(Index);
ee=[];
for jj=1:lenn
ii=Index(jj);
d=(s(:,(ii-1)*M+1:ii*M));
%     ee=[ee abs(det(d))];
ee=[ee (abs((1+det(d'*d))*log(1+det(d'*d))))];
end
f=mean(ee.^power);

%==========================================================================

function []=DispRes(F)
global M
global Qpre
ddd=[];
ccc=[];
ppp=[];
for ii=1:Qpre
for jj=ii:Qpre
if ii==jj
continue
end
f1=F(:,(ii-1)*M+1:ii*M);
f2=F(:,(jj-1)*M+1:jj*M);
d=abs(det(f2'*f1));
d=acos(d);
ddd=[ddd d];
ccc=[ccc 1/sqrt(2)*norm((f1*f1'-f2*f2'),'fro')];
ppp=[ppp norm((f1*f1'-f2*f2'),2)];

end
end
fsmin=min(ddd);
chmin=min(ccc);
pnmin=min(ppp);
disp(['Minimum Fubini Study distance = ',num2str(fsmin)])
disp(['Minimum Chordal distance = ',num2str(chmin)])
disp(['Minimum Projection two norm distance = ',num2str(pnmin)])
figure
subplot(3,1,1)
hist((ddd))
xlabel('Fubini Study distances ')
ylabel('Distance Distribution')
subplot(3,1,2)
hist((ccc))
xlabel('Chordal distances ')
ylabel('Distance Distribution')
subplot(3,1,3)
hist((ppp))
xlabel('Projection two norm distances ')
ylabel('Distance Distribution')

%==========================================================================
%==
function q=genorth(Mt,M,Qpre)
q=[];
for ii=1:Qpre
q=[q orth(randn(Mt,M)+1i*randn(Mt,M))];
end
%==========================================================================
function f=objectivefn2(F)
global M power
global Index
Dn=F'*F;
lenn=size(Dn,2)/M;
s=Dn(:,1:M);
for ii=2:lenn
s=[s;Dn(:,(ii-1)*M+1:ii*M)];
end
s=s.';
lenn=length(Index);
ee=[];
for jj=1:lenn
ii=Index(jj);
d=(s(:,(ii-1)*M+1:ii*M));
ee=[ee exp((1+det(d'*d))*log(1+det(d'*d)))];
end
f=norm(ee,'inf');
%==========================================================================
function index=genindex(Mt)
d=zeros(Mt,Mt);
for ii=1:Mt*Mt
d(ii)=ii;
end
WM=genweight(Mt);
d=d.*WM;
index=reshape(d,numel(d),1);
index=unique(index);
index=index(2:end);

%==========================================================================
function WM=genweight(Mt)
WM=zeros(Mt,Mt);
Qpre=Mt;
for ii=1:Qpre
for jj=ii:Qpre
WM(ii,jj)=1;
end
end
WM=1-WM;

```