[37] | 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 | |
---|