Logo Search packages:      
Sourcecode: semidef-oct version File versions  Download package

phase1.m

function [x,Z,z,ul,iters] = phase1(F,blck_szs,M,nu,abstol,maxiters);

% [x,Z,z,ul,iters] = phase1(F,blck_szs,M,nu,abstol,maxiters);
%
% Find an x s.t. F(x) > 0 and Tr F(x) < M 
% or prove that no such x exists.
%
% phase1.m determines whether the optimal value of the semidefinite
% program
%  
%  minimize    t 
%  subject to  F(x) + t*I = F_0 + x_1*F_1 + ... + x_m*F_m + t*I >= 0 
%              Tr F(x) <= M
% 
% is positive or negative.
% 
% It produces either x and t with t<0.0 or a solution Z,z to 
% the dual problem
%
%  maximize    -Tr F_0*(Z-z*I) - M*z
%  subject to  0 = Tr F_i*(Z - z*I), i=1,...,m
%              1 = Tr Z 
%              Z >= 0, z >= 0
% 
% with objective value positive.
%
% Convergence criteria:
% (1) maxiters is exceeded
% (2) duality gap is less than abstol
% (3) primal objective is less than zero or dual objective is greater
%     than zero
%  
% Input arguments:
% F:         (sum_i n_i^2) x (m+1) matrix
%            [ F_0^1(:) F_1^1(:) ... F_m^1(:) ]
%            [ F_0^2(:) F_1^2(:) ... F_m^2(:) ]
%                ...      ...          ...
%            [ F_0^L(:) F_1^L(:) ... F_m^L(:) ]
%            where F_i^j denotes the jth block of F_i, 
%            F_i^j has size n_i x n_i.
% blck_szs:  L-vector [n_1 ... n_L], contains dimensions of blocks.
% M:         scalar.
% nu:        >= 1.0.  Controls the rate of convergence.
% abstol:    absolute tolerance.
% maxiters:  maximum number of iterations.
%
% Output arguments:
% x:         m-vector; last primal iterate
% Z:         last dual iterate;  block-diagonal matrix stored as 
%            [ Z^1(:); Z^2(:); ... ; Z^L(:) ]
%            where Z^i is the ith diagonal block.
% z:         scalar part of last dual iterate. 
% ul:        2-vector, primal and dual objective.
%            If ul(1) < 0.0, the problem is feasible and x is a 
%            solution.  If ul(2) > 0.0, the problem is infeasible
%            and Z,z form a proof of infeasibility.  If ul(1) > 0.0
%            and ul(2) < 0.0, no conclusion can be made.
% iters:     scalar.

if (size(blck_szs,2) < size(blck_szs,1)), blck_szs = blck_szs'; end;
if (size(blck_szs,1) ~= 1), error('blck_szs must be a vector');  end;

m = size(F,2)-1;
if (size(F,1) ~= sum(blck_szs.^2)) 
   error('Dimensions of F do not match blck_szs.');
end;

% mineigF is the smallest eigenvalue of F_0
mineigF = 0.0;
k=0; for n=blck_szs,
   mineigF = min(mineigF, min(eig(reshape(F(k+[1:n*n],1),n,n))));
   k=k+n*n;   % k = sum n_i*n_i 
end;

% I is the identity
I = zeros(size(F,1),1);  
k=0; for n=blck_szs,
   I(k+[1:n*n]) = reshape(eye(n),n*n,1);   % identity
   k=k+n*n;   % k = sum n_i*n_i 
end;

if (M < I'*F(:,1)+1e-5), 
   error('M must be strictly greater than the trace of F_0.'); 
end;

% initial x0 
x0 = [zeros(m,1); max(-1.1*mineigF, 1e-5)];

% Z0 is the projection of I on the space Tr F_i*Z = 0
Z0 = I - F(:,2:m+1) * ( F(:,2:m+1) \ I );

% mineigZ is the smallest eigenvalue of Z0
mineigZ = 0.0;
k=0; for n=blck_szs,
   mineigZ = min( mineigZ, min(eig(reshape(Z0(k+[1:n*n]),n,n))) );
   k=k+n*n;   % k = sum n_i*n_i 
end;

% z = max ( 1e-5, -1.1*mineigZ)
Z0(k+1) = max( -1.1 * mineigZ, 1e-5 );  
Z0(1:k) = Z0(1:k) + Z0(k+1)*I; 
Z0 = Z0 / (I'*Z0(1:k));    % make Tr Z0 = 1

% add scalar block Tr F(x) <= M
F = ; 
blck_szs = ;
c = ;

[x,Z,ul,infostr,time] = ...
           sp(F,blck_szs,c,x0,Z0,nu,abstol,-1.0,0.0,maxiters);

z = Z(k+1);
Z = Z(1:k);
x = x(1:m);
iters = time(3);

Generated by  Doxygen 1.6.0   Back to index