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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
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.
// Copyright INRIA
// 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
//
// maximize -Tr F_0*(Z-z*I) - M*z
// subject to 0 = Tr F_i*(Z - z*I)
// 1 = Tr Z
// Z >= 0, z >= 0
//
// Convergence criterion:
// (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
// [ vec(F_0(1)) vec(F_1(1)) ... vec(F_m(1)) ]
// [ vec(F_0(2)) vec(F_1(2)) ... vec(F_m(2)) ]
// ... ... ...
// [ vec(F_0(L)) vec(F_1(L)) ... vec(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
// abstol: absolute tolerance
// maxiters: maximum number of iterations
//
// Output arguments:
// x: m-vector:
// Z: block-diagonal matrix stored as
// [ vec(Z(1)); vec(Z(2)); ... ; vec(Z(L)) ]
// where Z(i) is the ith block
// z: scalar part of last dual iterate,
// ul: 2-vector, primal and dual objective
// iters: scalar
if type(F)==5 then F=full(F);end
[mf,nf]=size(F)
m = nf-1;
if (mf ~= 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;
blck_szs=matrix(blck_szs,1,prod(size(blck_szs)));
k=0; for n=blck_szs,
mineigF = min(mineigF, min(real(spec(matrix(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(mf,1);
k=0; for n=blck_szs,
I(k+[1:n*n]) = matrix(eye(n,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(real(spec(matrix(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 ); // z
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 = [F, I; M-I'*F(:,1), -I'*F(:,2:m+1), 0];
blck_szs = [blck_szs,1];
c = [zeros(m,1); 1];
//Pack Z0 and F
Z0=pack(Z0,blck_szs);F=pack(F,blck_szs);
pause;
[x,Z,ul,info]=semidef(x0,Z0,F,blck_szs,c,[nu,abstol,-1,0,maxiters]);
nz=prod(size(Z))
z=Z(nz)
Z=unpack(Z(1:nz-1),blck_szs(1:prod(size(blck_szs))-1))
Z = Z(1:k);
x = x(1:m);
iters = info(2);
select info(1)
case 1
error('Max. iters. exeeded');
case 2 then
disp('Absolute accuracy reached');
case 3 then
disp('relative accuracy reached');
case 4 then
disp('target value reached');
case 5 then
error('target value not achievable');
else
error('semidef fails');
end
|
