[37] | 1 | function [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,candidateMonomials) |
---|
| 2 | %SOLVESOS Sum of squares decomposition |
---|
| 3 | % |
---|
| 4 | % [sol,m,B,residual] = solvesos(F,h,options,params,monomials) is used |
---|
| 5 | % for finding SOS decompositions of polynomials. |
---|
| 6 | % |
---|
| 7 | % The coefficients of the polynomials are assumed linear w.r.t a set of |
---|
| 8 | % decision variables params and polynomial with respect to a variable x. |
---|
| 9 | % |
---|
| 10 | % An extension with a nonlinear parameterization in params is possible. |
---|
| 11 | % Note though that this gives BMIs or PMIs, solvable (locally) only if |
---|
| 12 | % PENBMI is installed, or by specifying 'moment' as solver to try to |
---|
| 13 | % solve the nonconvex semidefinite programming problem using a |
---|
| 14 | % semidefinite relaxation based on moments. |
---|
| 15 | % |
---|
| 16 | % The SOS problem can be formulated as |
---|
| 17 | % |
---|
| 18 | % min h(params) |
---|
| 19 | % |
---|
| 20 | % subject to F(i) >(=) 0 or F(i) is SOS w.r.t x |
---|
| 21 | % |
---|
| 22 | % INPUT |
---|
| 23 | % F : SET object with SOS constrained polynomials and constraints on variables params |
---|
| 24 | % h : scalar SDPVAR object (can be []) |
---|
| 25 | % options : options structure obtained from SDPSETTINGS (can be []) |
---|
| 26 | % params : SDPVAR object defining parametric variables (can be []) |
---|
| 27 | % monomials : SDPVAR object with user-specified monomials for decomposition (can be []) |
---|
| 28 | % |
---|
| 29 | % OUTPUT |
---|
| 30 | % sol : Solution diagnostic from SDP problem |
---|
| 31 | % v : Cell with monomials used in decompositions |
---|
| 32 | % Q : Cell with Gram matrices, p = v{i}'*Q{i}*v{i}, where p is the ith SOS polynomial in your model. |
---|
| 33 | % residuals : Mismatch between p and decompositions. Same values (modulo numerical issue) as checkset(find(is(F,'sos'))) |
---|
| 34 | % Warning, these residuals are not computed on matrix sos-of-squares |
---|
| 35 | % |
---|
| 36 | % EXAMPLE |
---|
| 37 | % x = sdpvar(1);solvesos(set(sos(x^4+x^3+1))); % Simple decompositions |
---|
| 38 | % x = sdpvar(1);t = sdpvar(1);solvesos(set(sos(x^4+x^3+1-t)),-t); % Lower bound by maximizing t |
---|
| 39 | % |
---|
| 40 | % NOTES |
---|
| 41 | % |
---|
| 42 | % Variables not part of params, but part of non-SOS constraints in F |
---|
| 43 | % or objective h will automatically be appended to the params list. |
---|
| 44 | % |
---|
| 45 | % To extract SOS decomposition, use the command SOSD (or compute from use v and Q) |
---|
| 46 | % |
---|
| 47 | % If the 5th input argument is used, no additional monomial reduction is |
---|
| 48 | % performed (Newton, inconstency, congruence). It is thus assumed that |
---|
| 49 | % the supplied candidate monomials constitute a sufficient basis. |
---|
| 50 | % |
---|
| 51 | % The field options.sos can be used to tune the SOS-calculations. See HTML help for details |
---|
| 52 | % |
---|
| 53 | % sos.model - Kernel or image representation of SOS problem [0|1|2 (0)] |
---|
| 54 | % sos.newton - Use Newton polytope to reduce size [0|1 (1)] |
---|
| 55 | % sos.congruence - Block-diagonalize using congruence classes [0|1|2 (2)] |
---|
| 56 | % sos.scale - Scale polynomial [0|1 (1)] |
---|
| 57 | % sos.numblkdg - Try to perform a-posteriori block-diagonalization [real (0)] |
---|
| 58 | % sos.inconsistent - Remove diagonal-inconsistent monomials [0|1|2 (0)] |
---|
| 59 | % sos.clean - Remove monomials with coefficients < clean [real > 0 (1e-4)] |
---|
| 60 | % sos.traceobj - Minimize trace of Gram matrix in problems without objective function [0|1 (0)] |
---|
| 61 | % sos.extlp - Extract simple translated LP cones when performing dualization [0|1 (1)] |
---|
| 62 | % |
---|
| 63 | % See also SOS, SOSD, SDPSETTINGS, SOLVEMOMENT, SDPVAR, SDISPLAY |
---|
| 64 | |
---|
| 65 | %% Time YALMIP |
---|
| 66 | yalmip_time = clock; |
---|
| 67 | |
---|
| 68 | % ************************************************ |
---|
| 69 | %% Check #inputs |
---|
| 70 | % ************************************************ |
---|
| 71 | if nargin<5 |
---|
| 72 | candidateMonomials = []; |
---|
| 73 | if nargin<4 |
---|
| 74 | params = []; |
---|
| 75 | if nargin<3 |
---|
| 76 | options = sdpsettings; |
---|
| 77 | if nargin<2 |
---|
| 78 | obj = []; |
---|
| 79 | if nargin<1 |
---|
| 80 | help solvesos |
---|
| 81 | return |
---|
| 82 | end |
---|
| 83 | end |
---|
| 84 | end |
---|
| 85 | end |
---|
| 86 | end |
---|
| 87 | |
---|
| 88 | if isempty(options) |
---|
| 89 | options = sdpsettings; |
---|
| 90 | end |
---|
| 91 | |
---|
| 92 | % Lazy syntax (not official...) |
---|
| 93 | if nargin==1 & isa(F,'sdpvar') |
---|
| 94 | F = set(sos(F)); |
---|
| 95 | end |
---|
| 96 | |
---|
| 97 | if ~isempty(options) |
---|
| 98 | if options.sos.numblkdg |
---|
| 99 | [sol,m,Q,residuals,everything] = solvesos_find_blocks(F,obj,options,params,candidateMonomials); |
---|
| 100 | return |
---|
| 101 | end |
---|
| 102 | end |
---|
| 103 | |
---|
| 104 | % ************************************************************************* |
---|
| 105 | %% Extract all SOS constraints and candidate monomials |
---|
| 106 | % ************************************************************************* |
---|
| 107 | if ~any(is(F,'sos')) |
---|
| 108 | error('At-least one constraint should be an SOS constraints!'); |
---|
| 109 | end |
---|
| 110 | p = []; |
---|
| 111 | ranks = []; |
---|
| 112 | for i = 1:length(F) |
---|
| 113 | if is(F(i),'sos') |
---|
| 114 | pi = sdpvar(F(i)); |
---|
| 115 | p{end+1} = pi; |
---|
| 116 | ranks(end+1) = getsosrank(pi); % Desired rank : Experimental code |
---|
| 117 | end |
---|
| 118 | end |
---|
| 119 | if isempty(candidateMonomials) |
---|
| 120 | for i = 1:length(F) |
---|
| 121 | candidateMonomials{i}=[]; |
---|
| 122 | end |
---|
| 123 | elseif isa(candidateMonomials,'sdpvar') |
---|
| 124 | cM=candidateMonomials; |
---|
| 125 | candidateMonomials={}; |
---|
| 126 | for i = 1:length(p) |
---|
| 127 | candidateMonomials{i}=cM; |
---|
| 128 | end |
---|
| 129 | elseif isa(candidateMonomials,'cell') |
---|
| 130 | if length(p)~=length(candidateMonomials) |
---|
| 131 | error('Dimension mismatch between the candidate monomials and the number of SOS constraints'); |
---|
| 132 | end |
---|
| 133 | end |
---|
| 134 | |
---|
| 135 | % ************************************************************************* |
---|
| 136 | %% Get the parametric constraints |
---|
| 137 | % ************************************************************************* |
---|
| 138 | F_original = F; |
---|
| 139 | F_parametric = F(find(~is(F,'sos'))); |
---|
| 140 | if isempty(F_parametric) |
---|
| 141 | F_parametric = set([]); |
---|
| 142 | end |
---|
| 143 | |
---|
| 144 | % ************************************************************************* |
---|
| 145 | %% Expand the parametric constraints |
---|
| 146 | % ************************************************************************* |
---|
| 147 | if ~isempty(yalmip('extvariables')) |
---|
| 148 | [F_parametric,failure] = expandmodel(F_parametric,obj,options); |
---|
| 149 | F_parametric = expanded(F_parametric,1); |
---|
| 150 | obj = expanded(obj,1); |
---|
| 151 | if failure |
---|
| 152 | Q{1} = [];m{1} = [];residuals = [];everything = []; |
---|
| 153 | sol.yalmiptime = etime(clock,yalmip_time); |
---|
| 154 | sol.solvertime = 0; |
---|
| 155 | sol.info = yalmiperror(14,'YALMIP'); |
---|
| 156 | sol.problem = 14; |
---|
| 157 | end |
---|
| 158 | end |
---|
| 159 | |
---|
| 160 | if ~isempty(params) |
---|
| 161 | if ~isa(params,'sdpvar') |
---|
| 162 | error('Fourth argment should be a SDPVAR variable or empty') |
---|
| 163 | end |
---|
| 164 | end |
---|
| 165 | |
---|
| 166 | % ************************************************************************* |
---|
| 167 | % Collect all possible parametric variables |
---|
| 168 | % ************************************************************************* |
---|
| 169 | ParametricVariables = uniquestripped([depends(obj) depends(F_parametric) depends(params)]); |
---|
| 170 | |
---|
| 171 | if options.verbose>0; |
---|
| 172 | disp('-------------------------------------------------------------------------'); |
---|
| 173 | disp('YALMIP SOS module started...'); |
---|
| 174 | disp('-------------------------------------------------------------------------'); |
---|
| 175 | end |
---|
| 176 | |
---|
| 177 | % ************************************************************************* |
---|
| 178 | %% INITIALIZE SOS-DECOMPOSITIONS SDP CONSTRAINTS |
---|
| 179 | % ************************************************************************* |
---|
| 180 | F_sos = set([]); |
---|
| 181 | |
---|
| 182 | % ************************************************************************* |
---|
| 183 | %% FIGURE OUT ALL USED PARAMETRIC VARIABLES |
---|
| 184 | % ************************************************************************* |
---|
| 185 | AllVariables = uniquestripped([depends(obj) depends(F_original) depends(F_parametric)]); |
---|
| 186 | ParametricVariables = intersect(ParametricVariables,AllVariables); |
---|
| 187 | MonomVariables = setdiff(AllVariables,ParametricVariables); |
---|
| 188 | params = recover(ParametricVariables); |
---|
| 189 | if isempty(MonomVariables) |
---|
| 190 | error('No independent variables? Perhaps you added a constraint set(p(x)) when you meant set(sos(p(x)))'); |
---|
| 191 | end |
---|
| 192 | if options.verbose>0;disp(['Detected ' num2str(length(ParametricVariables)) ' parametric variables and ' num2str(length(MonomVariables)) ' independent variables.']);end |
---|
| 193 | |
---|
| 194 | % ************************************************ |
---|
| 195 | %% ANY BMI STUFF |
---|
| 196 | % ************************************************ |
---|
| 197 | NonLinearParameterization = 0; |
---|
| 198 | if ~isempty(ParametricVariables) |
---|
| 199 | monomtable = yalmip('monomtable'); |
---|
| 200 | ParametricMonomials = monomtable(uniquestripped([getvariables(obj) getvariables(F_original)]),ParametricVariables); |
---|
| 201 | if any(sum(abs(ParametricMonomials),2)>1) |
---|
| 202 | NonLinearParameterization = 1; |
---|
| 203 | end |
---|
| 204 | end |
---|
| 205 | |
---|
| 206 | % ************************************************ |
---|
| 207 | %% ANY INTEGER DATA |
---|
| 208 | % ************************************************ |
---|
| 209 | IntegerData = 0; |
---|
| 210 | if ~isempty(ParametricVariables) |
---|
| 211 | globalInteger = [yalmip('binvariables') yalmip('intvariables')]; |
---|
| 212 | integerVariables = getvariables(F_parametric(find(is(F_parametric,'binary') | is(F_parametric,'integer')))); |
---|
| 213 | integerVariables = [integerVariables intersect(ParametricVariables,globalInteger)]; |
---|
| 214 | integerVariables = intersect(integerVariables,ParametricVariables); |
---|
| 215 | IntegerData = ~isempty(integerVariables); |
---|
| 216 | end |
---|
| 217 | |
---|
| 218 | % ************************************************ |
---|
| 219 | %% ANY UNCERTAIN DATA |
---|
| 220 | % ************************************************ |
---|
| 221 | UncertainData = 0; |
---|
| 222 | if ~isempty(ParametricVariables) |
---|
| 223 | UncertainData = any(is(F_parametric,'uncertain')); |
---|
| 224 | end |
---|
| 225 | |
---|
| 226 | % ************************************************ |
---|
| 227 | %% DISPLAY WHAT WE FOUND |
---|
| 228 | % ************************************************ |
---|
| 229 | if options.verbose>0 & ~isempty(F_parametric) |
---|
| 230 | nLP = 0; |
---|
| 231 | nEQ = 0; |
---|
| 232 | nLMI = sum(full(is(F_parametric,'lmi')) & full(~is(F_parametric,'element-wise'))); %FULL due to bug in ML 7.0.1 |
---|
| 233 | for i = 1:length(F_parametric) |
---|
| 234 | if is(F_parametric,'element-wise') |
---|
| 235 | nLP = nLP + prod(size(F_parametric(i))); |
---|
| 236 | end |
---|
| 237 | if is(F_parametric,'equality') |
---|
| 238 | nEQ = nEQ + prod(size(F_parametric(i))); |
---|
| 239 | end |
---|
| 240 | end |
---|
| 241 | disp(['Detected ' num2str(full(nLP)) ' linear inequalities, ' num2str(full(nEQ)) ' equality constraints and ' num2str(full(nLMI)) ' LMIs.']); |
---|
| 242 | end |
---|
| 243 | |
---|
| 244 | % ************************************************ |
---|
| 245 | %% IMAGE OR KERNEL REPRESENTATION? |
---|
| 246 | % ************************************************ |
---|
| 247 | noRANK = all(isinf(ranks)); |
---|
| 248 | switch options.sos.model |
---|
| 249 | case 0 |
---|
| 250 | constraint_classes = constraintclass(F); |
---|
| 251 | noCOMPLICATING = ~any(ismember([7 8 9 10 12 13 14 15],constraint_classes)); |
---|
| 252 | if noCOMPLICATING & ~NonLinearParameterization & noRANK & ~IntegerData |
---|
| 253 | options.sos.model = 1; |
---|
| 254 | if options.verbose>0;disp('Using kernel representation (options.sos.model=1).');end |
---|
| 255 | else |
---|
| 256 | if NonLinearParameterization |
---|
| 257 | if options.verbose>0;disp('Using image representation (options.sos.model=2). Nonlinear parameterization found');end |
---|
| 258 | elseif ~noRANK |
---|
| 259 | if options.verbose>0;disp('Using image representation (options.sos.model=2). SOS-rank constraint was found.');end |
---|
| 260 | elseif IntegerData |
---|
| 261 | if options.verbose>0;disp('Using image representation (options.sos.model=2). Integrality constraint was found.');end |
---|
| 262 | elseif UncertainData |
---|
| 263 | if options.verbose>0;disp('Using image representation (options.sos.model=2). Uncertain data was found.');end |
---|
| 264 | else |
---|
| 265 | if options.verbose>0;disp('Using image representation (options.sos.model=2). Integer data, KYPs or similar was found.');end |
---|
| 266 | end |
---|
| 267 | options.sos.model = 2; |
---|
| 268 | end |
---|
| 269 | case 1 |
---|
| 270 | if NonLinearParameterization |
---|
| 271 | if options.verbose>0;disp('Switching to image model due to nonlinear parameterization (not supported in kernel model).');end |
---|
| 272 | options.sos.model = 2; |
---|
| 273 | end |
---|
| 274 | if ~noRANK |
---|
| 275 | if options.verbose>0;disp('Switching to image model due to SOS-rank constraints (not supported in kernel model).');end |
---|
| 276 | options.sos.model = 2; |
---|
| 277 | end |
---|
| 278 | if IntegerData |
---|
| 279 | if options.verbose>0;disp('Switching to image model due to integrality constraints (not supported in kernel model).');end |
---|
| 280 | options.sos.model = 2; |
---|
| 281 | end |
---|
| 282 | case 3 |
---|
| 283 | otherwise |
---|
| 284 | end |
---|
| 285 | |
---|
| 286 | if ~isempty(yalmip('extvariables')) & options.sos.model == 2 & nargin<4 |
---|
| 287 | disp(' ') |
---|
| 288 | disp('**Using nonlinear operators in SOS problems can cause problems.') |
---|
| 289 | disp('**Please specify all parametric variables using the fourth argument'); |
---|
| 290 | disp(' '); |
---|
| 291 | end |
---|
| 292 | |
---|
| 293 | % ************************************************ |
---|
| 294 | %% SKIP DIAGONAL INCONSISTENCY FOR PARAMETRIC MODEL |
---|
| 295 | % ************************************************ |
---|
| 296 | if ~isempty(params) & options.sos.inconsistent |
---|
| 297 | if options.verbose>0;disp('Turning off inconsistency based reduction (not supported in parametric models).');end |
---|
| 298 | options.sos.inconsistent = 0; |
---|
| 299 | end |
---|
| 300 | |
---|
| 301 | % ************************************************ |
---|
| 302 | %% INITIALIZE OBJECTIVE |
---|
| 303 | % ************************************************ |
---|
| 304 | if ~isempty(obj) |
---|
| 305 | options.sos.traceobj = 0; |
---|
| 306 | end |
---|
| 307 | parobj = obj; |
---|
| 308 | obj = []; |
---|
| 309 | |
---|
| 310 | % ************************************************ |
---|
| 311 | %% SCALE SOS CONSTRAINTS |
---|
| 312 | % ************************************************ |
---|
| 313 | if options.sos.scale |
---|
| 314 | for constraint = 1:length(p) |
---|
| 315 | normp(constraint) = sqrt(norm(full(getbase(p{constraint})))); |
---|
| 316 | p{constraint} = p{constraint}/normp(constraint); |
---|
| 317 | sizep(constraint) = size(p{constraint},1); |
---|
| 318 | end |
---|
| 319 | else |
---|
| 320 | normp = ones(length(p),1); |
---|
| 321 | end |
---|
| 322 | |
---|
| 323 | % ************************************************ |
---|
| 324 | %% Some stuff not supported for |
---|
| 325 | % matrix valued SOS yet, turn off for safety |
---|
| 326 | % ************************************************ |
---|
| 327 | for constraint = 1:length(p) |
---|
| 328 | sizep(constraint) = size(p{constraint},1); |
---|
| 329 | end |
---|
| 330 | if any(sizep>1) |
---|
| 331 | options.sos.postprocess = 0; |
---|
| 332 | options.sos.reuse = 0; |
---|
| 333 | end |
---|
| 334 | |
---|
| 335 | % ************************************************ |
---|
| 336 | %% SKIP CONGRUENCE REDUCTION WHEN SOS-RANK |
---|
| 337 | % ************************************************ |
---|
| 338 | if ~all(isinf(ranks)) |
---|
| 339 | options.sos.congruence = 0; |
---|
| 340 | end |
---|
| 341 | |
---|
| 342 | % ************************************************ |
---|
| 343 | %% Create an LP model to speed up things in Newton |
---|
| 344 | % polytope reduction |
---|
| 345 | % ************************************************ |
---|
| 346 | if options.sos.newton |
---|
| 347 | temp=sdpvar(1,1); |
---|
| 348 | tempops = options; |
---|
| 349 | tempops.solver = 'cdd,glpk,*'; % CDD is generally robust on these problems |
---|
| 350 | tempops.verbose = 0; |
---|
| 351 | tempops.saveduals = 0; |
---|
| 352 | [aux1,aux2,aux3,LPmodel] = export(set(temp>0),temp,tempops); |
---|
| 353 | else |
---|
| 354 | LPmodel = []; |
---|
| 355 | end |
---|
| 356 | |
---|
| 357 | |
---|
| 358 | % ************************************************ |
---|
| 359 | %% LOOP THROUGH ALL SOS CONSTRAINTS |
---|
| 360 | % ************************************************ |
---|
| 361 | for constraint = 1:length(p) |
---|
| 362 | % ********************************************* |
---|
| 363 | %% FIND THE VARIABLES IN p, SORT, UNIQUE ETC |
---|
| 364 | % ********************************************* |
---|
| 365 | if options.verbose>1;disp(['Creating SOS-description ' num2str(constraint) '/' num2str(length(p)) ]);end |
---|
| 366 | |
---|
| 367 | pVariables = depends(p{constraint}); |
---|
| 368 | AllVariables = uniquestripped([pVariables ParametricVariables]); |
---|
| 369 | MonomVariables = setdiff1D(pVariables,ParametricVariables); |
---|
| 370 | x = recover(MonomVariables); |
---|
| 371 | z = recover(AllVariables); |
---|
| 372 | MonomIndicies = find(ismember(AllVariables,MonomVariables)); |
---|
| 373 | ParametricIndicies = find(ismember(AllVariables,ParametricVariables)); |
---|
| 374 | |
---|
| 375 | if isempty(MonomIndicies) |
---|
| 376 | % This is the case set(sos(t)) where t is a parametric (matrix) variable |
---|
| 377 | % This used to create an error message befgore to avoid some silly |
---|
| 378 | % bug in the model generation. Creating this error message is |
---|
| 379 | % stupid, but at the same time I can not remember where the bug was |
---|
| 380 | % and I have no regression test for this case. To avoid |
---|
| 381 | % introducing same bug again by mistake, I create all data |
---|
| 382 | % specifically for this case |
---|
| 383 | previous_exponent_p_monoms = [];%exponent_p_monoms; |
---|
| 384 | n = length(p{constraint}); |
---|
| 385 | A_basis = getbase(sdpvar(n,n,'full'));d = find(triu(ones(n)));A_basis = A_basis(d,2:end); |
---|
| 386 | BlockedA{constraint} = {A_basis}; |
---|
| 387 | Blockedb{constraint} = p{constraint}(d); |
---|
| 388 | BlockedN{constraint} = {zeros(1,0)}; |
---|
| 389 | Blockedx{constraint} = x; |
---|
| 390 | Blockedvarchange{constraint}=zeros(1,0); |
---|
| 391 | continue |
---|
| 392 | % error('You have constraints of the type set(sos(f(parametric_variables))). Please use set(f(parametric_variables) > 0) instead') |
---|
| 393 | end |
---|
| 394 | |
---|
| 395 | % ********************************************* |
---|
| 396 | %% Express p in monimials and coefficients |
---|
| 397 | % ********************************************* |
---|
| 398 | [exponent_p,p_base] = getexponentbase(p{constraint},z); |
---|
| 399 | |
---|
| 400 | % ********************************************* |
---|
| 401 | %% Powers for user defined candidate monomials |
---|
| 402 | % (still experimental) |
---|
| 403 | % ********************************************* |
---|
| 404 | if ~all(cellfun('isempty',candidateMonomials)) |
---|
| 405 | exponent_c = []; |
---|
| 406 | if isa(candidateMonomials{constraint},'cell') |
---|
| 407 | for i = 1:length(candidateMonomials{constraint}) |
---|
| 408 | exponent_c{i} = getexponentbase(candidateMonomials{constraint}{i},z); |
---|
| 409 | exponent_c{i} = exponent_c{i}(:,MonomIndicies); |
---|
| 410 | end |
---|
| 411 | else |
---|
| 412 | exponent_c{1} = getexponentbase(candidateMonomials{constraint},z); |
---|
| 413 | exponent_c{1} = exponent_c{1}(:,MonomIndicies); |
---|
| 414 | end |
---|
| 415 | else |
---|
| 416 | exponent_c = []; |
---|
| 417 | end |
---|
| 418 | |
---|
| 419 | % ********************************************* |
---|
| 420 | %% STUPID PROBLEM WITH ODD HIGHEST POWER?... |
---|
| 421 | % ********************************************* |
---|
| 422 | if isempty(ParametricIndicies) |
---|
| 423 | max_degrees = max(exponent_p(:,MonomIndicies),[],1); |
---|
| 424 | bad_max = any(max_degrees-fix((max_degrees/2))*2); |
---|
| 425 | if bad_max |
---|
| 426 | for i = 1:length(p) |
---|
| 427 | Q{i}=[]; |
---|
| 428 | m{i}=[]; |
---|
| 429 | end |
---|
| 430 | residuals=[]; |
---|
| 431 | everything = []; |
---|
| 432 | sol.yalmiptime = etime(clock,yalmip_time); |
---|
| 433 | sol.solvertime = 0; |
---|
| 434 | sol.info = yalmiperror(1,'YALMIP'); |
---|
| 435 | sol.problem = 2; |
---|
| 436 | return |
---|
| 437 | end |
---|
| 438 | end |
---|
| 439 | |
---|
| 440 | % ********************************************* |
---|
| 441 | %% Can we make a smart variable change (no code) |
---|
| 442 | % ********************************************* |
---|
| 443 | exponent_p_monoms = exponent_p(:,MonomIndicies); |
---|
| 444 | varchange = ones(1,size(MonomIndicies,2)); |
---|
| 445 | |
---|
| 446 | % ********************************************* |
---|
| 447 | %% Unique monoms (copies due to parametric terms) |
---|
| 448 | % ********************************************* |
---|
| 449 | exponent_p_monoms = uniquesafe(exponent_p_monoms,'rows'); |
---|
| 450 | |
---|
| 451 | if options.sos.reuse & constraint > 1 & isequal(previous_exponent_p_monoms,exponent_p_monoms) |
---|
| 452 | % We don't have to do anything, candidate monomials can be-used |
---|
| 453 | if options.verbose>1;disp(['Re-using all candidate monomials (same problem structure)']);end |
---|
| 454 | else |
---|
| 455 | |
---|
| 456 | |
---|
| 457 | % ********************************************* |
---|
| 458 | %% PRUNE W.R.T USER DEFINED |
---|
| 459 | %********************************************* |
---|
| 460 | % if ~isempty(exponent_c) |
---|
| 461 | % hash = randn(size(exponent_m{1},2),1); |
---|
| 462 | % for i = 1:length(exponent_m) |
---|
| 463 | % hash = randn(size(exponent_m{i},2),1); |
---|
| 464 | % keep = find(ismember(exponent_m{i}*hash,exponent_c{1}*hash)); |
---|
| 465 | % if length(keep)>0 |
---|
| 466 | % exponent_m{i} = exponent_m{i}(keep,:); |
---|
| 467 | % else |
---|
| 468 | % exponent_m{i}=[]; |
---|
| 469 | % end |
---|
| 470 | % end |
---|
| 471 | % end |
---|
| 472 | |
---|
| 473 | % ********************************************* |
---|
| 474 | % User has supplied the whole candidate structure |
---|
| 475 | % Don't process this |
---|
| 476 | % ********************************************* |
---|
| 477 | if ~isempty(exponent_c) |
---|
| 478 | exponent_m{1} = []; |
---|
| 479 | N = {}; |
---|
| 480 | for i = 1:length(exponent_c) |
---|
| 481 | exponent_m{i} = [exponent_m{1};exponent_c{i}]; |
---|
| 482 | N{i,1} = exponent_c{i}; |
---|
| 483 | end |
---|
| 484 | else |
---|
| 485 | % ********************************************* |
---|
| 486 | %% CORRELATIVE SPARSITY PATTERN |
---|
| 487 | % ********************************************* |
---|
| 488 | [C,csclasses] = corrsparsity(exponent_p_monoms,options); |
---|
| 489 | |
---|
| 490 | % ********************************************* |
---|
| 491 | %% GENERATE MONOMIALS |
---|
| 492 | % ********************************************* |
---|
| 493 | exponent_m = monomialgeneration(exponent_p_monoms,csclasses); |
---|
| 494 | |
---|
| 495 | % ********************************************* |
---|
| 496 | %% REDUCE #of MONOMIALS |
---|
| 497 | % ********************************************* |
---|
| 498 | % Fix for matrix case, perform newton w.r.t |
---|
| 499 | % diagonal polynomials only. This can be |
---|
| 500 | % improved, but for now, keep it simple... |
---|
| 501 | n = length(p{constraint});diag_elements = 1:(n+1):n^2;used_diagonal = find(any(p_base(diag_elements,:),1)); |
---|
| 502 | exponent_p_monoms_diag = exponent_p(used_diagonal,MonomIndicies); |
---|
| 503 | exponent_m = monomialreduction(exponent_m,exponent_p_monoms_diag,options,csclasses,LPmodel); |
---|
| 504 | |
---|
| 505 | % ********************************************* |
---|
| 506 | %% BLOCK PARTITION THE MONOMIALS BY CONGRUENCE |
---|
| 507 | % ********************************************* |
---|
| 508 | N = congruenceblocks(exponent_m,exponent_p_monoms,options,csclasses); |
---|
| 509 | % ********************************************* |
---|
| 510 | %% REDUCE FURTHER BY EXPLOITING BLOCK-STRUCTURE |
---|
| 511 | % ********************************************* |
---|
| 512 | N = blockmonomialreduction(exponent_p_monoms_diag,N,options); |
---|
| 513 | |
---|
| 514 | end |
---|
| 515 | |
---|
| 516 | |
---|
| 517 | % ********************************************* |
---|
| 518 | %% PREPARE FOR SDP FORMULATION BY CALCULATING ALL |
---|
| 519 | % POSSIBLE MONOMIAL PRODUCS IN EACH BLOCK |
---|
| 520 | % ********************************************* |
---|
| 521 | [exponent_m2,N_unique] = monomialproducts(N); |
---|
| 522 | |
---|
| 523 | % ********************************************* |
---|
| 524 | %% CHECK FOR BUG/IDIOT PROBLEMS IN FIXED PROBLEM |
---|
| 525 | % ********************************************* |
---|
| 526 | if isempty(ParametricIndicies) |
---|
| 527 | if ~isempty(setdiff(exponent_p_monoms,N_unique(:,3:end),'rows')) |
---|
| 528 | for i = 1:length(p) |
---|
| 529 | Q{i} = []; |
---|
| 530 | m{i} = []; |
---|
| 531 | end |
---|
| 532 | residuals = [];everything = []; |
---|
| 533 | sol.problem = 2; |
---|
| 534 | sol.info = yalmiperror(1,'YALMIP'); |
---|
| 535 | warning('Problem is trivially infeasible (odd highest power?)'); |
---|
| 536 | return |
---|
| 537 | end |
---|
| 538 | end |
---|
| 539 | end |
---|
| 540 | |
---|
| 541 | previous_exponent_p_monoms = exponent_p_monoms; |
---|
| 542 | |
---|
| 543 | % ********************************************* |
---|
| 544 | %% GENERATE DATA FOR SDP FORMULATIONS |
---|
| 545 | % ********************************************* |
---|
| 546 | p_base_parametric = []; |
---|
| 547 | n = length(p{constraint}); |
---|
| 548 | for i=1:length(p{constraint}) |
---|
| 549 | for j = 1:length(p{constraint}) |
---|
| 550 | p_base_parametric = [p_base_parametric parameterizedbase(p{constraint}(i,j),z,params,ParametricIndicies,exponent_p,p_base((i-1)*n+j,:))]; |
---|
| 551 | end |
---|
| 552 | end |
---|
| 553 | [BlockedA{constraint},Blockedb{constraint}] = generate_kernel_representation_data(N,N_unique,exponent_m2,exponent_p,p{constraint},options,p_base_parametric,ParametricIndicies,MonomIndicies,constraint==1); |
---|
| 554 | |
---|
| 555 | % SAVE FOR LATER |
---|
| 556 | BlockedN{constraint} = N; |
---|
| 557 | Blockedx{constraint} = x; |
---|
| 558 | Blockedvarchange{constraint}=varchange; |
---|
| 559 | end |
---|
| 560 | |
---|
| 561 | % ********************************************* |
---|
| 562 | %% And now get the SDP formulations |
---|
| 563 | % |
---|
| 564 | % The code above has generated matrices A and b |
---|
| 565 | % in AQ == b(parametric) |
---|
| 566 | % |
---|
| 567 | % We use these to generate kernel or image models |
---|
| 568 | % ********************************************* |
---|
| 569 | sol.problem = 0; |
---|
| 570 | switch options.sos.model |
---|
| 571 | case 1 |
---|
| 572 | % Kernel model |
---|
| 573 | [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,[]); |
---|
| 574 | case 2 |
---|
| 575 | % Image model |
---|
| 576 | [F,obj,BlockedQ,sol] = create_imagemodel(BlockedA,Blockedb,F_parametric,parobj,options); |
---|
| 577 | |
---|
| 578 | case 3 |
---|
| 579 | % Un-official model to solve bilinearly parameterized SOS using SDPLR |
---|
| 580 | [F,obj,options] = create_lrmodel(BlockedA,Blockedb,F_parametric,parobj,options,ParametricVariables); |
---|
| 581 | |
---|
| 582 | otherwise |
---|
| 583 | end |
---|
| 584 | |
---|
| 585 | % Unofficial fifth output with pseudo-numerical data |
---|
| 586 | everything.BlockedA = BlockedA; |
---|
| 587 | everything.Blockedb = Blockedb; |
---|
| 588 | everything.F = F; |
---|
| 589 | everything.obj = obj; |
---|
| 590 | |
---|
| 591 | % % ********************************************** |
---|
| 592 | % %% SOLVE SDP |
---|
| 593 | % % ********************************************** |
---|
| 594 | if sol.problem == 0 |
---|
| 595 | if options.verbose > 0 |
---|
| 596 | disp(' '); |
---|
| 597 | end |
---|
| 598 | if all(isinf(ranks)) |
---|
| 599 | % Standard case |
---|
| 600 | %sol = solvesdp(F+set(-10<recover(depends(F)) < 10),obj,sdpsettings('bmibnb.maxit',200,'solver','bmibnb','bmibnb.lpred',1)) |
---|
| 601 | sol = solvesdp(F,obj,options); |
---|
| 602 | else |
---|
| 603 | % We have to alter the problem slightly if there are rank |
---|
| 604 | % constraints on the decompositions |
---|
| 605 | sol = solveranksos(F,obj,options,ranks,BlockedQ); |
---|
| 606 | end |
---|
| 607 | end |
---|
| 608 | |
---|
| 609 | % ********************************************** |
---|
| 610 | %% Recover solution (and rescale model+solution) |
---|
| 611 | % ********************************************** |
---|
| 612 | for constraint = 1:length(p) |
---|
| 613 | for i = 1:length(BlockedQ{constraint}) |
---|
| 614 | doubleQ{constraint}{i} = normp(constraint)*double(BlockedQ{constraint}{i}); |
---|
| 615 | end |
---|
| 616 | doubleb{constraint}=normp(constraint)*double(Blockedb{constraint}); |
---|
| 617 | end |
---|
| 618 | |
---|
| 619 | % ********************************************** |
---|
| 620 | %% Minor post-process |
---|
| 621 | % ********************************************** |
---|
| 622 | if all(sizep<=1) |
---|
| 623 | [doubleQ,residuals] = postprocesssos(BlockedA,doubleb,doubleQ,[],options); |
---|
| 624 | else |
---|
| 625 | residuals = 0; |
---|
| 626 | end |
---|
| 627 | |
---|
| 628 | % ********************************************** |
---|
| 629 | %% Safety check for bad solvers... |
---|
| 630 | % ********************************************** |
---|
| 631 | if any(residuals > 1e-3) & (sol.problem == 0) & options.verbose>0 |
---|
| 632 | disp(' '); |
---|
| 633 | disp('-> Although the solver indicates no problems,') |
---|
| 634 | disp('-> the residuals in the problem are really bad.') |
---|
| 635 | disp('-> My guess: the problem is probably infeasible.') |
---|
| 636 | disp('-> Make sure to check how well your decomposition') |
---|
| 637 | disp('-> matches your polynomial (see manual)') |
---|
| 638 | disp('-> You can also try to change the option sos.model') |
---|
| 639 | disp('-> or use another SDP solver.') |
---|
| 640 | disp(' '); |
---|
| 641 | end |
---|
| 642 | |
---|
| 643 | % ********************************************** |
---|
| 644 | %% Confused users. Primal dual kernel image?... |
---|
| 645 | % ********************************************** |
---|
| 646 | if options.verbose > 0 |
---|
| 647 | if sol.problem == 1 |
---|
| 648 | if options.sos.model == 1 |
---|
| 649 | disp(' ') |
---|
| 650 | disp('-> Solver reported infeasible dual problem.') |
---|
| 651 | disp('-> Your SOS problem is probably unbounded.') |
---|
| 652 | elseif options.sos.model == 2 |
---|
| 653 | disp(' ') |
---|
| 654 | disp('-> Solver reported infeasible primal problem.') |
---|
| 655 | disp('-> Your SOS problem is probably infeasible.') |
---|
| 656 | end |
---|
| 657 | elseif sol.problem == 2 |
---|
| 658 | if options.sos.model == 1 |
---|
| 659 | disp(' ') |
---|
| 660 | disp('-> Solver reported unboundness of the dual problem.') |
---|
| 661 | disp('-> Your SOS problem is probably infeasible.') |
---|
| 662 | elseif options.sos.model == 2 |
---|
| 663 | disp(' ') |
---|
| 664 | disp('-> Solver reported unboundness of the primal problem.') |
---|
| 665 | disp('-> Your SOS problem is probably unbounded.') |
---|
| 666 | end |
---|
| 667 | elseif sol.problem == 12 |
---|
| 668 | disp(' ') |
---|
| 669 | disp('-> Solver reported unboundness or infeasibility of the primal problem.') |
---|
| 670 | disp('-> Your SOS problem is probably unbounded.') |
---|
| 671 | end |
---|
| 672 | end |
---|
| 673 | |
---|
| 674 | % ********************************************** |
---|
| 675 | %% De-block |
---|
| 676 | % ********************************************** |
---|
| 677 | for constraint = 1:length(p) |
---|
| 678 | Qtemp = []; |
---|
| 679 | for i = 1:length(BlockedQ{constraint}) |
---|
| 680 | Qtemp = blkdiag(Qtemp,doubleQ{constraint}{i}); |
---|
| 681 | end |
---|
| 682 | Q{constraint} = full(Qtemp); |
---|
| 683 | end |
---|
| 684 | |
---|
| 685 | % ********************************************** |
---|
| 686 | %% Experimental code for declaring sparsity |
---|
| 687 | % ********************************************** |
---|
| 688 | if options.sos.impsparse == 1 |
---|
| 689 | somesmall = 0; |
---|
| 690 | for i = 1:length(BlockedQ) |
---|
| 691 | for j = 1:length(BlockedQ{i}) |
---|
| 692 | small = find(abs(double(BlockedQ{i}{j}))<options.sos.sparsetol); |
---|
| 693 | nullThese{i}{j} = small; |
---|
| 694 | if ~isempty(small) |
---|
| 695 | somesmall = 1; |
---|
| 696 | end |
---|
| 697 | end |
---|
| 698 | end |
---|
| 699 | if somesmall |
---|
| 700 | [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,nullThese); |
---|
| 701 | sol = solvesdp(F,obj,options); |
---|
| 702 | for constraint = 1:length(p) |
---|
| 703 | for i = 1:length(BlockedQ{constraint}) |
---|
| 704 | doubleQ{constraint}{i} = normp(constraint)*double(BlockedQ{constraint}{i}); |
---|
| 705 | end |
---|
| 706 | doubleb{constraint}=normp(constraint)*double(Blockedb{constraint}); |
---|
| 707 | end |
---|
| 708 | % for i = 1:length(doubleQ) |
---|
| 709 | % for j = 1:length(doubleQ{i}) |
---|
| 710 | % doubleQ{i}{j}(nullThese{i}{j}) = 0; |
---|
| 711 | % end |
---|
| 712 | % end |
---|
| 713 | % ********************************************** |
---|
| 714 | %% Post-process |
---|
| 715 | % ********************************************** |
---|
| 716 | [doubleQ,residuals] = postprocesssos(BlockedA,doubleb,doubleQ,nullThese,options); |
---|
| 717 | for constraint = 1:length(p) |
---|
| 718 | Qtemp = []; |
---|
| 719 | for i = 1:length(BlockedQ{constraint}) |
---|
| 720 | Qtemp = blkdiag(Qtemp,doubleQ{constraint}{i}); |
---|
| 721 | end |
---|
| 722 | Q{constraint} = Qtemp; |
---|
| 723 | end |
---|
| 724 | end |
---|
| 725 | end |
---|
| 726 | |
---|
| 727 | % ********************************************* |
---|
| 728 | %% EXTRACT DECOMPOSITION |
---|
| 729 | % ********************************************* |
---|
| 730 | switch sol.problem |
---|
| 731 | case {0,1,2,3,4,5} % Well, it didn't f**k up completely at-least |
---|
| 732 | |
---|
| 733 | % % If SDPLR was sucessful, the last dual should have rank 1 |
---|
| 734 | % Z = dual(F(end)); |
---|
| 735 | % i = 1e-16; |
---|
| 736 | % [R,fail] = chol(Z); |
---|
| 737 | % while fail & i<10 |
---|
| 738 | % i = i*100; |
---|
| 739 | % [R,fail] = chol(Z+i*eye(length(Z))); |
---|
| 740 | % end |
---|
| 741 | % Z = R(1,2:end)'; |
---|
| 742 | % assign(recover(ParametricVariables),Z); |
---|
| 743 | % start = 1; |
---|
| 744 | % for constraint = 1:length(p) |
---|
| 745 | % Qi = []; |
---|
| 746 | % for j = 1:length(BlockedA{constraint}) |
---|
| 747 | % Qi = blkdiag(Qi,dual(F(start))); |
---|
| 748 | % start = start + 1; |
---|
| 749 | % end |
---|
| 750 | % Q{constraint} = Qi; |
---|
| 751 | % end |
---|
| 752 | |
---|
| 753 | % ********************************************* |
---|
| 754 | %% GENERATE MONOMIALS IN SOS-DECOMPOSITION |
---|
| 755 | % ********************************************* |
---|
| 756 | for constraint = 1:length(p) |
---|
| 757 | |
---|
| 758 | if constraint > 1 & isequal(BlockedN{constraint},BlockedN{constraint-1}) & isequal(Blockedx{constraint},Blockedx{constraint-1}) & isequal(Blockedvarchange{constraint},Blockedvarchange{constraint-1}) & isequal(sizep(constraint),sizep(constraint-1)) |
---|
| 759 | monoms{constraint} = monoms{constraint-1}; |
---|
| 760 | else |
---|
| 761 | monoms{constraint} = []; |
---|
| 762 | totalN{constraint} = []; |
---|
| 763 | N = BlockedN{constraint}; |
---|
| 764 | x = Blockedx{constraint}; |
---|
| 765 | for i = 1:length(N) |
---|
| 766 | % Original variable |
---|
| 767 | for j = 1:size(N{i},1) |
---|
| 768 | N{i}(j,:)=N{i}(j,:).*Blockedvarchange{constraint}; |
---|
| 769 | end |
---|
| 770 | if isempty(N{i}) |
---|
| 771 | monoms{constraint} = [monoms{constraint};[]]; |
---|
| 772 | else |
---|
| 773 | mi = kron(eye(sizep(constraint)),recovermonoms(N{i},x)); |
---|
| 774 | monoms{constraint} = [monoms{constraint};mi]; |
---|
| 775 | end |
---|
| 776 | end |
---|
| 777 | if isempty(monoms{constraint}) |
---|
| 778 | monoms{constraint}=1; |
---|
| 779 | end |
---|
| 780 | end |
---|
| 781 | |
---|
| 782 | % For small negative eigenvalues |
---|
| 783 | % this is a good quick'n'dirty approximation |
---|
| 784 | % Improve...use shifted eigenvalues and chol or what ever... |
---|
| 785 | if ~any(any(isnan(Q{constraint}))) |
---|
| 786 | if isempty(Q{constraint}) |
---|
| 787 | Q{constraint}=0; |
---|
| 788 | h{constraint}=0; |
---|
| 789 | else |
---|
| 790 | usedVariables = find(any(Q{constraint},2)); |
---|
| 791 | if length(usedVariables)<length(Q{constraint}) |
---|
| 792 | Qpart = Q{constraint}(usedVariables,usedVariables); |
---|
| 793 | [U,S,V]=svd(Qpart); |
---|
| 794 | R = sqrt(S)*V'; |
---|
| 795 | h0 = R*monoms{constraint}(usedVariables); |
---|
| 796 | if isa(h0,'sdpvar') |
---|
| 797 | h{constraint} = clean(R*monoms{constraint}(usedVariables),options.sos.clean); |
---|
| 798 | h{constraint} = h{constraint}(find(h{constraint})); |
---|
| 799 | else |
---|
| 800 | h{constraint} = h0; |
---|
| 801 | end |
---|
| 802 | else |
---|
| 803 | [U,S,V]=svd(Q{constraint}); |
---|
| 804 | R = sqrt(S)*V'; |
---|
| 805 | h0 = R*monoms{constraint}; |
---|
| 806 | |
---|
| 807 | if isa(h0,'sdpvar') |
---|
| 808 | h{constraint} = clean(R*monoms{constraint},options.sos.clean); |
---|
| 809 | h{constraint} = h{constraint}(find(sum(h{constraint},2)),:); |
---|
| 810 | else |
---|
| 811 | h{constraint} = h0; |
---|
| 812 | end |
---|
| 813 | end |
---|
| 814 | end |
---|
| 815 | if isempty(ParametricVariables) |
---|
| 816 | ParametricVariables = []; |
---|
| 817 | end |
---|
| 818 | setsos(p{constraint},h{constraint},ParametricVariables,Q{constraint},monoms{constraint}); |
---|
| 819 | else |
---|
| 820 | if options.verbose>0; |
---|
| 821 | if UncertainData |
---|
| 822 | disp(' '); |
---|
| 823 | disp('-> Only partial decomposition is returned (since you have uncertain data).'); |
---|
| 824 | disp(''); |
---|
| 825 | else |
---|
| 826 | disp(' '); |
---|
| 827 | disp('-> FAILURE : SOS decomposition not available.'); |
---|
| 828 | disp('-> The reason is probably that you are using a solver that does not deliver a dual (LMILAB)'); |
---|
| 829 | disp('-> Use sdsettings(''sos.model'',2) to circumvent this, or use another solver (SDPT3, SEDUMI,...)'); |
---|
| 830 | disp(''); |
---|
| 831 | disp('-> An alternative reason is that YALMIP detected infeasibility during the compilation phase.'); |
---|
| 832 | end |
---|
| 833 | end |
---|
| 834 | end |
---|
| 835 | end |
---|
| 836 | |
---|
| 837 | m = monoms; |
---|
| 838 | |
---|
| 839 | otherwise |
---|
| 840 | Q = []; |
---|
| 841 | m = []; |
---|
| 842 | end |
---|
| 843 | |
---|
| 844 | % Don't need these outside |
---|
| 845 | yalmip('cleardual') |
---|
| 846 | |
---|
| 847 | % Done with YALMIP, this is the time it took, minus solver |
---|
| 848 | if ~isfield(sol,'solvertime') |
---|
| 849 | sol.solvertime = 0; |
---|
| 850 | end |
---|
| 851 | |
---|
| 852 | sol.yalmiptime = etime(clock,yalmip_time)-sol.solvertime; |
---|
| 853 | |
---|
| 854 | |
---|
| 855 | function p_base_parametric = parameterizedbase(p,z, params,ParametricIndicies,exponent_p,p_base) |
---|
| 856 | |
---|
| 857 | % Check for linear parameterization |
---|
| 858 | parametric_basis = exponent_p(:,ParametricIndicies); |
---|
| 859 | if all(sum(parametric_basis,2)==0) |
---|
| 860 | p_base_parametric = full(p_base(:)); |
---|
| 861 | return |
---|
| 862 | end |
---|
| 863 | if all(sum(parametric_basis,2)<=1) |
---|
| 864 | p_base_parametric = full(p_base(:)); |
---|
| 865 | n = length(p_base_parametric); |
---|
| 866 | |
---|
| 867 | |
---|
| 868 | if 1 |
---|
| 869 | [ii,vars] = find(parametric_basis); |
---|
| 870 | ii = ii(:)'; |
---|
| 871 | vars = vars(:)'; |
---|
| 872 | else |
---|
| 873 | ii = []; |
---|
| 874 | vars = []; |
---|
| 875 | js = sum(parametric_basis,1); |
---|
| 876 | indicies = find(js); |
---|
| 877 | for i = indicies |
---|
| 878 | if js(i) |
---|
| 879 | j = find(parametric_basis(:,i)); |
---|
| 880 | ii = [ii j(:)']; |
---|
| 881 | vars = [vars repmat(i,1,js(i))]; |
---|
| 882 | end |
---|
| 883 | end |
---|
| 884 | end |
---|
| 885 | |
---|
| 886 | k = setdiff1D(1:n,ii); |
---|
| 887 | if isempty(k) |
---|
| 888 | p_base_parametric = p_base_parametric.*sparse(ii,repmat(1,1,n),params(vars)); |
---|
| 889 | else |
---|
| 890 | pp = params(vars); % Must do this, bug in ML 6.1 (x=sparse(1);x([1 1]) gives different result in 6.1 and 7.0!) |
---|
| 891 | p_base_parametric = p_base_parametric.*sparse([ii k(:)'],repmat(1,1,n),[pp(:)' ones(1,1,length(k))]); |
---|
| 892 | end |
---|
| 893 | else |
---|
| 894 | % Bummer, nonlinear parameterization sucks... |
---|
| 895 | for i = 1:length(p_base) |
---|
| 896 | j = find(exponent_p(i,ParametricIndicies)); |
---|
| 897 | if ~isempty(j) |
---|
| 898 | temp = p_base(i); |
---|
| 899 | |
---|
| 900 | for k = 1:length(j) |
---|
| 901 | if exponent_p(i,ParametricIndicies(j(k)))==1 |
---|
| 902 | temp = temp*params(j(k)); |
---|
| 903 | else |
---|
| 904 | temp = temp*params(j(k))^exponent_p(i,ParametricIndicies(j(k))); |
---|
| 905 | end |
---|
| 906 | end |
---|
| 907 | xx{i} = temp; |
---|
| 908 | else |
---|
| 909 | xx{i} = p_base(i); |
---|
| 910 | end |
---|
| 911 | end |
---|
| 912 | p_base_parametric = stackcell(sdpvar(1,1),xx)'; |
---|
| 913 | end |
---|
| 914 | |
---|
| 915 | |
---|
| 916 | |
---|
| 917 | function [A,b] = generate_kernel_representation_data(N,N_unique,exponent_m2,exponent_p,p,options,p_base_parametric,ParametricIndicies,MonomIndicies,FirstRun) |
---|
| 918 | |
---|
| 919 | persistent saveData |
---|
| 920 | |
---|
| 921 | exponent_p_parametric = exponent_p(:,ParametricIndicies); |
---|
| 922 | exponent_p_monoms = exponent_p(:,MonomIndicies); |
---|
| 923 | pcoeffs = getbase(p); |
---|
| 924 | if any(exponent_p_monoms(1,:)) |
---|
| 925 | pcoeffs=pcoeffs(:,2:end); % No constant term in p |
---|
| 926 | end |
---|
| 927 | b = []; |
---|
| 928 | |
---|
| 929 | parametric = full((~isempty(ParametricIndicies) & any(any(exponent_p_parametric)))); |
---|
| 930 | |
---|
| 931 | % For problems with a lot of similar cones, this saves some time |
---|
| 932 | reuse = 0; |
---|
| 933 | if ~isempty(saveData) & isequal(saveData.N,N) & ~FirstRun |
---|
| 934 | n = saveData.n; |
---|
| 935 | ind = saveData.ind; |
---|
| 936 | if isequal(saveData.N_unique,N_unique) & isequal(saveData.exponent_m2,exponent_m2)% & isequal(saveData.epm,exponent_p_monoms) |
---|
| 937 | reuse = 1; |
---|
| 938 | end |
---|
| 939 | else |
---|
| 940 | % Congruence partition sizes |
---|
| 941 | for k = 1:size(N,1) |
---|
| 942 | n(k) = size(N{k},1); |
---|
| 943 | end |
---|
| 944 | % Save old SOS definition |
---|
| 945 | saveData.N = N; |
---|
| 946 | saveData.n = n; |
---|
| 947 | saveData.N_unique = N_unique; |
---|
| 948 | saveData.exponent_m2 = exponent_m2; |
---|
| 949 | saveData.N_unique = N_unique; |
---|
| 950 | end |
---|
| 951 | |
---|
| 952 | if reuse & options.sos.reuse |
---|
| 953 | % Get old stuff |
---|
| 954 | if size(exponent_m2{1},2)==2 % Stupid set(sos(parametric)) case |
---|
| 955 | ind = spalloc(1,1,0); |
---|
| 956 | ind(1)=1; |
---|
| 957 | allj = 1:size(exponent_p_monoms,1); |
---|
| 958 | used_in_p = ones(size(exponent_p_monoms,1),1); |
---|
| 959 | else |
---|
| 960 | allj = []; |
---|
| 961 | used_in_p = zeros(size(exponent_p_monoms,1),1); |
---|
| 962 | hash = randn(size(exponent_p_monoms,2),1); |
---|
| 963 | exponent_p_monoms_hash = exponent_p_monoms*hash; |
---|
| 964 | for i = 1:size(N_unique,1) |
---|
| 965 | monom = sparse(N_unique(i,3:end)); |
---|
| 966 | j = find(exponent_p_monoms_hash == (monom*hash)); |
---|
| 967 | |
---|
| 968 | if isempty(j) |
---|
| 969 | b = [b 0]; |
---|
| 970 | allj(end+1,1) = 0; |
---|
| 971 | else |
---|
| 972 | used_in_p(j) = 1; |
---|
| 973 | allj(end+1,1:length(j)) = j(:)'; |
---|
| 974 | end |
---|
| 975 | end |
---|
| 976 | ind = saveData.ind; |
---|
| 977 | end |
---|
| 978 | else |
---|
| 979 | allj = []; |
---|
| 980 | used_in_p = zeros(size(exponent_p_monoms,1),1); |
---|
| 981 | if size(exponent_m2{1},2)==2 % Stupid set(sos(parametric)) case |
---|
| 982 | ind = spalloc(1,1,0); |
---|
| 983 | ind(1)=1; |
---|
| 984 | allj = 1:size(exponent_p_monoms,1); |
---|
| 985 | used_in_p = ones(size(exponent_p_monoms,1),1); |
---|
| 986 | else |
---|
| 987 | % To speed up some searching, we random-hash data |
---|
| 988 | hash = randn(size(exponent_p_monoms,2),1); |
---|
| 989 | for k = 1:length(exponent_m2) |
---|
| 990 | if isempty(exponent_m2{k}) |
---|
| 991 | exp_hash{k}=[]; |
---|
| 992 | else |
---|
| 993 | exp_hash{k} = sparse((exponent_m2{k}(:,3:end)))*hash; % SPARSE NEEDED DUE TO STRANGE NUMERICS IN MATLAB ON 0s (the stuff will differ on last bit in hex format) |
---|
| 994 | end |
---|
| 995 | end |
---|
| 996 | |
---|
| 997 | p_hash = exponent_p_monoms*hash; |
---|
| 998 | ind = spalloc(size(N_unique,1),sum(n.^2),0); |
---|
| 999 | |
---|
| 1000 | for i = 1:size(N_unique,1) |
---|
| 1001 | monom = N_unique(i,3:end); |
---|
| 1002 | monom_hash = sparse(monom)*hash; |
---|
| 1003 | LHS = 0; |
---|
| 1004 | start = 0; |
---|
| 1005 | for k = 1:size(N,1) |
---|
| 1006 | j = find(exp_hash{k} == monom_hash); |
---|
| 1007 | if ~isempty(j) |
---|
| 1008 | pos=exponent_m2{k}(j,1:2); |
---|
| 1009 | nss = pos(:,1); |
---|
| 1010 | mss = pos(:,2); |
---|
| 1011 | indicies = nss+(mss-1)*n(k); |
---|
| 1012 | ind(i,indicies+start) = ind(i,indicies+start) + 1; |
---|
| 1013 | end |
---|
| 1014 | start = start + (n(k))^2; |
---|
| 1015 | % start = start + (matrixSOSsize*n(k))^2; |
---|
| 1016 | end |
---|
| 1017 | |
---|
| 1018 | j = find(p_hash == monom_hash); |
---|
| 1019 | |
---|
| 1020 | if isempty(j) |
---|
| 1021 | allj(end+1,1) = 0; |
---|
| 1022 | else |
---|
| 1023 | used_in_p(j) = 1; |
---|
| 1024 | allj(end+1,1:length(j)) = j(:)'; |
---|
| 1025 | end |
---|
| 1026 | end |
---|
| 1027 | end |
---|
| 1028 | end |
---|
| 1029 | saveData.ind = ind; |
---|
| 1030 | |
---|
| 1031 | % Some parametric terms in p(x,t) do not appear in v'Qv |
---|
| 1032 | % So these have to be added 0*Q = b |
---|
| 1033 | not_dealt_with = find(used_in_p==0); |
---|
| 1034 | while ~isempty(not_dealt_with) |
---|
| 1035 | j = findrows(exponent_p_monoms,exponent_p_monoms(not_dealt_with(1),:)); |
---|
| 1036 | allj(end+1,1:length(j)) = j(:)'; |
---|
| 1037 | used_in_p(j) = 1; |
---|
| 1038 | not_dealt_with = find(used_in_p==0); |
---|
| 1039 | ind(end+1,1)=0; |
---|
| 1040 | end |
---|
| 1041 | |
---|
| 1042 | matrixSOSsize = length(p); |
---|
| 1043 | if parametric |
---|
| 1044 | % Inconsistent behaviour in MATLAB |
---|
| 1045 | if size(allj,1)==1 |
---|
| 1046 | uu = [0;p_base_parametric]; |
---|
| 1047 | b = sum(uu(allj+1))'; |
---|
| 1048 | else |
---|
| 1049 | b = []; |
---|
| 1050 | for i = 1:matrixSOSsize |
---|
| 1051 | for j = i:matrixSOSsize |
---|
| 1052 | if i~=j |
---|
| 1053 | uu = [0;2*p_base_parametric(:,(i-1)*matrixSOSsize+j)]; |
---|
| 1054 | else |
---|
| 1055 | uu = [0;p_base_parametric(:,(i-1)*matrixSOSsize+j)]; |
---|
| 1056 | end |
---|
| 1057 | b = [b sum(uu(allj+1),2)']; |
---|
| 1058 | end |
---|
| 1059 | end |
---|
| 1060 | end |
---|
| 1061 | else |
---|
| 1062 | if matrixSOSsize == 1 |
---|
| 1063 | uu = [zeros(size(pcoeffs,1),1) pcoeffs]'; |
---|
| 1064 | b = sum(uu(allj+1,:),2)'; |
---|
| 1065 | else |
---|
| 1066 | b = []; |
---|
| 1067 | for i = 1:matrixSOSsize |
---|
| 1068 | for j = i:matrixSOSsize |
---|
| 1069 | if i~=j |
---|
| 1070 | uu = [0;2*pcoeffs((i-1)*matrixSOSsize+j,:)']; |
---|
| 1071 | else |
---|
| 1072 | uu = [0;pcoeffs((i-1)*matrixSOSsize+j,:)']; |
---|
| 1073 | end |
---|
| 1074 | b = [b;sum(uu(allj+1,:),2)']; |
---|
| 1075 | end |
---|
| 1076 | end |
---|
| 1077 | end |
---|
| 1078 | % uu = [0;pcoeffs(:)]; |
---|
| 1079 | % b = sum(uu(allj+1),2)'; |
---|
| 1080 | end |
---|
| 1081 | |
---|
| 1082 | b = b'; |
---|
| 1083 | dualbase = ind; |
---|
| 1084 | |
---|
| 1085 | j = 1; |
---|
| 1086 | A = cell(size(N,1),1); |
---|
| 1087 | for k = 1:size(N,1) |
---|
| 1088 | if matrixSOSsize==1 |
---|
| 1089 | A{k} = dualbase(:,j:j+n(k)^2-1); |
---|
| 1090 | else |
---|
| 1091 | % Quick fix for matrix SOS case, should be optimized |
---|
| 1092 | A{k} = inflate(dualbase(:,j:j+n(k)^2-1),matrixSOSsize,n(k)); |
---|
| 1093 | end |
---|
| 1094 | j = j + n(k)^2; |
---|
| 1095 | end |
---|
| 1096 | b = b(:); |
---|
| 1097 | |
---|
| 1098 | |
---|
| 1099 | |
---|
| 1100 | function newAi = inflate(Ai,matrixSOSsize,n); |
---|
| 1101 | % Quick fix for matrix SOS case, should be optimized |
---|
| 1102 | newAi = []; |
---|
| 1103 | for i = 1:matrixSOSsize |
---|
| 1104 | for r = i:matrixSOSsize |
---|
| 1105 | for m = 1:size(Ai,1) |
---|
| 1106 | ai = reshape(Ai(m,:),n,n); |
---|
| 1107 | V = zeros(matrixSOSsize,matrixSOSsize); |
---|
| 1108 | V(i,r)=1; |
---|
| 1109 | V(r,i)=1; |
---|
| 1110 | ai = kron(V,ai); |
---|
| 1111 | newAi = [newAi;ai(:)']; |
---|
| 1112 | end |
---|
| 1113 | end |
---|
| 1114 | end |
---|
| 1115 | |
---|
| 1116 | |
---|
| 1117 | function [Z,Q1,R] = sparsenull(A) |
---|
| 1118 | |
---|
| 1119 | [Q,R] = qr(A'); |
---|
| 1120 | n = max(find(sum(abs(R),2))); |
---|
| 1121 | Q1 = Q(:,1:n); |
---|
| 1122 | R = R(1:n,:); |
---|
| 1123 | Z = Q(:,n+1:end); % New basis |
---|
| 1124 | |
---|
| 1125 | |
---|
| 1126 | function [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,sparsityPattern); |
---|
| 1127 | |
---|
| 1128 | % To get the primal kernel representation, we simply use |
---|
| 1129 | % the built-in dualization module! |
---|
| 1130 | % First, write the problem in primal kernel format |
---|
| 1131 | traceobj = 0; |
---|
| 1132 | dotraceobj = options.sos.traceobj; |
---|
| 1133 | F = F_parametric; |
---|
| 1134 | for i = 1:length(Blockedb) |
---|
| 1135 | |
---|
| 1136 | |
---|
| 1137 | sizematrixSOS = sqrt(size(Blockedb{i},2)); |
---|
| 1138 | for k = 1:sizematrixSOS |
---|
| 1139 | for r = k:sizematrixSOS |
---|
| 1140 | res{(k-1)*sizematrixSOS+r} = 0; |
---|
| 1141 | end |
---|
| 1142 | end |
---|
| 1143 | |
---|
| 1144 | for j = 1:length(BlockedA{i}) |
---|
| 1145 | n = sqrt(size(BlockedA{i}{j},2)); |
---|
| 1146 | BlockedQ{i}{j} = sdpvar(n*sizematrixSOS,n*sizematrixSOS); |
---|
| 1147 | F = F + set(BlockedQ{i}{j}); |
---|
| 1148 | if sizematrixSOS>0 |
---|
| 1149 | % Matrix valued sum of sqaures |
---|
| 1150 | % Loop over all elements |
---|
| 1151 | starttop = 1; |
---|
| 1152 | for k = 1:sizematrixSOS |
---|
| 1153 | startleft = 1; |
---|
| 1154 | for r = 1:sizematrixSOS |
---|
| 1155 | if k<=r |
---|
| 1156 | Qkr = BlockedQ{i}{j}(starttop:starttop+n-1,startleft:startleft+n-1); |
---|
| 1157 | res{(k-1)*sizematrixSOS+r} = res{(k-1)*sizematrixSOS+r} + BlockedA{i}{j}*reshape(Qkr,n^2,1); |
---|
| 1158 | end |
---|
| 1159 | startleft = startleft + n; |
---|
| 1160 | end |
---|
| 1161 | starttop = starttop + n; |
---|
| 1162 | end |
---|
| 1163 | else |
---|
| 1164 | % Standard case |
---|
| 1165 | res{1} = res{1} + BlockedA{i}{j}*reshape(BlockedQ{i}{j},n^2,1); |
---|
| 1166 | end |
---|
| 1167 | if dotraceobj |
---|
| 1168 | traceobj = traceobj + trace(BlockedQ{i}{j}); |
---|
| 1169 | end |
---|
| 1170 | end |
---|
| 1171 | for k = 1:sizematrixSOS |
---|
| 1172 | for r = k:sizematrixSOS |
---|
| 1173 | F = F + set(res{(k-1)*sizematrixSOS+r} == Blockedb{i}(:,(k-1)*sizematrixSOS+r)); |
---|
| 1174 | end |
---|
| 1175 | end |
---|
| 1176 | end |
---|
| 1177 | |
---|
| 1178 | % % Constrain elements according to a desired sparsity |
---|
| 1179 | if ~isempty(sparsityPattern) |
---|
| 1180 | res = []; |
---|
| 1181 | for i = 1:length(BlockedQ) |
---|
| 1182 | for j = 1:length(BlockedQ{i}) |
---|
| 1183 | if ~isempty(sparsityPattern{i}{j}) |
---|
| 1184 | H = spalloc(length(BlockedQ{i}{j}),length(BlockedQ{i}{j}),length(sparsityPattern{i}{j})); |
---|
| 1185 | H(sparsityPattern{i}{j}) = 1; |
---|
| 1186 | k = find(triu(H)); |
---|
| 1187 | res = [res;BlockedQ{i}{j}(k)]; |
---|
| 1188 | end |
---|
| 1189 | end |
---|
| 1190 | end |
---|
| 1191 | F = F + set(0 == res); |
---|
| 1192 | end |
---|
| 1193 | |
---|
| 1194 | % And get the primal model of this |
---|
| 1195 | if isempty(parobj) |
---|
| 1196 | if options.sos.traceobj |
---|
| 1197 | [F,obj,Primal_matrices,Free_variables] = dualize(F,traceobj,1,options.sos.extlp); |
---|
| 1198 | else |
---|
| 1199 | [F,obj,Primal_matrices,Free_variables] = dualize(F,[],1,options.sos.extlp); |
---|
| 1200 | end |
---|
| 1201 | else |
---|
| 1202 | [F,obj,Primal_matrices,Free_variables] = dualize(F,parobj,1,options.sos.extlp); |
---|
| 1203 | end |
---|
| 1204 | % In dual mode, we maximize |
---|
| 1205 | obj = -obj; |
---|
| 1206 | |
---|
| 1207 | |
---|
| 1208 | function [F,obj,BlockedQ,sol] = create_imagemodel(BlockedA,Blockedb,F_parametric,parobj,options); |
---|
| 1209 | |
---|
| 1210 | |
---|
| 1211 | % ********************************* |
---|
| 1212 | % IMAGE REPRESENTATION |
---|
| 1213 | % Needed for nonlinearly parameterized problems |
---|
| 1214 | % More efficient in some cases |
---|
| 1215 | % ********************************* |
---|
| 1216 | g = []; |
---|
| 1217 | vecQ = []; |
---|
| 1218 | sol.problem = 0; |
---|
| 1219 | for i = 1:length(BlockedA) |
---|
| 1220 | q = []; |
---|
| 1221 | A = []; |
---|
| 1222 | for j = 1:length(BlockedA{i}) |
---|
| 1223 | n = sqrt(size(BlockedA{i}{j},2)); |
---|
| 1224 | BlockedQ{i}{j} = sdpvar(n,n); |
---|
| 1225 | q = [q;reshape(BlockedQ{i}{j},n^2,1)]; |
---|
| 1226 | A = [A BlockedA{i}{j}]; |
---|
| 1227 | end |
---|
| 1228 | vecQ = [vecQ;q]; |
---|
| 1229 | g = [g;Blockedb{i}-A*q]; |
---|
| 1230 | end |
---|
| 1231 | |
---|
| 1232 | g_vars = getvariables(g); |
---|
| 1233 | q_vars = getvariables(vecQ); |
---|
| 1234 | x_vars = setdiff(g_vars,q_vars); |
---|
| 1235 | |
---|
| 1236 | base = getbase(g); |
---|
| 1237 | if isempty(x_vars) |
---|
| 1238 | A = base(:,1);base = base(:,2:end); |
---|
| 1239 | B = (base(:,ismember(g_vars,q_vars))); |
---|
| 1240 | Bnull = sparsenull(B); |
---|
| 1241 | t = sdpvar(size(Bnull,2),1); |
---|
| 1242 | imQ = -B\A+Bnull*t; |
---|
| 1243 | else |
---|
| 1244 | A = base(:,1);base = base(:,2:end); |
---|
| 1245 | C = base(:,ismember(g_vars,x_vars)); |
---|
| 1246 | B = (base(:,ismember(g_vars,q_vars))); |
---|
| 1247 | [Bnull,Q1,R1] = sparsenull(B);Bnull(abs(Bnull) < 1e-12) = 0; |
---|
| 1248 | t = sdpvar(size(Bnull,2),1); |
---|
| 1249 | imQ = -Q1*(R1'\(A+C*recover(x_vars)))+Bnull*t; |
---|
| 1250 | end |
---|
| 1251 | notUsed = find(sum(abs(B),2)==0); |
---|
| 1252 | if ~isempty(notUsed) |
---|
| 1253 | ff=g(notUsed); |
---|
| 1254 | if isa(ff,'double') |
---|
| 1255 | if nnz(ff)>0 |
---|
| 1256 | sol.yalmiptime = 0; % FIX |
---|
| 1257 | sol.solvertime = 0; |
---|
| 1258 | sol.problem = 2; |
---|
| 1259 | sol.info = yalmiperror(1,'YALMIP'); |
---|
| 1260 | F = []; |
---|
| 1261 | obj = []; |
---|
| 1262 | if options.verbose > 0 |
---|
| 1263 | disp(' '); |
---|
| 1264 | disp('-> During construction of data, I encountered a situation') |
---|
| 1265 | disp('-> situation that tells me that the problem is trivially infeasible!') |
---|
| 1266 | disp('-> Have you forgotten to define some parametric variables?,') |
---|
| 1267 | disp('-> or perhaps you have a parametric problem where the highest') |
---|
| 1268 | disp('-> power in some of the independent variables is odd for any choice') |
---|
| 1269 | disp('-> of parametric variables, such as x^8+x^7+t*y^3') |
---|
| 1270 | disp('-> Anyway, take a good look at your model again...'); |
---|
| 1271 | end |
---|
| 1272 | return |
---|
| 1273 | % error('You seem to have a strange model. Have you forgotten to define some parametric variable?'); |
---|
| 1274 | end |
---|
| 1275 | else |
---|
| 1276 | F_parametric = F_parametric + set(g(notUsed)==0); |
---|
| 1277 | end |
---|
| 1278 | end |
---|
| 1279 | F_sos = set([]); |
---|
| 1280 | obj = 0; |
---|
| 1281 | for i = 1:length(BlockedQ) |
---|
| 1282 | for j = 1:size(BlockedQ{i},2) |
---|
| 1283 | Q_old = BlockedQ{i}{j}; |
---|
| 1284 | Q_old_vars = getvariables(Q_old); |
---|
| 1285 | Q_old_base = getbase(Q_old); |
---|
| 1286 | in_this = []; |
---|
| 1287 | for k = 1:length(Q_old_vars) |
---|
| 1288 | in_this = [in_this;find(Q_old_vars(k)==q_vars)]; |
---|
| 1289 | end |
---|
| 1290 | Q_new = Q_old_base(:,1) + Q_old_base(:,2:end)*imQ(in_this); |
---|
| 1291 | Q_new = reshape(Q_new,length(BlockedQ{i}{j}),length(BlockedQ{i}{j})); |
---|
| 1292 | obj = obj+trace(Q_new); |
---|
| 1293 | if ~isa(Q_new,'double') |
---|
| 1294 | F_sos = F_sos + set(Q_new); |
---|
| 1295 | elseif min(eig(Q_new))<-1e-8 |
---|
| 1296 | sol.yalmiptime = 0; % FIX |
---|
| 1297 | sol.solvertime = 0; |
---|
| 1298 | sol.problem = 2; |
---|
| 1299 | sol.info = yalmiperror(1,'YALMIP'); |
---|
| 1300 | F = []; |
---|
| 1301 | obj = []; |
---|
| 1302 | error('Problem is trivially infeasible. After block-diagonalizing, I found constant negative definite blocks!'); |
---|
| 1303 | return |
---|
| 1304 | end |
---|
| 1305 | BlockedQ{i}{j} = Q_new; |
---|
| 1306 | end |
---|
| 1307 | end |
---|
| 1308 | |
---|
| 1309 | F = F_parametric + F_sos; |
---|
| 1310 | |
---|
| 1311 | if isa(obj,'double') | (options.sos.traceobj == 0) |
---|
| 1312 | obj = []; |
---|
| 1313 | end |
---|
| 1314 | |
---|
| 1315 | if ~isempty(parobj) |
---|
| 1316 | obj = parobj; |
---|
| 1317 | end |
---|
| 1318 | |
---|
| 1319 | |
---|
| 1320 | function [F,obj,options] = create_lrmodel(BlockedA,Blockedb,F_parametric,parobj,options,ParametricVariables) |
---|
| 1321 | % Some special code for ther low-rank model in SDPLR |
---|
| 1322 | % Experimental code, not official yet |
---|
| 1323 | allb = []; |
---|
| 1324 | allA = []; |
---|
| 1325 | K.s = []; |
---|
| 1326 | for i = 1:length(Blockedb) |
---|
| 1327 | allb = [allb;Blockedb{i}]; |
---|
| 1328 | Ai = []; |
---|
| 1329 | for j = 1:size(BlockedA{i},2) |
---|
| 1330 | Ai = [Ai BlockedA{i}{j}]; |
---|
| 1331 | K.s = [K.s sqrt(size(BlockedA{i}{j},2))]; |
---|
| 1332 | end |
---|
| 1333 | %blkdiag bug in 7.0... |
---|
| 1334 | [n1,m1] = size(allA); |
---|
| 1335 | [n2,m2] = size(Ai); |
---|
| 1336 | allA = [allA spalloc(n1,m2,0);spalloc(n2,m1,0) Ai]; |
---|
| 1337 | end |
---|
| 1338 | options.solver = 'sdplr'; |
---|
| 1339 | z = recover(ParametricVariables) |
---|
| 1340 | start = size(BlockedA,2)+1; |
---|
| 1341 | Mi = []; |
---|
| 1342 | for i = 1:length(allb) |
---|
| 1343 | if isa(allb(i),'sdpvar') |
---|
| 1344 | [Qi,ci,fi,xi,infoi] = quaddecomp(allb(i),z); |
---|
| 1345 | else |
---|
| 1346 | Qi = zeros(length(z)); |
---|
| 1347 | ci = zeros(length(z),1); |
---|
| 1348 | fi = allb(i); |
---|
| 1349 | end |
---|
| 1350 | Z = -[fi ci'/2;ci/2 Qi]; |
---|
| 1351 | Mi = [Mi;Z(:)']; |
---|
| 1352 | end |
---|
| 1353 | K.s = [K.s length(z)+1]; |
---|
| 1354 | zeroRow = zeros(1,size(allA,2)); |
---|
| 1355 | allA = [allA Mi;zeroRow 1 zeros(1,K.s(end)^2-1)]; |
---|
| 1356 | b = zeros(size(allA,1),1);b(end) = 1; |
---|
| 1357 | y = sdpvar(length(b),1); |
---|
| 1358 | CminusAy = -allA'*y; |
---|
| 1359 | start = 1; |
---|
| 1360 | |
---|
| 1361 | % Get the cost, expressed in Z |
---|
| 1362 | [Qi,ci,fi,xi,infoi] = quaddecomp(parobj,z); |
---|
| 1363 | C = [fi ci'/2;ci/2 Qi]; |
---|
| 1364 | F = set([]); |
---|
| 1365 | for i = 1:length(K.s) |
---|
| 1366 | if i<length(K.s) |
---|
| 1367 | F = F + set(reshape(CminusAy(start:start+K.s(i)^2-1),K.s(i),K.s(i))); |
---|
| 1368 | else |
---|
| 1369 | F = F + set(reshape(C(:) + CminusAy(start:start+K.s(i)^2-1),K.s(i),K.s(i))); |
---|
| 1370 | end |
---|
| 1371 | start = start + K.s(i)^2; |
---|
| 1372 | end |
---|
| 1373 | obj = -b'*y; |
---|
| 1374 | |
---|
| 1375 | options.sdplr.forcerank = ceil(K.s/2); |
---|
| 1376 | options.sdplr.forcerank(end) = 1; |
---|
| 1377 | options.sdplr.feastol = 1e-7; |
---|
| 1378 | |
---|
| 1379 | |
---|
| 1380 | |
---|
| 1381 | function [exponent_m2,N_unique] = expandmatrix(exponent_m2,N_unique,n); |
---|
| 1382 | |
---|
| 1383 | |
---|
| 1384 | function [sol,m,Q,residuals,everything] = solvesos_find_blocks(F,obj,options,params,candidateMonomials) |
---|
| 1385 | |
---|
| 1386 | tol = options.sos.numblkdg; |
---|
| 1387 | if tol > 1e-2 |
---|
| 1388 | disp(' '); |
---|
| 1389 | disp('-> Are you sure you meant to have a tolerance in numblk that big!') |
---|
| 1390 | disp('-> The options numblkdiag controls the tolerance, it is not a 0/1 switch.') |
---|
| 1391 | disp(' '); |
---|
| 1392 | end |
---|
| 1393 | options.sos.numblkdg = 0; |
---|
| 1394 | [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,candidateMonomials); |
---|
| 1395 | |
---|
| 1396 | % Save old structure to find out when we have stalled |
---|
| 1397 | for i = 1:length(Q) |
---|
| 1398 | oldlengths{i} = length(Q{i}); |
---|
| 1399 | end |
---|
| 1400 | |
---|
| 1401 | go_on = (sol.problem == 0 | sol.problem == 4); |
---|
| 1402 | while go_on |
---|
| 1403 | |
---|
| 1404 | for sosfun = 1:length(Q) |
---|
| 1405 | Qtemp = Q{sosfun}; |
---|
| 1406 | keep = diag(Qtemp)>tol; |
---|
| 1407 | Qtemp(:,find(~keep)) = []; |
---|
| 1408 | Qtemp(find(~keep),:) = []; |
---|
| 1409 | |
---|
| 1410 | m{sosfun} = m{sosfun}(find(keep)); |
---|
| 1411 | |
---|
| 1412 | Qtemp(abs(Qtemp) < tol) = 0; |
---|
| 1413 | [v1,dummy1,r1,dummy3]=dmperm(Qtemp+eye(length(Qtemp))); |
---|
| 1414 | lengths{sosfun} = []; |
---|
| 1415 | n{sosfun} = {}; |
---|
| 1416 | for blocks = 1:length(r1)-1 |
---|
| 1417 | i1 = r1(blocks); |
---|
| 1418 | i2 = r1(blocks+1)-1; |
---|
| 1419 | if i2>i1 |
---|
| 1420 | n{sosfun}{blocks} = m{sosfun}(v1(i1:i2)); |
---|
| 1421 | else |
---|
| 1422 | n{sosfun}{blocks} = m{sosfun}(v1(i1)); |
---|
| 1423 | end |
---|
| 1424 | lengths{sosfun} = [lengths{sosfun} length(n{sosfun}{blocks})]; |
---|
| 1425 | end |
---|
| 1426 | lengths{sosfun} = sort(lengths{sosfun}); |
---|
| 1427 | end |
---|
| 1428 | go_on = ~isequal(lengths,oldlengths); |
---|
| 1429 | oldlengths = lengths; |
---|
| 1430 | if go_on |
---|
| 1431 | [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,n); |
---|
| 1432 | go_on = go_on & (sol.problem == 0 | sol.problem == 4); |
---|
| 1433 | if sol.problem == 1 |
---|
| 1434 | disp('-> Feasibility was lost during the numerical block-diagonalization.') |
---|
| 1435 | disp('-> The setting sos.numblkdiag is probably too big') |
---|
| 1436 | end |
---|
| 1437 | end |
---|
| 1438 | end |
---|
| 1439 | |
---|
| 1440 | |
---|
| 1441 | function sol = solveranksos(F,obj,options,ranks,BlockedQ) |
---|
| 1442 | |
---|
| 1443 | Frank = set([]); |
---|
| 1444 | for i = 1:length(ranks) |
---|
| 1445 | if ~isinf(ranks(i)) |
---|
| 1446 | Frank = Frank + set(rank(BlockedQ{i}{1}) <= ranks(i)); |
---|
| 1447 | end |
---|
| 1448 | end |
---|
| 1449 | % rank adds the pos.def constraints again!!, so we remove them |
---|
| 1450 | check = ones(length(F),1); |
---|
| 1451 | keep = ones(length(F),1); |
---|
| 1452 | for i = 1:length(BlockedQ) |
---|
| 1453 | for j = 1:length(BlockedQ{i}) |
---|
| 1454 | Qij = BlockedQ{i}{j}; |
---|
| 1455 | for k = find(check)' |
---|
| 1456 | if isequal(Qij,sdpvar(F(k))) |
---|
| 1457 | keep(k) = 0; |
---|
| 1458 | check(k) = 0; |
---|
| 1459 | end |
---|
| 1460 | end |
---|
| 1461 | end |
---|
| 1462 | end |
---|
| 1463 | % Let's hope LMIRANK is there |
---|
| 1464 | sol = solvesdp(F(find(keep)) + Frank,[],sdpsettings(options,'solver','')); |
---|
| 1465 | |
---|
| 1466 | |
---|