1 | function blks(p,Q_in,v_in) |
---|
2 | |
---|
3 | % Start by placing blocks in separate matrices |
---|
4 | [v1,dummy1,r1,dummy3]=dmperm(Q_in{1}+eye(length(Q_in{1}))); |
---|
5 | for blocks = 1:length(r1)-1 |
---|
6 | i1 = r1(blocks); |
---|
7 | i2 = r1(blocks+1)-1; |
---|
8 | Q{blocks} = Q_in{1}(i1:i2,i1:i2); |
---|
9 | v{blocks} = v_in{1}(i1:i2); |
---|
10 | S{blocks} = eye(length(v{blocks})); |
---|
11 | end |
---|
12 | |
---|
13 | pdecomp = 0; |
---|
14 | F = set([]); |
---|
15 | obj = 0; |
---|
16 | |
---|
17 | for blocks = 1:length(Q) |
---|
18 | if length(Q{blocks})> 1 |
---|
19 | if blocks <= 100 |
---|
20 | [Si,index] = blkOne(Q{blocks}); |
---|
21 | else |
---|
22 | Si = eye(length(Q{blocks})); |
---|
23 | end |
---|
24 | else |
---|
25 | Si = 1; |
---|
26 | index = 1; |
---|
27 | end |
---|
28 | Qreduced{blocks} = sdpvar(size(Si,2)); |
---|
29 | F = F + set(Qreduced{blocks} >= 0); |
---|
30 | S{blocks} = S{blocks}*Si; |
---|
31 | vi = S{blocks}'*v{blocks}; |
---|
32 | pdecomp = pdecomp + vi'*Qreduced{blocks}*vi; |
---|
33 | obj = obj + trace(Qreduced{blocks}); |
---|
34 | end |
---|
35 | |
---|
36 | solvesdp(F+set(coefficients(p-pdecomp,recover(depends(p))) == 0),obj,sdpsettings('dualize',1)); |
---|
37 | |
---|
38 | for i = 1:length(Q) |
---|
39 | Q{i} = double(Qreduced{i}); |
---|
40 | end |
---|
41 | |
---|
42 | clear ss |
---|
43 | for i = 1:length(Q) |
---|
44 | ss(i) = norm(Q{i}) > 1e-4; |
---|
45 | end |
---|
46 | Q = {Q{ss}}; |
---|
47 | S = {S{ss}}; |
---|
48 | v = {v{ss}}; |
---|
49 | |
---|
50 | F |
---|
51 | |
---|
52 | |
---|
53 | function [S,index] = blkOne(T) |
---|
54 | |
---|
55 | s = 1; |
---|
56 | S = []; |
---|
57 | clear v |
---|
58 | for i = 1:size(T,1) |
---|
59 | for j = i+1:size(T,1) |
---|
60 | v(i,j) = (T(i,:)*T(j,:)'./(norm(T(i,:))*norm(T(j,:))));%var(T(i,:)./T(j,:))/abs(mean(T(i,:)./T(j,:))); |
---|
61 | end |
---|
62 | end |
---|
63 | v = abs((abs(v)-1)) < 1e-6 ; |
---|
64 | |
---|
65 | top = 1; |
---|
66 | S = []; |
---|
67 | remove = []; |
---|
68 | for i = 1:size(T,1) |
---|
69 | if any(v(:,i)) |
---|
70 | j = find(v(:,i)); |
---|
71 | for k = 1:length(j) |
---|
72 | nn = (norm(T(i,:))/norm(T(j(k),:))); |
---|
73 | s = (T(i,:)*T(j(k),:)'./(norm(T(i,:))*norm(T(j(k),:)))); |
---|
74 | if j(k) <= size(S,2) |
---|
75 | if S(j(k),j(k)) |
---|
76 | S(i,j(k)) = s*(round(10*nn)/10); |
---|
77 | S(i,j(k)) = round(S(i,j(k))*10)/10; |
---|
78 | remove = [remove i]; |
---|
79 | end |
---|
80 | end |
---|
81 | end |
---|
82 | else |
---|
83 | S(i,i) = 1; |
---|
84 | end |
---|
85 | end |
---|
86 | |
---|
87 | S = S(:,any(S,1)); |
---|
88 | index = setdiff(1:size(T,2),remove); |
---|
89 | |
---|
90 | TT = T; |
---|
91 | TT(remove,:) = []; |
---|
92 | TT(:,remove) = []; |
---|
93 | S*TT*S' - T |
---|
94 | |
---|
95 | |
---|
96 | function [Q,v] = prune(Q,v,Qtemp) |
---|
97 | |
---|
98 | for sosfun = 1:length(Q) |
---|
99 | keep = diag(Qtemp{sosfun})>tol; |
---|
100 | Qtemp(:,find(~keep)) = []; |
---|
101 | Qtemp(find(~keep),:) = []; |
---|
102 | |
---|
103 | m{sosfun} = m{sosfun}(find(keep)); |
---|
104 | |
---|
105 | Qtemp(abs(Qtemp) < tol) = 0; |
---|
106 | [v1,dummy1,r1,dummy3]=dmperm(Qtemp+eye(length(Qtemp))); |
---|
107 | lengths{sosfun} = []; |
---|
108 | n{sosfun} = {}; |
---|
109 | for blocks = 1:length(r1)-1 |
---|
110 | i1 = r1(blocks); |
---|
111 | i2 = r1(blocks+1)-1; |
---|
112 | if i2>i1 |
---|
113 | n{sosfun}{blocks} = m{sosfun}(v1(i1:i2)); |
---|
114 | else |
---|
115 | n{sosfun}{blocks} = m{sosfun}(v1(i1)); |
---|
116 | end |
---|
117 | lengths{sosfun} = [lengths{sosfun} length(n{sosfun}{blocks})]; |
---|
118 | end |
---|
119 | lengths{sosfun} = sort(lengths{sosfun}); |
---|
120 | end |
---|
121 | |
---|