%A must be of the form [A11,A12;rxN,rxr] %A11 must be an NxN diagonalizable upper triangular matrix %A12 must be an rxr matrix with entries<1 %A11=[a_1I_(n_1), * , * ,...,*; zeros(n_2,n_1),a_2I_(n_2),*,...,*; ... ;zeros(n_k,n_1+...+n_(k-1)),a_kI_(n_k)] %n=[n_1,n_2,...,n_k]; %diagonal entries (a1,...,ak) of A11 must be less than 1 function [P,Q]=Twoposcon(A11,A12,k,n) format short %obtaining diagonal D=a_1 I_{n1}\oplus ... \oplus a_k I_{nk} %creating invertible V such that A11V=DV with unit columns and %col_i^{T}col_j=0 whenever d_{i}=d_{j} if min(eig(A11*A11'+A12*A12'))>1 fprintf('[A11,A12;0,0] is not a product of two positive projections.\n'); else V=[]; [W,D]=eig(A11); for ii=1:k [Q,~]=qr(W(:,size(V,2)+1:size(V,2)+n(ii)),0); V=[V,Q]; end %APPLYING ALTERNATING PROJECTION METHOD to obtain Gamma X=sqrt(D)*V'*inv(A11*A11'+A12*A12')*V*sqrt(D); Y=V'*V; if min(eig(X-Y))<0 fprintf('[A11,A12;0,0] is not a product of two positive projections.\n'); else err=10^(-7); %tolerance (set this to a smaller value for more accuracy) iterlimit=100000; %max number of iterations (set this to a higher value for more iteration) iter=1; %initial Gamma Gamma=[]; for ii=1:k indi=size(Gamma,1)+1:size(Gamma,1)+n(ii); Gamma=blkdiag(Gamma,(X(indi,indi)+Y(indi,indi))/2); end tic while iter<=iterlimit %projection onto X>=Gamma constraint (Omega1) [W,S]=eig(X-Gamma); S(S<0)=0; GammaT=X-W*S*W'; %projection onto Omega0 Gamma=[]; for jj=1:k indi=size(Gamma,1)+1:size(Gamma,1)+n(jj); [W,S]=eig(GammaT(indi,indi)); S(S<0)=0; Gamma=blkdiag(Gamma,W*S*W'); end %stopping criterion err1=abs(min(0,min(eig(Gamma-Y)))); err2=abs(min(0,min(eig(X-Gamma)))); error=err1+err2; fprintf('iter %d error %f \n',iter,error); if error=Y constraint (Omega2) [W,S]=eig(Gamma-Y); S(S<0)=0; GammaT=W*S*W'+Y; %projection onto Omega0 Gamma=[]; for jj=1:k indi=size(Gamma,1)+1:size(Gamma,1)+n(jj); [W,S]=eig(GammaT(indi,indi)); S(S<0)=0; Gamma=blkdiag(Gamma,W*S*W'); end %stopping criterion err1=abs(min(0,min(eig(Gamma-Y)))); err2=abs(min(0,min(eig(X-Gamma)))); error=err1+err2; fprintf('iter %d error %2.10f \n',iter,error); if error