leaders.zip_desk1oo_matlab_zip

% Created 2012-Feb-9 % Written by Fu Lin % A customized interior point method to solve the relaxed convex problem of % the leader selection problem. % Input: L -> Laplacian matrix % k -> number of leaders % Da -> absolute feedback gain matrix % Output: Jcr -> objective function of the convex relaxation % dcr -> relaxed vector from [0,1] % % Convex relaxation % minimize J(d) = Tr( (L + Da * diag(d))^(-1) ) % subject to 0 <= d_i <= 1 % sum(d) = k % where Da is the diagonal matrix of the absolute feedback gains of the % nodes, % k is the number of leaders. % % We use interior point method to eliminate the inequality constraints % to have an approximate problem % minimize J_gam(d) = gam*J(d) - ( sum( log(d_i) + log(1-d_i) ) ) % subject to sum(d) = k. % % We solve a sequence of these approximate problems parameterized by gam. As % gam increases to infinity, we have good approximation of the solution for % the convex relaxation. % % To solve the approximate problem, we use Newton's method subject to a % linear constraint function [Jcr,dcr] = cvxrelax(L,k,Da) n = size(L,1); % absolute feedback gain alpha Delta = Da; % approximation accuracy gam = 1; mu = 5; gamtol = 1.e+5; % initial condition for d d = k/n*ones(n,1); if abs(sum(d)-k) > 1.e-5 error('Initial condition does not choose k leaders!') end % vector of all ones One = ones(n,1); % backtracking line search parameters alpha = 0.3; beta = 0.5; % tolerance for the stopping criterion of the Newton's method tol = 1.e-6; % step number stepno = 0; % display message control quiet = 0; if ~quiet disp('StepNo Objective value NTdecre stepsize gamma') end while 1 while 1 invLD = inv( L + Delta*diag(d) ); Jgam = gam * trace(invLD) - ( sum( log(d) + log(1-d) ) ); invLD2 = invLD * invLD; %#ok<*MINV> gradJ = - gam * diag(invLD2*Delta) - ( 1./d + 1./(d-1) ); hessJ = gam * 2 * (Delta * invLD2 * Delta) .* invLD + diag( 1./(d.^2) + 1./((1-d).^2) ); % Newton's method to solve the gam-parameterized problem v1 = hessJ\gradJ; v2 = hessJ\One; lam = - One'*v1/(One'*v2); dnt = - v1 - lam*v2; % test if Newton's direction satisfies the linear constraint if abs(sum(dnt)) > 1.e-5 error('Newton direction does not satisfy linear constraint!') end % if the Newton decrement is small enough, break decre = dnt'*hessJ*dnt/2; if decre < tol break; end % use backtracking line search to find the step-size stepsize = 1; while 1 dtemp = d + stepsize*dnt; invLDtemp = inv(L + Delta*diag(dtemp) ); Jgam_temp = gam * trace(invLDtemp) - ( sum( log(dtemp) + log(1-dtemp) ) ); % L + Dtemp needs to be positive definite minEigLDtemp = min(eig(L + Delta*diag(dtemp))); if minEigLDtemp > 0 && Jgam - Jgam_temp > stepsize*alpha*(-gradJ'*dnt) && min(dtemp) > 0 && min(1-dtemp) > 0 break; end stepsize = stepsize*beta; if stepsize < 1.e-16 error('extremly small stepsize!'); end end d = dtemp; Jgam = Jgam_temp; stepno = stepno + 1; if ~quiet disp([num2str(stepno),' ',num2str(Jgam,'%11.4E'),' ', num2str(decre,'%11.4E'),' ',num2str(stepsize),' ',num2str(gam)]) end end gam = gam*mu; if gam > gamtol break; end end dcr = d; Jcr = trace(invLD);