diff --git a/benchmark.m b/benchmark.m new file mode 100644 index 0000000..ba6d484 --- /dev/null +++ b/benchmark.m @@ -0,0 +1,235 @@ +function [s] = benchmark(path, filename, file, tfile, solverlist) +%% BECNHMARKING - a routine to perform a benchmark +% s = benchmarking(path, filename, file, tfile) +% +%% >> Description: +% A benchmark test will be performed on problems found on path 'path'. +% The results (optimal values ,rigorous lower and upper bounds and times) +% are saved in the textfiles 'filename' reps. 'filename'_timings or +% in the files given by handels 'file' resp. 'tfile'. Some Tex sequences +% for tables in Latex will be created. +% +%% >> Input: +% path: path of test problems +% filename: name of the file for the results +% file, tfile: file handels for results if not given by 'filename' +% +%% >> Output: +% s: 0 if benchmark test completed, otherwise an error will be thrown +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 12/08/12 M. Lange, small adaptions for overworked VSDP +% + +%% preparation +sysstr = computer; +if isempty(strfind(sysstr,'WIN')) + nlstr= '\n'; +else + nlstr = '\r\n'; +end + +% list of solvers that shall be benchmarked +if nargin<5 || ~iscell(solverlist) + solverlist = {'sdpt3','sedumi','sdpa','csdp'}; +end +ml = length(solverlist); + +if nargin<3 || isempty(file) || file<=0 + file = fopen(filename,'w'); + if file==0 + error('VSDP:BENCHMARK','Could not open the file %s',filename); + end +end + +if nargin<4 || isempty(tfile) || tfile<=0 + tfilename = [filename(1:end-4),'_timings.txt']; + tfile = fopen(tfilename,'w'); + if tfile==0 + error('VSDP:BENCHMARK','Could not open the file %s',tfilename); + end +end + +% root contains files and sub-folders in path folder +root = dir(path); + +if length(root)==2 % empty folder + warning('VSDP:BENCHMARK',['"',path,'" is empty']); + return; +elseif (length(root)>2) % if folder: exclude './' and '../' in root + root([1 2]) = []; +end +% folders first than files +[tmp idx] = sort([root(:).isdir],'descend'); +root = root(idx); + +% write table head if test subdirectory reached +if ~root(1).isdir + % table for results + fprintf(file,['\\begin{center} ',nlstr]); + fprintf(file,'\\begin{table} {Test Problems from %s} \\\\[1ex]',path); + fprintf(file,[nlstr,'\\begin{tiny}',nlstr]); + fprintf(file,['\\begin{tabular}[b]{|l|c|c|c|c|c|c|l|} \\hline',... + '& & & & & & & \\\\[-2ex]',nlstr]); + fprintf(file,['\\textbf{Problem} & $p^*$ & $d^*$ & $\\bar{d}$ & ',... + '$\\underline{p}$ & $\\mu(p^*,d^*)$ & $\\mu(',... + '\\underline{p},\\bar{d})$ & Solver \\\\ \\hline',nlstr]); + % table for timings + fprintf(tfile,['\\begin{center} ',nlstr]); + fprintf(tfile,'\\begin{table} {Timings for %s} \\\\[1ex]',path); + fprintf(tfile,[nlstr,'\\begin{tiny}',nlstr]); + fprintf(tfile,['\\begin{tabular}[b]{|l|c|c|c|l|} \\hline',... + 'Problem & $t_s$ & $t_u$ & $t_l$ & ',... + 'Solver \\\\ \\hline',nlstr]); +end + +for j = 1:length(root) % start from third entry to exclude './' and '../' in mr + if root(j).isdir + newpath = fullfile(path,root(j).name); + benchmark(newpath,'',file,tfile,solverlist); + else + % cleaning/preparing memory + clear A b c K x0 y0 z0 objt info; + pack; + % go the next problem + probname = root(j).name; + fullname = fullfile(path,probname); + try + % load matlab data + if length(probname)>=5 && strcmpi(fullname(end-3:end),'.mat') + load(fullname); + ww = whos; + for jw=1:length(ww) + if strcmp(ww(jw).name,'At') + A = At; clear At; + end + if strcmp(ww(jw).name,'Amatrix') + A = Amatrix; clear Amatrix; + end + if strcmp(ww(jw).name,'C') + c = C; clear C ww; + break; + end + end + % [A,b,c] = sedumi2vsdp(A,b,c,K); + % load sparse sdpa data + elseif strcmpi(probname(end-min(4,length(probname)-1):end),'dat-s') + [A,b,c,K] = sdpa2vsdp(fullname); + % load nothing, continue to next file + elseif strcmpi(probname(end-min(2,length(probname)-1):end),'SIF') + [A,b,c,K] = mps2vsdp(fullname); + else + warning('VSDP:BENCHMARK','not supported file found'); + continue; + end + catch loaderror + msgstr = loaderror.message; + msglen = min(length(msgstr),20); + msgstr = msgstr(1:msglen); + fprintf(file,['%s & \\multicolumn{6}{|l|}{',... + msgstr,'} & \\\\ \\hline',nlstr],... + probname); + fprintf(tfile,['%s & \\multicolumn{3}{|l|}{ & ',... + msgstr,'} & \\\\ \\hline',nlstr],... + probname); + continue; + end + % some preprocessing of probleme names + sidx = (probname=='_'); + probname(sidx) = ''; + sidx = find(probname=='.') - 1; + if isempty(sidx) + sidx = length(probname); + end + probname = probname(1:min(sidx,12)); + %try to solve the problem + for i = 1:ml + opts.SOLVER = solverlist{i}; + % try to solve with an approximate solver + try + tic; [objt,x0,y0,z0] = mysdps(A,b,c,K,[],[],[],opts); ts = toc; + ps = objt(1); ds = objt(2); + fprintf(file,'%s & %1.8e & %1.8e & ',probname,ps,ds); + catch solvererror + msgstr = solvererror.message; + msglen = min(length(msgstr),20); + msgstr = msgstr(1:msglen); + fprintf(file,['%s & \\multicolumn{6}{|l|}{',... + msgstr,'} & %s\\\\ \\hline',nlstr],... + probname,solverlist{i}); + fprintf(tfile,['%s & \\multicolumn{3}{|l|}{',... + msgstr,'} & %s\\\\ \\hline',nlstr],... + probname,solverlist{i}); + continue; + end + % try to find upper bound + try + tic; fU = vsdpup(A,b,c,K,x0,y0,z0,[],opts); tu = toc; + fprintf(file,'%1.8e & ',fU); + catch vsdpuperror + msgstr = vsdpuperror.message; + msglen = min(length(msgstr),20); + msgstr = msgstr(1:msglen); + fprintf(file,[' & \\multicolumn{4}{|l|}{',... + msgstr,'} & %s\\\\ \\hline',nlstr],... + 'VSDP:VSDPUP'); + fprintf(tfile,[' & \\multicolumn{2}{|l|}{',... + msgstr,'} & %s\\\\ \\hline',nlstr],... + 'VSDP:VSDPUP'); + continue; + end + % try to find lower bound + try + tic; fL = vsdplow(A,b,c,K,x0,y0,z0,[],opts); tl = toc; + fprintf(file,'%1.8e & ',fL); + catch dummy + msgstr = dummy.message; + msglen = min(length(msgstr),20); + msgstr = msgstr(1:msglen); + fprintf(file,[' & \\multicolumn{1}{|l|}{',... + msgstr,'} & %s\\\\ \\hline',nlstr],... + 'VSDP:VSDPLOW'); + fprintf(tfile,[' & \\multicolumn{1}{|l|}{',... + msgstr,'} & %s\\\\ \\hline',nlstr],... + 'VSDP:VSDPLOW'); + continue; + end + mup = (ps - ds)/max(1,(abs(ps) + abs(ds))/2); + muv = (fU - fL)/max(1,(abs(fL) + abs(fU))/2); + % write rest results + fprintf(file,['%1.8e & %1.8e & %s \\\\ \\hline',nlstr],... + mup, muv, solverlist{i}); + fprintf(tfile,['%s & %6.3d & %6.3d & %6.3d & %s \\\\ ',... + '\\hline',nlstr],probname,ts,tu,tl,solverlist{i}); + clear objt x0 y0 z0; pack; + end + end +end + +% write table closing sequences +if ~root(1).isdir + % end of results table + fprintf(file,['\\end{tabular}',nlstr]); + fprintf(file,['\\end{tiny}',nlstr]); + fprintf(file,['\\end{table}',nlstr]); + fprintf(file,['\\end{center}',nlstr,nlstr]); + % end of timings table + fprintf(tfile,['\\end{tabular}',nlstr]); + fprintf(tfile,['\\end{tiny}',nlstr]); + fprintf(tfile,['\\end{table}',nlstr]); + fprintf(tfile,['\\end{center}',nlstr]); +end + +s = 0; diff --git a/bnd4sd.m b/bnd4sd.m new file mode 100644 index 0000000..26357fd --- /dev/null +++ b/bnd4sd.m @@ -0,0 +1,291 @@ +function [lambdaMin,trN,dshift] = bnd4sd(A,sym,mode) +% NBND4SDP calculates a lower bound for the smallest eigenvalue and the +% sum over all negative eigenvalues of A. +% The matrix A is assumed to be nearly positive semidefinite. +% BND4SD can be used to find rigorous bounds for non-positive +% definiteness due to rounding errors. The function is designed +% for the use in VSDP. +% +% A may also be an interval matrix, full or sparse. This interval +% matrix can be defined as: +% 1. a structure with the fields 'mid' and 'rad', +% 2. a cell array, where the first entry will be interpreted +% as 'mid' and the second as 'rad', +% 3. a class of type 'intval' [see INTLAB toolbox]. +% +% @params +% A: Real symmetric input matrix. +% +% sym: The function requires a symmetric input matrix. +% In case that A is not symmetric (for memory resons etc.) +% set 'sym=false'. Then S = 0.5*(A+A') will be used instead. +% -> default: true +% +% mode: inclusion mode: +% 0 - fast verification by applying cholesky +% factorization, if cholesky fails, lower bounds is +% -Inf +% 1 - if verification by cholesky decomposition failed +% function will call 'vsdpneig' for a stronger enclosure +% of all eigenvalues +% default: 1 +% +% @output +% lambdaN: A lower bound for the minimum eigenvalue of A. +% +% trN: A lower bound for the sum of all negative eigenvalues of A. +% +% dshift: Shift of diagonal elements. Possible permutation of the +% matrix: If positive semidefiniteness could not be +% verified, it may be possible for the slightly +% perturbated matrix A+diag(dshift). +% +% +% @dependencies: 'setround'-function [see Intlab] +% + +%% +% written 05/04/12 Marko Lange +% modified 05/07/12 -> additional parameter MINEIG +% modified 12/06/12 -> using eigs instead of eig +% modified 17/06/12 -> completely rewritten for adaptive shifting +% modified 21/08/12 -> added 'mode' parameter to enable call of vsdpneig +% + +%% +% ToDo +% - support for complex Hermitian matrices +% + + +%% check input data %% + +% input matrix +if nargin<1 || isempty(A) + error('BND4SD: No input matrix set.'); +elseif ~(isnumeric(A) || isa(A,'intval') || all(isfield(A,{'mid','rad'}))) + error('BND4SD: Datatype of input matrix is not accepted.'); +end % size & dimension is not checked. Error will occur in the code. + + +%% prepare input data %% + +% get rounding mode and set round to nearest +rnd = getround(); +setround(0); % for preparation only + +% prepare interval input +if ~isnumeric(A) % A is an interval structure or of type 'intval' + Arad = A.rad; + A = A.mid; +else + Arad = sparse(size(A,1),size(A,2)); +end + + +%% prepare A +if ~isreal(A) + error('BND4SD: Complex input not yet supported') +elseif size(A,1)~=size(A,2) && size(A,2)~=1 + A = A(:); + Arad = Arad(:); +end +n = size(A,1); + +if size(A,2)==1 % compact vectorized input + n = sqrt(2*n+0.25) - 0.5; + if n~=round(n) + error('BND4SD: Wrong dimension of vectorized input'); + end + E(triu(true(n))) = full(A); % E is place-holder, will be overwritten later + A = reshape(E,n,n); + A = A + triu(A,1)'; + if nnz(A)=2 && ~isempty(sym) && ~sym % non-symmetric input + setround(1); + E = 0.5 * A; % E and L are used as place-holder + L = E'; + A = E + L; % sup(0.5*(A+A')) + Arad = 0.5 * (Arad+Arad') + (A + ((-E) - L)); % Arad + (sup(A)-inf(A)) + setround(0); +end + + +%% calculate expected error for precondition shift +E = abs(A); +d = sqrt(diag(E)); +setround(-1); +L = E~=d*d'; +setround(1); +L = L | E~=d*d'; +E = E.*L; +dshift = 8*eps*sqrt(sum(max(E))*sum(E)) + sum(Arad) + realmin; +setround(0); + + +%% remove zero rows/columns +index = any(A-diag(sparse(diag(A)))) | any(Arad - diag(sparse(diag(Arad)))); +lambdaMin = inf; +trN = 0; + +if ~all(index) + setround(-1); + lambdaMin = min(diag(A)-diag(Arad)); + trN = sum(min(diag(A)-diag(Arad),0)); + setround(0); + A = A(index,index); + Arad = Arad(index,index); + n = length(A); +end + +if isempty(A) % A was pure diagonal matrix + dshift = min(lambdaMin/mean(dshift),0) * dshift; + setround(rnd); + return; +end + + +%% approximative cholesky factorization with adaptive error shift %% + +% approximate cholesky +[L,p] = chol(A-diag(sparse(dshift(index))),'lower'); +if p~=0 % not semidefinite + [L,p] = chol(A-diag(sparse(0.4*dshift(index))),'lower'); +end +if p~=0 + [L,p] = chol(A,'lower'); +end +if p==0 && ~isempty(find(isnan(L),1)) + p = n; +end + + +%% cholesky decomposition did not work +if nargin==2 && mode==0 && p~=0 + lambdaMin = -Inf; + trN = -Inf; + return; +elseif p~=0 % call vsdpneig + A = struct('mid',A,'rad',Arad); + clear Arad D E L; % not needed anymore + lambda = vsdpneig(A,1); % full enclosure + setround(-1); + lambda = lambda.mid - lambda.rad; % inf(lambda) + lambdaMin = min(min(lambda),lambdaMin); + trN = sum(min(lambda,0)) + trN; + dshift = min((-lambdaMin)/mean(-dshift),0) * dshift; + setround(rnd); + return; +end + + +%% calculate error matrix %% +setround(-1); +E = -(L*L'-A) - diag(sparse(inf(n,1))); % -inf(L*L'-A) = sup(A-L*L'), without diagonal +setround(1); +E = max(E,L*L'-A)+Arad; % max(sup(A-L*L')-inf*I,sup(L*L'-A)) + + +if find(isnan(E) | tril(E,-1)<0,1) + disp('break'); +end + +% clear variables to save memory +clear L D; + +% verify that E is symmetric +E = min(E,E'); + +% keep setround(1) as default rounding for verified norm bounds + + +%% find lower bound for minimum eigenvalue %% + +% apply adaptive shift of diagonal elements +d = -E(1:n+1:end); % diagonal of E had opposite sign +maxd = max(d); +E(1:n+1:end) = maxd - d; + +% we have Enew = abs(Eold - D) with D = maxd*I +% hence lambda_i(Eold)>= lambda_i(D)-norm(Enew,2) + +% using collatz theorem for upper bound of 2-norm +y = sum(E,2); +norm2bnd = max(y); + +while norm2bnd + oldnorm = norm2bnd; + x = y / norm2bnd; + y = E*x; + norm2bnd = max(y./x); + if norm2bnd*1.005>=oldnorm + break; + end +end + +% lambda_min(D+E) >= lambda_min(D) - norm(E,2) +% => lambdaMin >= inf(maxd-norm(E,2)) = -sup(norm(E,2)-maxd) +lambdaMinO = lambdaMin; % old upper bound for the eigenvalues +lambdaMin = min(-(norm2bnd-maxd), lambdaMinO); + +if lambdaMin>=0 + trN = 0; + return; +elseif nargin>2 && mode>0 % do full eigenvalue enclosure + A = struct('mid',A,'rad',Arad); + clear Arad E; + lambda = vsdpneig(A,1); + lambda = lambda.rad - lambda.mid; % -inf(lambda) + lambdaMin = max(min(-max(lambda),lambdaMinO),lambdaMin); + trN = - (sum(max(lambda,0)) - trN); + dshift = min(lambdaMin,0) * (dshift/mean(dshift)); + setround(rnd); + return; +end + + +%% calculate trace bound %% + +% trN := trace of negative semidefinite part of A +% -> sum(|lambda(E)|) <= sqrt(n) * norm(E,'fro') +% -> sum(lambda_pos(E))+sum(lamda_neg(E)) = trace(E); +% => 2*sum(|lamda_neg(E)|) <= sqrt(n) * norm(E,'fro') - trace(E) +lambda_sn = 0.5 * ( sqrtsup(n*(E(:)'*E(:))) + sum(-diag(E)) ); + +if maxd<0 % if maxd<0: trN = trace(D) - sum(|lamda_neg(E)|) + trN = -((-n)*maxd + lambda_sn); +else + % for a given lambda_sn = sum(|lamda_neg(E)|), with maxd>0 and + % nneg := number of negative eigenvalues it holds: + % trN = - lambda_sn - nneg*maxd + % => trN is maximized when nneg is minimized + nneg = floor(lambda_sn/norm2bnd); % number of negative eigenvalues + trN = -max((-nneg)*lambdaMin,lambda_sn+(-nneg-1)*maxd); +end + + +% reset rounding mode +setround(rnd); + +%_____________________________END OF NBND4SDP____________________________ + + + +% the following is a small function to calculate a verified upper bound +% for the square root -> setround(1) is assumed +% only used in case that the sqrt function does not regard the rounding mode +function res = sqrtsup(a) +res = sqrt(a); +res = res + realmin * (res.*(-res) > -a); \ No newline at end of file diff --git a/conversion/lp2vsdp.m b/conversion/lp2vsdp.m new file mode 100644 index 0000000..c6fa854 --- /dev/null +++ b/conversion/lp2vsdp.m @@ -0,0 +1,117 @@ +function [A,b,c,K,pd] = lp2vsdp(A,b,c,e,lb,ub) +%% LP2VSDP: transforms problem data from LP_SOLVE format to VSDP3 format +% [At,b,c,K] = lp2vsdp(A,b,c,e,lb,ub) +% +%% >> Description: +% The LP_SOLVE solver solves MILP problems of the form: +% max c'*x +% s.t. A*x <> b +% lb <= x <= ub +% Function converts this form into the standard primal/dual form of conic +% problems treated by VSDP3. Note, that if the problem gets dualized the +% optimal value must be negated to get the right sign. +% +%% >> Input: +% A: m x n matrix representing linear constraints +% b: m x 1 vector of right sides for the inequality constraints +% c: n x 1 vector of coefficients for primal objective function +% e: m x 1 flag vector that determines the sense of the inequalities +% -1 - less than +% 0 - equality +% 1 - greater than +% +%% >> Output: +% A: Matrix of linear equations +% b, c: - coefficients of primal/dual objective functions +% K: K.f, K.l, K.q, K.s converted to column vectors +% pd: a scalar with one of 'p' or 'd'. If 'p' problem was not dualized. +% If 'd' problem was dualized and optimal value must be negated. +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/08/10 V. Haerter, comments added +% 09/08/12 M. Lange, speed improvement and code reduction +% 10/08/12 M. Lange, extended primal form conversion +% +%% +% TODO: improvements: - computations in primal form with interval +% inclusion +% - mirroring of variables [-inf 0] to R+ +% + +%% check input +b = b(:); +c = c(:); +[m n] = size(A); +if m~=length(b) || n~=length(c) + error('VSDP:LP2VSDP','c or b not compatible with A'); +elseif m~=length(e) + error('VSDP:LP2VSDP','dimension of flag vector "e" does not match'); +end +if nargin<5 + lb = -inf(n,1); + ub = inf(n,1); +elseif nargin<6 + ub = inf(n,1); +end + +% resolve fixed bounds +ind = find(lb==ub); % index for fixed bounds +nfx = length(ind); % number of fixed bounds +if nfx>0 + A = [A; sparse(1:nfx,ind,1,nfx,n)]; + e = [e; zeros(nfx,1)]; + b = [b; lb(ind)]; + lb(ind) = -inf; + ub(ind) = inf; + m = m + nfx; +end + +% indices +indl = find(e<0); % inequalities less than (lower) +indu = find(e>0); % inequalities greater than (upper) + +Ai = [A(indl,:); -A(indu,:)]; % part of coefficient matrix for inequalities +A([indl; indu],:) = []; % part of coefficient matrix for equality constraints + +bi = [b(indl); -b(indu)]; +b([indl; indu]) = []; +mi = length(bi); + +indl = find(lb>-inf); % index for (applied) lower bounds +nl = length(indl); % number of lower bounds + +if all(lb(indl)==0) && min(ub)==inf && n+mi> Input: +% problem: can be the filename of a text file in MPS format or a problem +% structure as created by the read_mps function [see VSDP/read_mps] +% +%% >> Output: +% A: matrix of linear equations +% b, c: - coefficients of primal/dual objective functions +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% pd: a scalar with one of 'p' or 'd'. If 'p' problem was not dualized. +% If 'd' problem was dualized and optimal value must be negated. +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 29/08/12 M. Lange, written +% +%% ToDo: +% - support ranges +% - check of unsupported flags etc +% + +%% check input +if ischar(problem) + problem = read_mps(problem); +elseif isempty(problem) || ~all(isfield(problem,{'A','rowtypes'})) + error('VSDP:MPS2VSDP','insufficient input data'); +end + + +%% read A and c +idx = problem.rowtypes~='N'; +c = -sum(problem.A(~idx,:),1)'; % (-) because of maximization +A = problem.A(idx,:); + + +%% read e +e = problem.rowtypes(idx); +e = (e=='G') - (e=='L'); + + +%% read b +if isfield(problem,'rhs') && ~isempty(problem.rhs) + b = problem.rhs(idx,1); +else + b = zeros(size(A,1),1); +end +if isfield(problem,'ranges') && ~isempty(problem.ranges) + ranges = problem.ranges(idx); + % set equalities: ranges==0 + e(ranges==0) = 0; + % find two-sided inequalities + idx = find(ranges~=0 & ~isinf(ranges)); + ranges = ranges(idx); + rhsH = b(idx); + eR = e(idx); + rhsL = rhsH - (eR<0).*abs(ranges) + (eR==0).*min(ranges,0); + rhsH = rhsH + (eR>0).*abs(ranges) + (eR==0).*max(ranges,0); + % extend problem by ranges + b(idx) = rhsL; + b = [b; rhsH]; + e(idx) = 1; + e = [e; -ones(size(rhsH))]; + A = [A; A(idx,:)]; +end + + +%% read lb and ub +if isfield(problem,'lowerbounds') && ~isempty(problem.lowerbounds) + lb = problem.lowerbounds(:,1); +else + lb = zeros(size(A,2),1); +end +if isfield(problem,'upperbounds') && ~isempty(problem.upperbounds) + ub = problem.upperbounds(:,1); +else + ub = inf(size(A,2),1); +end + + +%% transform LP format to VSDP format +[A,b,c,K,pd] = lp2vsdp(A,b,c,e,lb,ub); + +%___________________________End of MPS2VSDP____________________________ \ No newline at end of file diff --git a/conversion/read_mps.m b/conversion/read_mps.m new file mode 100644 index 0000000..5ec80bc --- /dev/null +++ b/conversion/read_mps.m @@ -0,0 +1,472 @@ +function problem = read_mps(filename) +%% READ_MPS reads a file in MPS format and creates a problem structure +% +% @output problem with the fields: +% - name: name of the problem +% - objsense: 'MIN','MINIMIZE','MAX','MAXIMIZE' +% - objname: name of the objective function rw +% - refrow: name of reference row for SOS +% - rownames: rowcnt x 1 cell array of row names +% - rowtypes: rowcnt x 1 vector of row types ('L','G','N','E') +% - columnnames: colcnt x 1 cell array of column names +% - rhsnames: cell array of names of right hand sides +% - rangenames: cell array of names of ranges +% - boundnames: cell array of names of bounds +% - rhs: rowcnt x length(rhsnames) matrix of right hand sides +% - ranges: rowcnt x length(rangenames) matrix of ranges +% - lowerbounds: colcnt x length(boundnames) matrix of lower bounds +% - upperbounds: colcnt x length(boundnames) matrix of upper bounds +% - bintflags: colcnt x length(boundnames) matrix of bintflags +% bintflag(i,j)=1 for every integer column i in boundset j +% - intflags: numbers of interger columns +% - sosflags: cell of SOS sets, sosflags{1}=[2 i j k] means that +% columns i, j and k are in an SOS2 set +% + +% initiailize default output +problem.name = ''; +problem.objsense = 'MINIMIZE'; +problem.objname = ''; +problem.refrow = ''; + +% initialize variables +[acnt,rowcnt,colcnt,bndcnt,rngcnt,intflag,setsos] = ... + deal(0,0,0,0,0,0,0); +premem = 100; + +% prepare internal variables for output data, preallocate memory +A = zeros(3,2*premem); +colnames = cell(premem,1); +rownames = cell(premem,1); +rowtypes = zeros(premem,1); +rhsnames = cell(0); +rngnames = cell(0); +bndnames = cell(0); +rhs = []; +ranges = []; +intflags = []; +sosflags = cell(0); +lbnds = sparse([]); +ubnds = []; +bintflags = sparse([]); + +% open file +fid = fopen(filename,'r'); +if fid==-1 + error('READ_MPS: cannot open file'); +end + +% first line. +[line,f,lf] = getfields(fid); + + +%% main loop +while 1 + + % check for missing section card. + if line(1)==' ' + warning('VSDP:READ_MPS','Missing Section Card'); + [line,f,lf] = getfields(fid); + continue; + end + + %% handle simple section cards + section = find(strcmp(f{1},{'ENDATA','NAME','OBJSENSE','REFROW','OBJNAME'}),1); + + if ~isempty(section) + switch section + case 1 % ENDATA + problem.rhsnames = rhsnames; + problem.rhs = rhs; + problem.rangenames = rngnames; + problem.ranges = ranges; + problem.intflags = intflags; + problem.sosflags = sosflags; + problem.boundnames = bndnames; + problem.lowerbounds = lbnds; + problem.upperbounds = ubnds; + problem.bintflags = bintflags; + return; + case 2 % NAME + if lf>1 + problem.name = f{2}; + end + case 3 % OBJSENSE + [line,f] = getfields(fid); + problem.objsense = f{1}; + case 4 % REFROW + [line,f] = getfields(fid); + problem.refrow = f{1}; + case 5 % OBJNAME + [line,f] = getfields(fid); + problem.objname = f{1}; + end + + [line,f,lf] = getfields(fid); + continue; + end + + + %% Handle the ROWS Section. + if strcmp(f{1},'ROWS') + [line,f,lf] = getfields(fid); + precnt = premem; % preallocation count + while line(1)==' ' + rowcnt = rowcnt + 1; + rownames{rowcnt} = f{2}; + rowtypes(rowcnt) = upper(f{1}(1)); + if rowcnt==precnt % preallocate memory + precnt = 2 * precnt; + rownames{precnt} = ''; + rowtypes(precnt) = 0; + end + [line,f,lf] = getfields(fid); + end + problem.rownames = rownames(1:rowcnt); + problem.rowtypes = char(rowtypes(1:rowcnt)); + + % create row hash-table + rown = 2 * rowcnt; + rowtbl = cell(rown,1); + for i = 1:rowcnt + rowtbl{hash(rownames{i},rown)}(end+1) = i; + end + + continue; % move on to the next section + end + + + %% handle COLUMNS section. + if strcmp(f{1},'COLUMNS') + colname = ''; + [line,f,lf] = getfields(fid); + precnt = 1; % preallocation count + preacnt = 1; % preallocation count for A + while line(1)==' ' + olf = lf ~= 2*floor(0.5*lf); + + % check number of fields + if lf<2 || lf>5 + error('READ_MPS: wrong number of fields in record'); + end + + % check for 'MARKER' records + idx = 0; + if lf>2 && strcmp('''MARKER''',f{2}) + idx = 3; + elseif lf>3 && strcmp('''MARKER''',f{3}) + idx = 4; + end + if idx>0 + intflag = (intflag + ... + strcmp(f{idx},{'''INTORG''','''INTEND'''})*[1;-1]) > 0; + if strcmp(f{idx},'''SOSORG''') % SOSORG -> find sosflag + setsos = true; + sosflags{end+1}(1) = 1; + if f{1}(1)=='S' + sosflags{end}(1) = str2double(f{1}(2:end)); + end + elseif strcmp(f{idx},'''SOSEND''') + setsos = false; + end + + [line,f,lf] = getfields(fid); + continue; + end + + % new column + if olf && ~strcmp(f{1},colname) + colname = f{1}; + colcnt = colcnt + 1; + colnames{colcnt} = f{1}; + if intflag + intflags{end}(end+1) = colcnt; + end + if setsos + sosflags{end}(end+1) = colcnt; + end + if colcnt==precnt % preallocate memory + precnt = 2 * precnt; + colnames{precnt} = ''; + end + end + + % read entries in record + if strcmp(f{1},'DPLWRB25') + tt = 0; + end + for i = 1+olf:2:lf-1 + rowno = rowtbl{hash(f{i},rown)}; + if length(rowno)>1 + rowno = rowno(strcmp(f{i},rownames(rowno))); + end + if length(rowno)~=1 + error('READ_MPS: COLUMNS entry specifies no unique row'); + end + acnt = acnt + 1; + A(:,acnt) = [rowno; colcnt; str2double(f{i+1})]; + if acnt==preacnt % preallocation + preacnt = 2 * preacnt; + A(3,preacnt) = 0; + end + end + + [line,f,lf] = getfields(fid); % move on to the next record + end + + % all columns are read, add to problem + problem.columnnames = colnames(1:colcnt); + A = A(:,1:acnt)'; + problem.A = sparse(A(:,1),A(:,2),A(:,3),rowcnt,colcnt); + clear A; + + % create column hash-table + coln = 2 * colcnt; + coltbl = cell(coln,1); + for i = 1:colcnt + coltbl{hash(colnames{i},coln)}(end+1) = i; + end + + continue; % to the next section + end + + + %% handle RHS section + if strcmp(f{1},'RHS') + [line,f,lf] = getfields(fid); + + % initial rhs + rhsname = 'DRHS'; % a default RHS name + if lf ~= 2*floor(0.5*lf) + rhsname = f{1}; + end + rhscnt = 1; + rhsno = 1; + rhsnames{rhscnt} = rhsname; + rhs = zeros(rowcnt,rhscnt); + + % read rhs + while line(1)==' ' + olf = lf ~= 2*floor(0.5*lf); + + if olf && ~strcmp(f{1},rhsname) % new rhs + rhsname = f{1}; + rhsno = find(strcmp(rhsname,rhsnames),1); + if isempty(rhsno) + rhscnt = rhscnt + 1; + rhsnames{rhscnt} = rhsname; + rhsno = rhscnt; + rhs(rowcnt,rhscnt) = 0; + end + end + + for i = 1+olf:2:lf-1 % read entries + rowno = rowtbl{hash(f{i},rown)}; + if length(rowno)>1 + rowno = rowno(strcmp(f{i},rownames(rowno))); + end + if length(rowno)~=1 + error('READ_MPS: RHS entry specifies row that does not exist'); + end + rhs(rowno,rhsno) = str2double(f{i+1}); + end + + [line,f,lf] = getfields(fid); % move on to the next record + end + + continue; % to the next section + end + + + %% handle the RANGES section. + if strcmp(f{1},'RANGES') + [line,f,lf] = getfields(fid); + rngname = 'DRNG'; + % default range values for N, E, L, G constraints + rngVAL = zeros(rowcnt,1); + rngVAL(rowtypes=='L' | rowtypes=='G') = Inf; + while line(1)==' ' + olf = lf ~= 2*floor(0.5*lf); + + if isempty(rngnames) || (olf && ~strcmp(f{1},rngname)) % new range + rngname = f{1}; + rngno = find(strcmp(rngname,rngnames),1); + if isempty(rngno) + rngcnt = rngcnt + 1; + rngnames{rngcnt} = rngname; + rngno = rngcnt; + ranges(:,rngno) = rngVAL; % initial range values + end + end + + for i = 1+olf:2:lf-1 + rowno = rowtbl{hash(f{i},rown)}; + if length(rowno)>1 + rowno = rowno(strcmp(f{i},rownames(rowno))); + end + if length(rowno)~=1 + error('READ_MPS: RANGES entry specifies row that does not exist'); + end + ranges(rowno,rngno) = str2double(f{i+1}); + end + + [line,f,lf] = getfields(fid); % move on to the next record + end + + continue; % to the next section + end + + + %% handle the BOUNDS section. + if strcmp(f{1},'BOUNDS') + [line,f,lf] = getfields(fid); + bndname = 'DBOUND'; + lastbndname = ''; + while line(1)==' ' + + switch lf + case 4 + bndname = f{2}; + colname = f{3}; + entry = str2double(f{4}); + case 2 + % only legal if no boundname is given and type is FR or BV + if any(strcmpi(f{1},{'BV','FR'})) + colname = f{2}; + entry = 0; + else + disp(line); + error('READ_MPS: invalid entry in BOUNDS section.'); + end + case 3 + % a BV, FR, MI or PL bound with no entry + if any(strcmpi(f{1},{'BV','FR','MI','PL'})) + bndname = f{2}; + colname = f{3}; + entry = 0; + else % no bound name given + colname = f{2}; + entry = str2double(f{3}); + end + otherwise + disp(line); + error('READ_MPS: wrong number of fields in BOUNDS record'); + end + + if ~strcmp(bndname,lastbndname) % new bounds + lastbndname = bndname; + bndno = find(strcmp(bndname,bndnames),1); + if isempty(bndno) + bndcnt = bndcnt + 1; + bndnames{bndcnt} = bndname; + bndno = bndcnt; + % default bounds + lbnds(colcnt,bndno) = 0; + ubnds(1:colcnt,bndno) = Inf; + ubnds(intflags,bndno) = 1; + bintflags(colcnt,bndno) = 0; + end + end + + colno = coltbl{hash(colname,coln)}; + if length(colno)>1 + colno = colno(strcmp(colname,colnames(colno))); + end + try + if length(colno)~=1 + disp(line); + error('READ_MPS: invalid column in BOUNDS section'); + end + catch dummy + disp('hallo'); + end + + % handle bound types + btype = strfind(['LO','UP','FX','FR','MI','PL','LI',... + 'UI','BV'],upper(f{1})) / 2 - 0.5; + + if isempty(btype) || btype~=round(btype) || length(f{1})~=2 + disp(line); + error('unhandled bound type'); + end + + switch btype + case 0 % LO + lbnds(colno,bndno) = entry; + case 1 % UP + ubnds(colno,bndno) = entry; + case 2 % FX + lbnds(colno,bndno) = entry; + ubnds(colno,bndno) = entry; + case 3 % FR + lbnds(colno,bndno) = -Inf; + ubnds(colno,bndno) = +Inf; + case 4 % MI + lbnds(colno,bndno) = -Inf; + case 5 % PL + ubnds(colno,bndno) = +Inf; + case 6 % LI + lbnds(colno,bndno) = entry; + bintflags(colno,bndno) = 1; + case 7 % UI + ubnds(colno,bndno) = entry; + bintflags(colno,bndno) = 1; + case 8 % BV + ubnds(colno,bndno) = 1; + bintflags(colno,bndno) = 1; + end + + [line,f,lf] = getfields(fid); % move to next record + end + + continue % to the next section + end + + % section that cannot be handled + disp('unhandled section'); + disp(line); + [line,f,lf] = getfields(fid); + while line(1)==' ' + [line,f,lf] = getfields(fid); + end +end + +%_________________________ END OF READ_MPS _____________________________ + + + +%% HELPER FUNCTIONS %% + + +function [line,fields,lf] = getfields(fid) +% GETFIELDS +% [line,fields,lf] = getfields(fid) +% reads the next nonempty, non-comment line of the MPS input file and +% returns the fields of this line in a cell array + +while 1 + line = fgets(fid); + if line==-1 + error('READ_MPS:GETFIELDS','unexpected end of file'); + elseif ~isempty(line) && line(1)~='*' + fields = textscan(line,'%s'); + fields = fields{1}; + lf = length(fields); + if lf>0 + break; + end + end +end + +%_________________________ END OF GETFIELDS ____________________________ + + +function h = hash(key,n) +% HASH +% h = hash(key,n) +% computes a hash value for the 'key' string between 1 and n + +h = sum(key./(27.3:27+length(key))); +h = round((h-floor(h))*(n-1)) + 1; + +%__________________________ END OF HASH _____________________________ \ No newline at end of file diff --git a/conversion/sdpa2vsdp.m b/conversion/sdpa2vsdp.m new file mode 100644 index 0000000..e16f434 --- /dev/null +++ b/conversion/sdpa2vsdp.m @@ -0,0 +1,278 @@ +function [A,b,c,K] = sdpa2vsdp(filename,blksize) +%% SDPA2VSDP: reads data from a SDPA file and converts them to VSDP +% [A,b,c,K] = sdpa2vsdp(sparse-problem-data.dat-s) +% [A,b,c,K] = sdpa2vsdp(dense-problem-data.dat) +% [x,y,INFO] = sdpa2vsdp(optimization-result.out) +% [x0,y0,z0,K] = sdpa2vsdp(sparse-initial-data.ini-s,blksize) +% [x0,y0,z0,K] = sdpa2vsdp(dense-initial-data.ini,blksize) +% +%% >> Input: +% filename: string with the path + filname to the problem +% +%% >> Output: +% A: a nA3 x M Matrix, +% whereas nA3 = dimf+diml+dimq+dims3 +% dimf: number of free variables: dimf = sum(K.f>0) +% diml: number of nonnegative variables: diml = sum(K.l>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims3: sum of all sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% b: a M x 1 vector +% c: a nC x 1 vector, +% whereas nC = dimf+diml+dimq+dims +% dimf, diml, dimq are just the same as for nA3 +% dims = sum_i(K.s(i)^2) +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 25/07/12 M. Lange, rewrite for improved performance +% 04/08/12 M. Lange, add support for initial and full data +% +%% +% TODO: (error when reading large scale problems) +% + +%% open file and import data +if nargin~=1 || ~ischar(filename) || length(filename)<4 + error('VSDP:SDPA2VSDP','input argument must be a filename with extension'); +end + +% type of sdpa-data +extLIST = {'at-s','.dat','.out','ni-s','.ini'}; +dtype = find(strncmpi(filename(end-3:end),extLIST,4),1); +if isempty(dtype) + error('VSDP:SDPA2VSDP','extension "%s" is not supported',filename(end-3:end)); +end + +% open file +fid = fopen(filename,'r'); +if fid==-1 + error('VSDP:SDPA2VSDP','Cannot open %s',filename); +end + +switch dtype + case 1 % '.dat-s' sparse problem data + + % skip comments + str = fgetl(fid); + while str(1)=='*' || str(1)=='"' + str = fgetl(fid); + end + + m = sscanf(str,'%d',1); % number of decision variables + nblocks = fscanf(fid,'%d',1); % number of blocks + blksize = fscanf(fid,'%*[^0-9+-]%d',nblocks); % size of each block + + % right hand side vector (b) + b = -fscanf(fid,'%*[^0-9+-]%lf',m); + + %% read data of A and c + [data cnt] = fscanf(fid,'%*[^0-9+-]%d %d %d %d %lf',[5 inf]); + fclose(fid); + if cnt==0 || mod(cnt,5)~=0 + error('VSDP:SDPA2VSDP','Could not read SDPA data'); + end + [col blk i j data] = deal(data(1,:), data(2,:), data(3,:), data(4,:), -data(5,:)); + + % any i>j ? + idx = find(i>j); + if ~isempty(idx) + warning('VSDP:SDPA2VSDP','Lower triangular elements will be ignored.'); + col(idx) = []; blk(idx) = []; i(idx) = []; j(idx) = []; data(idx) = []; + end + + %% block information + idx = blksize>1; % index vector for sdp blocks + dims = abs(blksize(~idx)); % dimensions of lp blocks + + % cone structure + K.l = sum(dims); % number of all linear variables + K.s = reshape(blksize(idx),[],1); % dimensions of sdp blocks + + %% calculate offset positions + ofs = zeros(1,length(blksize)); % block offsets + ofs(~idx) = cumsum([0 dims(1:end-1)]); % linear block offsets + ofs(idx) = cumsum([K.l; K.s(1:end-1).*(K.s(1:end-1)+1)/2]); + % offsets = i + block-offsets + column-offsets + idx = 0.5*idx'; % factor 0.5 for sdp part column offset: + i = i + ofs(blk) + idx(blk).*j.*(j-1); % overwrite i by offsets + + %% create output + dim3 = K.l + (K.s'*(K.s+1))/2; + j = col>0; % overwrite j by index vector for A + A = sparse(i(j),col(j),data(j),dim3,m); + c = sparse(i(~j),1,data(~j),dim3,1); + + case 2 % '.dat' dense problem data + + % skip comments + str = fgetl(fid); + while str(1)=='*' || str(1)=='"' + str = fgetl(fid); + end + + m = sscanf(str,'%d',1); % number of decision variables + nblocks = fscanf(fid,'%d',1); % number of blocks + blksize = fscanf(fid,'%*[^0-9+-]%d',nblocks); % size of each block + + % right hand side vector (b) + b = -fscanf(fid,'%*[^0-9+-]%lf',m); + + %% block information + idx = blksize>1; % index vector for sdp blocks + dims = abs(blksize(~idx)); % dimensions of lp blocks + + % cone structure + K.l = sum(dims); % number of all linear variables + K.s = reshape(blksize(idx),[],1); % dimensions of sdp blocks + + dim = K.l + K.s'*K.s; + + %% read data of A and c + [c cnt] = fscanf(fid,'%*[^0-9+-]%f',[dim 1]); + if cnt~=dim + fclose(fid); + error('VSDP:SDPA2VSDP','Could not read SDPA data'); + end + [A cnt] = fscanf(fid,'%*[^0-9+-]%f',[dim m]); + if cnt~=dim*m + fclose(fid); + error('VSDP:SDPA2VSDP','Could not read SDPA data'); + end + fclose(fid); + + %% sort data: lp first + if find(idx,1)<=find(~idx,1,'last') + idxs = cell(length(idx),1); + for j = find(idx)' % sdp part + idxs{j} = true(blksize(j)^2,1); + end + for j = find(~idx)' % lp part + idxs{j} = false(abs(blksize(j)),1); + end + idxs = vertcat(idxs{:}); % sdp index + idx = ~idxs; % lp index + c = -[c(idx); c(idxs)]; + A = -[A(idx,:); A(idxs,:)]; + else + c = -c; + A = -A; + end + + case 3 % '.out' solver output + + error('output data not yet supported'); + + case 4 % '.ini-s' sparse initial data + + % cone information from blksize + if nargin<2 || isempty(blksize) + error('VSDP:SDPA2VSDP','block sizes "blksize" has to be set to read initial data'); + end + + % dual initial vector y0 + str = fgetl(fid); + b = sscanf(str,'%*[^0-9+-]%lf',[inf 1]); + + %% read data of x0 and z0 + [data cnt] = fscanf(fid,'%*[^0-9+-]%d %d %d %d %lf',[5 inf]); + fclose(fid); + if cnt==0 || mod(cnt,5)~=0 + error('VSDP:SDPA2VSDP','Could not read SDPA data'); + end + [col blk i j data] = deal(data(1,:), data(2,:), data(3,:), data(4,:), data(5,:)); + + % any i>j ? + idx = i>j; + if any(idx) + warning('VSDP:SDPA2VSDP','Lower triangular elements will be ignored.'); + col(idx) = []; blk(idx) = []; i(idx) = []; j(idx) = []; data(idx) = []; + end + + %% block information + idx = blksize>1; % index vector for sdp blocks + dims = abs(blksize(~idx)); % dimensions of lp blocks + + % cone structure + K.l = sum(dims); % number of all linear variables + K.s = reshape(blksize(idx),[],1); % dimensions of sdp blocks + + %% calculate offset positions + ofs = zeros(1,length(blksize)); % block offsets + ofs(~idx) = cumsum([0 dims(1:end-1)]); % linear block offsets + ofs(idx) = cumsum([K.l; K.s(1:end-1).*(K.s(1:end-1)+1)/2]); + % offsets = i + block-offsets + column-offsets + idx = 0.5*idx'; % factor 0.5 for sdp part column offset: + blk = i + ofs(blk) + idx(blk).*j.*(j-1); % overwrite blk by offsets + + %% create output + dim3 = K.l + (K.s'*(K.s+1))/2; + col = col<2; % overwrite col by index vector for z0 + A = sparse(blk(col),1,data(col),dim3,1); % z0 + data = data .* (1 + (i1; % index vector for sdp blocks + dims = abs(blksize(~idx)); % dimensions of lp blocks + + % cone structure + K.l = sum(dims); % number of all linear variables + K.s = reshape(blksize(idx),[],1); % dimensions of sdp blocks + + dim = K.l + K.s'*K.s; + + %% read data of z0 and x0 + [c cnt] = fscanf(fid,'%*[^0-9+-]%f',[dim 1]); + if (cnt~=dim) + fclose(fid); + error('VSDP:SDPA2VSDP','Could not read SDPA data'); + end + [A cnt] = fscanf(fid,'%*[^0-9+-]%f',[dim 1]); + if (cnt~=dim) + fclose(fid); + error('VSDP:SDPA2VSDP','Could not read SDPA data'); + end + fclose(fid); + + %% sort data: lp first + if (find(idx,1)<=find(~idx,1,'last')) + idxs = cell(length(idx),1); + for j = find(idx)' % sdp part + idxs{j} = true(blksize(j)^2,1); + end + for j = find(~idx)' % lp part + idxs{j} = false(abs(blksize(j)),1); + end + idxs = vertcat(idxs{:}); % sdp index + idx = ~idxs; % lp index + c = [c(idx); c(idxs)]; + A = [A(idx,:); A(idxs,:)]; + end +end diff --git a/conversion/sdpam2vsdp.m b/conversion/sdpam2vsdp.m new file mode 100644 index 0000000..713c903 --- /dev/null +++ b/conversion/sdpam2vsdp.m @@ -0,0 +1,131 @@ +function [A,b,c,K,x,y,z] = sdpam2vsdp(bLOCKsTRUCT,c,F,x0,X0,Y0) +%% VSDP2SDPT3: transforms problem data from SDPAM to VSDP3 format +% [A,b,c,K,x,y,z] = sdpam2vsdp(bLOCKsTRUCT,c,F,x0,X0,Y0) +% +%% >> Input: +% mDIM: scalar - number of primal variables +% nBLOCK: scalar - number of blocks of F +% bLOCKsTRUCT: scalar - represents the block structure of F +% c: vector - coefficient vector +% F: cell array of coefficient matrices +% x0: vector for initial solution +% X0,Y0: cell arrays of initial points +% +%% >> Output: +% A: a nA3 x M Matrix, +% whereas nA = dimf+diml+dimq+dims3 +% dimf: number of free variables: dimf = sum(K.f>0) +% diml: number of nonnegative variables: diml = sum(K.l>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims3: sum of all sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% b: M x 1 vector - right hand side of linear constraints +% c: nA3 x 1 vector - primal objective function +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% x: a nA3 x 1 vector - approx. primal optimal solution - svec(mu=2) +% y: a M x 1 vector - approx. dual optimal solution +% z: a nA3 x 1 vector - approx. dual optimal solution (slack vars) +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 22/08/12 M. Lange, written +% + + +%% prepare data +if nargin<1 || isempty(bLOCKsTRUCT) + error('VSDP:SDPAM2VSDP','"bLOCKsTRUCT" has to be set'); +end + +% initial output +A = []; b = []; x = []; y = []; z = []; + +% indices of diagonal blocks - linear part +indl = find(bLOCKsTRUCT<2); + +% indices of semidefinite blocks +inds = find(bLOCKsTRUCT>1); + +% index vector to speed up svec +Ivec = []; + + +%% create K +K = struct('l',sum(abs(bLOCKsTRUCT(indl))),'s',reshape(bLOCKsTRUCT(inds),[],1)); + + +%% convert c to b +if nargin>1 + b = -c; +end + + +%% convert F to c & A +if nargin>2 && ~isempty(F) + % vectorize linear part, sort blocks + F = [ cellfun(@(x) x(linspace(1,numel(x),length(x))),F(indl,:),... + 'UniformOutput',false); F(inds,:) ]; + + % concat + mdim = size(F,2); + for i = 1:size(F,1) + F{i,1} = reshape(cat(2,F{i,:}),[],mdim); + end + F = -cat(1,F{:,1}); % [c A] in SeDuMi format + + % compact VSDP format + [F Ivec] = vsvec(F,K,1,1); + + % split into c and A + c = F(:,1); + A = F(:,2:end); + clear F; +else + c = []; +end + + +%% write y +if nargin>4 + y = x0; +end + + +%% convert X0 to z +if nargin>3 && ~isempty(X0) && nargout>6 + % vectorize + sort blocks + X0 = [ cellfun(@(x) reshape(x(linspace(1,numel(x),length(x))),[],1),... + X0(indl),'UniformOutput',false) ; ... + cellfun(@(x) x(:),X0(inds),'UniformOutput',false) ]; + + % concat + svec + [z Ivec] = vsvec(cat(1,X0{:}),K,1,1,Ivec); +end + + +%% convert Y0 to x +if nargin>5 && ~isempty(Y0) && nargout>4 + % vectorize + sort blocks + Y0 = [ cellfun(@(x) reshape(x(linspace(1,numel(x),length(x))),[],1),... + Y0(indl),'UniformOutput',false) ; ... + cellfun(@(x) x(:),Y0(inds),'UniformOutput',false) ]; + + % concat + svec + x = vsvec(cat(1,Y0{:}),K,2,1,Ivec); % mu=2 +end + + +%__________________________End of SDPAM2VSDP___________________________ \ No newline at end of file diff --git a/conversion/sdpt2vsdp.m b/conversion/sdpt2vsdp.m new file mode 100644 index 0000000..10204f8 --- /dev/null +++ b/conversion/sdpt2vsdp.m @@ -0,0 +1,148 @@ +function [K,A,c,x,z] = sdpt2vsdp(blk,A,c,x,z) +%% SDPT2VSDP - converts data in SDPT3 input format to VSDP format +% [K,A,c,x,z] = sdpt2vsdp(blk,At,c,x,z) +% +% except for blk all input parameter are optional +% +%% >> Input: +% blk: a cell array describing the block diagonal structure of SQL data +% At: a cell array with At{1} = A(1:dimf,:), At{2} = A(dimf+1:dimf+diml,:) +% and the same for socp cone, and At{i} = [svec(Ai1) ... svec(Aim)] +% for positive semidefinite cone constraint matrices, where Aij is the +% matrix of j-th block and i-th constraint. +% c: a cell array of matrices of dimensions given by K +% -> correspond to primal objective function +% x: a cell array of matrices of dimensions given by K +% -> approximate primal optimal solution +% z: a cell array of matrices of dimensions given by K +% +%% >> Output: +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% A: a nA3 x M Matrix, +% whereas nA = dimf+diml+dimq+dims +% dimf: number of free variables: dimf = sum(K.f>0) +% diml: number of nonnegative variables: diml = sum(K.l>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims3: sum of sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% c: an nA3 x 1 vector - svec(mu=1) +% x: an nA3 x 1 vector -> approximate optimal solution - svec(mu=2) +% z: an nA3 x 1 vector -> approximate LMI dual vector - svec(mu=1) +% +%% +% Note that the right hand side of the linear constraints (b) and the +% dual optimal solution vector (y) have the same format in VSDP and SDPT3. +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + + +%% check input +if nargin<1 || isempty(blk) + error('VSDP:SDPT2VSDP','not enough input parameter'); +elseif nargin>1 && iscell(A) && all(min(size(A))~=[0 1]) + error('VSDP:SDPT2VSDP',['this is an old SDPT3 input format '... + 'not supported by this function']); +elseif size(blk,2)>2 && any([blk{:,3}]) + error('VSDP:SDPT2VSDP','unsupported SDPT3 low rank structure input'); +end + + +%% permutation index for cells if necessary + +% declare blk-vector +blkvec = cat(1,blk{:,1}); + +% indices to sort variables f -> l -> q -> s +indf = find(blkvec=='u'); +indl = find(blkvec=='l'); +indq = find(blkvec=='q'); +inds = find(blkvec=='s'); + +ind = [indf; indl; indq; inds]; + +if size(ind,1)~=size(blk,1) + error('VSDP:SDPT2VSDP','unsupported cell format or empty cells'); +end + + +%% create K +K = struct('f',sum([blk{indf,2}]),'l',sum([blk{indl,2}]),... + 'q',[blk{indq,2}]','s',[blk{inds,2}]'); + +% what has to be processed +doA = nargin>1 && iscell(A) && ~isempty(A{1}) && nargout>1; +doC = nargin>2 && iscell(c) && ~isempty(c{1}) && nargout>2; +doX = nargin>3 && iscell(x) && ~isempty(x{1}) && nargout>3; +doZ = nargin>4 && iscell(z) && ~isempty(z{1}) && nargout>4; + + +%% convert c,x,z +if doC || doX || doZ + % vectorize sdp blocks + for i = 1:length(inds) + j = inds(i); % index of current sdp block group + js = blk{j,2}(:); % index set for sdp blocks of j-th sdp group + if length(js)==1 + blkI = triu(true(js)); + else % extract block diagonal + % block diagonal index computation to avoid additional for-loop + blkI = ones((js'*(js+1))/2,1); + ofI = ones(sum(js),1); + ofV = -ofI; + js(end) = []; + ofI(cumsum([2; js])) = [0; 1-js]; + ofV(cumsum([1; js])) = [1; js-1]; + ofV(2) = length(ofV) - (js(1)>1); + blkI(cumsum(cumsum(ofI(1:length(ofV))))) = cumsum(ofV); + blkI = cumsum(blkI); + end + if doC + c{j} = c{j}(blkI); + end + if doX + x{j} = x{j}(blkI); + end + if doZ + z{j} = z{j}(blkI); + end + end +end + +% cell to mat +if doC + c = cat(1,c{ind}); +else + c = []; +end +if doX + x = sscale(cat(1,x{ind}),K,2); +else + x = []; +end +if doZ + z = cat(1,z{ind}); +else + z = []; +end + + +%% convert A +if doA + A = sscale(cat(1,A{ind}),K,sqrt(0.5)); +else + A = []; +end + +%_____________________________End SDPT2VSDP___________________________ \ No newline at end of file diff --git a/conversion/sscale.m b/conversion/sscale.m new file mode 100644 index 0000000..297347f --- /dev/null +++ b/conversion/sscale.m @@ -0,0 +1,126 @@ +function vA = sscale(vA,K,mu) +%% SSCALE: scale operator to scale the off-diagonal elements of sdp blocks +% vA = vscale(vA,K,mu) +% +%% >> Description: +% For a symmetric matrix X it applies: +% vx = vscale(X(:),K,mu) => vx = vX(:), where vX = mu*X + (1-mu)*I +% The same holds valid for matrices in the SDPT3 format. (-> see svec) +% +%% >> Input: +% vA: an nA x M, M x nA, nA3 x M, or M x nA3 matrix, alternatively, +% whereas nA = dimf+diml+dimq+dims, nA3 = dimf+diml+dimq+dims3 +% dimf: number of free variables: dimf = sum_i(K.f(i)>0) +% diml: number of nonnegative variables: diml = sum_i(K.l(i)>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims: sum of all sdp variables: dims = sum_i(K.s(i)^2) +% dims3: sum of sdp variables: dims = sum_i(K.s(i)*(K.s(i)+1)/2) +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% mu: scaling factor for off-diagonal elements +% for performance reasons and reversibilty mu=0 is not allowed +% +%% >> Output: +% vA: matrix of same dimension as input vA +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 10/07/12 M. Lange, rewrite for faster indexing and both formats +% 28/07/12 M. Lange, speed improvements + allow scaling along rows +% + +%% check input parameter +if nargin~=3 || ~isstruct(K) + error('VSDP:VSCALE','all input parameters have to be set\n'); +end + +% get problem data dimensions +nos = 0; % number of variables that are not in SDP-cone +fields = isfield(K,{'f','l','q','s'}); +if fields(1) + nos = sum(K.f); +end +if fields(2) + nos = nos + sum(K.l); +end +if fields(3) + nos = nos + sum(K.q); +end +if fields(4) + K.s = reshape(K.s(K.s>0),[],1); + dim = nos + sum(K.s.*K.s); + dim3 = nos + sum(K.s.*(K.s+1))/2; +else + dim = nos; + dim3 = nos; +end + +% svec or vec formated, column or row-wise +[isize idim] = max(size(vA)); +if all(isize~=[0 dim3 dim]) + error('VSDP:SSCALE','cone dimension does not fit to size of input'); +end + +% anything to do ? +if isize==0 || isempty(mu) || mu==1 || isempty(K.s) + % warning('VSDP:SSCALE','nothing to do - no sdp block found or mu=1\n'); + return; +elseif mu==0 + error('VSDP:SSCALE','Input mu=0 is not allowed'); +end + + +%% create index vector for diagonal entries +if isize==dim % full matrix format + I = zeros(sum(K.s),1); + I(cumsum([2; K.s(1:end-1)])) = K.s; + I(1) = 1; +else % reduced vectorized format + I = ones(sum(K.s),1); +end +I(cumsum(K.s(1:end-1))+1) = (isize==dim3) - K.s(1:end-1); +I = cumsum(I); +I(1) = nos + 1; +I = cumsum(I); + + +%% scale off-diagonal elements +transInt = false; % transposed interval input +if isa(vA,'intval') && idim==1 + vA = vA'; + idim = 2; + transInt = true; +end +switch 2*(nos~=0) + (idim==2) + case 0 % nos==0 && idim==1 + vA = mu*vA + bsxfun(@times,vA,sparse(I,1,1-mu,isize,1)); + case 1 % nos==0 && idim==2 + tmp = vA(:,I); + vA = mu*vA; + vA(:,I) = tmp; + case 2 % nos>0 && idim==1 + vA = [vA(1:nos,:); mu*vA(nos+1:end,:)] + ... + bsxfun(@times,vA,sparse(I,1,1-mu,isize,1)); + otherwise % nos>0 && idim==2 + tmp = vA(:,I); + vA = [vA(:,1:nos) mu*vA(:,nos+1:end)]; + vA(:,I) = tmp; +end +if transInt + vA = vA'; +end + +%________________________________End SSCALE______________________________ \ No newline at end of file diff --git a/conversion/vsdp12vsdp.m b/conversion/vsdp12vsdp.m new file mode 100644 index 0000000..256cfb8 --- /dev/null +++ b/conversion/vsdp12vsdp.m @@ -0,0 +1,89 @@ +function [K,A,c,x,z] = vsdp12vsdp(blk,At,C,Xt,Zt) +%% VSDP2SEDUMI: transforms problem data from VSDP 0.1 to new VSDP format +% [K,A,c,x,z] = vsdp12vsdp(blk,At,C,Xt,Zt) +% +%% >> Input: +% blk: a cell array describing the block diagonal structure of problem data +% At: a cell array with At{i,1:m} = [svec(Ai1) ... svec(Aim)] +% for positive semidefinite cone constraint matrices, where Aij is the +% matrix of j-th block and i-th constraint +% C: a cell array of matrices of dimensions given by blk +% Xt: a cell array of matrices, initial primal solution +% Zt: a cell array of matrices, initial dual solution +% +%% >> Output: +% A: a dims3 x M Matrix, +% whereas dims3: sum of all sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% c: dims3 x 1 vector for primal objective function +% x: dims3 x 1 vector, initial primal solution +% z: dims3 x 1 vector, initial dual solution (slack vars) +% K: a structure with following fields +% - K.s lists the dimensions of semidefinite blocks +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 13/04/11 V. Haerter, comments added +% 24/07/12 M. Lange, rewritten +% 11/14/12 M. Lange, vectorization for speed-up +%% + +% create K structure +if nargin<1 || ~iscell(blk) || isempty(blk{1,2}) + error('VSDP:VSDP12VSDP','cannot read "blk"'); +else + K.s = cat(1,blk{:,2}); +end + +% initial output +A = []; +c = []; +x = []; +z = []; + +% index vector to speed up svec +Ivec = []; + +% coefficient matrix +if nargin>1 && iscell(At) && ~isempty(At{1}) + [mblk nblk] = size(At); + A = cell(nblk,1); + for i = nblk:-1:1 + % not memory but runtime efficient + A{i} = reshape(cat(2,At{:,i}),[],mblk); + At(:,i) = []; + end + [A Ivec] = vsvec(cat(1,A{:}),K,1,1); +end + +% primal objective vector +if nargin>2 && iscell(C) && ~isempty(C{1}) + C = cellfun(@(x) x(:),C,'UniformOutput',false); + [c Ivec] = vsvec(cat(1,C{:}),K,1,1,Ivec); + clear C; +end + +% initial primal solution vector +if nargin>3 && iscell(Xt) && ~isempty(Xt{1}) + Xt = cellfun(@(x) x(:),Xt,'UniformOutput',false); + [x Ivec] = vsvec(cat(1,Xt{:}),K,2,1,Ivec); % mu=2 + clear Xt; +end + +% initial dual solution vector (slack variables) +if nargin>4 && iscell(Zt) && ~isempty(Zt{1}) + Zt = cellfun(@(x) x(:),Zt,'UniformOutput',false); + z = vsvec(cat(1,Zt{:}),K,1,1,Ivec); + clear Zt; +end + +%_____________________________END OF VSDP12VSDP_________________________ \ No newline at end of file diff --git a/conversion/vsdp2sdpam.m b/conversion/vsdp2sdpam.m new file mode 100644 index 0000000..63c76e7 --- /dev/null +++ b/conversion/vsdp2sdpam.m @@ -0,0 +1,153 @@ +function [mDIM,nBLOCK,bLOCKsTRUCT,c,F,x0,X0,Y0] = vsdp2sdpam(A,b,c,K,x,y,z) +%% VSDP2SDPT3: transforms problem data from VSDP3 to SDPAM format +% [mDIM,nBLOCK,bLOCKsTRUCT,c,F,x0,X0,Y0] = vsdp2sdpt3(At,b,c,K) +% +%% >> Input: +% A: a nA3 x M Matrix, +% whereas nA = dimf+diml+dimq+dims3 +% dimf: number of free variables: dimf = sum(K.f>0) +% diml: number of nonnegative variables: diml = sum(K.l>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims3: sum of all sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% b: M x 1 vector - right hand side of linear constraints +% c: nA3 x 1 vector - primal objective function +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% x: a nA3 x 1 vector - approx. primal optimal solution - svec(mu=2) +% y: a M x 1 vector - approx. dual optimal solution +% z: a nA3 x 1 vector - approx. dual optimal solution (slack vars) +% +%% >> Output: +% mDIM: scalar - number of primal variables +% nBLOCK: scalar - number of blocks of F +% bLOCKsTRUCT: scalar - represents the block structure of F +% c: vector - coefficient vector +% F: cell array of coefficient matrices +% x0: vector for initial solution +% X0,Y0: cell arrays of initial points +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 08/08/11 V. Haerter, corrected errors in transforming linear blocks +% 14/08/12 M. Lange, rewrite for improved performance and new functions +% + + +%% check input +if nargin<4 || ~isstruct(K) || isempty(A) || isempty(b) || isempty(c) + error('VSDP:VSDP2SDPAM','the given problem is incomplete'); +elseif nargin<5 + x = []; y = []; z = []; +elseif nargin<6 + y = []; z = []; +elseif nargin<7 + z = []; +end + +c = c(:); x = x(:); y = y(:); z = z(:); + + +%% prepare data + +% prepare K +fields = isfield(K,{'f','l','q','s'}); +if fields(1) && sum(K.f)>0 + error('VSDP:VSDP2SDPAM','free variables are not supported by SDPA'); +elseif fields(3) && sum(K.q)>0 + error('VSDP:VSDP2SDPAM','SOCP cone is not supported by SDPA'); +elseif fields(2) + K.l = sum(K.l); +else + K.l = 0; +end +if fields(4) + K.s = reshape(K.s(K.s>0),[],1); +else + K.s = []; +end + +% non-compact vectorized format +xdim = ~isempty(x); +zdim = ~isempty(z); +[c Imat] = vsmat(c,K,1,1); +if zdim + c = [vsmat(z,K,1,1,Imat) c]; +end +if xdim + c = [vsmat(x,K,0.5,1,Imat) c]; +end + +% A = [x z c A] +if size(A,1)0)+1:nBLOCK + nj = K.s(j-(K.l>0)); + At = reshape(A(blkj+1:blkj+nj*nj,:),nj,[]); + blkj = blkj + nj*nj; + Y0{j} = At(:,1:xdim*nj); + X0{j} = At(:,xdim*nj+1:(xdim+zdim)*nj); + blki = xdim + zdim; + for i = 1:mDIM+1; + F{j,i} = At(:,blki+1:blki+nj); + blki = blki + nj; + end +end + +% empty unused cell arrays +if ~zdim + X0 = []; +end +if ~xdim + Y0 = []; +end + +%__________________________End of VSDP2SDPAM___________________________ \ No newline at end of file diff --git a/conversion/vsdp2sdpt3.m b/conversion/vsdp2sdpt3.m new file mode 100644 index 0000000..0fdbe37 --- /dev/null +++ b/conversion/vsdp2sdpt3.m @@ -0,0 +1,430 @@ +function [blk, At, ct, xt, zt] = vsdp2sdpt3(K,A,c,x0,z0,opts) +%% VSDP2SDPT3: transforms problem data from VSDP to SDPT3 format +% [blk,At,ct,xt,zt] = vsdp2sdpt3(K,A,c,x0,z0,opts) +% +% except for K all input parameter are optional +% +%% >> Input: +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% A: a nA3 x M Matrix, +% whereas nA = dimf+diml+dimq+dims +% dimf: number of free variables: dimf = sum(K.f>0) +% diml: number of nonnegative variables: diml = sum(K.l>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims3: sum of sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% c: nC x 1 vector - primal objective function +% x0: a nC x 1 vector - approx. primal optimal solution +% z0: a nC x 1 vector - approx. dual optimal solution (slack vars) +% opts: structure for additional parameter settings: +% regarded fields: +% 'SDPT3_VERSION' Version number of SDPT3 solver. +% - default: 4.0 +% 'MIN_SDPBLK_SIZE' Minimum size of an sdp block. For efficiency +% reasons SDPT3 allows to group smaller sdp +% blocks. Every sdp block with a dimension that is +% smaller MIN_SDPBLK_SIZE is considered as a small +% block that should be grouped. +% - default: 50 +% +%% >> Output: +% blk: a cell array describing the block diagonal structure of problem data +% At: a cell array with At{1} = A(1:dimf,:), At{2} = A(dimf+1:dimf+diml,:) +% and the same for socp cone, and At{i} = [svec(Ai1) ... svec(Aim)] +% for positive semidefinite cone constraint matrices, where Aij is the +% matrix of j-th block and i-th constraint +% ct: a cell array of matrices of dimensions given by K +% xt: a cell array of matrices of dimensions given by K +% zt: a cell array of matrices of dimensions given by K +% Is: index vector which describes the reorganization of the sdp blocks to +% group small sdp blocks +% +%% +% Note that the right hand side of the linear constraints (b) and the +% dual optimal solution vector (y) have the same format in VSDP and SDPT3. +% +% This is a private conversion function. The dimension of the input +% parameter is not checked. This has to be done in the calling function. +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 30/07/12 M. Lange, rewrite for new parameter and lower memory usage +% +%% ToDo +% - code reduction +% + + +%% input parameter + +% check input +if nargin<1 || ~isstruct(K) + error('VSDP:VSDP2SDPT3','not enough input parameter'); +elseif nargin<2 + A = []; c = []; x0 = []; z0 = []; opts = []; +elseif nargin<3 + c = []; x0 = []; z0 = []; opts = []; +elseif nargin<4 + x0 = []; z0 = []; opts = []; +elseif nargin<5 + z0 = []; opts = []; +elseif nargin<6 + opts = []; +end + +% global VSDP setting +global VSDP_OPTIONS; + +% version number +version = 4.0; % assumed version number +if isfield(opts,'SDPT3_VERSION') + version = opts.SDPT3_VERSION; +elseif isfield(VSDP_OPTIONS,'SDPT3_VERSION') + version = VSDP_OPTIONS.SDPT3_VERSION; +end + +% max blocksize for 'small' sdp blocks +blksize = []; +if version<4.0 % ignore MIN_SDPBLK_SIZE if SOLVER_VERSION < 4.0 + blksize = 0; +elseif isfield(opts,'MIN_SDPBLK_SIZE') + blksize = opts.MIN_SDPBLK_SIZE; +elseif isfield(VSDP_OPTIONS,'MIN_SDPBLK_SIZE') + blksize = VSDP_OPTIONS.MIN_SDPBLK_SIZE; +end + +% create sdp block group if blk3s > blk3size = blksize*(blksize+1)/2 +if isempty(blksize) + blksize = 2500; +else + blksize = blksize * (blksize+1) / 2; +end + + +%% preparation + +% get problem data dimensions +fields = isfield(K,{'f','l','q','s'}); +if fields(1) + K.f = sum(K.f); +else + K.f = 0; +end +if fields(2) + K.l = sum(K.l); +else + K.l = 0; +end +if fields(3) + K.q = reshape(K.q(K.q>0),[],1); +else + K.q = []; +end +if fields(4) + K.s = reshape(K.s(K.s>0),[],1); +else + K.s = 0; +end +dim = K.f + K.l + sum(K.q) + sum(K.s.*K.s); +dim3 = K.f + K.l + sum(K.q) + sum(K.s.*(K.s+1))/2; +n = (K.f>0) + (K.l>0) + length(K.q) + length(K.s); + +% convert to appropriate format +if length(A)~=dim3 + A = vsvec(A,K,sqrt(2),0); +else + A = sscale(A,K,sqrt(2)); +end +Imat = []; +if length(c)~=dim + [c Imat] = vsmat(c,K,2,0); +end +if all(length(x0)~=[0 dim]) % vx0 = svec(x0(:),K,2), !! mu=2 !! + [x0 Imat] = vsmat(x0,K,1,0,Imat); +end +if all(length(z0)~=[0 dim]) + z0 = vsmat(z0,K,2,0,Imat); +end + +% row- or column-wise +idim = 1 + (size(A,2)==dim3); + +% error if free variables are used with version < 4.0 +if version<4.0 && K.f>0 + error('VSDP:VSDP2SDPT3', ... + 'SDPT3 version < 4.0 does not support unconstrained variables'); +end + +% create cell arrays for output +blk = cell(n,2); +At = cell(n,1); +ct = cell(n,1); +xt = cell(n,1); +zt = cell(n,1); + +% reducing the size of the input data after building the corresponding +% cells saves memory and increases the speed for the remaining data +% access +% the reduction, however, can have a significant effect on cpu-time +% the following is a tradeoff for cpu and memory efficiency +mem3tol = dim3 / min(n+1,3+sqrt(n)); + +% cell index +bli = 1; + +% position indices +blk2e = 0; +blk3e = 0; + + +%% transform unconstrained variables +if K.f>0 + blk{bli,1} = 'u'; + blk{bli,2} = K.f; + if ~isempty(A) && idim==1 + At{bli} = A(blk3e+1:blk3e+K.f,:); + blk3e = blk3e+K.f; + elseif ~isempty(A) + At{bli} = A(:,blk3e+1:blk3e+K.f)'; + blk3e = blk3e+K.f; + end + if ~isempty(c) + ct{bli} = c(blk2e+1:blk2e+K.f); + end + if ~isempty(x0) + xt{bli} = x0(blk2e+1:blk2e+K.f); + end + if ~isempty(z0) + zt{bli} = z0(blk2e+1:blk2e+K.f); + end + blk2e = blk2e + K.f; + bli = bli + 1; +end + + +%% transform linear variables +if K.l>0 + blk{bli,1} = 'l'; + blk{bli,2} = K.l; + if ~isempty(A) && idim==1 + At{bli} = A(blk3e+1:blk3e+K.l,:); + blk3e = blk3e+K.l; + elseif ~isempty(A) + At{bli} = A(:,blk3e+1:blk3e+K.l)'; + blk3e = blk3e+K.l; + end + if ~isempty(c) + ct{bli} = c(blk2e+1:blk2e+K.l); + end + if ~isempty(x0) + xt{bli} = x0(blk2e+1:blk2e+K.l); + end + if ~isempty(z0) + zt{bli} = z0(blk2e+1:blk2e+K.l); + end + blk2e = blk2e + K.l; + bli = bli + 1; +end + + +%% transform socp part +if ~isempty(K.q) + dimq = sum(K.q); + blk{bli,1} = 'q'; + blk{bli,2} = K.q(:)'; + if ~isempty(A) && idim + At{bli} = A(blk3e+1:blk3e+dimq,:); + blk3e = blk3e+dimq; + elseif ~isempty(A) + At{bli} = A(:,blk3e+1:blk3e+dimq)'; + blk3e = blk3e+dimq; + end + if ~isempty(c) + ct{bli} = c(blk2e+1:blk2e+dimq); + end + if ~isempty(x0) + xt{bli} = x0(blk2e+1:blk2e+dimq); + end + if ~isempty(z0) + zt{bli} = z0(blk2e+1:blk2e+dimq); + end + blk2e = blk2e + dimq; + bli = bli + 1; +end + + +%% transform sdp part, consider small blocks +if ~isempty(K.s) + + % save index positions before processing sdp cells + bls = bli; + blk2s = blk2e; + + % starting index for current group + kstart = 1; + + %% create blk + blke = 0; + for k = 1:length(K.s) + blke = blke + K.s(k)*(K.s(k)+1)/2; + + % group great enough or end of sdp data + if blke>blksize || k==length(K.s) + blk{bli,1} = 's'; + blk{bli,2} = reshape(K.s(kstart:k),1,[]); + blke = 0; + kstart = k + 1; + bli = bli + 1; + end + end + + + %% create At cells for sdp part + if ~isempty(A) + % keep sdp part only + if blk3e>0 && idim + A(1:blk3e,:) = []; + blk3e = 0; + elseif blk3e>0 + A(:,1:blk3e) = []; + blk3e = 0; + end + % create cell blocks + if idim==1 + for k = bls:bli-1 + nk3 = sum(blk{k,2}.*(blk{k,2}+1))/2; + At{k} = A(blk3e+1:blk3e+nk3,:); + blk3e = blk3e + nk3; + if blk3e>mem3tol + A = A(blk3e+1:end,:); + blk3e = 0; + end + end + else + for k = bls:bli-1 + nk3 = sum(blk{k,2}.*(blk{k,2}+1))/2; + At{k} = A(:,blk3e+1:blk3e+nk3)'; + blk3e = blk3e + nk3; + if blk3e>mem3tol + A = A(:,blk3e+1:end); + blk3e = 0; + end + end + end + end % transformation of A + + + %% block - cell transformation for sdp part of c + if ~isempty(c) + blk2e = blk2s; % reset index position + for k = bls:bli-1 + nk = blk{k,2}; + lk = length(nk); + if lk==1 % single sdp block + ct{k} = reshape(c(blk2e+1:blk2e+nk*nk),nk,nk); + ct{k} = 0.5 * (ct{k} + ct{k}'); + blk2e = blk2e + nk*nk; + else % block diagonal form of sdp group + ct{k} = cell(1,lk); + for i = 1:lk + nki = nk(i); + ct{k}{i} = reshape(c(blk2e+1:blk2e+nki*nki),nki,nki); + ct{k}{i} = 0.5 * sparse(ct{k}{i} + ct{k}{i}'); + blk2e = blk2e + nki*nki; + end + nki = nk(end); + ct{k} = blkdiag(ct{k}{:}); + end + end + end % transformation of c + + + %% block - cell transformation for sdp part of x0 + if ~isempty(x0) + blk2e = blk2s; % reset index position + for k = bls:bli-1 + nk = blk{k,2}; + lk = length(nk); + if lk==1 % single sdp block + xt{k} = reshape(x0(blk2e+1:blk2e+nk*nk),nk,nk); + xt{k} = 0.5 * (xt{k} + xt{k}'); + blk2e = blk2e + nk*nk; + else % block diagonal form of sdp group + xt{k} = cell(1,lk); + for i = 1:lk + nki = nk(i); + xt{k}{i} = reshape(x0(blk2e+1:blk2e+nki*nki),nki,nki); + xt{k}{i} = 0.5 * sparse(xt{k}{i} + xt{k}{i}'); + blk2e = blk2e + nki*nki; + end + nki = nk(end); + xt{k} = blkdiag(xt{k}{:}); + blk2e = blk2e + nki*nki; + end + end + end % transformation of x0 + + + %% block - cell transformation for sdp part of z0 + if ~isempty(z0) + blk2e = blk2s; % reset index position + for k = bls:bli-1 + nk = blk{k,2}; + lk = length(nk); + if lk==1 % single sdp block + zt{k} = reshape(z0(blk2e+1:blk2e+nk*nk),nk,nk); + zt{k} = 0.5 * (zt{k} + zt{k}'); + blk2e = blk2e + nk*nk; + else % block diagonal form of sdp group + zt{k} = cell(1,lk); + for i = 1:lk + nki = nk(i); + zt{k}{i} = reshape(z0(blk2e+1:blk2e+nki*nki),nki,nki); + zt{k}{i} = 0.5 * sparse(zt{k}{i} + zt{k}{i}'); + blk2e = blk2e + nki*nki; + end + nki = nk(end); + zt{k} = blkdiag(zt{k}{:}); + blk2e = blk2e + nki*nki; + end + end + end % transformation of z0 + +end % transform sdp part + + +%% remove unused cells +blk(bli:end,:) = []; +At(bli:end) = []; +ct(bli:end) = []; +xt(bli:end) = []; +zt(bli:end) = []; + +if isempty(At{1}) + At = []; +end +if isempty(c) + ct = []; +end +if isempty(x0) + xt = []; +end +if isempty(z0) + zt = []; +end + +%_____________________________End VSDP2SDPT3___________________________ \ No newline at end of file diff --git a/conversion/vsdp2sedumi.m b/conversion/vsdp2sedumi.m new file mode 100644 index 0000000..7405b69 --- /dev/null +++ b/conversion/vsdp2sedumi.m @@ -0,0 +1,118 @@ +function [A,c,x,z] = vsdp2sedumi(A,c,x,z,K,opts) +%% VSDP2SEDUMI: transforms problem data from VSDP (intern) format to SeDuMi format +% [A c x z] = vsdp2sedumi(A,c,x,z,K,opts) +% +%% >> Input: +% A: an M x nA3 matrix, or nA3 x M matrix in VSDP format +% whereas nA3 = dimf+diml+dimq+dims3 +% dimf: number of free variables: dimf = sum(K.f>0) +% diml: number of nonnegative variables: diml = sum(K.l>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims3: sum of sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% c,x,z: nA3 x 1 vectors, whereas x has been created by svec(mu=2) +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% opts: structure for additional parameter settings: +% regarded fields: +% - 'ALLOW_TRIANGULAR' Some solvers like SeDuMi allow non-symmetric +% input. If 'ALLOW_TRIANGULAR' is set true, VSDP +% is increasing the speed of data conversion by +% processing only a triangular part of the +% matrices. +% -> default: false +% +%% >> Output: +% A: linear constraints matrix (nA x M or M x nA) in SeDuMi format +% c,x,z: nA x 1 vectors in SeDuMi format. The primal ojective vector may be +% created such that lower triangulars of sdp parts are zero. +% +%% +% Note that the right hand side of the linear constraints (b) and the +% vector for primal objective function (c) as well as the cone structure +% (K) have the same format in VSDP and SeDuMi. Thus only At has to be +% converted. +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 11/07/12 M. Lange, added opts parameter and adapted for new functions +% 27/07/12 M. Lange, adapted for new format and modified smat +% + + +%% input parameter +if nargin<5 || ~isstruct(K) + error('VSDP:VSDP2SEDUMI','not enough input parameters\n'); +elseif nargin<6 + opts = []; +end + +global VSDP_OPTIONS; + +sflag = 1; % symmetry flag -> default: symmetric output +if isfield(opts,'ALLOW_TRIANGULAR') + sflag = (opts.ALLOW_TRIANGULAR==false); +elseif isfield(VSDP_OPTIONS,'ALLOW_TRIANGULAR') + sflag = (VSDP_OPTIONS.ALLOW_TRIANGULAR==false); +end + +% get problem data dimensions +fields = isfield(K,{'f','l','q','s'}); +dim = 0; +if fields(1) + dim = sum(K.f); +end +if fields(2) + dim = dim + sum(K.l); +end +if fields(3) + dim = dim + sum(K.q); +end +if fields(4) + K.s = K.s(K.s>0); + dim3 = dim + sum(K.s.*(K.s+1))/2; + dim = dim + sum(K.s.*K.s); +else + dim3 = dim; +end + +% check dimension of input +adim = [0 dim3 dim]; % allowed dimensions +if all(length(A)~=adim) || all(length(c)~=adim) || all(length(x)~=adim) ... + || all(length(z)~=adim) + error('VSDP:VSDP2SEDUMI','cone structure does not match to matrix dimension'); +end + + +%% transform data +% use mu=2 for triangular output else mu=1 +mu = 2-sflag; +Imat = []; +adim = [0 dim]; +if all(length(A)~=adim) + [A Imat] = vsmat(A,K,mu,sflag); +end +if all(length(c)~=adim) + [c Imat] = vsmat(c,K,mu,sflag,Imat); +end +if all(length(x)~=adim) + [x Imat] = vsmat(x,K,0.5,1,Imat); % x has been created by svec(mu=2) +end +if all(length(z)~=adim) + z = vsmat(z,K,1,1,Imat); +end + +%_____________________________End VSDP2SEDUMI___________________________ \ No newline at end of file diff --git a/conversion/vsmat.m b/conversion/vsmat.m new file mode 100644 index 0000000..6ffa1fc --- /dev/null +++ b/conversion/vsmat.m @@ -0,0 +1,153 @@ +function [A, I] = vsmat(vA,K,mu,sflag,I) +%% VSMAT: smat operator for VSDP3 +% A = vsmat(vA,K,mu,sflag) +% [A, I] = vsmat(vA,K,mu,sflag,I) +% +%% >> Description: +% For symmetric matrices X and Y it applies: +% := trace(X*Y) = svec(X,K,sqrt(2))'*svec(Y,K,sqrt(2)). +% The 'smat' operator is the converse operator of 'svec'. +% For vX = svec(X,mu=2) the vector mX = smat(vX,mu=0.5) equals the matrix X. +% +% For a symmetric matrix X let vA=X(:). Then for K.s=n (n=length(X)) +% svec(vA,K,mu) = [x11 mu*x12 ... mu*x1n x22 mu*x23 ...mu*xn3 ... xnn] +% +% it is: vA = smat(svec(vA,K,mu),K,1/mu) +% +%% >> Input: +% vA: an M x nA3 matrix, or nA3 x M matrix, alternatively +% whereas nA3 = dimf+diml+dimq+dims3 +% dimf: number of free variables: dimf = sum_i(K.f(i)>0) +% diml: number of nonnegative variables: diml = sum_i(K.l(i)>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims3: sum of all sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% mu: scaling factor for off-diagonal elements +% for performance reasons and reversibilty mu=0 is not allowed +% sflag: symmetry flag: how to treat coefficient matrices +% - 0 -> do not care about lower triangular part of matrix +% the corresponding elements are filled with zeros +% - 1 -> output will be symmetric +% - default: 1 +% I: index vector for matrix conversion. Setting the index vector "I" avoids +% creation of the index vector and saves some computation time. If +% "sflag" is changed the index vector has to be computed again. +% +%% >> Output: +% A: a M x nA matrix, or nA x M matrix, depending on the input +% whereas nA = dimf+diml+dimq+dims +% dims: sum of all sdp variables: dims = sum_i(K.s(i)^2) +% I: index vector for applied conversion. Can be used to speed up +% following conversions of similar cone variables by smat. +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 09/07/12 M. Lange, rewrite for faster indexing and non-symmetric input +% 28/07/12 M. Lange, speed improvements + adaption for new format +% + +%% check input parameter +if nargin<2 + error('VSDP:VSMAT','Not enough input parameter'); +end +if nargin<3 || isempty(mu) + mu = 1; +end +if nargin<4 || isempty(sflag) + sflag = 1; % assume symmetric case +end +if nargin<5 + I = []; +end + +% get problem data dimensions +nos = 0; % number of variables that are not in SDP-cone +fields = isfield(K,{'f','l','q','s'}); +if fields(1) + nos = sum(K.f); +end +if fields(2) + nos = nos + sum(K.l); +end +if fields(3) + nos = nos + sum(K.q); +end +if fields(4) + K.s = K.s(K.s>0); + dim = nos + sum(K.s.*K.s); + dim3 = nos + sum(K.s.*(K.s+1))/2; +else + dim = nos; + dim3 = nos; +end + + +% column or row-wise +[isize idim] = max(size(vA)); +if isize>3 && all(isize~=[dim3 dim]) + error('VSDP:VSMAT','Cone dimension does not fit to size of input'); +elseif isize<=3 || isize==dim + % nothing to do + A = vA; + return; +end + + +%% mat-transformation index vector +ns = length(K.s); +if length(I)~=dim || (sflag==1 && ~isnumeric(I)) || (sflag~=1 && ~islogical(I)) + I = cell(ns+1,1); % cell to hold index vectors for every sdp block + if sflag==1 % symmetric output + blks = nos + 1; + I{1} = (1:nos)'; + for k = 2:ns+1 + nk = K.s(k-1); + I{k} = tril(repmat((-1:nk-2)',1,nk),-1)+1; + I{k}(1,:) = cumsum([blks 1:nk-1]); + I{k} = reshape(cumsum(I{k},1),[],1); + blks = blks + nk*(nk+1)/2; + end + else % only upper triangular output + I{1} = true(nos,1); + for k = 1:ns + I{k+1} = reshape(triu(true(K.s(k))),[],1); + end + end + % create index vector from cell + I = vertcat(I{:}); +end + + +%% transform matrix +if mu~=1 + vA = sscale(vA,K,mu); +end + +switch 2*(sflag==1)+idim + case 1 % sflag==0 && idim==1 + A(I,:) = vA; + case 2 % sflag==0 && idim==2 + A(:,I) = vA; + case 3 % sflag==1 && idim==1 + A = vA'; % faster than A = vA(I,:) + A = A(:,I)'; + otherwise % sflag==1 && idim==2 + A = vA(:,I); +end + +%________________________________End VSMAT______________________________ \ No newline at end of file diff --git a/conversion/vsvec.m b/conversion/vsvec.m new file mode 100644 index 0000000..78ddb78 --- /dev/null +++ b/conversion/vsvec.m @@ -0,0 +1,162 @@ +function [vA, I] = vsvec(A,K,mu,sflag,I) +%% VSVEC: svec operator for VSDP3 +% vA = vsvec(A,K,mu,sflag) +% [vA, I] = vsvec(A,K,mu,sflag,I) +% +%% >> Description: +% For symmetric matrices X and Y it applies: +% := trace(X*Y) = svec(X,K,sqrt(2))'*svec(Y,K,sqrt(2)). +% Don't use mu = sqrt(2) for verified computations. Instead, +% to avoid rounding errors use: svec(X,K,1)'*svec(Y,K,2). +% +% For a symmetric matrix X let A=X(:). Then for K.s=n (n=length(X)) +% svec(A,K,mu) = [x11 mu*x12 x22 mu*x13 mu*x23 x33 ... mu*x1n ... xnn] +% +%% >> Input: +% A: a M x nA matrix, or nA x M matrix, alternatively +% whereas nA = dimf+diml+dimq+dims +% dimf: number of free variables: dimf = sum_i(K.f(i)>0) +% diml: number of nonnegative variables: diml = sum_i(K.l(i)>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims: sum of all sdp variables: dims = sum_i(K.s(i)^2) +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% mu: scaling factor for off-diagonal elements +% for performance reasons and reversibilty mu=0 is not allowed +% sflag: symmetry flag: how to treat coefficient matrices +% - 0 -> matrices have no symmetry and both triangular parts has to be +% processed - this is allowed for real matrices only +% - 1 -> default symmetric case +% I: index cell for matrix conversion. Setting the index vector "I" avoids +% creation of the index vector and saves some computation time. +% +%% >> Output: +% vA: an M x nA3 matrix, or nA3 x M matrix (depends on input) +% whereas nA3 = dimf+diml+dimq+dims3 +% dims3: sum of all sdp variables: dims3 = sum_i(K.s(i)*(K.s(i)+1)/2) +% I: index vector for applied conversion. Can be used to speed up +% following conversions of similar cone variables by svec. +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 09/07/12 M. Lange, rewrite for faster indexing and non-symmetric input +% 30/07/12 M. Lange, speed improvement + new input format +% + +%% preparation + +% check input parameter +if nargin<2 || ~isstruct(K) + error('VSDP:VSVEC','not enough input parameters\n'); +end +if nargin<3 || isempty(mu) + mu = 1; +end +if nargin<4 || isempty(sflag) + sflag = 1; % assume symmetric case +end +if nargin<5 + I = []; +end + +% get problem data dimensions +nos = 0; % number of variables that are not in SDP-cone +fields = isfield(K,{'f','l','q','s'}); +if fields(1) + nos = sum(K.f); +end +if fields(2) + nos = nos + sum(K.l); +end +if fields(3) + nos = nos + sum(K.q); +end +if fields(4) + K.s = K.s(K.s>0); + dim = nos + sum(K.s.*K.s); + dim3 = nos + sum(K.s.*(K.s+1))/2; +else + dim = nos; + dim3 = nos; +end + +% column or row-wise +[isize idim] = max(size(A)); +if isize>3 && all(isize~=[dim3 dim]) + error('VSDP:VSVEC','Cone dimension does not fit to size of input'); +elseif isize<=3 || isize==dim3 + % nothing to do + vA = A; + return; +end + +% input index +if iscell(I) && length(I)==2 && length(I{1})==dim + Ilow = I{2}; I = I{1}; +else + Ilow = []; I = []; +end + + +%% index creation - upper triangular part +ns = length(K.s); +if isempty(I) + I = cell(ns+1,1); + I{1} = true(nos,1); + for k = 1:ns + I{k+1} = reshape(triu(true(K.s(k))),[],1); + end + I = vertcat(I{:}); +end + + +%% index vector for lower triangular parts of sdp blocks +if sflag==0 && isempty(Ilow) + Ilow = cell(ns+1,1); + Ilow{1} = (1:nos)'; + blks = nos + 1; + for k = 2:ns+1 + nk = K.s(k-1); + Ilow{k} = nk * ones(nk*(nk+1)/2,1); + Ilow{k}(cumsum([1 1:nk-1])) = [blks 1:-nk:1-nk*(nk-2)]; + Ilow{k} = cumsum(Ilow{k}); + blks = blks + nk*nk; + end + Ilow = vertcat(Ilow{:}); +end + + +%% remove all rows that are not indexed, regard mu +switch 2*(sflag==1)+idim + case 1 % sflag==0 && idim==1 + vA = 0.5 * (A(I,:) + A(Ilow,:)); + case 2 % sflag==0 && idim==2 + vA = 0.5 * (A(:,I) + A(:,Ilow)); + case 3 % sflag==1 && idim==1 + vA = A(I,:); + otherwise % sflag==1 && idim==2 + vA = A(:,I); +end + +if mu~=1 + vA = sscale(vA,K,mu); +end + +% write index cell +I = {I Ilow}; + +%________________________________End VSVEC______________________________ \ No newline at end of file diff --git a/examples/examNeumaier.mat b/examples/examNeumaier.mat new file mode 100644 index 0000000..b61fb04 Binary files /dev/null and b/examples/examNeumaier.mat differ diff --git a/examples/nb_L1.mat b/examples/nb_L1.mat new file mode 100644 index 0000000..139f1bd Binary files /dev/null and b/examples/nb_L1.mat differ diff --git a/examples/nb_L1free.mat b/examples/nb_L1free.mat new file mode 100644 index 0000000..491c377 Binary files /dev/null and b/examples/nb_L1free.mat differ diff --git a/mysdps.m b/mysdps.m new file mode 100644 index 0000000..05c4d6a --- /dev/null +++ b/mysdps.m @@ -0,0 +1,356 @@ +function [obj,x,y,z,info] = mysdps(A,b,c,K,x0,y0,z0,opts) +%% MYSDPS - interface to several conic solvers for VSDP3 +% [obj,x,y,z,info] = mysdps(A,b,c,K,x0,y0,z0) +% +%% >> Description: +% uses a preselected solver to get an approx. solution of a conic +% problem in the standard primal-dual form: +% +% (P) min c'*x (D) max b'*y +% s.t. A*x = b s.t. z := c - A'*y +% x in K z in K* +% +% where K is a cartesian product of the cones R+, SOCP, PSD. +% For theoretical introduction into verified conic programming see: +% C. Jansson. On Verified Numerical Computations in Convex Programming. +% Japan J. Indust. Appl. Math., 26:337–363, 2009 +% +%% >> Input: +% A: nA x m coefficient matrix in SeDuMi or VSDP internal format +% b: a M x 1 vector +% c: a nA x 1 vector, primal objective +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% x0: a nA x 1 vector - a primal feasible (eps-optimal) solution +% y0: a M x 1 vector - a dual feasible (eps-optimal) solution +% z0: a nA x 1 vector - a dual feasible (eps-optimal) solution (slack vars) +% opts: structure for additional parameter settings: +% regarded fields: +% 'SOLVER' to select one of the supported solvers: +% 'sedumi','sdpt3','sdpa','csdp','sdplr', 'lp_solve','linprog' +% 'USE_STARTING_POINT' decides whether initial starting point shall +% be used or not +% -> default: true +% +%% >> Output: +% obj: obj(1) - primal approx. optimal value +% obj(2) - dual approx. optimal value +% x: a nC x 1 vector - approx. primal optimal solution +% y: a M x 1 vector - approx. dual optimal solution +% z: a nC x 1 vector - approx. dual optimal solution (slack vars) +% info: exit code of conic solver: +% info = 0: normal termination +% info = 1: problem indicated to be primal infeasible +% info = 2: problem indicated to be dual infeasible +% info = 3: problem indicated to be primal and dual infeasible +% info = -1: an error occured +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 27/07/10 V. Haerter, 3 new interfaces for additional solvers, +% now we can use VSDP with overall 5 solvers: +% SDPT3, SDPA, SEDUMI, CSDP and SDPLR +% 28/07/10 V. Haerter, adopted to new VSDP version, 2 new +% functions added: sdpt3_to_vsdp, sdpa_to_vsdp +% 11/08/10 V. Haerter, new interfaces for LP_SOLVe and LINPROG +% +% 06/01/12 V. Haerter updated the lpsolve and linprog interfaces +% 07/01/12 V. Haerter corrected second order cone test/insert test of +% false orientation of the matrix A +% 30/03/12 V. Haerter adopted for use of plain Sedumi input data format. +% 20/07/12 M. Lange, rewrite for new internal format. +% +%% ToDo: +% - rotated quadratic cones +% - "info" for sdpa solver +% + + +%% check input parameter +if nargin<4 || isempty(A) || isempty(b) || isempty(c) || isempty(K) + error('VSDP:MYSDPS','not enpugh input parameter'); +elseif nargin<5 + x0 = []; y0 = []; z0 = []; opts = []; +elseif nargin<6 + y0 = []; z0 = []; opts = []; +elseif nargin<7 + z0 = []; opts = []; +elseif nargin<8 + opts = []; +end + +global VSDP_OPTIONS; % global option settings + +solver = []; % default: allow every solver +if isfield(opts,'SOLVER') + solver = opts.SOLVER; +elseif isfield(VSDP_OPTIONS,'SOLVER') + solver = VSDP_OPTIONS.SOLVER; +end + +useSTART = true; % use starting point +if isfield(opts,'USE_STARTING_POINT') + useSTART = logical(opts.USE_STARTING_POINT); +elseif isfield(VSDP_OPTIONS,'USE_STARTING_POINT') + useSTART = logical(VSDP_OPTIONS.USE_STARTING_POINT); +end + +% import data +[A,dum,b,dum,c,dum,K,x0,y0,z0,IF] = import_vsdp(A,b,c,K,x0,y0,z0); + + +%% get solver number +solverLIST = {'sedumi','sdpt3','sdpa','csdp','sdplr', ... + 'lp_solve','linprog'}; +nsolvers = length(solverLIST); +if ischar(solver) + solverID = find(strncmpi(solver,solverLIST,4)); +elseif isnumeric(solver) + solverID = solver; +else + solverID = []; +end + +% if no solver selected -> allow all solvers +if isempty(solverID) + solverID = 1:nsolvers; +end + +% initialization of default output +obj = [inf -inf]; +x = []; y = []; z = []; +info = -1; + + +%% call solver for problem + +if any(solverID==1) && exist('sedumi','file')==2 + %% call sedumi + + % options parameter for sedumi + OPTIONS = []; + if isfield(opts,'SEDUMI_OPTIONS') + OPTIONS = opts.SEDUMI_OPTIONS; + elseif isfield(VSDP_OPTIONS,'SEDUMI_OPTIONS') + OPTIONS = VSDP_OPTIONS.SEDUMI_OPTIONS; + end + % transform data into sedumi format + [A c] = vsdp2sedumi(A,c,[],[],K,opts); + % call sedumi + [x,y,INFO] = sedumi(A,b,c,K,OPTIONS); + % transform results to vsdp format + if ~isempty(x) + obj(1) = c'*x; + end + if ~isempty(y) + obj(2) = b'*y; + z = c - A*y; + end + [x z] = export_vsdp(IF,K,x,z); + info = INFO.pinf + 2*INFO.dinf; + +elseif any(solverID==2) && exist('sqlp','file')==2 + %% call sdpt3 + + % options parameter for sdpt3 + OPTIONS = []; + if isfield(opts,'SDPT3_OPTIONS') + OPTIONS = opts.SDPT3_OPTIONS; + elseif isfield(VSDP_OPTIONS,'SDPT3_OPTIONS') + OPTIONS = VSDP_OPTIONS.SDPT3_OPTIONS; + end + % use starting point ? + if ~useSTART || isempty(x0) || isempty(y0) || isempty(z0) + x0 = []; y0 = []; z0 = []; + end + % transform data into sdpt3 format + [blk,A,c,x0,z0] = vsdp2sdpt3(K,A,c,x0,z0,opts); + % call sdpt3 solver [with intial point] + [obj,x,y,z,INFO] = sqlp(blk,A,c,b,OPTIONS,x0,y0,z0); + % transform results to vsdp format + export data + [dum,dum,dum,x,z] = sdpt2vsdp(blk,[],[],x,z); + [x z] = export_vsdp(IF,K,x,z); + % save info codes + if isstruct(INFO) + info = INFO.termcode; % SDPT3-4.0 output + else + info = INFO(1); % SDPT3-3.x output + end + +elseif any(solverID==3) && exist('sdpam','file')==2 + %% call sdpa + + % check cones + if K.f>0 || sum(K.q)>0 + error('VSDP:MYSDPS','Second order cone and free variables are not supported by SDPAM'); + end + % use starting point ? + if ~useSTART || isempty(x0) || isempty(y0) || isempty(z0) + x0 = []; y0 = []; z0 = []; + end + % transform data into sdpam format + [mDIM,nBLOCK,bLOCKsTRUCT,ct,F,x0,X0,Y0] = vsdp2sdpam(A,b,c,K,x0,y0,z0); + % options parameter for sdpam + OPTIONS = []; + if isfield(opts,'SDPAM_OPTIONS') + OPTIONS = opts.SDPAM_OPTIONS; + elseif isfield(VSDP_OPTIONS,'SDPAM_OPTIONS') + OPTIONS = VSDP_OPTIONS.SDPAM_OPTIONS; + end + % call csdp solver [with intial] + [dum,x0,X0,Y0,INFO] = sdpam(mDIM,nBLOCK,bLOCKsTRUCT,ct,F,x0,X0,Y0,OPTIONS); + % transform results to vsdp format + export + [dum,dum,dum,dum,x,y,z] = sdpam2vsdp(bLOCKsTRUCT,[],[],x0,X0,Y0); + if ~isempty(x) && ~any(isnan(x)) + obj(1) = c'*x; + end + if ~isempty(y) && ~any(isnan(y)) + obj(2) = b'*y; + end + [x z] = export_vsdp(IF,K,x,z); + info = 0; + + +elseif any(solverID==4) && exist('csdp','file')==2 + %% call csdp + + % check cones + if sum(K.q)>0 + error('VSDP:MYSDPS','Second order cone is not supported by CSDP'); + elseif K.f>0 + disp(['Attention: CSDP supports free variables by converting ' ... + 'them to the difference of positive variables.\n' ... + 'The resulting problem is ill-posed.']); + end + % transform data into sedumi format which can be read by csdp + [A,c,x0,z0] = vsdp2sedumi(A,c,x0,z0,K,struct('ALLOW_TRIANGULAR',false)); + % options parameter for csdp + OPTIONS.check = 0; + if isfield(opts,'CSDP_OPTIONS') + OPTIONS = opts.CSDP_OPTIONS; + elseif isfield(VSDP_OPTIONS,'CSDP_OPTIONS') + OPTIONS = VSDP_OPTIONS.CSDP_OPTIONS; + end + % call csdp solver [with intial] + if useSTART && ~isempty(x0) && ~isempty(y0) && ~isempty(z0) + [x,y,z,INFO] = csdp(A,full(b),full(c),K,OPTIONS,full(x0),full(y0),full(z0)); + else + [x,y,z,INFO] = csdp(A,full(b),full(c),K,OPTIONS); + end + % transform results to vsdp format + export + if ~isempty(x) && ~any(isnan(x)) + obj(1) = c'*x; + end + if ~isempty(y) && ~any(isnan(y)) + obj(2) = b'*y; + end + [x z] = export_vsdp(IF,K,x,z); + if any(INFO==[0 1 2]) + info = INFO; + end + +elseif any(solverID==5) && exist('sdplr','file')==2 + %% call sdplr + + % check cones + if K.f>0 || sum(K.q)>0 + error('VSDP:MYSDPS','Second order cone and free variables are not supported by SDPLR'); + end + % options parameter for sdplr + OPTIONS = []; + if isfield(opts,'SDPLR_OPTIONS') + OPTIONS = opts.SDPLR_OPTIONS; + elseif isfield(VSDP_OPTIONS,'SDPLR_OPTIONS') + OPTIONS = VSDP_OPTIONS.SDPLR_OPTIONS; + end + % transform data into sedumi format, which can be read by sdplr + [A,c,x0] = vsdp2sedumi(A,c,x0,[],K,struct('ALLOW_TRIANGULAR',false)); + % call csdp solver [with intial] + if useSTART && ~isempty(x0) && ~isempty(y0) + [x,y] = sdplr(A,b,c,K,OPTIONS,[],x0,y0); + else + [x,y] = sdplr(A,b,c,K,OPTIONS); + end + % transform results to vsdp format + if ~isempty(x) + obj(1) = c'*x; + end + if ~isempty(y) + obj(2) = b'*y; + z = c - A*y; + end + [x z] = export_vsdp(IF,K,x,z); + info = 0; + +elseif any(solverID==6) && exist('lp_solve','file')==2 + %% call lp_solve + + % test of not supported cones, + if sum(K.q)>0 || sum(K.s)>0 + error('VSDP:MYSDPS','Second order and semidefinite cones are not supported by LPSOLVE'); + end + % call lp_solve solver + [objt x y stat] = lp_solve(-full(c),A',full(b),zeros(length(b),1),[-inf(K.f,1);zeros(K.l,1)]); + if ~isempty(x) + obj(1) = c'*x; + end + if ~isempty(y) + y = -y; + obj(2) = b'*y; + z = c - A*y; + end + info = (stat==2) + 2*(stat==3); + +elseif any(solverID==7) && exist('linprog','file')==2 + %% call linprog + + % options parameter for linprog + OPTIONS = optimset('LargeScale','on'); + if isfield(opts,'LINPROG_OPTIONS') + OPTIONS = opts.LINPROG_OPTIONS; + elseif isfield(VSDP_OPTIONS,'LINPROG_OPTIONS') + OPTIONS = VSDP_OPTIONS.LINPROG_OPTIONS; + end + + % test of not supported cones, + if sum(K.q)>0 || sum(K.s)>0 + error('VSDP:MYSDPS','Second order and semidefinite cones are not supported by LINPROG'); + end + % use starting point ? + if ~useSTART || isempty(x0) + x0 = []; + end + % call linprog solver [with intial point] + [x,tmp,flag,tmp,lambda] = ... + linprog(full(c),[],[],A',full(b),[-inf(K.f,1); zeros(K.l,1)],inf(length(c),1),x0,OPTIONS); + % transform results to vsdp format + if ~isempty(x) + y = -lambda.eqlin; + z = c - A*y; + obj(1) = c'*x; + obj(2) = b'*y; + end + info = [1 2] * (abs(flag)==[2; 3]); + +else + + % no solver called + warning('VSDP:MYSDPS','Solver not found!'); + disp('The selected solver could not be detected!'); + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \ No newline at end of file diff --git a/private/export_vsdp.m b/private/export_vsdp.m new file mode 100644 index 0000000..b02440c --- /dev/null +++ b/private/export_vsdp.m @@ -0,0 +1,106 @@ +function [x z] = export_vsdp(IF,K,x,z) +%% EXPORT_VSDP - transform internal format back to imported problem format +% [x z] = read_vsdp(IF,K,x,z) +% +%% >> Input: +% IF: 'SEDUMI', 'VSDP01' or [] for SeDuMi, old VSDP and new internal VSDP +% format, respectively +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% x,z: primal/dual solution vector in SeDuMi or VSDP internal format +% +%% >> Output: +% x,y: in chosen format +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 17/08/12 M. Lange, written for common data export +% +%% +% TODO: rotated quadratic cones +% + + +%% check input +if nargin<3 || ~isstruct(K) + error('VSDP:EXPORT_VSDP','not enough input parameter or wrong data format'); +elseif nargin<4 + z = []; +end + +%% conversion +if isempty(IF) + %% export internal VSDP format + % export x + if ~isempty(x) + x = vsvec(x,K,2,0); % mu=2, does nothing if already compact vec format + end + % export z + if ~isempty(z) + z = vsvec(z,K,1,0); + end +elseif strcmpi(IF,'VSDP01') + %% export old vsdp format + % check supported cones + if isfield(K,'f') && any(K.f>0) || isfield(K,'l') && any(K.l>0) || ... + isfield(K,'q') && any(K.q>0) + error('VSDP:EXPORT_VSDP','old vsdp format supports only semidefinite cone'); + end + % convert x + if ~isempty(x) + x = vsvec(x,K,2,0); % to verify that internal VSDP format is given + blke = length(x); + for j = length(K.s):-1:1 + nj = K.s(j); + blks = blke - nj*(nj+1)/2 + 1; + Xt{j}(triu(true(nj))) = 0.5 * x(blks:blke); % mu==2 + Xt{j} = reshape(Xt{j},nj,nj); + Xt{j} = Xt{j} + Xt{j}'; + blke = blks - 1; + end + x = Xt; + end + % convert z + if ~isempty(z) + z = vsvec(z,K,1,0); % to verify that internal VSDP format is given + blke = length(z); + for j = length(K.s):-1:1 + nj = K.s(j); + blks = blke - nj*(nj+1)/2 + 1; + Zt{j}(triu(true(nj))) = z(blks:blke); + Xt{j} = reshape(Xt{j},nj,nj); + Zt{j} = Xt{j} + Xt{j}' - diag(sparse(diag(Xt{j}))); + blke = blks - 1; + end + z = Zt; + end +elseif strcmpi(IF,'SEDUMI') + %% export sedumi format + Imat = []; % index vector to improve speed of smat + % convert x + if ~isempty(x) + [x Imat] = vsmat(x,K,0.5,1); + end + % convert z + if ~isempty(z) + z = vsmat(z,K,1,1,Imat); + end +else + %% format not supported or detected + error('VSDP:EXPORT_VSDP','tries to export unsupported format'); +end + +%_____________________________End EXPORT_VSDP____________________________ \ No newline at end of file diff --git a/private/import_vsdp.m b/private/import_vsdp.m new file mode 100644 index 0000000..882bdb1 --- /dev/null +++ b/private/import_vsdp.m @@ -0,0 +1,191 @@ +function [A Arad b brad c crad K x y z IF] = import_vsdp(A,b,c,K,x,y,z) +%% IMPORT_VSDP - allows to read different supported formats +% [A b c K] = read_vsdp(A,b,c,K) +% [A b c K x y z IF] = read_vsdp(A,b,c,K,x,y,z) +% +%% >> Input: +% problem data in SEDUMI, old VSDP, or new VSDP internal format +% +%% >> Output: +% A/Arad: an nA3 x M, +% whereas nA3 = dimf+diml+dimq+dims3 +% dimf: number of free variables: dimf = sum(K.f>0) +% diml: number of nonnegative variables: diml = sum(K.l>0) +% dimq: sum of all socp variables: dimq = sum_i(K.q(i)) +% dims3: sum of all sdp variables: dims3 = K.s'*K.s/2 +% b/brad: an M x 1 vector +% c/crad: an nA3 x 1 vector, +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% x: an nA3 x 1 vector - a primal feasible (eps-optimal) solution, mu=2 +% y: an M x 1 vector - a dual feasible (eps-optimal) solution +% z: an nA3 x 1 vector - a dual feasible (eps-optimal) solution (slack vars) +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 07/08/12 M. Lange, written for common data import +% +%% +% TODO: - rotated quadratic cones +% - import of cells containing mixed intervals +% + + +%% check input +if nargin<4 || isempty(A) || isempty(b) || isempty(c) || isempty(K) + error('VSDP:IMPORT_VSDP','not enough input parameter'); +elseif nargin<5 + x = []; y = []; z = []; +elseif nargin<6 + y = []; z = []; +elseif nargin<7 + z = []; +end + +% default input format +IF = []; % if empty: default internal format + +% import old vsdp format +if iscell(A) % (A,b,c,K,x,y,z) -> (blk,A,C,b,Xt,yt,Zt) + IF = 'VSDP01'; + tmp = K; % b + [K,A,c,x,z] = vsdp12vsdp(A,b,c,x,z); % (blk,A,C,Xt,Zt) + b = tmp; +end + +% index vector to speed-up svec +Ivec = []; + + +%% prepare cone structure +f = 0; l = 0; q = []; s = []; +fields = isfield(K,{'f','l','q','s'}); +if fields(1) + f = sum(K.f); +end +if fields(2) + l = sum(K.l); +end +if fields(3) + q = reshape(K.q(K.q>0),[],1); +end +if fields(4) + s = reshape(K.s(K.s>0),[],1); +end +K = struct('f',f,'l',l,'q',q,'s',s); +% appropriate dimensions +dim3 = f + l + sum(q) + sum(s.*(s+1))/2; + + +%% prepare interval input for b +if isnumeric(b) + b = b(:); brad = sparse(size(b,1),1); +elseif isa(b,'intval') || all(isfield(b,{'mid','rad'})) + brad = b.rad(:); b = b.mid(:); +else + error('VSDP:IMPORT_VSDP','cannot import dual objective "b"'); +end +% check size of b +m = size(b,1); +if size(brad,1)~=m + error('VSDP:IMPORT_VSDP','rad(b) has not the same dimension as mid(b)'); +end + + +%% prepare interval input for A +if isnumeric(A) + Arad = 0; +elseif isa(A,'intval') || all(isfield(A,{'mid','rad'})) + Arad = A.rad; A = A.mid; +else + error('VSDP:IMPORT_VSDP','cannot import coefficient matrix "A"'); +end +% compact vectorized format +if length(A)~=dim3 + IF = 'SEDUMI'; + [A Ivec] = vsvec(A,K,1,0,Ivec); +end +if length(Arad)~=dim3 && any(find(Arad,1)) + Arad = vsvec(Arad,K,1,0,Ivec); +elseif length(Arad)~=dim3 + Arad = sparse(dim3,m); +end +% transposed format +if size(A,2)==dim3 + A = A'; +end +if size(Arad,2)==dim3 + Arad = Arad'; +end +% check size of A +if any(size(A)~=[dim3 m] | size(Arad)~=[dim3 m]) + error('VSDP:IMPORT_VSDP','wrong dimension of coefficient matrix "A"'); +end + + +%% prepare interval input for c +if isnumeric(c) + c = c(:); crad = 0; +elseif isa(c,'intval') || all(isfield(c,{'mid','rad'})) + crad = c.rad(:); c = c.mid(:); +else + error('VSDP:IMPORT_VSDP','cannot import primal objective "c"'); +end +% compact vectorized format +if size(c,1)~=dim3 + [c Ivec] = vsvec(c,K,1,0,Ivec); +end +if size(crad,1)~=dim3 && any(find(crad,1)) + crad = vsvec(crad,K,1,0,Ivec); +elseif size(crad,1)~=dim3 + crad = sparse(dim3,1); +end + + +%% prepare x +if ~isempty(x) + if ~isnumeric(x) + error('VSDP:IMPORT_VSDP','primal solution vector "x" has to be numeric'); + end + x = x(:); + % compact vectorized format, mu=2 + if size(x,1)~=dim3 + [x Ivec] = vsvec(x,K,1,0,Ivec); % Ivec can only be used with mu=1 + x = sscale(x,K,2); + end +end + + +%% prepare y +if ~isnumeric(y) || all(length(y)~=[0 m]) + error('VSDP:IMPORT_VSDP','cannot import dual solution vector "y"'); +end +y = y(:); + + +%% prepare z +if ~isempty(z) + if ~isnumeric(z) + error('VSDP:IMPORT_VSDP','cannot import dual solution "z"'); + end + z = z(:); + % compact vectorized format + if size(z,1)~=dim3 + z = vsvec(z,K,1,0,Ivec); + end +end + +%_____________________________End IMPORT_VSDP____________________________ \ No newline at end of file diff --git a/private/prodsup.m b/private/prodsup.m new file mode 100644 index 0000000..c8d6eb7 --- /dev/null +++ b/private/prodsup.m @@ -0,0 +1,34 @@ +function R = prodsup(Amid,Bmid,Arad,Brad) +% helper function to compute a rigorous upper bound (supremum) for +% R = A*B +% +% R = prodsup(Amid,Bmid) +% R = prodsup(Amid,Bmid,Arad,Brad) +% +% this function is used within VSDP (for efficiency reasons) +% Important! setround(1) is assumed + +% sup(Amid*Bmid) +if ~isreal(Amid) && ~isreal(Bmid) + R = real(Amid)*Bmid + imag(Amid)*(1j*Bmid); +else + R = Amid*Bmid; +end + +% matrices are point-matrices -> finish! +if nargin<4 + return; +end + +% initial interval radius +Rrad = sparse(size(R,1),size(R,2)); + +if ~isempty(find(Arad,1)) % regard Arad + Rrad = Rrad + Arad*(abs(Bmid)+Brad); +end + +if ~isempty(find(Brad,1)) % regard Brad + Rrad = Rrad + abs(Amid)*Brad; +end + +R = R + Rrad; \ No newline at end of file diff --git a/private/resmag.m b/private/resmag.m new file mode 100644 index 0000000..4cbc796 --- /dev/null +++ b/private/resmag.m @@ -0,0 +1,63 @@ +function R = resmag(Amid,Bmid,Cmid,Dmid,Arad,Brad,Crad,Drad) +% helper function to compute a rigorous upper bound for +% R = magnitude(A*B-C*D), +% +% R = resmag(Amid,Bmid,Cmid,Dmid) +% R = resmag(Amid,Bmid,Cmid,Dmid,Arad,Brad,Crad,Drad) +% +% this function is used within VSDP (for efficiency reasons) +% Important! setround(1) is assumed + +if (isreal(Amid) && isreal(Bmid) && isreal(Cmid) && isreal(Dmid)) + % max(sup(A*B-C*D),-inf(A*B-C*D)) = abs(A*B-C*D) + R = Amid*Bmid + Cmid*(-Dmid); + R = max(Amid*(-Bmid) + Cmid*Dmid, R); +else + % extract imaginary part of A + if (isreal(Amid)) + iAmid = sparse(size(Amid,1),size(Amid,2)); + else + iAmid = imag(Amid); + end + % extract imaginary part of C + if (isreal(Cmid)) + iCmid = sparse(size(Cmid,1),size(Cmid,2)); + else + iCmid = imag(Cmid); + end + + % sup(A*B-C*D) - Rrad is used as placeholder + Rrad = real(Amid)*Bmid + iAmid*(1j*Bmid) + ... + (-real(Cmid))*Dmid + iCmid*(-1j*Dmid); + % -inf(A*B-C*D) + R = (-real(Amid))*Bmid + iAmid*(-1j*Bmid) + ... + real(Cmid)*Dmid + iCmid*(1j*Dmid); + % mag(A*B-C*D) + R = abs( max(real(Rrad),real(R)) + 1j*max(imag(Rrad),imag(R)) ); +end + +% matrices are point-matrices -> finish! +if (nargin<8) + return; +end + +% initial interval radius +Rrad = sparse(size(R,1),size(R,2)); + +if (any(any(Arad))) % regard Arad + Rrad = Rrad + Arad*(abs(Bmid)+Brad); +end + +if (any(any(Brad))) % regard Brad + Rrad = Rrad + abs(Amid)*Brad; +end + +if (any(any(Crad))) % regard Crad + Rrad = Rrad + Crad*(abs(Dmid)+Drad); +end + +if (any(any(Drad))) % regard Drad + Rrad = Rrad + abs(Cmid)*Drad; +end + +R = R + Rrad; \ No newline at end of file diff --git a/private/resmidrad.m b/private/resmidrad.m new file mode 100644 index 0000000..52ac161 --- /dev/null +++ b/private/resmidrad.m @@ -0,0 +1,69 @@ +function [R,Rrad] = resmidrad(Amid,Bmid,Cmid,Dmid,Arad,Brad,Crad,Drad) +% helper function to compute a rigorous midpoint radius enclosure for +% R = A*B-C*D +% +% [R,Rrad] = resmidrad(Amid,Bmid,Cmid,Dmid) +% [R,Rrad] = resmidrad(Amid,Bmid,Cmid,Dmid,Arad,Brad,Crad,Drad) +% +% this function is used within VSDP (for efficiency reasons) +% Important! setround(1) is assumed + + +%% calculate inclusion for point-matrices +if isreal(Amid) && isreal(Bmid) && isreal(Cmid) && isreal(Dmid) + % sup(A*B-C*D) + R = Amid*Bmid + Cmid*(-Dmid); + % -inf(A*B-C*D) - Rrad is used as placeholder + Rrad = Amid*(-Bmid) + Cmid*Dmid; + % 0.5*(sup(A*B-C*D)+inf(A*B-C*D)) = mid(A*B-C*D) + R = 0.5*(R-Rrad); + % Rrad = mid(A*B-C*D) - inf(A*B-C*D) + Rrad = R + Rrad; +else + % extract imaginary part of A + if isreal(Amid) + iAmid = sparse(size(Amid,1),size(Amid,2)); + else + iAmid = imag(Amid); + end + % extract imaginary part of C + if isreal(Cmid) + iCmid = sparse(size(Cmid,1),size(Cmid,2)); + else + iCmid = imag(Cmid); + end + + % sup(A*B-C*D) + R = real(Amid)*Bmid + iAmid*(1j*Bmid) + ... + (-real(Cmid))*Dmid + iCmid*(-1j*Dmid); + % -inf(A*B-C*D) - Rrad is used as placeholder + Rrad = (-real(Amid))*Bmid + iAmid*(-1j*Bmid) + ... + real(Cmid)*Dmid + iCmid*(1j*Dmid); + % 0.5*(sup(A*B-C*D)+inf(A*B-C*D)) = mid(A*B-C*D) + R = 0.5*(R-Rrad); + % Rrad = abs(mid(A*B-C*D)-inf(A*B-C*D)) + Rrad = abs(R+Rrad); +end + +% matrices are point-matrices -> finish! +if nargin<8 + return; +end + + +%% extend interval radius +if ~isempty(find(Arad,1)) % regard Arad + Rrad = Rrad + Arad*(abs(Bmid)+Brad); +end + +if ~isempty(find(Brad,1)) % regard Brad + Rrad = Rrad + abs(Amid)*Brad; +end + +if ~isempty(find(Crad,1)) % regard Crad + Rrad = Rrad + Crad*(abs(Dmid)+Drad); +end + +if ~isempty(find(Drad,1)) % regard Drad + Rrad = Rrad + abs(Cmid)*Drad; +end \ No newline at end of file diff --git a/private/spdotK.m b/private/spdotK.m new file mode 100644 index 0000000..d88f439 --- /dev/null +++ b/private/spdotK.m @@ -0,0 +1,315 @@ +function [res, resrad, q] = spdotK(varargin) +%% SPDOTK Dot products in K-fold precision +% +% [res, resrad] = spdotK(A1,x1,...,An,xn,K) +% +% where res is approximation of A1*x1 + A2*x2 + ... , +% or A1.'*x1 + A2.'*x2 + ... respectively (useful for sparse matrices) +% at least one of the factors has to be real +% +% the exact result is enclosured by [res-resrad, res+resrad] +% for K>2 the error vector of res(i) is saved in q{i} +% +% @parameters +% Ai - matrices, size(Ai,2) = length(res) +% xi - vectors, size(xi,1) = size(Ai,1) +% K - precision value: the greater K the more accurate the result +% K has to be an integer >1 +% default: 2 +% +% @dependencies +% for the computation of 'resrad' the setround function is needed, +% [for further information, see INTLAB] +% +% Note! +% this function implements algorithms from +% T. Ogita, S.M. Rump, S. Oishi: Accurate Sum and Dot Product, +% SIAM Journal on Scientific Computing (SISC), 26(6):1955-1988, 2005 +% +%% last modified +% 27/08/12 M.Lange, written +% 05/09/12 M.Lange, added code for full matrices +% + + +% fast check for rounding mode +rnd = getround(); +crnd = rnd; +if rnd~=0 + setround(0); + crnd = 0; +end + +% get precision value +if even(nargin) + K = 2; +else + K = varargin{end}; +end + +% check precision K +if K<2 + error('SPDOTK: invalid parameter K'); +end + +% initialization +factor = 134217729; % splitting factor 2^27+1 + + +%% prepare input + +% check whether Ai*xi or Ai.'*xi shall be computed +szs = zeros(floor(nargin/2),4); % [size(Ai) size(xi); ...] +for i = 2:2:nargin + szs(i/2,:) = [size(varargin{i-1}) size(varargin{i})]; +end +horz = all(szs(:,2)==szs(:,3)) && all(szs(:,1)==szs(1,1)); +vert = ~horz && all(szs(:,1)==szs(:,3)) && range(szs(:,2))==0; + +% check dimesion +if any(szs(:,4)~=1) + error('SPDOTK: second factor must be vector'); +elseif ~horz && ~vert + error('SPDOTK: dimension of factors do not match'); +end + +% concat Ai and xi +if nargin>3 && horz + A = cat(2,varargin{1:2:end-1}); + x = full(cat(1,varargin{2:2:end})); +elseif nargin>3 + A = cat(1,varargin{1:2:end-1}); + x = full(cat(1,varargin{2:2:end})); +else + A = varargin{1}; + x = full(varargin{2}); +end + + +if 10*nnz(A) setround(1) is assumed +% only used in case the sqrt function does not regard the rounding mode +res = sqrt(a); +res = res + realmin * (res.*(-res) > -a); \ No newline at end of file diff --git a/vsdpinfeas.m b/vsdpinfeas.m new file mode 100644 index 0000000..7d13b7b --- /dev/null +++ b/vsdpinfeas.m @@ -0,0 +1,303 @@ +function [infeas x y] = vsdpinfeas(A,b,c,K,choose,x0,y0,z0,opts) +%% VSDPINFEAS - infeasibility check for semidefinite-quadratic-linear programming +% [infeas x y] = vsdpinfeas(A,b,C,K,choose,x0,y0,z0) +% +%% >> Description: +% computes verified upper bound of primal optimal value and rigorous +% enclosure of primal strict feasible (near optimal) solution of a conic +% problem in the standard primal-dual form: +% +% (P) min c'*x (D) max b'*y +% s.t. A*x = b s.t. z := c - A'*y +% x in K z in K* +% +% where K is a cartesian product of the cones R+, SOCP, PSD. +% For theoretical introduction into verified conic programming see: +% C. Jansson. On Verified Numerical Computations in Convex Programming. +% Japan J. Indust. Appl. Math., 26:337–363, 2009 +% +%% >> Input: +% A: nA x m coefficient matrix in SeDuMi or VSDP internal format +% b: a M x 1 vector +% c: a nA x 1 vector, primal objective +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% choose: 'p' or 'd' charachter for primal reps. dual infeasibility check +% x0: a nA x 1 vector - a primal feasible (eps-optimal) solution +% y0: a M x 1 vector - a dual feasible (eps-optimal) solution +% z0: a nA x 1 vector - a dual feasible (eps-optimal) solution (slack vars) +% opts: structure for additional parameter settings: +% fields: +% 'FULL_EIGS_ENCLOSURE' if true code for stronger complete +% eigenvalue enclosure will be applied +% the code is a little bit slower +% -> default: true +% +%% >> Output: +% infeas: 1 if primal infeasibility is proved +% 0 if couldn't prove infeasibility +% -1 if dual infeasibility is proved +% x: a nC x 1 vector - a rigorous Farkas solution if dual +% infeasibility could be proved (primal unbounded ray) +% y: a M x 1 vector - a rigorous Farkas solution if primal +% infeasiblity could be proved (dual unbounded ray) +% +%% EXAMPLE: +% EPS = -0.01; DELTA = 0.1; +% min <[0 0; 0 0],X> +% s.t. <[1 0; 0 0],X> = [EPS]; +% <[0 1; 1 DELTA],X> = [1]; +% X is PSD +% +% b = [EPS; 1]; +% c = zeros(4,1); +% A = [1 0 0; 0 1 DELTA]'; +% choose = 'p'; +% +% [objt,x0,y0,z0,info] = mysdps(A,b,c,K); +% [isinfeas x y] = vsdpinfeas(blk,A,c,b,choose,x0,y0,z0); +% isinfeas = +% 1 +% y = +% -100.5998 +% -0.0060 +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 08/09/12 M. Lange, rewrite for new format +% +%% ToDo +% - rotated quadratic cones +% + + +%% check input +if nargin<7 + error('VSDP:VSDPINFEAS','to few input arguments'); +elseif nargin<9 + opts = []; +end + +% default output +infeas = 0; +x = NaN; y = NaN; + +% check if approximations are applicable +if choose=='p' && (isempty(y0) || any(isnan(y0))) + disp('VSDPINFEAS: not applicable approximations given (NaN)'); + return; +elseif choose=='d' && (isempty(x0) || any(isnan(x0))) + disp('VSDPINFEAS: not applicable approximations given (NaN)'); + return; +end + +% global options structure +global VSDP_OPTIONS; + +% read option parameter +FULL_EIGS_ENCLOSURE = true; +if isfield(opts,'FULL_EIGS_ENCLOSURE') + FULL_EIGS_ENCLOSURE = opts.FULL_EIGS_ENCLOSURE; +elseif isfield(VSDP_OPTIONS,'FULL_EIGS_ENCLOSURE') + FULL_EIGS_ENCLOSURE = VSDP_OPTIONS.FULL_EIGS_ENCLOSURE; +end + +% rounding mode +rnd = getround(); +setround(0); + +% import data +[A,Arad,b,brad,c,crad,K,x0,y0,z0,IF] = import_vsdp(A,b,c,K,x0,y0,[]); + + +%% check primal infeasibility +if choose=='p' + + %% 1.step: rigorous enclosure for z = -A*y + if K.f>0 + % solve dual linear constraints rigorously + y0 = vuls([],[],struct('mid',A(1:K.f,:),'rad',Arad(1:K.f,:)),... + zeros(K.f,1),[],[],y0,[],opts); + if ~isstruct(y0) + disp('VSDPINFEAS: could not verify primal infeasibility'); + setround(rnd); + return; + else + y0rad = y0.rad; + y0 = y0.mid; + end + else + y0rad = sparse(size(y0,1),1); + end + % default rounding mode for verification + setround(1); + % z = -A*y (with free variables) + [z,zrad] = resmidrad(0,0,A,y0,0,0,Arad,y0rad); + + %% 2.step: check positive cost + alpha = prodsup(-b',y0,brad',y0rad); + if alpha>=0 + disp('VSDPINFEAS: could not verify primal infeasibility'); + setround(rnd); + return; + end + + %% 3.step: verified lower bounds on cone eigenvalues + dl = inf; + % bound for linear variables + if K.l>0 + ind = K.f+1:K.f+K.l; + dl = - max(zrad(ind) - z(ind)); + end + % eigenvalue bound for second-order cone variables + ind = K.f + K.l; + for j = 1:length(K.q) + ind = ind(end)+2:ind(end)+K.q(j); + zj1 = zrad(ind(1)-1) - z(ind(1)-1); % -inf(zq(1)) + zjn = abs(z(ind)) + zrad(ind); % sup(abs(zq(2:end))) + zjn = sqrtsup(zjn'*zjn); % sup(||zq(2:end)||) + dl = min(dl,-(zj1+zjn)); % inf(zq(1)-||zq(2:end)||) + end + % finish if dl<0 + if dl<0 + disp('VSDPINFEAS: could not verify primal infeasibility'); + setround(rnd); + return; + end + % eigenvalue bound for semidefinite cone + blke = K.f + K.l + sum(K.q) + sum(K.s.*(K.s+1))/2; + for j = length(K.s):-1:1 + blks = blke - K.s(j)*(K.s(j)+1)/2 + 1; + lambdaMin = bnd4sd(struct('mid',z(blks:blke),... + 'rad',zrad(blks:blke)),1,FULL_EIGS_ENCLOSURE); + if lambdaMin<0 + disp('VSDPINFEAS: could not verify primal infeasibility'); + setround(rnd); + return; + end + dl = min(dl,lambdaMin); + blke = blks - 1; + end + + %% 4.step: final feasibility check, compute result + if dl>=0 + infeas = 1; + if any(y0rad) + y = midrad(full(y0),full(y0rad)); + else + y = y0; + end + end + setround(rnd); + return; + + +%% check dual infeasibility +elseif choose=='d' + + %% 1.step: check c'*x > 0 + [alpha alpharad] = spdotK(c,x0,3); + setround(1); % default rounding mode for verification + alpharad = alpharad + crad'*x0; + if alpharad>=-alpha + disp('VSDPINFEAS: could not verify dual infeasibility'); + setround(rnd); + return; + end + + %% 2.step: inclusion for x such that A*x = 0 + vx = vuls([],[],struct('mid',A,'rad',Arad),0*b,[],[],x0,[],opts); + if ~isstruct(vx) + disp('VSDPINFEAS: could not verify dual infeasibility'); + setround(rnd); + return; + elseif prodsup(vx.mid',c,full(vx.rad)',crad)>=0 % c'*x < 0 ? + % find inclusion such that c'x = alpha + A = cat(2,A,c); + Arad = cat(2,Arad,crad); + b(end+1) = alpha; + brad(end+1) = 0; + vx = vuls([],[],struct('mid',A,'rad',Arad),... + struct('mid',b,'rad',brad),[],[],x0,[],opts); + if ~isstruct(vx) + disp('VSDPINFEAS: could not verify dual infeasibility'); + setround(rnd); + return; + end + end + x0rad = sparse(vx.rad); + x0 = full(vx.mid); + + %% 3.step: verified lower bounds on cone eigenvalues + lb = inf; + % bound for linear variables + if K.l>0 + ind = K.f+1:K.f+K.l; + lb = - max(x0rad(ind) - x0(ind)); + end + % eigenvalue bound for second-order cone variables + ind = K.f + K.l; + for j = 1:length(K.q) + ind = ind(end)+2:ind(end)+K.q(j); + xj1 = x0rad(ind(1)-1) - x0(ind(1)-1); % -inf(xq(1)) + xjn = abs(x0(ind)) + x0rad(ind); % sup(abs(xq(2:end))) + xjn = sqrtsup(xjn'*xjn); % sup(||xq(2:end)||) + lb = min(lb,-(xj1+xjn)); % inf(xq(1)-||xq(2:end)||) + end + % finish if lb<0 + if lb<0 + disp('VSDPINFEAS: could not verify dual infeasibility'); + setround(rnd); + return; + end + % eigenvalue bound for semidefinite cone + blke = K.f + K.l + sum(K.q) + sum(K.s.*(K.s+1))/2; + for j = length(K.s):-1:1 + dind = cumsum(1:K.s(j)); % index for diagonal entries + blk3 = dind(end); + vx = struct('mid',0.5*x0(blke-blk3+1:blke),... + 'rad',0.5*x0rad(blke-blk3+1:blke)); + vx.mid(dind) = 2 * vx.mid(dind); % regard mu=2 for x + vx.rad = vx.rad + vx.rad.*sparse(dind,1,1,blk3,1); + lambdaMin = bnd4sd(vx,1,FULL_EIGS_ENCLOSURE); + if lambdaMin<0 + disp('VSDPINFEAS: could not verify dual infeasibility'); + setround(rnd); + return; + end + lb = min(lb,lambdaMin); + blke = blke - blk3; + end + + %% 4.step: final feasibility check, compute result + if lb>=0 + infeas = -1; + x = export_vsdp(IF,K,midrad(full(x0),full(x0rad))); + end + setround(rnd); + return; + + +%% wrong "choose" parameter +else + + error('VSDP:VSDPINFEAS','"choose" must be p or d'); + +end + + +%______________________________End VSDPINFEAS_____________________________ \ No newline at end of file diff --git a/vsdpinit.m b/vsdpinit.m new file mode 100644 index 0000000..6e22191 --- /dev/null +++ b/vsdpinit.m @@ -0,0 +1,285 @@ +function [stat opts] = vsdpinit(opts,display) +%% VSDPINIT - Initialization and defaults for VSDP3 +% [stat] = vsdpinit(par,val) +% +%% >> Description: +% VSDPINIT can be used to initialize VSDP or setup option settings. +% +% vsdpinit('clear') to clear all settings +% vsdpinit(opts) updates global setting structure with opts structur +% vsdpinit('sdpt3') same as opts.SOLVER = 'sdpt3'; vsdpinit(opts) +% +%% >> Input: +% +% - opts: option structure with the following fields +% +% 'SOLVER' select one of the supported solvers: +% 'sedumi','sdpt3','sdpa','csdp','sdplr', ... +% 'lp_solve','linprog' +% +% 'USE_STARTING_POINT' decides whether initial starting point shall be +% used or not +% +% 'ITER_MAX' maximum number of iterations that can be used +% to perturbate the problem and find a feasible +% solution +% +% 'ALPHA' growing factor for problem perturbation +% +% 'FULL_EIGS_ENCLOSURE' if true code for stronger complete eigenvalue +% enclosure will be applied +% +% 'VERIFY_FULL_LSS' if true the full non-symmetric matrix lss +% enclosure is applied within vuls-function +% +% 'ALLOW_TRIANGULAR' some solvers like SeDuMi allow non-symmetric +% input, if 'ALLOW_TRIANGULAR' is set true VSDP +% is increasing the speed of data conversion by +% processing only a triangular part of the sdp +% block matrices +% +% 'SDPT3_VERSION' version number of SDPT3 solver +% +% 'MIN_SDPBLK_SIZE' minimum size of an sdp block, if block size is +% smaller the blocks will be grouped when +% transforming into SDPT3 format +% +% '$SOLVER$_OPTIONS' option structure for approximate solver +% the field for sdpt3 options is 'SDPT3_OPTIONS' +% +% - display: -1 - display only errors +% 0 - display warnings +% 1 - display warnings and additional information +% +%% >> Output: +% - stat: number of option structure fields which have been detected in +% opts +% +% - opts: copy of global option structure +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/08/10 V. Haerter, comments added +% 09/09/10 V. Haerter, store Matlab version +% 09/09/12 M. Lange, rewrite for new options structure +% +%% ToDo: +% + +%% check input + +% global option settings +global VSDP_OPTIONS; + +% check input parameter +if nargin<1 + opts = []; display = 0; +elseif nargin<2 || isempty(display) + display = 0; +end + +% vsdpinit('solver') or vsdpinit('clear') call +if ischar(opts) + if strcmpi(opts,'clear') + clear GLOBAL VSDP_OPTIONS; + stat = 0; + opts = []; + return; + end + opts = struct('SOLVER',opts); +end + +wng = warning; % old warning status +warning off; + +% get current vsdp path +vsdppath = which('vsdpinit'); +vsdppath(end-9:end) = []; + +% check whether path already exists, add path if necessary +if isempty(strfind(path,vsdppath)) + % add vsdp to paths + addpath(vsdppath); + addpath([vsdppath 'conversion']); + addpath([vsdppath 'examples']); + % refresh path + path(path); +end + +% display warnings +if display<0 + warning off; +else + warning on; +end + + +%% update options +stat = 0; % number of read options fields +wfields = 0; % number of fields that could not be imported + +% solver +solverLIST = {'sedumi','sdpt3','sdpa','csdp','sdplr', ... + 'lp_solve','linprog'}; +if isfield(opts,'SOLVER') + stat = stat + 1; + ofield = opts.SOLVER; + if isnumeric(ofield) && ofield>0 && ofield<=length(solverLIST) + VSDP_OPTIONS.SOLVER = solverLIST(ofield); + elseif isempty(ofield) || any(strcmp(ofield,solverLIST)) + VSDP_OPTIONS.SOLVER = ofield; + else + warning('VSDP:VSDPINIT','selected solver is not supported'); + wfields = wfields + 1; + end +elseif ~isfield(VSDP_OPTIONS,'SOLVER') + warning('VSDP:VSDPINIT','no solver has been selected'); +end + +% use starting point +if isfield(opts,'USE_STARTING_POINT') + stat = stat + 1; + ofield = opts.USE_STARTING_POINT; + if length(ofield)==1 && any(ofield==[0 1]) + VSDP_OPTIONS.USE_STARTING_POINT = ofield; + else + warning('VSDP:VSDPINIT','USE_STARTING_POINT has to be logical'); + wfields = wfields + 1; + end +end + +% maximum iteration in vsdplow and vsdpup +if isfield(opts,'ITER_MAX') + stat = stat + 1; + ofield = opts.MAX_ITER; + if isnumeric(ofield) && length(ofield)==1 && ofield>=1 + VSDP_OPTIONS.MAX_ITER = ofield; + else + warning('VSDP:VSDPINIT','MAX_ITER not accepted'); + wfields = wfields + 1; + end +end + +% growing factor for perturbations +if isfield(opts,'ALPHA') + stat = stat + 1; + ofield = opts.ALPHA; + if isnumeric(ofield) && length(ofield)==1 && ofield~=0 + VSDP_OPTIONS.ALPHA = ofield; + else + warning('VSDP:VSDPINIT','ALPHA has wrong format or is zero'); + wfields = wfields + 1; + end +end + +% full eigenvalue enclosure +if isfield(opts,'FULL_EIGS_ENCLOSURE') + stat = stat + 1; + ofield = opts.FULL_EIGS_ENCLOSURE; + if length(ofield)==1 && any(ofield==[0 1]) + VSDP_OPTIONS.FULL_EIGS_ENCLOSURE = ofield; + else + warning('VSDP:VSDPINIT','FULL_EIGS_ENCLOSURE has to be logical'); + wfields = wfields + 1; + end +end + +% full lss verification +if isfield(opts,'VERIFY_FULL_LSS') + stat = stat + 1; + ofield = opts.VERIFY_FULL_LSS; + if length(ofield)==1 && any(ofield==[0 1]) + VSDP_OPTIONS.VERIFY_FULL_LSS = ofield; + else + warning('VSDP:VSDPINIT','VERIFY_FULL_LSS has to be logical'); + wfields = wfields + 1; + end +end + +% allow non-symmetric, upper triangular format conversion +if isfield(opts,'ALLOW_TRIANGULAR') + stat = stat + 1; + ofield = opts.ALLOW_TRIANGULAR; + if length(ofield)==1 && any(ofield==[0 1]) + VSDP_OPTIONS.ALLOW_TRIANGULAR = ofield; + else + warning('VSDP:VSDPINIT','ALLOW_TRIANGULAR has to be logical'); + wfields = wfields + 1; + end +end + +% SDPT3 version number +if isfield(opts,'SDPT3_VERSION') + stat = stat + 1; + ofield = opts.SDPT3_VERSION; + if isnumeric(ofield) && length(ofield)==1 + VSDP_OPTIONS.SDPT3_VERSION = ofield; + else + warning('VSDP:VSDPINIT','wrong SDPT3_VERSION format'); + wfields = wfields + 1; + end +end + +% blocksize for sdp block groups +if isfield(opts,'MIN_SDPBLK_SIZE') + stat = stat + 1; + ofield = opts.MIN_SDPBLK_SIZE; + if isnumeric(ofield) && length(ofield)==1 + VSDP_OPTIONS.MIN_SDPBLK_SIZE = ofield; + else + warning('VSDP:VSDPINIT','wrong MIN_SDPBLK_SIZE format'); + wfields = wfields + 1; + end +end + +% solver options - not checked +soLIST = cellfun(@(x) [upper(x) '_OPTIONS'],solverLIST,'UniformOutput',false); +for i = 1:length(soLIST) + if isfield(opts,soLIST{i}) + stat = stat + 1; + eval(['VSDP_OPTIONS.' soLIST{i} ' = opts.' soLIST{i} ';']); + end +end + +% output options + statistics +if display>=0 && ~isempty(opts) + fnum = numel(structfun(@(x) 1,opts)); + if stat0 + fprintf('\n___________________________________________________________________\n'); + fprintf('********************** Initialization of VSDP *********************\n'); + fprintf('* *\n'); + fprintf('* VSDP 2012: Verified Conic Programming *\n'); + fprintf('* Viktor Haerter, viktor.haerter@tu-harburg.de *\n'); + fprintf('* Christian Jansson, jansson@tu-harburg.de *\n'); + fprintf('* Marko Lange, marko.lange@tu-harburg.de *\n'); + fprintf('* Institute for Reliable Computing *\n'); + fprintf('* Hamburg University of Technology, *\n'); + fprintf('* Germany, October 2012 *\n'); + fprintf('* *\n'); + fprintf('******************* VSDP is now ready for use *********************'); + fprintf('\n*******************************************************************\n\n'); +end + +% reset warning state +warning(wng); + +%______________________________END OF VSDPINIT____________________________ diff --git a/vsdplow.m b/vsdplow.m new file mode 100644 index 0000000..c0d8a86 --- /dev/null +++ b/vsdplow.m @@ -0,0 +1,286 @@ +function [fL, y, dl, info] = vsdplow(A,b,c,K,x0,y,z0,xu,opts) +%% VSDPLOW - verified lower bound for an SQLP problem +% [fL y dl info] = vsdplow(A,b,c,K,[],y0) +% [fL y dl info] = vsdplow(A,b,c,K,x0,y0,z0,xu,opts) +% +%% >> Description: +% computes verified lower bound of primal optimal value and rigorous +% enclosure of dual strict feasible (near optimal) solution of a conic +% problem in the standard primal-dual form: +% +% (P) min c'*x (D) max b'*y +% s.t. A*x = b s.t. z := c - A'*y +% x in K z in K* +% +% where K is a cartesian product of the cones R+, SOCP, PSD. +% For theoretical introduction into verified conic programming see: +% C. Jansson. On Verified Numerical Computations in Convex Programming. +% Japan J. Indust. Appl. Math., 26:337–363, 2009 +% +%% >> Input: +% A: nA x m coefficient matrix in SeDuMi or VSDP internal format +% b: a M x 1 vector +% c: a nA x 1 vector, primal objective +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% x0: a nA x 1 vector - a primal feasible (eps-optimal) solution +% y0: a M x 1 vector - a dual feasible (eps-optimal) solution +% z0: a nA x 1 vector - a dual feasible (eps-optimal) solution (slack vars) +% xu: upper bounds of eigenvalues or spectral values of primal optimal +% solution x +% opts: structure for additional parameter settings: +% regarded fields: +% 'ITER_MAX' maximum number of iterations that can be used to +% perturbate the problem and find a feasible solution +% 'ALPHA' growing factor for problem perturbation +% -> default: 0.5 +% 'FULL_EIGS_ENCLOSURE' if true code for stronger complete +% eigenvalue enclosure will be applied +% the code is a little bit slower +% -> default: true +% 'SOLVER' to select one of the supported solvers: +% 'sedumi','sdpt3','sdpa','csdp','sdplr', 'lp_solve','linprog' +% 'USE_STARTING_POINT' decides whether initial starting point shall +% be used or not +% -> default: false +% +%% >> Output: +% fL: verified lower bound of the primal optimal value +% y: an M x 1 vector - rigorous enclosure of dual strict feasible solution +% dl: verified lower bounds of eigenvalues or spectral values of z=c-A'*y +% info.iter: number of iterations +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 16/06/12 M. Lange, adapted for new funtions and data format +% 06/08/12 M. Lange, new helper functions for faster computation +% +%% +% TODO: rotated quadratic cones +% + +%% input parameter + +% check number of input arguments +if nargin<6 || isempty(A) || isempty(b) || isempty(c) || isempty(K) + error('VSDP:VSDPLOW','more input arguments are required'); +elseif nargin<8 + xu = []; opts = []; +elseif nargin<9 + opts = []; +end + +global VSDP_OPTIONS; % global options structure + +% use starting point - default: false +if ~isfield(VSDP_OPTIONS,'USE_STARTING_POINT') + VSDP_OPTIONS.USE_STARTING_POINT = false; +end +if ~isfield(opts,'USE_STARTING_POINT') + opts.USE_STARTING_POINT = VSDP_OPTIONS.USE_STARTING_POINT; +end + +% parameter list +optLIST = {'ITER_MAX','ALPHA','FULL_EIGS_ENCLOSURE'}; +optfs = isfield(opts,optLIST); +goptfs = isfield(VSDP_OPTIONS,optLIST); + +% read parameters +[ITER_MAX, ALPHA, FULL_EIGS_ENCLOSURE] = deal(10, 0.5, true); +for i = 1:length(optLIST) + if optfs(i) + eval([optLIST{i},' = opts.',optLIST{i},';']); + elseif goptfs(i) + eval([optLIST{i},' = VSDP_OPTIONS.',optLIST{i},';']); + end +end + + +%% Preliminary steps / Prealocations + +% initial output +fL = -Inf; +dl = NaN; +info.iter = 0; + +% rounding mode +rnd = getround(); +setround(0); + +% import data +[A,Arad,b,brad,c,crad,K,x0,y,z0] = import_vsdp(A,b,c,K,x0,y,z0); + +% check if approximations are applicable +if isempty(y) || any(isnan(y)) + warning('VSDP:VSDPLOW','not applicable approximations given (NaN)'); + y = NaN; + return; +elseif nargin<7 || any(isnan(x0)) || any(isnan(z0)) + x0 = []; z0 = []; +end + +% get problem data dimensions +dim3 = length(c); % dimension +nc = K.l + length(K.q) + length(K.s); % number of cone constraints + +% variable declarations and allocations +I = []; % basic indices +xu = xu(:); +if isempty(xu) + xu = inf(K.f+nc,1); +elseif size(xu,1)~=K.f+nc + error('VSDP:VSDPLOW','upper bound vector has wrong dimension'); +end +xuf = xu(1:K.f); xu(1:K.f) = []; +yrad = sparse(size(y,1),1); % interval radius y + +dl = -inf(nc,1); % dual lower bounds / trace bounds +epsj = ones(nc,1); % factor for perturbation amount +ceps = sparse(dim3,1); % perturbation for c +pertS = cell(length(K.s),1); % for diagonal perturbations of sdp blocks + +% extract free part +Af = A(1:K.f,:); Afrad = Arad(1:K.f,:); +cf = c(1:K.f); cfrad = crad(1:K.f); + +% create index vector for perturbation entries +pertI = ones(sum(K.s),1); +pertI(cumsum(K.s(1:end-1))+1) = 1 - K.s(1:end-1); +pertI = [ones(K.l+(~isempty(K.q)),1); K.q(1:end-1); cumsum(pertI)]; +pertI(1) = K.f + 1; +pertI = cumsum(pertI); + + +%% Algorithm with finite/infinite primal bounds xu +% **** main loop **** +while info.iter<=ITER_MAX + info.iter = info.iter + 1; + setround(1); % default for rigorous computation in steps 1-3 + + %% 1.step: defect computation, free variables handling + if K.f>0 && max(xuf)==inf + % solve dual linear constraints rigorously + [y,I] = vuls([],[],struct('mid',Af,'rad',Afrad),... + struct('mid',cf,'rad',cfrad),[],[],y,I,opts); + if ~isstruct(y) + disp('VSDPLOW: could not find solution of dual equations'); + break; + else + yrad = y.rad; + y = y.mid; + end + end + + % compute rigorous enclosure for z = c - A*y (with free variables) + [z,zrad] = spdotK(c,1,A,-y,2); % for point matrices + zrad = zrad + crad; % regard radii of other parameters + if any(yrad) + zrad = zrad + abs(A)*yrad; + end + if ~isempty(find(Arad,1)) + zrad = zrad + Arad*(abs(y)+yrad); + end + + defect = 0; % defect by free variables + if K.f>0 && max(xuf)0 % bound for linear variables + ind = K.f+1:K.f+K.l; + zjn = zrad(ind) - z(ind); % -inf(zl) + % interprete all lp vars as one cone: + % if any element is negative all constraints close to zero will be + % perturbated + zj1 = max(zjn); + if zj1>0 % there is a negative bound + ind = zjn > -zj1; + zj1 = max((-1e-13)*min(zjn),zj1); + zjn(ind) = min(zjn(ind)+zj1,1.05*zj1); + end + dl(1:K.l) = - zjn; + end + % eigenvalue bound for second-order cone variables + ind = K.f + K.l; + for j = 1:length(K.q) + ind = ind(end)+2:ind(end)+K.q(j); + zj1 = zrad(ind(1)-1) - z(ind(1)-1); % -inf(zq(1)) + zjn = abs(z(ind)) + zrad(ind); % sup(abs(zq(2:end))) + zjn = sqrtsup(zjn'*zjn); % sup(||zq(2:end)||) + dl(K.l+j) = -(zj1+zjn); % inf(zq(1)-||zq(2:end)||) + end + % eigenvalue bound for semidefinite cone + blke = K.f + K.l + sum(K.q) + sum(K.s.*(K.s+1))/2; + for j = length(K.s):-1:1 + ofs = K.l + length(K.q) + j; + blks = blke - K.s(j)*(K.s(j)+1)/2 + 1; + [lmin,dl(ofs),pertS{j}] = bnd4sd(struct('mid',z(blks:blke),... + 'rad',zrad(blks:blke)),1,FULL_EIGS_ENCLOSURE); + if lmin>0 + dl(ofs) = lmin; + end + pertS{j} = epsj(ofs) * pertS{j}(:); + blke = blks - 1; + end + + %% 3.step: cone feasibility check, computing lower bound + dli = find(dl<0); + if ~any(isinf(xu(dli))) + % inf(min(dl,0)*xu + b'*y - defect) + fL = -(sum(dl(dli)'*(-xu(dli))) + prodsup(-b',y,brad',yrad) + defect); + y = midrad(full(y),full(yrad)); + setround(rnd); % reset rounding mode + return; + end + + %% 4.step: create some perturbed problem + setround(0); % no code for rigorous computations + ind = 1:K.l+length(K.q); + if isempty(ind) + ceps = ceps + sparse(pertI,1,cat(1,pertS{:}),dim3,1); + else + ceps = ceps + sparse(pertI,1,cat(1,epsj(ind).*min(dl(ind),0),... + pertS{:}),dim3,1); + end + if any(isinf(ceps)) || any(isnan(ceps)) + disp('VSDLOW: perturbation extended range'); + break; + end + epsj(dli) = epsj(dli) * (1+ALPHA); % update perturbation factor + + %% 5.step: solve the perturbed problem + clear dli ind z zrad; % free some memory before calling solver + [obj,x0,y,z0,INFO] = mysdps(A,b,c+ceps,K,x0,y,z0,opts); + % if could not found solution or dual infeasible, break + if isempty(y) || any(isnan(y)) || any(isinf(y)) || any(INFO(1)==[(1) 2 3]) + disp('VSDLOW: conic solver could not find solution for perturbed problem'); + break; + end +end +% **** end of main loop **** + +% reset rounding mode +setround(rnd); + +% write output +if info.iter==ITER_MAX + disp('VSDPLOW: maximum number of iterations reached'); +end +y = NaN; fL = -Inf; dl = NaN; + +%________________________________End VSDPLOW______________________________ \ No newline at end of file diff --git a/vsdpneig.m b/vsdpneig.m new file mode 100644 index 0000000..5456371 --- /dev/null +++ b/vsdpneig.m @@ -0,0 +1,345 @@ +function lambda = vsdpneig(A,hermit) +% VSDPNEIG rigorous enclosure of all eigenvalues of a normal matrix. +% Input matrix A has to be normal. +% A may also be an interval matrix, full or sparse. This interval +% matrix can be defined as: +% 1. a structure with the fields 'mid' and 'rad', +% 2. a cell array, where the first entry will be interpreted +% as 'mid' and the second as 'rad', +% 3. a class of type 'intval' [see INTLAB toolbox]. +% +% @params +% A: Normal input matrix. Normality is not checked. However, if the +% condition of the approximative eigen-vector matrix is too bad +% the function will throw a warning. +% +% hermit: Can be used to explicitly force the function to +% treat A as a hermitian matrix. If not set and A is not +% equal to A' the matrix will be treated as an +% arbitrary normal matrix. +% +% @output +% lambda: If A is of type 'intval' lambda will be of the same type, +% else lambda is a structure with the fields 'mid' and 'rad' +% which are describing the corresponding interval vector. +% +% @dependencies: 'setround'-function [see Intlab] +% +% Replaces "veigsym" by Christian Jansson +% + +% +% written 19/03/12 Marko Lange +% modified 20/04/12 - allow interval structure input +% - removed most dependencies, still need 'setround' +% function +% modified 07/08/12 - allow compact vectorized input +% + +% +% Example: +% +% A = [1 2 3; 2 1 4; 3 4 5]; +% A = midrad(A, 0.01*ones(3)); +% lambda = vsdpneig(A,true) +% intval lambda = +% [ -1.5125, -1.4615] +% [ -0.6172, -0.5679] +% [ 9.0506, 9.1084] +% +% or: +% +% A.mid = [1 2 3; 2 1 4; 3 4 5]; +% A.rad = 0.01*ones(3); +% lambda = vsdpneig(A,true) +% lambda = +% mid: [3x1 double] +% rad: [3x1 double] +% lambda.mid +% ans = +% -1.4870 +% -0.5925 +% 9.0795 +% lambda.rad +% ans = +% 0.0254 +% 0.0246 +% 0.0289 +% + + +%% check input data %% + +% input matrix +if nargin<1 || isempty(A) + error('VSDPNEIG: No input matrix set.'); +elseif ~(isnumeric(A) || isa(A,'intval') || all(isfield(A,{'mid','rad'}))) + error('VSDPNEIG: Datatype of input matrix is not accepted.'); +end % size & dimension is not checked. Error will occur when calling eig. + +if nargin<2 || ~islogical(hermit) + hermit = false; +end + + +%% prepare input data %% + +% get rounding mode and set round to nearest +rnd = getround(); +setround(0); % for preparation only + +% is user using INTLAB? +isintval = isa(A,'intval'); + +% prepare interval input +if ~isnumeric(A) % A is an interval structure or of type 'intval' + Arad = A.rad; + A = A.mid; +else + Arad = sparse(size(A,1),size(A,2)); +end + +%% prepare A +if size(A,1)~=size(A,2) && size(A,2)~=1 + A = A(:); + Arad = Arad(:); +end +n = size(A,1); + +if size(A,2)==1 % vectorized input + hermit = true; + n = sqrt(2*n+0.25)-0.5; + if n~=round(n) + error('VSDPNEIG: Wrong dimension of vectorized input'); + end + E(triu(true(n))) = full(A); % E is place-holder, will be overwritten later + A = reshape(E,n,n); + A = A + triu(A,1)'; + if any(Arad) + X(triu(true(n))) = full(Arad); % L is used as place-holder again + Arad = reshape(X,n,n); + Arad = Arad + triu(Arad,1)'; + if nnz(Arad) up +setround(1); + +% E = mag(I-X'*X) +if isreal(X) + setround(-1); + E = -(X'*X-speye(n)); % -inf(X'*X-I) + setround(1); + % E = max(-inf(X'*X-I),sup(X'*X-I)) >= ||I-X'*X|| + E = max(E,X'*X-speye(n)); +else + realX = real(X); % by using these variables the compiler will + imagX = imag(X); % recognize symmetry -> faster + + % sup(imag(X')*real(X)+real(X')*imag(X)) = sup(imag(X'*X)) + E = imagX'*(-realX) + realX'*imagX; + % abs(imag(X'*X)) <= max(sup(imag(X'*X)),-inf(imag(X'*X))) + % sup(imag(X'*X)) = -inf(imag(X'*X))' + E = max(E,E'); % overwrite E for memory reasons, here E >= abs(imag(X'*X)) + + % E >= abs(abs(real(X'*X)+j*abs(imag(X'*X))) + % abs(real(X'*X)) <= max(sup(real(X'*X)),-inf(real(X'*X))) + setround(-1); + INF = realX'*realX + imagX'*imagX - speye(size(X,2)); + setround(1); + SUP = realX'*realX + imagX'*imagX - speye(size(X,2)); + E = abs( max(SUP,-INF) + 1i*E); + + clear INF SUP realX imagX; + + % significantly easier code with Intlab, but slower + more dependencies + % E = mag(eye(n)-intval(X)*X'); +end + +% kappa <= 1 - ||X'*X-I|| <= omin(X)^2 +kappa = -(sqrtsup(norm22pos(E))-1); % avoid switch of rounding mode + +% warning if A may not be normal +if kappa= abs(abs(real(A*X-X*D)) + j*abs(imag(A*X-X*D))) >= abs(A*X-X*D) + R = abs( max(real(SUP),real(INF)) + 1i * max(imag(SUP),imag(INF)) ); +end + +if ~isempty(find(Arad,1)) % if interval radius~=0 + R = R + Arad*abs(X); +end + +% free some memory +clear INF SUP Arad X + + +% r >= ||A*X-X*D|| / omin(X) % +if hermit + r = sqrtsup( norm22pos(R) / kappa ); +else % else frobenius-norm + r = sqrtsup( (R(:)'*R(:)) / kappa ); +end + +drad = r*e; % radius of eigenvalue diagonal vector + + +%% using generalized perturbation theorem for tighter inclusion %% + +% split into clusters of eigenvalues and update residual norms +if hermit % Cao: sharp version of Kahan's perturbation theorem + + [d,index] = sort(real(d),1,'ascend'); + + k = 1; % index of first eigenvalue in current cluster + + for i = 1:n + if i>=n || d(i+1)>d(i)+r % if end of cluster + % upper bound for matrix norm + kappa = -( sqrtsup( norm22pos(E(index(k:i),index(k:i))) ) - 1 ); + drad(k:i) = sqrtsup( norm22pos(R(:,index(k:i))) / kappa ); + k = i+1; % set new starting position for next cluster + end + end + +else % matrix is normal: use frobenius-norm for inclusion, see Kahan + + % build cluster graph -> abs(e*d'-d*e')=l % last position of current cluster ? + + kappa = -( sqrtsup(norm22pos(E(index(k:l),index(k:l)))) - 1 ); + Rp = reshape(R(:,index(k:l)),[],1); % get residuum for current cluster + drad(index(k:l)) = sqrtsup((Rp'*Rp) / kappa); + k = i + 1; % set position of next cluster + l = k; + + end + + end + +end % is normal + + +% update interval vector for lambda inclusion +if isintval + lambda = midrad(d,drad); +else + lambda = struct('mid',d,'rad',drad); +end + + +% reset rounding mode +setround(rnd); + +%________________________________END OF VSDPNEIG____________________________ + + + +%%% modified normposA from INTLAB by Siegfried M. Rump %%% + +function normbnd = norm22pos(A) +% upper bound for squared spectral norm of nonnegative matrix A +% applying collatz theorem to power iteration +% +% note: this is a helper function used by vsdpneig, +% parameters are not checked + +% setround(1) % here not set, already done in calling function + +x = full(sum(A,2)'*A)'; +normbnd = max(x); +x = x / normbnd; + +while normbnd + oldnorm = normbnd; + y = ((A*x)'*A)'; + normbnd = max(y./x); + if normbnd*1.005>=oldnorm + break; + end + x = y / normbnd; +end + +%__________________________End of NORM22POS___________________________ + + +% the following is a small function to calculate a verified upper bound +% for the square root -> setround(1) is assumed +% only used in case the sqrt function does not regard the rounding mode +function res = sqrtsup(a) +res = sqrt(a); +res = res + realmin * (res*(-res) > -a); \ No newline at end of file diff --git a/vsdpup.m b/vsdpup.m new file mode 100644 index 0000000..fc7ad5d --- /dev/null +++ b/vsdpup.m @@ -0,0 +1,312 @@ +function [fU x lb info] = vsdpup(A,b,c,K,x,y0,z0,yu,opts) +%% VSDPLOW - verified upper bound for semidefinite-quadratic-linear programming +% [fU x lb info] = vsdpup(A,b,c,K,x0,y0,z0) +% +%% >> Description: +% computes verified upper bound of primal optimal value and rigorous +% enclosure of primal strict feasible (near optimal) solution of a conic +% problem in the standard primal-dual form: +% +% (P) min c'*x (D) max b'*y +% s.t. A*x = b s.t. z := c - A'*y +% x in K z in K* +% +% where K is a cartesian product of the cones R+, SOCP, PSD. +% For theoretical introduction into verified conic programming see: +% C. Jansson. On Verified Numerical Computations in Convex Programming. +% Japan J. Indust. Appl. Math., 26:337–363, 2009 +% +%% >> Input: +% A: nA x m coefficient matrix in SeDuMi or VSDP internal format +% b: a M x 1 vector +% c: a nA x 1 vector, primal objective +% K: a structure with following fields +% - K.f stores the number of free variables +% - K.l is the number of nonnegative components +% - K.q lists the lengths of socp blocks +% - K.s lists the dimensions of semidefinite blocks +% x0: a nA x 1 vector - a primal feasible (eps-optimal) solution +% y0: a M x 1 vector - a dual feasible (eps-optimal) solution +% z0: a nA x 1 vector - a dual feasible (eps-optimal) solution (slack vars) +% yu: upperbounds of the absolute value of dual optimal solution y +% opts: structure for additional parameter settings: +% fields: +% 'ITER_MAX' maximum number of iterations that can be used to +% perturbate the problem and find a feasible solution +% 'ALPHA' growing factor for problem perturbation +% -> default: 0.5 +% 'FULL_EIGS_ENCLOSURE' if true code for stronger complete +% eigenvalue enclosure will be applied +% the code is a little bit slower +% -> default: true +% 'SOLVER' to select one of the supported solvers: +% 'sedumi','sdpt3','sdpa','csdp','sdplr', 'lp_solve','linprog' +% 'USE_STARTING_POINT' decides whether initial starting point shall +% be used or not +% -> default: false +% +%% >> Output: +% fU: verified upper bound of the primal optimal value +% x: a nA x 1 vector - rigorous enclosure of primal strict feasible solution +% lb: verified lower bounds of eigenvalues or spectral values of X with +% respect to K +% info.iter: number of iterations +% + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + +%% Last modified: +% 31/07/10 V. Haerter, comments added +% 08/18/12 M. Lange, complete rewrite +% 18/08/12 M. Lange, preconditioning to find better basis indices +% +%% +% TODO: rotated quadratic cones +% + + +%% input parameter + +% check number of input arguments +if nargin<5 || isempty(A) || isempty(b) || isempty(c) || isempty(K) + error('VSDP:VSDPUP','more input arguments are required'); +elseif nargin<6 + y0 = []; z0 = []; yu = []; opts = []; +elseif nargin<7 + z0 = []; yu = []; opts = []; +elseif nargin<8 + yu = []; opts = []; +elseif nargin<9 + opts = []; +end + +global VSDP_OPTIONS; % global options structure + +% use starting point - default: false +if ~isfield(VSDP_OPTIONS,'USE_STARTING_POINT') + VSDP_OPTIONS.USE_STARTING_POINT = false; +end +if ~isfield(opts,'USE_STARTING_POINT') + opts.USE_STARTING_POINT = VSDP_OPTIONS.USE_STARTING_POINT; +end + +% parameter list +optLIST = {'ITER_MAX','ALPHA','FULL_EIGS_ENCLOSURE'}; +optfs = isfield(opts,optLIST); +goptfs = isfield(VSDP_OPTIONS,optLIST); + +% read parameters +[ITER_MAX, ALPHA, FULL_EIGS_ENCLOSURE] = deal(10, 0.5, true); +for i = 1:length(optLIST) + if optfs(i) + eval([optLIST{i},' = opts.',optLIST{i},';']); + elseif goptfs(i) + eval([optLIST{i},' = VSDP_OPTIONS.',optLIST{i},';']); + end +end + + +%% Preliminary steps / Prealocations + +% default output +fU = Inf; +lb = NaN; +info.iter = 0; + +% rounding mode +rnd = getround(); +setround(0); + +% import data +[A,Arad,b,brad,c,crad,K,x,y0,z0,IF] = import_vsdp(A,b,c,K,x,y0,z0); + +% check if approximations are applicable +if isempty(x) || any(isnan(x)) + warning('VSDP:VSDPUP','not applicable approximations given (NaN)'); + x = NaN; + return; +elseif any(isnan(y0)) || any(isnan(z0)) + y0 = []; z0 = []; +end + +% get problem data dimensions +m = length(b); % number of linear system constraints +dim3 = length(x); +nc = K.l + length(K.q) + length(K.s); % number of cone constraints + +% check upper dual bound +yu = yu(:); +if isempty(yu) + yu = inf(m,1); +elseif length(yu)~=m + error('VSDP:VSDPUP','upper bound vector has wrong dimension'); +end + + +%% Algorithm with finite dual bounds yu, projection into cone +if max(yu) K: now regard defect = |A*x-b| + defect = resmag(x',A,b',1,0,Arad,brad',0); + % fU = sup(x'*c + defect*yu) + fU = prodsup(x',c,0,crad) + defect*yu; + x = NaN; lb = NaN; + % reset rounding mode and return + setround(rnd); + return; +end + +%% Algorithm with infinite dual bounds yu +% variable declarations and allocations +x0 = x; % starting point +lb = -inf(nc,1); % dual lower bounds +epsj = ones(nc,1); % factor for perturbation amount +xeps = sparse(length(x),1); % wanted perturbation for x +pertS = cell(length(K.s),1); % for diagonal perturbations of sdp blocks +I = []; % basic indices + +% create index vector for perturbation entries +pertI = ones(sum(K.s),1); +pertI(cumsum(K.s(1:end-1))+1) = 1-K.s(1:end-1); +pertI = [ones(K.l+(~isempty(K.q)),1); K.q(1:end-1); cumsum(pertI)]; +pertI(1) = K.f + 1; +pertI = cumsum(pertI); + + +%% **** main loop **** +while info.iter<=ITER_MAX + info.iter = info.iter + 1; + setround(1); % default rounding for verification part + + %% 1.step: compute rigorous enclosure for x + [x, I] = vuls([],[],struct('mid',A,'rad',Arad),... + struct('mid',b,'rad',brad),[],[],x,I,opts); + if ~isstruct(x) + disp('VSDPUP: could not find solution of primal equations'); + break; + else + xrad = sparse(x.rad); + x = full(x.mid); + end + + %% 2.step: verified lower bounds on cone eigenvalues + if K.l>0 % bound for linear variables + ind = K.f+1:K.f+K.l; + xjn = xrad(ind) - x(ind); % -inf(xl) + % interprete all lp vars as one cone: + % if any element is negative all constraints close to zero will be + % perturbated + xj1 = max(xjn); + if xj1>0 % there is a negative bound + ind = xjn > -0.5*xj1; + xjn(ind) = xjn(ind) + xj1; + end + lb(1:K.l) = - xjn; + end + % eigenvalue bound for second-order cone variables + ind = K.f + K.l; + for j = 1:length(K.q) + ind = ind(end)+2:ind(end)+K.q(j); + xj1 = xrad(ind(1)-1) - x(ind(1)-1); % -inf(xq(1)) + xjn = abs(x(ind)) + xrad(ind); % sup(abs(xq(2:end))) + xjn = sqrtsup(xjn'*xjn); % sup(||xq(2:end)||) + lb(K.l+j) = -(xj1+xjn); % inf(xq(1)-||xq(2:end)||) + end + % eigenvalue bound for semidefinite cone + blke = K.f + K.l + sum(K.q) + sum(K.s.*(K.s+1))/2; + for j = length(K.s):-1:1 + ofs = K.l + length(K.q) + j; + dind = cumsum(1:K.s(j)); % index for diagonal entries + blk3 = dind(end); + vx = struct('mid',0.5*x(blke-blk3+1:blke),... + 'rad',0.5*xrad(blke-blk3+1:blke)); + vx.mid(dind) = 2 * vx.mid(dind); % regard mu=2 for x + vx.rad = vx.rad + vx.rad.*sparse(dind,1,1,blk3,1); + [lmin,lb(ofs),pertS{j}] = bnd4sd(vx,1,FULL_EIGS_ENCLOSURE); + if lmin>0 + lb(ofs) = lmin; + end + pertS{j} = epsj(ofs) * pertS{j}(:); + blke = blke - blk3; + end + + %% 3.step: cone feasibility check, computing upper bound + lbi = lb<0; + if ~any(lbi) + fU = prodsup(x',c,xrad',crad); % sup(x'*c) + x = export_vsdp(IF,K,midrad(full(x),full(xrad))); + setround(rnd); + return; + end + + %% 4.step: create some perturbed problem + setround(0); % no code for rigorous computations + ind = 1:K.l+length(K.q); + if isempty(ind) + xeps = xeps + sparse(pertI,1,cat(1,pertS{:}),dim3,1); + else + xeps = xeps + sparse(pertI,1,cat(1,epsj(ind).*min(lb(ind),0),... + pertS{:}),dim3,1); + end + if any(isinf(xeps)) || any(isnan(xeps)) + disp('VSDPUP: perturbation extended range'); + break; + end + epsj(lbi) = epsj(lbi) * (1+ALPHA); % update perturbation factor + + %% 5.step: solve the perturbed problem + clear lbi ind vx x xrad; % free some memory before calling solver + [objt,x0,y0,z0,INFO] = mysdps(A,b+(xeps'*A)',c,K,x0,y0,z0,opts); + % if could not find solution or primal infeasible: break + if isempty(x0) || any(isnan(x0) | isinf(x0)) || any(INFO(1)==[1 (2) 3]) + disp('VSDPUP: could not find solution for perturbed problem'); + break; + end + x = x0 - xeps; % undo perturbation +end +% **** end of main loop **** +setround(rnd); + +% write output +if info.iter==ITER_MAX + disp('VSDPUP: maximum number of iterations reached'); +end +fU = Inf; x = NaN; lb = NaN; + +%________________________________End VSDPUP_______________________________ \ No newline at end of file diff --git a/vuls.m b/vuls.m new file mode 100644 index 0000000..6b0ff23 --- /dev/null +++ b/vuls.m @@ -0,0 +1,228 @@ +function [X, I, J] = vuls(A,a,B,b,xl,xu,x0,I,opts) +%% VULS: Verification for Underdetermined Linear Systems +% of inequalities and equations: +% A * x <= a, +% B * x = b, or if B is transposed: x' * B = b' +% xl <= x <= xu +% The input A (m*n matrix), a (m-vector), B (p*n or n*p matrix), and b (p-vector) +% can be real or interval quantities; the simple bounds xl, xu must +% be real, but may be infinite; the approximate solution x0 must be real. +% The optional input vector I must contain p indices out of {1,...,n} such +% that the submatrix B(:,I) is nonsingular; if opts.VERIFY_FULL_LSS is true +% the full non-symmetric lss enclosure algorithm will be applied. +% +%% Output: +% X a box (n-interval vector), containing for every real +% input (A,a,B,b) within the interval input data a solution +% x of the above system, provided J is empty. Especially, +% existence of solutions is verified, and moreover +% X is computed close to x0 in a specified manner; for details see +% C. Jansson, Rigorous Lower and Upper Bounds in Linear +% Programming, SIAM J. OPTIM. Vol.14, No.3, pp. 914-935. +% +% I index vector such that the p*p submatrix B(:,I) (or in the +% transposed case: B(I,:)) is nonsingular. +% +% X := nan(n,1), I = [], J = []; +% if existence of solutions cannot be proved, and verified finite +% bounds cannot be computed. This is the case, if +% (i) B has no full rank, or (ii) the linear interval solver +% VERIFYLSS cannot compute rigorous bounds, or (iii) the box +% of simple bounds [xl,xu] has no appropriate interior. Otherwise, +% J returns the vector of indices of violated inequalities +% for X: +% J.ineqlin: violated row indices of A * X <= b, +% J.lower: violated indices of xl <= X , +% J.upper: violated indices of X <= xu . +% +% Used Files: verifylss (INTLAB) +% +% written 01/12/05 Christian Jansson +% modified 03/21/06 +% modified 27/07/10 Vikor Haerter +% modified 10/08/12 Marko Lange, parameters + pre-conditioning of I +% modified 17/09/12 Marko Lange, - removed pre-conditioning +% - added use of spdotK +% - stronger pivoting when using UMFPACK +% +% Reference: C. JANSSON, Termination and Verification for +% Ill-posed Semidefinite Programming Problems, +% to appear +% http://optimization-online.org/DBHTML/|2005/06/1150.html. + +%% ********************************************************************* %% +%% This file is part of VSDP by V. Haerter, C. Jansson and M. Lange %% +%% Copyright (c) 2012, C. Jansson %% +%% Technical University of Hamburg (TUHH) %% +%% Institute for Reliable Computing (IRC) %% +%% VSDP can be freely used for private and academic purposes. %% +%% Commercial use or use in conjunction with a commercial program which %% +%% requires VSDP or part of it to function properly is prohibited. %% +%% ********************************************************************* %% + + +%% check parameter +if nargin<7 || isempty(x0) || ~isnumeric(x0) + error('VULS: not enough input parameter or no initial solution set'); +end +x0 = x0(:); +if ~isnumeric(a) + error('VULS: inappropriate a'); +end +a = a(:); +isIntval = 1; +if isstruct(A) || isstruct(B) || isstruct(b) + isIntval = 0; +end +% read interval b +if isnumeric(b) + b = b(:); brad = sparse(size(b,1),1); +elseif isa(b,'intval') || all(isfield(b,{'mid','rad'})) + brad = b.rad(:); b = b.mid(:); +else + error('VULS: unexpected data format of b'); +end +% dimensions +m = length(a); +n = length(x0); +p = length(b); +% read interval A +if isnumeric(A) && ~isempty(A) + Arad = sparse(size(A,1),size(A,2)); +elseif isa(A,'intval') || all(isfield(A,{'mid','rad'})) + Arad = A.rad; A = A.mid; +elseif ~isempty(A) + error('VULS: unexpected data format of A'); +end +if all(size(A)==[m n]) + A = A'; Arad = Arad'; +elseif ~isempty(A) && (any(size(A)~=[n m]) || any(size(Arad)~=[n m])) + error('VULS: dimension of A does not match to a and x0'); +end +% read interval B +if isnumeric(B) && ~isempty(B) + Brad = sparse(max(size(B)),min(size(B))); +elseif isa(B,'intval') || all(isfield(B,{'mid','rad'})) + Brad = B.rad; B = B.mid; +else + error('VULS: cannot import data format of B'); +end +if all(size(B)==[p n]) + B = B'; +end +if all(size(Brad)==[p n]) + Brad = Brad'; +elseif ~isempty(B) && (any(size(B)~=[n p]) || any(size(Brad)~=[n p])) + error('VULS: dimension of B does not match to b and x0'); +end +% read bounds for x0 +xl = xl(:); +if length(xl)~=n + xl = -inf(n,1); +end +xu = xu(:); +if length(xu)~=n + xu = inf(n,1); +end +% check index and opts parameter +if nargin<8 || length(I)~=p + I = []; +end +if nargin<9 + opts = []; +end + + +%% import options +global VSDP_OPTIONS; + +% full non-symmetric matrix lss verification +full_lss = false; % function default +if isfield(opts,'VERIFY_FULL_LSS') + full_lss = opts.VERIFY_FULL_LSS; +elseif isfield(VSDP_OPTIONS,'VERIFY_FULL_LSS') + full_lss = VSDP_OPTIONS.VERIFY_FULL_LSS; +end + + +%% preparation +% initial output +X = nan(n,1); +J.ineqlin = []; +J.lower = []; +J.upper = []; + +% projection of x0 into the interior of [xl,xu] +if any(xu