Merge pull request #4 from OxfordControl/developer

Merge Developer into Master for release v1.1

Author: giofantuzzi

.gitignore

@@ -216,6 +216,13 @@ pythontex-files-*/
 TSWLatexianTemp*
 
 ## Editors:
+
+# MATLAB editor
+*.asv
+
+# gedit
+*.*~
+
 # WinEdt
 *.bak
 *.sav
@@ -228,5 +235,3 @@ TSWLatexianTemp*
 
 # KBibTeX
 *~[0-9]*
-admmdual/dataChorDec.m
-admmdual/dataChorDec.m

README.md

@@ -1,11 +1,15 @@
 # CDCS
 
 CDCS (Cone Decomposition Conic Solver) is an open-source MATLAB solver for sparse conic programs with partially decomposable conic constraints. CDCS implements the alternating direction method of multipliers (ADMM) 
-described in on our paper [_Fast ADMM for Semidefinite Programs with Chordal Sparsity_](https://arxiv.org/pdf/1609.06068v2.pdf) (included in the `doc/` folder).
+described in our papers [_Fast ADMM for Semidefinite Programs with Chordal Sparsity_](https://arxiv.org/pdf/1609.06068v2.pdf) and [_Fast ADMM for homogeneous self-dual embeddings of sparse SDPs_](https://arxiv.org/pdf/1611.01828.pdf) (included in the `doc/` folder)
 
-**Current version:** 1.0
+**Current version:** 1.1
 
-**Release notes:** CDCS is based on a temporary research code called ADMM-PDCP, which is no longer maintained. If you downloaded ADMM-PDCP, please replace it with CDCS.
+**Release notes:** 
+
+* Homogeneous self-dual embedding is the new default method 
+* The termination codes have changed. This means that if you use CDCS from YALMIP, the termination code returned by YALMIP will be incorrect. This should be fixed in the next YALMIP release!
+* CDCS is based on a temporary research code called ADMM-PDCP, which is no longer maintained. If you downloaded ADMM-PDCP, please replace it with CDCS.
 
 
 ## Contents
@@ -21,7 +25,7 @@ described in on our paper [_Fast ADMM for Semidefinite Programs with Chordal Spa
 CDCS solves in the standard primal and dual vectorized forms
 
 		minimize 	c'x						maximize 	b'y
-	(1)	subject to	Ax = b,				(2)	subject to	c - A'y = z,	
+	(1)	subject to	Ax = b,				(2)	subject to	A'y + z = c,	
 					x \in K								z \in K*
 
 where the conic constraint `x \in K` are partially decomposable. This means that
@@ -41,6 +45,9 @@ sparsity pattern. The other supported cone types are not decomposed.
 This means that CDCS is most suitable for large sparse semidefinite programs (SDPs),
 although it can be used for any conic program over the supported cones.
 
+CDCS offers a choice to solve the primal problem (1) only, the dual problem (2) only,
+or the homogeneous self-dual embedding of the two problems. _From version 1.1, the 
+homogeneous self-dual embedding is the default method._
 
 
 ## Quick start<a name="QuickStart"></a>
@@ -72,26 +79,36 @@ you have any suggestions for improvement, or find any bugs, feel free to [contac
 
 ## How to cite<a name="References"></a>
 
-If you find CDCS useful, please cite:
+If you find CDCS useful, please cite at least one of the following papers as appropriate:
 
 ```
-@article{ZFPGW2016,
-	archivePrefix = {arXiv},
-	eprint = {1609.06068v2},
-	primaryClass = "math-OC",
-	author = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
-	title = {{Fast ADMM for Semidefinite Programs with Chordal Sparsity}}
-	}
+@article{{ZFPGWhsde2016,
+    archivePrefix= {arXiv},
+    eprint       = {1611.01828},
+    primaryClass = "math-OC",
+    author       = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
+    title        = {{Fast ADMM for homogeneous self-dual embeddings of sparse SDPs}}
+    }
+
+@article{ZFPGWpd2016,
+    archivePrefix= {arXiv},
+    eprint       = {1609.06068v2},
+    primaryClass = "math-OC",
+    author       = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
+    title        = {{Fast ADMM for Semidefinite Programs with Chordal Sparsity}}
+    }
 	
 @misc{CDCS,
     author       = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
-    title        = {{CDCS}: Cone Decomposition Conic Solver, version 1.0},
+    title        = {{CDCS}: Cone Decomposition Conic Solver, version 1.1},
     howpublished = {\url{https://github.com/giofantuzzi/CDCS}},
     month        = Sep,
     year         = 2016
     }
 ```
 A selection of BibTex styles that support arXiv preprints can be found [here](http://arxiv.org/hypertex/bibstyles/).
+
+
 ## Contact us<a name="Contacts"></a>
 To contact us about CDCS, suggest improvements and report bugs, email either [Giovanni Fantuzzi](mailto:gf910@ic.ac.uk?Subject=CDCS) or [Yang Zheng](mailto:yang.zheng@eng.ox.ac.uk?Subject=CDCS).
 

README.txt

@@ -4,11 +4,15 @@
 
 CDCS (Cone Decomposition Conic Solver) is an open-source MATLAB solver for sparse
 conic programs with partially decomposable conic constraints. CDCS implements the
-alternating direction method of multipliers (ADMM) described in on our paper: 
-"Fast ADMM for Semidefinite Programs with Chordal Sparsity" available from
-(https://arxiv.org/pdf/1609.06068v2.pdf).
+alternating direction method of multipliers (ADMM) described in on our papers
+ 
+* "Fast ADMM for Semidefinite Programs with Chordal Sparsity" available from
+(https://arxiv.org/pdf/1609.06068v2.pdf)
 
-Current version: 1.0
+* "Fast ADMM for homogeneous self-dual embeddings of sparse SDPs" available from
+(https://arxiv.org/abs/1611.01828)
+
+Current version: 1.1
 Release notes: CDCS is based on a temporary research code called ADMM-PDCP, which
                is no longer maintained. If you downloaded ADMM-PDCP, please 
                replace it with CDCS.
@@ -32,7 +36,7 @@ Release notes: CDCS is based on a temporary research code called ADMM-PDCP, whic
 CDCS solves in the standard primal and dual vectorized forms
 
 		minimize 	c'x						maximize 	b'y
-	(1)	subject to	Ax = b,				(2)	subject to	c - A'y = z,	
+	(1)	subject to	Ax = b,				(2)	subject to	A'y + z = c,	
 					x \in K								z \in K*
 
 where the conic constraint `x \in K` are partially decomposable. This means that
@@ -87,9 +91,17 @@ us (see the Contact Us section below).
                                 HOW TO CITE
 ================================================================================
 
-If you find CDCS useful, please cite it with
+If you find CDCS useful, please cite
+
+@article{{ZFPGWhsde2016,
+    archivePrefix = {arXiv},
+	eprint = {1611.01828},
+	primaryClass = "math-OC",
+	author = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
+	title = {{Fast ADMM for homogeneous self-dual embeddings of sparse SDPs}}
+	}
 
-@article{ZFPGW2016,
+@article{ZFPGWpd2016,
 	archivePrefix = {arXiv},
 	eprint = {1609.06068v2},
 	primaryClass = "math-OC",

cdcs.m

@@ -1,44 +1,50 @@
 function [x,y,z,info] = cdcs(At,b,c,K,userOpts,initVars)
-
 % CDCS
 %
 % Syntax:
 %
-% [x,y,z,info] = CDCS(At,b,c,K,opts)
+% [x,y,z,info] = CDCS(At,b,c,K,options)
 %
-% Solve a sparse conic program using chordal decomposition for the semidefinite
+% Solve a sparse conic program using chordal decomposition for the positive semidefinite
 % cones and ADMM. CDCS solves the primal (P) or dual (D) standard forms of
 % the conic problem,
 %
 %         min <c,x>                             max <b,y>
-%   (P)   s.t. Ax = b,                  (D)     s.t. c - A^Ty = z
+%   (P)   s.t. Ax = b,                  (D)     s.t. A^Ty + z = c
 %         x \in K                               z \in K*
 %
-% where A,b and c are the problem date and K is the cone (K* is the dual cone).
+% where A,b and c are the problem data and K is the cone (K* is the dual cone).
+% CDCS supports the following cones: Free, Linear, second-order,
+% Semi-definite, called as called K.f, K.l, K.q, and K.s.
+%
 % The standard form to be solved is specified by the "solver" field of the
 % options structure:
 %
-%   opts.solver = 'primal' (default): solve the problem in primal standard form
-%   opts.solver = 'dual'            : solve the problem in dual standard form
+%   options.solver = 'hsde' (default): solve the problem in homogeneous self-dual embedding form
+%   options.solver = 'primal'        : solve the problem in primal standard form
+%   options.solver = 'dual'          : solve the problem in dual standard form
 %
 % The chordal decomposition can be carried out in two ways, specified by the
 % "chordalize" option:
 %
-%   opts.chordalize = 1 (default): split the data equally between the cliques
-%   opts.chordalize = 2          : assign data to one clique only
+%   options.chordalize = 1 (default): split the data equally between the cliques
+%   options.chordalize = 2          : assign data to one clique only
 %
 % <a href="matlab:help('cdcsOpts')">Click here for a complete list of options</a>.
 %
 % The output structure 'info' contains the following information:
 %
 % info.problem: - 0: CDCS terminated succesfully 
-%               - 1: the maximum number of iterations was reached
-%               - 2: the ADMM iterations terminated succesfully, but the positive 
-%                 matrix completion algorithm threw an error
+%               - 1: primal infeasibility detected
+%               - 2: dual infeasibility detected
+%               - 3: maximum number of iterations reached
+%               - 4: the ADMM iterations terminated succesfully, but the positive 
+%                    matrix completion algorithm threw an error
 % info.iter: number of iterations
 % info.cost: terminal cost
 % info.pres: terminal primal ADMM residual
 % info.dres: terminal dual ADMM residual
+% info.log : history log of the ADMM residuals, cost, etc.
 % info.time: some timing information (setup, ADMM iterations, cleanup, total)
 %
 % See also CDCSOPTS
@@ -61,10 +67,11 @@
 
 
 %============================================
-% Solver options
+% Solver options & import cdcs_utils
 %============================================
 tstart = tic;
 opts = cdcsOpts;
+import cdcs_utils.*
 
 
 %============================================
@@ -76,13 +83,13 @@
 end
 
 % Checks on specified solver type and method
-if ~any(strcmpi(opts.solver,{'primal','dual'}))
-    error('Unknown opts.solver. Please use "primal" or "dual".')
+if ~any(strcmpi(opts.solver,{'primal','dual','hsde'}))
+    error('Unknown opts.solver. Please use "primal", "dual" or "hsde".')
 end
 
 % Print nice welcoming header
 if opts.verbose
-    myline1 = [repmat('=',1,64),'\n'];
+    [header,myline1,myline2] = printHeader(opts);
     fprintf(myline1)
     fprintf('CDCS by G. Fantuzzi, Y. Zheng -- v1.0\n')
     fprintf(myline1)
@@ -93,15 +100,14 @@
 % start timing
 proctime = tic;
 
-% sparsify everything, check cone constraints, rescale
+% sparsify everything, check cone constraints
 [At,b,c,K,opts] = checkInputs(At,b,c,K,opts);
-[At,b,c,K,opts] = rescaleData(At,b,c,K,opts);
 [At,b,c,K,opts] = splitBlocks(At,b,c,K,opts);
 [opts.n,opts.m] = size(At);
 
-% chordal decomposition
+% rescale & chordal decomposition
 Kold = K;
-[At,b,c,K,Ech,cd,chstuff] = chordalize(At,b,c,K,opts);
+[At,b,c,K,Ech,chstuff,opts] = preprocess(At,b,c,K,opts);
 
 % basic decomposed problem dimensions:  no. of cones, no. of vectorized conic
 % variables, and no. of free primal variables
@@ -114,80 +120,77 @@
 [X,Y,Z,others] = makeVariables(K,initVars,opts);
 
 % Make operators for ADMM
-[step1,step2,step3,checkConv] = makeADMM(At,b,c,K,cd,Ech,opts);
+[updateX,updateY,updateZ,checkConvergence] = makeADMM(At,b,c,K,Ech,opts);
 
 % Time setup and display
 proctime = toc(proctime);
 if opts.verbose
-    myline2 = [repmat('-',1,64),'\n'];
+    % Set method to display
+    if strcmpi(opts.solver,'hsde')
+        method = 'homogeneous self-dual embedding';
+    else
+        method = opts.solver;
+    end
     fprintf('done in %.4f seconds.      \n',proctime);
-    fprintf('Standard form          : %s\n',opts.solver);
+    fprintf('Algorithm              : %s\n',method);
     fprintf('Chordalization method  : %i\n',opts.chordalize);
     fprintf('Adaptive penalty       : %i\n',opts.adaptive);
     fprintf('Scale data             : %i\n',opts.rescale);
     fprintf('Free variables         : %i                \n',K.f);
     fprintf('Non-negative variables : %i                \n',K.l);
-    fprintf('Second-order cones     : %i (max. size: %i)\n',length(K.q),max(K.q));
-    fprintf('Semidefinite cones     : %i (max. size: %i)\n',length(K.s),max(K.s));
+    fprintf('Second-order cones     : %i (max. size: %i)\n',length(find(K.q ~=0)),max(K.q));
+    fprintf('Semidefinite cones     : %i (max. size: %i)\n',length(find(K.s ~=0)),max(K.s));
     fprintf('Affine constraints     : %i                \n',opts.m);
     fprintf('Consensus constraints  : %i                \n',sum(accumarray(Ech,1)));
-    fprintf(myline1)
+    fprintf(myline1);
+    fprintf(header);
+    fprintf(myline2);
 end
 
 %============================================
 % Run ADMM
 %============================================
-% Display
-if opts.verbose
-    fprintf(' iter |   pres   |   dres   |    cost    |   rho    | time (s) |\n')
-    fprintf(myline2)
-end
-
 admmtime = tic;
-opts.feasCode = 1;
 for iter = 1:opts.maxIter
 
     % Save current iterate for convergence test
     YOld = Y;
 
     % Update block variables
-    [X,others] = step1(X,Y,Z,opts.rho,others);
-    [Y,others] = step2(X,Y,Z,opts.rho,others);
-    [Z,others] = step3(X,Y,Z,opts.rho,others);
+    [X,others] = updateX(X,Y,Z,opts.rho,others);
+    [Y,others] = updateY(X,Y,Z,opts.rho,others);
+    [Z,others] = updateZ(X,Y,Z,opts.rho,others);
 
     % log errors / check for convergence
-    [isConverged,log,opts] = checkConv(X,Y,Z,YOld,others,iter,admmtime,opts);
-    if isConverged
-        opts.feasCode = 0;
+    [stop,info,log,opts] = checkConvergence(X,Y,Z,YOld,others,iter,admmtime,opts);
+    if stop
         break;
     end
 end
 admmtime = toc(admmtime);
-if opts.verbose
-    fprintf(myline1)
-end
 
 %============================================
 % Outputs
 %============================================
 % Variables in sedumi format
 posttime = tic;
-[x,y,z,opts] = setOutputs(X,Y,Z,others,Kold,cd,Ech,chstuff,opts);
+[x,y,z,info,opts] = setOutputs(X,Y,Z,others,Kold,c,Ech,chstuff,info,opts);
 posttime = toc(posttime);
 
 % Info
-info.problem = opts.feasCode;              % diagnostic code
 info.iter    = iter;                       % # of iterations
-info.cost    = log.cost;                   % terminal cost
-info.pres    = log.pres;                   % terminal primal ADMM res
-info.dres    = log.dres;                   % terminal dual ADMM res
+info.cost    = log(iter).cost;             % terminal cost
+info.pres    = log(iter).pres;             % terminal primal ADMM res
+info.dres    = log(iter).dres;             % terminal dual ADMM res
+info.log     = log;                        % log of residuals etc
 info.time.setup   = proctime;              % setup time
 info.time.admm    = admmtime;              % ADMM time
 info.time.cleanup = posttime;              % post-processing time
 info.time.total   = toc(tstart);           % total CPU time
 
 % Print summary
 if opts.verbose
+    fprintf(myline1)
     fprintf(' SOLUTION SUMMARY:\n')
     fprintf('------------------\n')
     fprintf(' Termination code     : %11.1d\n',info.problem)

cdcsInstall.m

@@ -12,11 +12,13 @@
 cd('include')
 cs_install
 cd(here);
-movefile(['include',filesep,'cs_lsolve.mex*'],[here,filesep,'private']);
-movefile(['include',filesep,'cs_ltsolve.mex*'],[here,filesep,'private']);
+movefile(['include',filesep,'cs_lsolve.mex*'], ...
+         [here,filesep,'packages',filesep,'+cdcs_utils']);
+movefile(['include',filesep,'cs_ltsolve.mex*'], ...
+         [here,filesep,'packages',filesep,'+cdcs_utils']);
 
 % Then compile some mex files from this package
-cd('private')
+cd(['packages',filesep,'+cdcs_utils',filesep,'private'])
 if (~isempty (strfind (computer, '64')))
     mexcmd = 'mex -largeArrayDims' ;
 else
@@ -27,6 +29,7 @@
 cd(here)
 
 % Finally add to path and save
+addpath([here,filesep,'packages']);
 addpath(here);
 savepath