Code covered by the BSD License  

Highlights from
Grassmannian Design Package

Grassmannian Design Package

by

 

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);
Dx=mag*NumGradI(objectivefn,x);
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    
%=========================================================================
function [Dx] = NumGradI(fcn,x)
f0 = feval(fcn,x);
t=zeros(size(x));
delta = 1e-4;
gradR=t;
gradI=t;
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);
        gradR(ii,kk)=(f1-f0)/delta;
    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);
        gradI(ii,kk)=(f1-f0)/delta;
    end
end
Dx=gradR+1i*gradI;
%==========================================================================    

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;

    
    


Contact us