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

log_cheb.m

function [x, ul] = log_cheb(A)

% [x, ul] = log_cheb(A);
%
% solves the logarithmic Chebychev approximation problem
%
% (primal) minimize    t 
%          subject to  (1/t)*e <= A*x <= t*e
%
% Input arguments:
% - A:    pxk matrix.
%
% Output arguments:
% - x:    k-vector.  Optimal solution.
% - ul:   2-vector.  Upper and lower bound on optimal t.


[p,k] = size(A);

%
% build F and c
%

% F = 
F = zeros(5*p,k+2);
blck_szs = ;
F(p+[2:4:4*p],1) = ones(p,1);
F(p+[3:4:4*p],1) = ones(p,1);
F(1:p,2:k+1) = -A;
F(p+[1:4:4*p],2:k+1) = A;
F(1:p,k+2) = ones(p,1);
F(p+[4:4:4*p],k+2) = ones(p,1);

c = ;


%
% primal feasible point 
%

% find initial point with A*x > 0 using phase1
disp(' '); disp(' PHASE 1.');
nu=20.0; abstol=1e-8; maxiters=100;
[x0,Z,z,ul,iters] = ...
   phase1([zeros(p,1) A], ones(p,1), p,nu,abstol,maxiters);
if iters == maxiters
   error('Maximum number of iterations exceeded.');
end;
if (ul(1) > 0.0),  
   error('The inequalities A*x > 0 are not feasible.');  
end;

% t0 > a_i'*x0  and  t0 > 1.0/a_i'*x0 
x0 = ;


%
% Dual points have the form: Z = diag ( diag(z), Z^1, ... Z^p );
% with z a p-vector, and Z^i a 2x2 - matrix.
% Conditions for feasibility:
% 1. sum a_i * ( z_i  - Z_{11}^i ) = 0
% 2. sum (z_i + Z_{22}^i) = 1
% 3. Z >= 0
% Objective function: -2 * sum Z^i_{12}
% Initial dual feasible point (with obj. value zero):
% - solve A'*w = 0;  \| w \|_\infty = 0.5
% - split each w_i in two positive parts z_i and Z_{11}^i
% - scale everything to make e'*z < 1. Z_{22}^i = (1-e'*z)/p 
% - Z_{12}^i = 0
%   

Z0 = zeros(5*p,1);
[q,r] = qr(A);  
w = q(:,k+1);  w = (0.5/sum(abs(w)))*w;   % A'*w = 0; ||w||_1 = 0.5
Z0(1:p) = max(w,0) + 0.001/p;             % z
Z0(p+[1:4:4*p]) = -min(w,0) + 0.001/p;    % Z_{11}^i
Z0(p+[4:4:4*p]) = ((1.0 - sum(Z0(1:p)))/p) * ones(p,1);  % Z_{22}^i

disp(' '); disp(' PHASE 2.');
nu=20.0; abstol=1e-8; reltol=1e-3;  tv=0.0;  maxiters=100;
[x,Z,ul,info,time] = ...
   sp(F,blck_szs,c,x0,Z0,nu,abstol,reltol,tv,maxiters);
if time(3) == maxiters,
   error('Maximum number of iterations exceeded.');
end;
x = x(1:k);

Generated by  Doxygen 1.6.0   Back to index