1 | function x_extract = extractsolution(momentdata,options) |
---|
2 | %EXTRACTSOLUTIONS Tries to extract solutions from moment matrices |
---|
3 | % |
---|
4 | % |
---|
5 | % xoptimal = extractsolution(momentstructure) |
---|
6 | % |
---|
7 | % xoptimal : Extracted solutions |
---|
8 | % |
---|
9 | % momentdata : Problem data, obtained from SOLVEMOMENT (or SOLVESDP) |
---|
10 | % options : Options structure from SDPSETTINGS |
---|
11 | % |
---|
12 | % See also SOLVEMOMENT, SDPSETTINGS |
---|
13 | |
---|
14 | % Author Johan Löfberg |
---|
15 | % $Id: extractsolution.m,v 1.1 2006/03/30 13:36:38 joloef Exp $ |
---|
16 | |
---|
17 | moment = momentdata.moment; |
---|
18 | x = momentdata.x; |
---|
19 | monomials = momentdata.monomials; |
---|
20 | n = momentdata.n; |
---|
21 | d = momentdata.d; |
---|
22 | |
---|
23 | [U,S,V,ranks] = numranks(moment); |
---|
24 | |
---|
25 | if options.moment.extractrank>0 |
---|
26 | % We try a extraction from highest order moment no matter what |
---|
27 | flat = length(moment); |
---|
28 | ranks(flat) = options.moment.extractrank; |
---|
29 | else |
---|
30 | % Find a flat extension |
---|
31 | flat = d+min(find(ranks(1+d:end)-ranks(1:end-d)==0)); |
---|
32 | end |
---|
33 | |
---|
34 | if ~isempty(flat) |
---|
35 | |
---|
36 | % Find a basis |
---|
37 | r = ranks(flat); |
---|
38 | V = U{flat}(:,1:r)*sqrt(diag(diag(S{flat}(1:r,1:r)))); |
---|
39 | if options.moment.rceftol >= 0 |
---|
40 | [R,pivot] = rref(V',options.moment.rceftol);R = R'; |
---|
41 | else |
---|
42 | % Try to find a reasonable tolerance by avoiding severly badly |
---|
43 | % conditioned R. Hack.., but seem to behave rather robustly. |
---|
44 | cV = cond(V); |
---|
45 | tol = 1e-10; |
---|
46 | [R,pivot] = rref(V',tol);R = R'; |
---|
47 | while tol<1 & (cond(R)/cond(cV)>1e4) |
---|
48 | tol = tol*5; |
---|
49 | [R,pivot] = rref(V',tol);R = R'; |
---|
50 | end |
---|
51 | end |
---|
52 | |
---|
53 | % Figure out multiplying matrices using YALMIP code |
---|
54 | w = monomials(pivot); |
---|
55 | for i = 1:n |
---|
56 | xw = x(i)*w; |
---|
57 | k = []; |
---|
58 | for j = 1:length(xw) |
---|
59 | k = [k;find(ismember(xw(j),monomials))]; |
---|
60 | end |
---|
61 | N{i} = R(k,:); |
---|
62 | end |
---|
63 | |
---|
64 | % Things missing in the basis... |
---|
65 | if ~all(cellfun('prodofsize',N)==length(w)^2) |
---|
66 | x_extract = {[]}; |
---|
67 | return; |
---|
68 | end |
---|
69 | |
---|
70 | % Create random convex combination |
---|
71 | rands = rand(n,1);rands = rands/sum(rands); |
---|
72 | M = 0; |
---|
73 | for i = 1:n |
---|
74 | M = M + rands(i)*N{i}; |
---|
75 | end |
---|
76 | |
---|
77 | [Q,T] = schur(M); |
---|
78 | % Extract solution |
---|
79 | for i = 1:r |
---|
80 | for j = 1:n |
---|
81 | x_extract{i}(j,1) = Q(:,i)'*N{j}*Q(:,i); |
---|
82 | end |
---|
83 | end |
---|
84 | |
---|
85 | % Refine solutions v-Rw = e(x)=0 |
---|
86 | if options.moment.refine>0 |
---|
87 | xtemp = double(x); |
---|
88 | e = monomials(1:size(R,1))-R*w; |
---|
89 | dedx = jacobian(e,x); |
---|
90 | for j = 1:r |
---|
91 | assign(x,x_extract{j}); |
---|
92 | for i = 1:options.moment.refine |
---|
93 | assign(x,double(x)-double(dedx)\double(e)); |
---|
94 | end |
---|
95 | x_extract{j} = double(x); |
---|
96 | end |
---|
97 | assign(x,xtemp); |
---|
98 | end |
---|
99 | |
---|
100 | else |
---|
101 | x_extract = {}; |
---|
102 | end |
---|
103 | |
---|
104 | % Somewhat more stable rank detection |
---|
105 | % looking for sharp drops in singular value |
---|
106 | function [U,S,V,ranks] = numranks(moment) |
---|
107 | for i = 1:length(moment) |
---|
108 | [U{i},S{i},V{i}] = svd(moment{i}); |
---|
109 | s = diag(S{i}); |
---|
110 | decay = s(2:end)./(eps+s(1:end-1)); |
---|
111 | r = min(find(decay<1e-3)); |
---|
112 | if isempty(r) |
---|
113 | ranks(i) = rank(moment{i},1e-8); |
---|
114 | else |
---|
115 | ranks(i) = min(r,rank(moment{i},1e-8)); |
---|
116 | end |
---|
117 | end |
---|
118 | |
---|