1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
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 = [ diag(t*e - A*x) 0 ; 0 diag([a_i'*x 1.0; 1.0; t]) ]
F = zeros(5*p,k+2);
blck_szs = [ ones(p,1); 2*ones(p,1) ];
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 = [zeros(k,1); 1];
%
% 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 = [x0; 0.001 + max( max(A*x0), max(1./(A*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);
|
