function [R,Z,VT,t] = mainHKM(A,B,C); % min tr AXBX^t + CX^t % input: A = A^t ... n*n % B = B^t ... n*n % C ... n*n % matrix R: ((n-1)^2+1)*((n-1)^2+1) % matrix Z: ((n-1)^2+1)*((n-1)^2+1) % (P) min LR s.t. G_{P,Q}(VT*R*VT') = E_oo; R >= 0 % (D) max -t_0 % s.t. L+ VT'*G(t)*VT - Z = 0; Z >= 0, % 14.6.2002, wolkowicz & sotirov, waterloo % gangster indices: [Jb,i,j] = Jbar(n); % call: [R,Z,VT,t] = mainHKM(A,B,C); n = size(A,1); %size of the problem stoptoler=1e-5; [Jb,P,Q] = Jbar(n); P = [1;P]; Q = [1;Q]; eq = length(P); % initialization t = zeros(eq,1); % belongs to gangster part a = zeros(n*n+1,1); a(1,1) = -1; b = zeros(eq,1); b(1,1) = -1; % old definition of V V = [ eye( n-1 ); - ones( n-1,1)']; e = ones(n,1); VT = [1, zeros(1,(n-1)*(n-1));kron(e,e)/n, kron(V,V) ]; VT = sparse(VT); % new definition of V to guarantee zeros in R: nm1=(n-1)^2; MD=[1 spalloc(1,nm1,0);-ones(nm1,1)/n speye(nm1)]; VT=VT*MD; R = [ 1, zeros(1,(n-1)*(n-1)); zeros((n-1)*(n-1),1),kron(n*eye(n-1)-ones(n-1),n*eye(n-1) - ones(n-1))/(n*n*(n-1))]; %clanka Di=inv(MD); R = Di*R*Di'; c = C(:); LT = [ 0, c'/2; c/2, kron(B,A)]; %dimension (n^2+1)*(n^2+1) L = VT'*LT*VT; %dimension ((n - 1)^2 + 1)*((n - 1)^2 + 1) [Y,t] = initial(L,VT,P,Q,eq,n,Jb); Z = L + VT'*Y*VT; t(2:end)=t(2:end)*2; % just fixing for the first step done = 0; k1 = 0; disp(' '); disp('------------------------------------------------------------------------- '); disp(' min tr AXBX^t + CX^t, G_g( VT*R*VT^t)= E_oo, R \succ 0 '); disp('------------------------------------------------------------------------- '); disp(' num L(:)^t*R(:) dual mu t1 t2 relres') disp('------------------------------------------------------------------------- '); RVT = R*VT'; RT = VT*RVT; % RT = VT*R*VT' % calculating initial mu: mu = 1.5*Z(:)'*R(:)/((n-1)*(n -1) + 1); while done < 1; % while done < 1 k1 = k1 + 1; % left side invZ = inv(Z); VTiZ = VT*invZ; ZT = VTiZ * VT'; [M] = left(ZT,RT,P,Q,eq,n); % forming the matrix of op.al % right side [D,GTt,U1] = right(mu,n,VT,Z,invZ,R,L,VTiZ,t,RT,a,P,Q,eq,b); [Mpom,p] = chol(M); if p> 0, disp('problem with cho'); keyboard end; wttemp = (Mpom')\ D; rez = Mpom \ wttemp; % delta omega dt = rez; if norm(M*rez-D)> stoptoler; disp(['relaccuracy: ',num2str(norm(M*rez-D)/norm(D))]); end; [GTt] = GT(dt,P,Q,eq,n); U2 = VT'*GTt*VT; dZ = U1 + U2; dZ = (dZ + dZ')/2; dR = invZ*(mu*eye((n-1)*(n-1)+1) - dZ*R ) - R; VdRV = VT*dR*VT'; % FOR TESTING: % pom = 2:n*n+1; % L1 = diag(VdRV) - [0; (VdRV(pom,1) + VdRV(1,pom)')/2]; % R1 = - (diag(RT) - [0; (RT(pom,1) + RT(1,pom)')/2]) - a; nn1=norm(L1-R1); %%-------------------------------------------------------------------- % L2 = VT'*(ATw+GTt)*VT-dZ; R2 = - U1; nn2=norm(L2 - R2); %%-------------------------------------------------------------------- % L3= G(VdRV,P,Q,eq,n); R3 = G(RT,P,Q,eq,n); nn3=norm(L3-R3); %%-------------------------------------------------------------------- % L4 = dZ*R + Z*dR; R4 = mu*eye((n-1)*(n-1) + 1) - Z*R; nn4=norm(R4-L4); %nnn=[nn1 nn2 nn3 nn4]; %if max(nnn) >1e-3; disp('error'); keyboard; end dR = ( dR + dR' )/2; % delta R VdRV = VT*dR*VT'; t1 = 1; stop1 = 0; while stop1 < 1 Ztmp = Z + t1*dZ; [Z1,p] = chol(Ztmp); if p == 0, stop1 = 1; else t1 = t1 * 0.8; end; if t1 < 1e-7 stop1 = 1; disp('problem with stepsize...'); end; end; if t1 <1; t1 = t1 * 0.95; end; % to make sure the next point is not to close to boundary Z = Z + t1*dZ; t = t + t1*dt; [GTt] = GT(t,P,Q,eq,n); t2 = 1; stop2 = 0; Rpom = VT*R*VT'; while stop2 < 1 Rtmp = R + t2*dR; [R1,p] = chol(Rtmp); % first condition RTtmp = Rpom + t2*VdRV; if p == 0, stop2 = 1; else t2 = t2 * 0.8; end; end; if t2 <1; t2 = t2 * 0.95; end; %to make sure the next point is not to close to baundary R = R + t2*dR; RT = Rpom + t2*VdRV; if k1 >= 100; done = 1; disp(' max iter.'); end; pom = 2:n*n+1; arrow1 = diag(RT) - [0; (RT(pom,1) + RT(1,pom)')/2]; proba = norm(L + VT'*GTt*VT - Z,1) + norm( arrow1+a,1); if proba >0.1; mu_fact = .95; elseif min(t1,t2) < .2; mu_fact=1.05; elseif min(t1,t2) > .6; mu_fact = .25; else mu_fact = .8; end; [GRT] = G(RT,P,Q,eq,n); mu = mu * mu_fact; dual = -t(1); %**************************************** primal = L(:)'*R(:); fprintf(' %3.0d %12.5f %12.5f %12.7f %8.4f %8.4f %9.4e \n',[k1, L(:)'*R(:) dual mu t1 t2 (norm(M*rez-D)/norm(D))]); if (proba < stoptoler) & (abs(primal-dual)/max(1,abs(primal))