Domain-Driven Solver (DDS) Version 2.0: a MATLAB-based Software Package for Convex Optimization Problems in Domain-Driven Form

TL;DR

Domain-Driven Solver 2.0 is a MATLAB software package that efficiently solves a wide range of convex optimization problems using an infeasible-start primal-dual algorithm, with extensive examples and implementation insights.

Contribution

This paper introduces DDS 2.0, expanding the types of convex problems it can handle and providing detailed implementation and efficiency improvements.

Findings

Supports diverse convex constraints including symmetric cones and quantum entropy

Provides practical implementation details and efficiency techniques

Includes numerous examples demonstrating solver capabilities

Abstract

Domain-Driven Solver (DDS) is a MATLAB-based software package for convex optimization problems in Domain-Driven form [Karimi and Tun\c{c}el, arXiv:1804.06925]. The current version of DDS accepts every combination of the following function/set constraints: (1) symmetric cones (LP, SOCP, and SDP); (2) quadratic constraints that are SOCP representable; (3) direct sums of an arbitrary collection of 2-dimensional convex sets defined as the epigraphs of univariate convex functions (including as special cases geometric programming and entropy programming); (4) generalized power cone; (5) epigraphs of matrix norms (including as a special case minimization of nuclear norm over a linear subspace); (6) vector relative entropy; (7) epigraphs of quantum entropy and quantum relative entropy; and (8) constraints involving hyperbolic polynomials. DDS is a practical implementation of the infeasible-start primal-dual algorithm designed and analyzed in [Karimi and Tun\c{c}el, arXiv:1804.06925]. This manuscript contains the users' guide, as well as theoretical results needed for the implementation of the algorithms. To help the users, we included many examples. We also discussed some implementation details and techniques we used to improve the efficiency and further expansion of the software to cover the emerging classes of convex optimization problems.

Tables

Table 1. Table 1. Function/set constraints the current version of DDS accepts, and their abbreviations.

function/set constraint	abbreviation
LP	LP
SOCP	SOCP
Rotated SOCP	SOCPR
SDP	SDP
Generalized Power Cone	GPC
Quadratic Constraints	QC
Epigraph of a Matrix Norm	MN
Direct sum of 2-dimensional sets
(geometric, entropy, and $p$-norm	TD
programming)
Quantum Entropy	QE
Quantum Relative Entropy	QRE
Relative Entropy	RE
Hyperbolic Polynomials	HB
Equality Constraints	EQ

Table 2. Table 2. Some 2-dimensional convex sets and their s.c. barriers.

	set $(z,t)$	s.c. barrier $\Phi (z,t)$
1	$-\text{ln}\left(z\right)\le t,z>0$	$-\text{ln}\left(t+\text{ln}\left(z\right)\right)-\text{ln}\left(z\right)$
2	$e^z\le t$	$-\text{ln}\left(\text{ln}\left(t\right)-z\right)-\text{ln}\left(t\right)$
3	$z\text{ln}\left(z\right)\le t,z>0$	$-\text{ln}\left(t-z\text{ln}\left(z\right)\right)-\text{ln}\left(z\right)$
4	${\left|z\right|}^p\le t,p\ge 1$	$-\text{ln}\left(t^{\frac{2}{p}}-z^2\right)-2\text{ln}\left(t\right)$
5	$-z^p\le t,z>0,0\le p\le 1$	$-\text{ln}\left(z^p+t\right)-\text{ln}\left(z\right)$
6	$\frac{1}{z}\le t,z>0$	$-\text{ln}\left(zt-1\right)$

Table 3. Table 3. Numerical results for some problem from the Dimacs library for t ​ o ​ l = 10 − 8 𝑡 𝑜 𝑙 superscript 10 8 tol=10^{-8} .

Problem	size of $A$	Type of Constraints	Iterations
nb	$2383\ast 123$	SOCP-LP	42
nb_L1	$3176\ast 915$	SOCP-LP	37
nb_L2	$4195\ast 123$	SOCP-LP	23
nb_L2_bessel	$2641\ast 123$	SOCP-LP	27
filter48_socp	$3284\ast 969$	SDP-SOCP-LP	80
filtinf1	$3395\ast 983$	SDP-SOCP-LP	26
truss5	$3301\ast 208$	SDP	66
truss8	$11914\ast 496$	SDP	76
copo14	$3108\ast 1275$	SDP-LP	24
copo23	$3108\ast 1275$	SDP-LP	32
toruspm3-8-50	$262144\ast 512$	SDP	20
torusg3-8	$262144\ast 512$	SDP	24
sched_50_50_scaled	$4977\ast 2526$	SOCP-LP	81
mater-3	$39448\ast 1439$	SDP-LP	124
cnhil8	$14400\ast 1716$	SDP	31
cnhil10	$48400\ast 5005$	SDP	35
cphil10	$48400\ast 5005$	SDP	12
ros_500	$17944\ast 4988$	SDP	41
sensor_500	$245601\ast 3540$	SDP	65
taha1a	$231672\ast 3002$	SDP	23
taha1a	$231672\ast 3002$	SDP	42
G40mc	$4000000\ast 2000$	SDP	33
1tc.1024	$1048576\ast 7937$	SDP	41
yalsdp	$30000\ast 5051$	SDP	17

Table 4. Table 4. Numerical results for some EXP cone problems from CBLIB.

Problem	size of $A$	size of $A_{EQ}$	Type of Constraints	Iterations
syn30m	$324\ast 121$	$19\ast 121$	EXP-LP	36
syn30h	$528\ast 249$	$161\ast 249$	EXP-LP	51
syn30m02h	$1326\ast 617$	$381\ast 617$	EXP-LP	46
syn05m02m	$177\ast 67$	$11\ast 67$	EXP-LP	36
syn10h	$171\ast 84$	$55\ast 84$	EXP-LP	37
syn15m	$177\ast 67$	$12\ast 67$	EXP-LP	31
syn40h	$715\ast 331$	$213\ast 331$	EXP-LP	48
synthes1	$30\ast 13$	$1\ast 13$	EXP-LP	13
synthes2	$41\ast 16$	$1\ast 12$	EXP-LP	26
synthes3	$71\ast 26$	$3\ast 26$	EXP-LP	20
ex1223a	$59\ast 17$	$1\ast 17$	EXP-LP	11
gp_dave_1	$985\ast 705$	$453\ast 705$	EXP-LP	37
fiac81a	$163\ast 191$	$126\ast 191$	EXP-LP	15
fiac81b	$65\ast 87$	$57\ast 87$	EXP-LP	15
batch	$175\ast 58$	$13\ast 58$	EXP-LP	160
demb761	$93\ast 131$	$90\ast 131$	EXP-LP	24
demb762	$93\ast 131$	$90\ast 131$	EXP-LP	27
demb781	$14\ast 19$	$12\ast 19$	EXP-LP	8
enpro48	$506\ast 167$	$33\ast 167$	EXP-LP	216
isi101	$467\ast 393$	$261\ast 33$	EXP-LP	52 (infeas)
jha88	$866\ast 1131$	$825\ast 1131$	EXP-LP	33
LogExpCR-n20-m400	$2023\ast 2022$	$1200\ast 2022$	EXP-LP	40
LogExpCR-n100-m400	$2103\ast 2102$	$1200\ast 2102$	EXP-LP	46
LogExpCR-n20-m1200	$6023\ast 6022$	$3600\ast 6022$	EXP-LP	44
LogExpCR-n500-m400	$2503\ast 2502$	$1200\ast 2502$	EXP-LP	68
LogExpCR-n500-m1600	$8503\ast 8502$	$4800\ast 8502$	EXP-LP	63

Table 5. Table 5. Results for problems involving Quantum Entropy optimization problems.

Problem	size of $A$	Type of Constraints	Itr/  time(s)	CVXQUAD
QuanEntr-10	$101\ast 21$	QE	11/ 0.6	12/  1.3
QuanEntr-20	$401\ast 41$	QE	13/ 0.8	13/ 2.4
QuanEntr-50	$2501\ast 101$	QE	18/ 1.8	12/ 23.2
QuanEntr-100	$10001\ast 201$	QE	24/ 6.7	array size error
QuanEntr-200	$40001\ast 401$	QE	32/ 53.9	array size error
QuanEntr-LP-10	$111\ast 21$	QE-LP	19/ 1.1	16/ 3.4
QuanEntr-LP-20	$421\ast 41$	QE-LP	23 / 1.6	20/ 9.1
QuanEntr-LP-50	$2551\ast 101$	QE-LP	34/ 2.3	20/ 103.4
QuanEntr-LP-100	$10101\ast 201$	QE-LP	44/ 14.5	array size error
QuanEntr-LP-SOCP-10-infea	$122\ast 21$	QE-LP-SOCP	12/ 0.6 (infea)	16/ 5.1
QuanEntr-LP-SOCP-10	$122\ast 21$	QE-LP-SOCP	26/ 1.0	24/ 2.0
QuanEntr-LP-SOCP-20-infea	$441\ast 41$	QE-LP-SOCP	14/ 1.1 (infea)	14/ 5.8
QuanEntr-LP-SOCP-20	$441\ast 41$	QE-LP-SOCP	33/ 1.8	22/ 4.1
QuanEntr-LP-SOCP-100	$10201\ast 201$	QE-LP-SOCP	51/ 20.8	array size error
QuanEntr-LP-SOCP-200	$40402\ast 401$	QE-LP-SOCP	69/ 196	array size error
QuanEntr-LP-3norm-10-infea	$121\ast 21$	QE-LP-pnorm	16/ 0.6 (infea)	22/ 2.0
QuanEntr-LP-3norm-10	$121\ast 21$	QE-LP-pnorm	21/ 1.2	19/ 2.0
QuanEntr-LP-3norm-20-infea	$441\ast 41$	QE-LP-pnorm	16/ 0.9 (infea)	18/ 6.5
QuanEntr-LP-3norm-20	$441\ast 41$	QE-LP-pnorm	29/ 1.1	20/ 7.3
QuanEntr-LP-3norm-100-infea	$10201\ast 201$	QE-LP-pnorm	25/ 15.3 (infea)	array size error
QuanEntr-LP-3norm-100	$10201\ast 201$	QE-LP-pnorm	66/ 40.5	array size error

Table 6. Table 6. Results for problems involving Quantum Relative Entropy optimization problems.

Problem	size of $A$	Type of Constraints	Itr/  time(s)	CVXQUAD
QuanReEntr-6	$73\ast 13$	QRE	12/ 6.3	12/  18.1
QuanReEntr-10	$201\ast 21$	QRE	12/  21.2	25/  202
QuanReEntr-20	$801\ast 41$	QRE	17/  282	array size error
QuanReEntr-LP-6	$79\ast 13$	QRE-LP	25/ 23.2	21/  19.6
QuanReEntr-LP-6-infea	$79\ast 13$	QRE-LP	30/ 11.8	28/  21.5
QuanReEntr-LP-10	$101\ast 21$	QRE-LP	26/ 49.3	24/  223

Table 7. Table 7. Results for problems involving hyperbolic polynomials.

Problem	size of $A$	var/deg of $p\left(z\right)$	Type of Constraints	Itr/  time(s)
Vamos-1	$8\ast 4$	8/4	HB	7/  1.5
Vamos-LP-1	$12\ast 4$	8/4	HB-LP	8/  2
Vamos-SOCP-1	$17\ast 5$	8/4	HB-LP-SOCP	12/  2
Vamos-Entr-1	$17\ast 5$	8/4	HB-LP-Entropy	9/  1.2
VL-Entr-1	$41\ast 11$	20/4	HB-LP-Entropy	9/  24
VL-Entr-1	$61\ast 16$	30/4	HB-LP-Entropy	12/  211

Table 8. Table 8. Results for problems involving hyperbolic polynomials as product of linear forms and their derivatives.

Problem	size of $A$	var/deg of $p\left(z\right)$	Type of Constraints	Itr
HPLin-LP-1	$55\ast 50$	5/3	HB-LP	8
HPLinDer-LP-1	$55\ast 50$	5/3	HB-LP	7
HPLin-LP-2	$55\ast 50$	5/7	HB-LP	19
HPLinDer-LP-2	$55\ast 50$	5/7	HB-LP	22
HPLin-Entr-1	$111\ast 51$	10/15	HB-LP-Entropy	19
HPLinDer-Entr-1	$111\ast 51$	10/15	HB-LP-Entropy	18
HPLin-pn-1	$111\ast 51$	10/15	HB-LP-pnorm	30
HPLinDer-pn-1	$111\ast 51$	10/15	HB-LP-pnorm	32
Elem-unb-10	$10\ast 5$	10/3	HB	12 (unb)
Elem-LP-10	$15\ast 5$	10/3	HB-LP	9
Elem-LP-10-2	$1010\ast 1000$	10/3	HB-LP	30
Elem-Entr-10-1	$211\ast 101$	10/3	HB-LP-Entropy	17
Elem-Entr-10-2	$221\ast 101$	10/4	HB-LP-Entropy	18

Table 9. Table 9. LF conjugates for the first three s.c. barriers in Table 2 .

	$\Phi (z,t)$	${\Phi }_{\ast }(y,\eta )$
1	$-\text{ln}\left(t+\text{ln}\left(z\right)\right)-\text{ln}\left(z\right)$	$-1+\left(-\eta +1\right)\left[-1+\text{ln}\frac{-\left(-\eta +1\right)}{y}\right]-\text{ln}\left(-\eta \right)$
2	$-\text{ln}\left(\text{ln}\left(t\right)-z\right)-\text{ln}\left(t\right)$	$-1+\left(y+1\right)\left[-1+\text{ln}\frac{-\left(y+1\right)}{\eta }\right]-\text{ln}\left(y\right)$
3	$-\text{ln}\left(t-z\text{ln}\left(z\right)\right)-\text{ln}\left(z\right)$	$-\text{ln}\left(-\eta \right)+\theta \left(1+\frac{y}{\eta }-\text{ln}\left(-\eta \right)\right)-\frac{y}{\eta }+\frac{1}{\theta \left(1+\frac{y}{\eta }-\text{ln}\left(-\eta \right)\right)}-3$

Table 10. Table 10. s.c. barriers and their LF conjugate for rows 4 and 5 of Table 2

	s.c. barrier $\Phi (z,t)$	${\Phi }_{\ast }(y,\eta )$
4	$-\text{ln}\left(t^{\frac{2}{p}}-z^2\right)-2\text{ln}\left(t\right)$	$-\left(\frac{2}{p}+\left(\frac{1}{p}-1\right)yz(y,\eta )\right)-2+2\text{ln}\left(\frac{2\left(\frac{1}{p}+1\right)+\frac{1}{p}yz(y,\eta )}{-\eta }\right)$
		$+\text{ln}\left({\left(\frac{2\left(\frac{1}{p}+1\right)+\frac{1}{p}yz(y,\eta )}{-\eta }\right)}^{\frac{2}{p}}-z^2(y,\eta )\right)$
5	$-\text{ln}\left(z^p+t\right)-\text{ln}\left(z\right)$	$\eta \left(p-1\right)z^p(y,\eta )-2-\text{ln}\left(-\eta \right)+\text{ln}\left(z(y,\eta )\right)$

Equations

$$\underset{x}{\inf }\{⟨c,x⟩:Ax\in D\},$$

$$\underset{y}{\inf }\{{\delta }_{\ast }(y|D):A^{\top }y=-c,y\in D_{\ast }\},$$

$$D_{\ast }:=\{y:⟨y,h⟩\le 0,\forall h\in \text{rec}(D)\},$$

$$⟨c,x⟩+{\delta }_{\ast }(y|D).$$

$$A_L^{i}x+b_L^i\ge 0,i\in \{1,\dots ,\ell \},$$

$$\displaystyle\texttt{A\{k,1\}}=\left[\begin{array}[]{c}A_{L}^{1}\\ \vdots\\ A_{L}^{\ell}\end{array}\right],\ \ \ \texttt{b\{k,1\}}=\left[\begin{array}[]{c}b_{L}^{1}\\ \vdots\\ b_{L}^{\ell}\end{array}\right]$$

$$\text{cons{k,1}=’LP’},\text{cons{k,2}}=[m_L^1,\dots ,m_L^{\ell }].$$

$$\Vert A_S^{i}x+b_S^i\Vert \le (g_S^i{)}^{\top }x+d_S^i,i\in \{1,\dots ,\ell \},$$

$$\displaystyle\texttt{A\{k,1\}}=\left[\begin{array}[]{c}(g_{S}^{1})^{\top}\\ A_{S}^{1}\\ \vdots\\ (g_{S}^{\ell})^{\top}\\ A_{S}^{\ell}\end{array}\right],\ \ \ \texttt{b\{k,1\}}=\left[\begin{array}[]{c}d_{S}^{1}\\ b_{S}^{1}\\ \vdots\\ d_{S}^{\ell}\\ b_{S}^{\ell}\end{array}\right]$$

$$\text{cons{k,1}=’SOCP’},\text{cons{k,2}}=[m_S^1,\dots ,m_S^{\ell }].$$

$$\min$$

$$\displaystyle\texttt{cons\{1,1\}='LP'},\ \ \texttt{cons\{1,2\}=[1]},\ \ \texttt{A\{1,1\}}=\left[\begin{array}[]{cc}2&-1\end{array}\right],\ \ \ \texttt{b\{1,1\}}=\left[\begin{array}[]{c}1\end{array}\right],$$

$$\displaystyle\texttt{cons\{2,1\}='SOCP'},\ \ \texttt{cons\{2,2\}=[2]},\ \ \texttt{A\{2,1\}}=\left[\begin{array}[]{cc}0&0\\ 2&1\\ 1&3\end{array}\right],\ \ \ \texttt{b\{2,2\}}=\left[\begin{array}[]{c}2\\ 3\\ 4\end{array}\right].$$

$$\Phi (z)=-\text{ln}(z),z\in {\mathbb{R}}_{+},{\Phi }_{\ast }(\eta )=-1-\text{ln}(-\eta ),$$

$$\Phi (t,z)=-\text{ln}(t^2-z^{\top }z),{\Phi }_{\ast }(\eta ,w)=-2+\text{ln}(4)-\text{ln}({\eta }^2-w^{\top }w).$$

$$\{(z,t,s)\in {\mathbb{R}}^n\oplus \mathbb{R}\oplus \mathbb{R}:\Vert z{\Vert }^2\le ts,t\ge 0,s\ge 0\},$$

$$\Vert A_S^{i}x+b_S^i{\Vert }_2^2\le ((g_S^i{)}^{\top }x+d_S^i)(({\bar{g}}_S^i{)}^{\top }x+{\bar{d}}_S^i),i\in \{1,\dots ,\ell \},$$

$$(g_S^i{)}^{\top }x+d_S^i\ge 0,({\bar{g}}_S^i{)}^{\top }x+{\bar{d}}_S^i\ge 0,$$

$$\displaystyle\texttt{A\{k,1\}}=\left[\begin{array}[]{c}(g_{S}^{1})^{\top}\\ (\bar{g}_{S}^{1})^{\top}\\ A_{S}^{1}\\ \vdots\\ (g_{S}^{\ell})^{\top}\\ (\bar{g}_{S}^{\ell})^{\top}\\ A_{S}^{\ell}\end{array}\right],\ \ \ \texttt{b\{k,1\}}=\left[\begin{array}[]{c}d_{S}^{1}\\ \bar{d}_{S}^{1}\\ b_{S}^{1}\\ \vdots\\ d_{S}^{\ell}\\ \bar{d}_{S}^{\ell}\\ b_{S}^{\ell}\end{array}\right]$$

$$\text{cons{k,1}=’SOCPR’},\text{cons{k,2}}=[m_S^1,\dots ,m_S^{\ell }].$$

$$F_0^i+x_1{F}_1^i+\cdots +x_n{F}_n^i⪰0,i\in \{1,\dots ,\ell \}.$$

$$\text{vec2sm(sm2vec(X))=X}.$$

$$⟨X,Y⟩={\text{sm2vec(X)}}^{\top }\text{sm2vec(Y)}.$$

$$\displaystyle\texttt{A\{k,1\}}:=\left[\begin{array}[]{c}\texttt{sm2vec}(F^{1}_{1}),\cdots,\texttt{sm2vec}(F^{1}_{n})\\ \vdots\\ \texttt{sm2vec}(F^{\ell}_{1}),\cdots,\texttt{sm2vec}(F^{\ell}_{n})\end{array}\right],\ \ \ b\{k,1\}:=\left[\begin{array}[]{c}\texttt{sm2vec}(F^{1}_{0})\\ \vdots\\ \texttt{sm2vec}(F^{\ell}_{0})\end{array}\right],$$

$$\text{cons{k,1}=’SDP’}\text{cons{k,2}}=[n^1,\dots ,n^{\ell }].$$

$$\displaystyle A_{0}=\left[\begin{array}[]{ccc}2&-0.5&-0.6\\ -0.5&2&0.4\\ -0.6&0.4&3\end{array}\right],\ A_{1}=\left[\begin{array}[]{ccc}0&1&0\\ 1&0&0\\ 0&0&0\end{array}\right],\ A_{2}=\left[\begin{array}[]{ccc}0&0&1\\ 0&0&0\\ 1&0&0\end{array}\right],\ A_{3}=\left[\begin{array}[]{ccc}0&0&0\\ 0&0&1\\ 0&1&0\end{array}\right].$$

$$\min$$

$$\text{cons{1,1}=’LP’},\text{cons{1,2}}=[1],\text{cons{2,1}=’SDP’},\text{cons{2,2}}=[3],$$

$$\displaystyle\texttt{A\{1,1\}}=\left[\begin{array}[]{cccc}1&1&1&0\end{array}\right],\ \ \ \texttt{b\{1,1\}}=\left[\begin{array}[]{c}-1\end{array}\right],$$

$$\displaystyle\texttt{A\{2,1\}}=\left[\begin{array}[]{cccc}0&0&0&1\\ -1&0&0&0\\ 0&-1&0&0\\ -1&0&0&0\\ 0&0&0&1\\ 0&0&-1&0\\ 0&-1&0&0\\ 0&0&-1&0\\ 0&0&0&1\end{array}\right],\ \ \ \texttt{b\{2,1\}}=\left[\begin{array}[]{c}-2\\ 0.5\\ 0.6\\ 0.5\\ -2\\ -0.4\\ 0.6\\ -0.4\\ -3\end{array}\right],$$

$$\text{c}=(0,0,0,1{)}^{\top }.$$

$$x^{\top }{A}_i^{\top }{Q}_i{A}_{i}x+b_i^{\top }x+d_i\le 0,i\in \{1,\dots ,\ell \},$$

$$\displaystyle\texttt{A\{k,1\}}=\left[\begin{array}[]{c}b_{1}^{\top}\\ A_{1}\\ \vdots\\ b_{l}^{\top}\\ A_{\ell}\end{array}\right],\ \ \ \texttt{b\{k,1\}}=\left[\begin{array}[]{c}d_{1}\\ 0\\ \vdots\\ d_{\ell}\\ 0\end{array}\right]$$

$$\text{cons{k,1}=’QC’}\text{cons{k,2}}=[m_1,\dots ,m_{\ell }],$$

$$\text{cons{k,3,i}}=Q_i,i\in \{1,\dots ,\ell \}.$$

$$u^{\top }u+w+d\le 0$$

Full text

Mehdi Karimi (m7karimi@uwaterloo.ca) and Levent Tunçel (ltuncel@math.uwaterloo.ca): Department of Combinatorics and Optimization, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada.

Research of the authors was supported in part by Discovery Grants from NSERC and by U.S. Office of Naval Research under award numbers: N00014-12-1-0049, N00014-15-1-2171 and N00014-18-1-2078.

1. 1 Introduction

2. 1.1 Installation

3. 2 How to use the DDS code

4. 3 LP, SOCP, SDP

5. 3.1 Linear programming (LP) and second-order cone programming (SOCP)

6. 3.2 Semidefinite programming (SDP)

7. 4 Quadratic constraints

8. 5 Generalized Power Cone

9. 6 Epigraphs of matrix norms

10. 6.1 Minimizing nuclear norm

11. 7 Epigraphs of convex univariate functions (geometric, entropy, and $p$-norm programming)

12. 7.1 Constraints involving power functions

13. 8 Vector Relative Entropy

14. 9 Quantum entropy and Quantum relative entropy

15. 9.1 Adding quantum entropy based constraints

16. 9.2 Adding quantum relative entropy based constraints

17. 10 Hyperbolic polynomials

18. 10.1 Different formats for inputting multivariate polynomials

19. 10.2 Adding constraints involving hyperbolic polynomials

20. 11 Equality constraints

21. 12 Primal-heavy version of the algorithm

22. 13 Numerical Results

23. 13.1 LP-SOCP-SDP

24. 13.2 EXP cone optimization problems from CBLIB

25. 13.3 Minimizing Nuclear Norm

26. 13.4 Quantum Entropy and Quantum Relative Entropy

27. 13.5 Hyperbolic Polynomials

28. A Calculating the predictor and corrector steps

29. B Implementation details for SDP and generalized epigraphs of matrix norms

30. B.1 SDP

31. B.2 Generalized epigraphs of matrix norms

32. C Quantum entropy and quantum relative entropy

33. D Univariate convex functions

34. D.1 Legendre-Fenchel conjugates of univariate convex functions

35. D.2 Gradient and Hessian

36. E A s.c. barrier for vector relative entropy

37. F Comparison with some related solvers and modeling systems

38. F.1 MOSEK [23]

39. F.2 SDPT3 [36, 33]

40. F.3 SeDuMi [32]

41. F.4 CVX [16]

1. Introduction

The code DDS (Domain-Driven Solver) solves convex optimization problems of the form

[TABLE]

where $x\mapsto Ax:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is a linear embedding, $A$ and $c\in\mathbb{\mathbb{R}}^{n}$ are given, and $D\subset\mathbb{\mathbb{R}}^{m}$ is a closed convex set defined as the closure of the domain of a $\vartheta$-self-concordant (s.c.) barrier $\Phi$ [27, 26]. In practice, the set $D$ is typically formulated as $D=D_{1}\oplus\cdots\oplus D_{\ell}$, where $D_{i}$ is associated with a s.c. barrier $\Phi_{i}$, for $i\in\{1,\ldots,\ell\}$. Every input constraint for DDS may be thought of as either the convex set it defines or the corresponding s.c. barrier.

The current version of DDS accepts many functions and set constraints as we explain in this article. If a user has a nonlinear convex objective function $f(x)$ to minimize, one can introduce a new variable $x_{n+1}$ and minimize a linear function $x_{n+1}$ subject to the convex set constraint $f(x)\leq x_{n+1}$ (and other convex constraints in the original optimization problem). As a result, in this article we will talk about representing functions and convex set constraints interchangeably. The algorithm underlying the code also uses the Legendre-Fenchel (LF) conjugate $\Phi_{*}$ of $\Phi$ if it is computationally efficient to evaluate $\Phi_{*}$ and its first (and hopefully second) derivative. Any new discovery of a s.c. barrier allows DDS to expand the classes of convex optimization problems it can solve as any new s.c. barrier $\Phi$ with a computable LF conjugate can be easily added to the code. DDS is a practical implementation of the primal-dual algorithm designed and analyzed in [19], which has the current best iteration complexity bound available for conic formulations. Stopping criteria for DDS and the way DDS suggests the status (“has an approximately optimal solution”, “is infeasible”, “is unbounded”, etc.) is based on the analyses in [20].

Even though there are similarities between DDS and some modeling systems such as CVX [16] (such as the variety input constraints), there are major differences, including:

• •

DDS is not just a modeling system and it uses its own algorithm. The algorithm used in DDS is an infeasible-start primal-dual path-following algorithm, and is of predictor corrector type [19].

• •

The modeling systems for convex optimization that rely on SDP solvers have to use approximation for set constraints which are not efficiently representable by spectrahedra (for example epigraph of functions involving $exp$ or ln functions). However, DDS uses a s.c. barrier specifically well-suited to each set constraint without such approximations. This enables DDS to return proper certificates for all the input problems.

• •

As far as we know, some set constraints such as hyperbolic ones, are not accepted by other modeling systems. Some other constraints such as epigraphs of matrix norms, and those involving quantum entropy and quantum relative entropy are handled more efficiently than other existing options.

The main part of the article is organized to be more as a users’ guide, and contains the minimum theoretical content required by the users. In Section 2, we explain how to install and use DDS. Sections 3-11 explain how to input different types of set/function constraints into DDS. Section 13 contains several numerical results and tables. Many theoretical foundations needed by DDS are included in the appendices. The LHS matrix of the linear systems of equations determining the predictor and corrector steps have a similar form. In Appendix A, we explain how such linear systems are being solved for DDS. Some of the constraints accepted by DDS are involving matrix functions and many techniques are used in DDS to efficiently evaluate these functions and their derivatives, see Appendices B and C. DDS heavily relies on efficient evaluation of LF conjugates of s.c. barrier functions. For some of the functions, we show in the main text how to efficiently calculate the LF conjugate. For some others, the calculations are shifted to the appendices, for example Appendix D. DDS uses interior-point methods, and gradient and Hessian of the s.c. barriers are building blocks of the underlying linear systems. Explicit formulas for the gradient and Hessian of some s.c. barriers and their LF conjugates are given in Appendix D.2. Appendix F contains a brief comparison of DDS with some popular solvers and modeling systems.

1.1. Installation

The current version of DDS is written in MATLAB. This version is available from the websites:

http://www.math.uwaterloo.ca/~m7karimi/DDS.html

https://github.com/mehdi-karimi-math/DDS

To use DDS, the user can follow these steps:

• •

unzip DDS.zip;

• •

run MATLAB in the directory DDS;

• •

run the m-file DDS_startup.m.

• •

(optional) run the prepared small examples DDS_example_1.m and DDS_example_2.m.

The prepared examples contain many set constraints accepted by DDS and running them without error indicates that DDS is ready to use. There is a directory Text_Examples in the DDS package which includes many examples on different classes of convex optimization problems.

2. How to use the DDS code

In this section, we explain the format of the input for many popular classes of optimization problems. In practice, we typically have $D=\bar{D}-b$, where $\textup{int}\bar{D}$ is the domain of a “canonical” s.c. barrier and $b\in\mathbb{R}^{m}$. For example, for LP, we typically have $D=\mathbb{R}^{n}_{+}+b$, where $b\in\mathbb{R}^{m}$ is given as part of the input data, and $-\sum_{i=1}^{n}\textup{ln}(x_{i})$ is a s.c. barrier for $\mathbb{R}^{n}_{+}$. The command in MATLAB that calls DDS is

[x,y,info]=DDS(c,A,b,cons,OPTIONS);

Input Arguments:

cons: A cell array that contains the information about the type of constraints.

c,A,b: Input data for DDS: $A$ is the coefficient matrix, $c$ is the objective vector, $b$ is the RHS vector (i.e., the shift in the definition of the convex domain $D$).

OPTIONS (optional): An array which contains information about the tolerance and initial points.

Output Arguments:

x: Primal point.

y: Dual point which is a cell array. Each cell contains the dual solution for the constraints in the corresponding cell in A.

info: A structure array containing performance information such as info.time, which returns the CPU time for solving the problem.

Note that in the Domain-Driven setup, the primal problem is the main problem, and the dual problem is implicit for the user. This implicit dual problem is:

[TABLE]

where $\delta_{*}(y|D):=\sup\{\langle y,z\rangle:z\in D\},$ is the support function of $D$, and $D_{*}$ is defined as

[TABLE]

where $\text{rec}(D)$ is the recession cone of $D$. For the convenience of users, there is a built-in function in DDS package to calculate the dual objective value of the returned $y$ vector:

z=dual_obj_value(y,b,cons);

For a primal feasible point $x\in\mathbb{R}^{n}$ which satisfies $Ax\in D$ and a dual feasible point $y\in D_{*}$, the duality gap is defined in [19] as

[TABLE]

It is proved in [19] that the duality gap is well-defined and zero duality gap implies optimality. If DDS returns status “solved” for a problem (info.status=1), it means (x,y) is a pair of approximately feasible primal and dual points, with duality gap close to zero (based on tolerance). If info.status=2, the problem is suspected to be unbounded and the returned x is a point, approximately primal feasible with very small objective value ($\langle c,x\rangle\leq-1/tol$). If info.status=3, problem is suspected to be infeasible, and the returned y in $D_{*}$ approximately satisfies $A^{\top}y=0$ with $\delta_{*}(y|D)<0$. If info.status=4, problem is suspected to be ill-conditioned.

The user is not required to input any part of the OPTIONS array. The default settings are:

• •

$tol=10^{-8}$.

• •

The initial points $x^{0}$ and $z^{0}$ for the infeasible-start algorithm are chosen such that, assuming $D=D_{1}\oplus\cdots\oplus D_{\ell}$, the $i$th part of $Ax^{0}+z^{0}$ is a canonical point in $\textup{int}D_{i}$.

However, if a user chooses to provide OPTIONS as an input, here is how to define the desired parts: OPTIONS.tol may be given as the desired tolerance, otherwise the default $tol:=10^{-8}$ is used. OPTIONS.x0 and OPTIONS.z0 may be defined as the initial points as any pair of points $x^{0}\in\mathbb{R}^{n}$ and $z^{0}\in\mathbb{R}^{m}$ that satisfy $Ax^{0}+z^{0}\in\textup{int}D$. If only OPTIONS.x0 is given, then $x^{0}$ must satisfy $Ax^{0}\in\textup{int}D$. In other words, OPTIONS.x0 is a point that strictly satisfies all the constraints.

In the following sections, we discuss the format of each input function/set constraint. Table 1 shows the classes of function/set constraints the current version of DDS accepts, plus the abbreviation we use to represent the constraint. From now on, we assume that the objective function is “$\inf\ \langle c,x\rangle$”, and we show how to add various function/set constraints. Note that A, b, and cons are cell arrays in MATLAB. cons(k,1) represents type of the $k$th block of constraints by using the abbreviations of Table 1. For example, cons(2,1)=’LP’ means that the second block of constraints are linear inequalities. It is advisable to group the constraints of the same type in one block, but not necessary.

3. LP, SOCP, SDP

3.1. Linear programming (LP) and second-order cone programming (SOCP)

Suppose we want to add $\ell$ LP constraints of the form

[TABLE]

where $A_{L}^{i}$ is an $m_{L}^{i}$-by-$n$ matrix, as the $k$th block of constraints. Then, we define

[TABLE]

Similarly to add $\ell$ SOCP constraints of the form

[TABLE]

where $A_{S}^{i}$ is an $m_{S}^{i}$-by-$n$ matrix for $i=\in\{1,\ldots,\ell\}$, as the $k$th block, we define

[TABLE]

Let us see an example:

Example 3.1**.**

Suppose we are given the problem:

[TABLE]

Then we define

[TABLE]

The s.c. barriers and their LF conjugates being used in DDS for these constraints are

[TABLE]

DDS also accepts constraints defined by the rotated second order cones:

[TABLE]

which is handled by the s.c. barrier $-\textup{ln}(ts-z^{\top}z)$. The abbreviation we use is ’SOCPR’. To add $\ell$ rotated SOCP constraints of the form

[TABLE]

where $A_{S}^{i}$ is an $m_{S}^{i}$-by-$n$ matrix for $i\in\{1,\ldots,\ell\}$, as the $k$th block, we define

[TABLE]

3.2. Semidefinite programming (SDP)

Consider $\ell$ SDP constraints in standard inequality (linear matrix inequality (LMI)) form:

[TABLE]

$F^{i}_{j}$’s are $n_{i}$-by-$n_{i}$ symmetric matrices. The above optimization problem is in the matrix form. To formulate it in our setup, we need to write it in the vector form. DDS has two internal functions sm2vec and vec2sm. sm2vec takes an $n$-by-$n$ symmetric matrix and changes it into a vector in $\mathbb{R}^{n^{2}}$ by stacking the columns of it on top of one another in order. vec2sm changes a vector into a symmetric matrix such that

[TABLE]

By this definition, it is easy to check that for any pair of $n$-by-$n$ symmetric matrices $X$ and $Y$ we have

[TABLE]

To give (53) to DDS as the $k$th input block, we define:

[TABLE]

The s.c. barrier used in DDS for SDP is the well-known function $-\textup{ln}(\det(X))$ defined on the convex cone of symmetric positive definite matrices.

Example 3.2**.**

Assume that we want to find scalars $x_{1}$, $x_{2}$, and $x_{3}$ such that $x_{1}+x_{2}+x_{3}\geq 1$ and the maximum eigenvalue of $A_{0}+x_{1}A_{1}+x_{2}A_{2}+x_{3}A_{3}$ is minimized, where

[TABLE]

We can write this problem as

[TABLE]

To solve this problem, we define:

[TABLE]

Then DDS(c,A,b,cons) gives the answer $x=(1.1265,0.6,-0.4,3)^{\top}$, which means the minimum largest eigenvalue is $3$.

4. Quadratic constraints

Suppose we want to add the following constraints to DDS:

[TABLE]

where each $A_{i}$ is $m_{i}$-by-$n$ with rank $n$, and $Q_{i}\in\mathbb{S}^{m_{i}}$. In general, this type of constraints may be non-convex and difficult to handle. Currently, DDS handles two cases:

• •

$Q_{i}$ is positive semidefinite,

• •

$Q_{i}$ has exactly one negative eigenvalue. In this case, DDS considers the intersection of the set of points satisfying (80) and a shifted hyperbolicity cone defined by the quadratic inequality $y^{\top}Q_{i}y\leq 0$.

To give constraints in (80) as input to DDS as the $k$th block, we define

[TABLE]

If cons{k,3} is not given as the input, DDS takes all $Q_{i}$’s to be identity matrices.

If $Q_{i}$ is positive semidefinite, then the corresponding constraint in (80) can be written as

[TABLE]

where $Q_{i}=R_{i}^{\top}R_{i}$ is a Cholesky factorization of $Q_{i}$. We associate the following s.c. barrier and its LF conjugate to such quadratic constraints:

[TABLE]

If $Q_{i}$ has exactly one negative eigenvalue with eigenvector $v$, then $-y^{\top}Q_{i}y$ is a hyperbolic polynomial with respect to $v$. The hyperbolicity cone is the connected component of $y^{\top}Q_{i}y\leq 0$ which contains $v$ and $-\textup{ln}(-y^{\top}Q_{i}y)$ is a s.c. barrier for this cone.

If for any of the inequalities in (80), $Q_{i}$ has exactly one negative eigenvalue while $b_{i}=0$ and $d_{i}=0$, DDS considers the hyperbolicity cone defined by the inequality as the set constraint.

5. Generalized Power Cone

We define the $(m,n)$-generalized power cone with parameter $\alpha$ as

[TABLE]

where $\alpha$ belongs to the simplex $\{\alpha\in\mathbb{R}^{m}_{+}:\sum_{i=1}^{m}\alpha_{i}=1\}$. Note that the rotated second order cone is a special case where $m=2$ and $\alpha_{1}=\alpha_{2}=\frac{1}{2}$. Different s.c. barriers for this cone or special cases of it have been considered [8, 35]. Chares conjectured a s.c. barrier for the cone in (95), which is proved in [31]. This function that we use in DDS is:

[TABLE]

To add generalized power cone constraints to DDS, we use the abbreviation ’GPC’. Therefore, if the $k$th block of constraints is GNC, we define cons{k,1}=’GPC’. Assume that we want to input the following $\ell$ constraints to DDS:

[TABLE]

where $A_{s}^{i}$, $b_{s}^{i}$, $A_{u}^{i}$, and $b_{u}^{i}$ are matrices and vectors of proper size. Then, to input these constraints as the $k$th block, we define cons{k,2} as a MATLAB cell array of size $\ell$-by-$2$, each row represents one constraint. We then define:

[TABLE]

For matrices $A$ and $b$, we define:

[TABLE]

Example 5.1**.**

Consider the following optimization problem with ’LP’ and ’GPC’ constraints:

[TABLE]

Then we define:

[TABLE]

6. Epigraphs of matrix norms

Assume that we have constraints of the form

[TABLE]

where $A_{i}$, $i\in\{1,\ldots,\ell\}$, are $m$-by-$m$ symmetric matrices, and $B_{i}$, $i\in\{1,\ldots,\ell\}$, are $m$-by-$n$ matrices. The set $\{(Z,U)\in\mathbb{S}^{m}\oplus\mathbb{R}^{m\times n}:Z-UU^{\top}\succeq 0\}$ is handled by the following s.c. barrier:

[TABLE]

with LF conjugate

[TABLE]

where $Y\in\mathbb{R}^{m\times m}$ and $V\in\mathbb{R}^{m\times n}$ [25]. This constraint can be reformulated as an SDP constraint using a Schur complement. However, $\Phi(Z,U)$ is a $m$-s.c. barrier while the size of SDP reformulation is $m+n$. For the cases that $m\ll n$, using the Domain-Driven form may be advantageous. A special but very important application is minimizing the nuclear norm of a matrix, which we describe in a separate subsection in the following.

DDS has two internal functions m2vec and vec2m for converting matrices (not necessarily symmetric) to vectors and vice versa. For an $m$-by-$n$ matrix $Z$, m2vec(Z,n) change the matrix into a vector. vec2m(v,m) reshapes a vector $v$ of proper size to a matrix with $m$ rows. The abbreviation we use for epigraph of a matrix norm is MN. If the $k$th input block is of this type, cons{k,2} is a $\ell$-by-$2$ matrix, where $\ell$ is the number of constraints of this type, and each row is of the form $[m\ \ n]$. For each constraint of the form (6), the corresponding parts in $A$ and $b$ are defined as

[TABLE]

For implementation details involving epigraph of matrix norms, see Appendix B.2. Here is an example:

Example 6.1**.**

Assume that we have matrices

[TABLE]

and our goal is to solve

[TABLE]

Then the input to DDS is defined as

[TABLE]

6.1. Minimizing nuclear norm

The nuclear norm of a matrix $Z$ is $\|Z\|_{*}:=\textup{Tr}\left((ZZ^{\top})^{1/2}\right)$. The dual norm of $\|\cdot\|_{*}$ is the operator 2-norm $\|\cdot\|$ of a matrix. Minimization of nuclear norm has application in machine learning and matrix sparsification. The following optimization problems are a primal-dual pair [29].

[TABLE]

where $A$ is a linear transformation on matrices and $A^{*}$ is its adjoint. $(P_{N})$ is a very popular relaxation of the problem of minimizing $\textup{rank}(X)$ subject to $A(X)=b$, with applications in machine learning and compressed sensing. The dual problem $(D_{N})$ is a special case of (6) where $Z=I$ and $U=A^{*}(z)$. As we will show on an example, solving $(D_{N})$ by [x,y]=DDS(c,A,b,Z) leads us to $y$, which gives a solution for $(P_{N})$.

More specifically, consider the optimization problem

[TABLE]

where $X$ is $n$-by-$m$. DDS can solve the dual problem by defining

[TABLE]

Then, if we run [x,y]=DDS(c,A,b,cons) and define V:=(vec2m(y{1}(1:m*n),m))⊤, then $V$ is an optimal solution for (134). In subsection 13.3, we present numerical results for solving problem (134) and show that for the cases that $n\gg m$, DDS can be more efficient than SDP based solvers. Here is an example:

Example 6.2**.**

We consider minimizing the nuclear norm over a subspace. Consider the following optimization problem:

[TABLE]

where

[TABLE]

By using (133), the dual of this problem is

[TABLE]

To solve this problem with our code, we define

[TABLE]

If we solve the problem using [x,y]=DDS(c,A,b,cons), the optimal value is $-2.2360$. Now V:=(vec2m(y{1}(1:8),2))⊤ is the solution of (141) with objective value $2.2360$. We have

[TABLE]

7. Epigraphs of convex univariate functions (geometric, entropy, and $p$-norm programming)

DDS accepts constraints of the form

[TABLE]

where $\alpha_{i}\geq 0$ and $f_{i}(x)$, $i\in\{1,\ldots,\ell\}$, can be any function from Table 2. Note that every univariate convex function can be added to this table in the same fashion. By using this simple structure, we can model many interesting optimization problems. Geometric programming (GP) [2] and entropy programming (EP) [12] with many applications in engineering are constructed with constraints of the form (154) when $f_{i}(z)=e^{z}$ for $i\in\{1,\ldots,\ell\}$ and $f_{i}(z)=z\textup{ln}(z)$ for $i\in\{1,\ldots,\ell\}$, respectively. The other functions with $p$ powers let us solve optimization problems related to $p$-norm minimization.

The corresponding s.c. barriers used in DDS are shown in Table 2. There is a closed form expression for the LF conjugate of the first two functions. For the last four, the LF conjugate can be calculated to high accuracy efficiently. In Appendix D, we show how to calculate the LF conjugates for the functions in Table 2 and the internal functions we have in DDS.

To represent a constraint of the from (154), for given $\gamma\in\mathbb{R}$ and $\beta_{i}\in\mathbb{R}$, $i\in\{1,\ldots,\ell\}$, we can define the corresponding convex set $D$ as

[TABLE]

and our matrix $A$ represents $w=\sum_{i=1}^{\ell}\alpha_{i}u_{i}+g^{\top}x$ and $s_{i}=a_{i}^{T}x$, $i\in\{1,\ldots,\ell\}$. As can be seen, to represent our set as above, we introduce some auxiliary variables $u_{i}$’s to our formulations. DDS code does this internally. Let us assume that we want to add the following $s$ constraints to our model

[TABLE]

From now on, $type$ indexes the rows of Table 2. The abbreviation we use for these constraints is TD. Hence, if the $k$th input block are the constraints in (156), then we have cons{k,1}=’TD’. cons{k,2} is a cell array of MATLAB with $s$ rows, each row represents one constraint. For the $j$th constraint we have:

• •

cons{k,2}{j,1} is a matrix with two columns: the first column shows the type of a function from Table 2 and the second column shows the number of that function in the constraint. Let us say that in the $j$th constraint, we have $l_{2}^{j}$ functions of type 2 and $l_{3}^{j}$ functions of type 3, then we have

[TABLE]

The types can be in any order, but the functions of the same type are consecutive and the order must match with the rows of $A$ and $b$.

• •

cons{k,2}{j,2} is a vector with the coefficients of the functions in a constraint, i.e., $\alpha_{i}^{j,type}$ in (156). Note that the coefficients must be in the same order as their corresponding rows in $A$ and $b$. If in the $j$th constraint we have 2 functions of type 2 and 1 function of type 3, it starts as

[TABLE]

To add the rows to $A$, for each constraint $j$, we first add $g_{j}$, then $a_{i}^{j,type}$’s in the order that matches cons{k,2}. We do the same thing for vector $b$ (first $\gamma_{j}$, then $\beta_{i}^{j,type}$’s). The part of $A$ and $b$ corresponding to the $j$th constraint is as follows if we have for example five types

[TABLE]

Let us see an example:

Example 7.1**.**

Assume that we want to solve

[TABLE]

For this problem, we define:

[TABLE]

Note:* As we mentioned, modeling systems for convex optimization that are based on SDP solvers, such as CVX, have to use approximation for functions involving $exp$ and ln. Approximation makes it hard to return dual certificates, specifically when the problem is infeasible or unbounded.*

7.1. Constraints involving power functions

The difference between these two types (4 and 5) and the others is that we also need to give the value of $p$ for each function. To do that, we add another column to cons{k,2}.

Note: For TD constraints, cons{k,2} can have two or three columns. If we do not use types 4 and 5, it has two, otherwise three columns. cons{k,2}{j,3} is a vector which contains the powers $p$ for functions of types 4 and 5. The powers are given in the same order as the coefficients in cons{k,2}{j,2}. If the constraint also has functions of other types, we must put 0 in place of the power.

Let us see an example:

Example 7.2**.**

[TABLE]

For this problem, we define:

[TABLE]

8. Vector Relative Entropy

Consider the relative entropy function $f:\mathbb{R}_{++}^{\ell}\oplus\mathbb{R}_{++}^{\ell}\rightarrow\mathbb{R}$ defined as

[TABLE]

The epigraph of this function accepts the following $(2\ell+1)$-s.c. barrier (see Appendix E):

[TABLE]

This function covers many applications as discussed in [10]. The special case of $\ell=1$ is used for handling EXP cone in the commercial solver MOSEK. We claim that the LF of $\Phi$ can be efficiently calculated using $\ell$ one-dimensional minimizations, i.e., having the function $\theta(r)$ as the solution of $\frac{1}{\theta}-\textup{ln}(\theta)=r$. We consider the case of $\ell=1$ and the generalization is straightforward. We have

[TABLE]

Writing the optimality conditions, we get (defining $T:=t-u\textup{ln}(u)+u\textup{ln}(z)$)

[TABLE]

We can see that at the optimal solution, $T=-1/\alpha$, and if we calculate $u$, we can easily get $z$. $u$ is calculated by the following equation:

[TABLE]

Therefore, to calculate $\Phi_{*}$ at every point, we need to evaluate $\theta$ once. For the general $\ell$, we evaluate $\theta$ by $\ell$ times.

The abbreviation we use for relative entropy is RE. So, for the $k$th block of $s$ constraints of the form

[TABLE]

we define cons{k,1} = ’RE’ and cons{k,2} is a vector of length $s$ with the $i$th element equal to $2\ell+1$. We also define:

[TABLE]

Example 8.1**.**

Assume that we want to minimize a relative entropy function under a linear constraint:

[TABLE]

We add an auxiliary variable $x_{3}$ to model the objective function as a constraint. For this problem we define:

[TABLE]

If we solve this problem by DDS, for $\beta=2$ the problem is infeasible, and for $\beta=7$ it returns an optimal solution $x^{*}:=(5.93,1.06)^{\top}$ with the minimum of function equal to $-7.259$.

9. Quantum entropy and Quantum relative entropy

Quantum entropy and quantum relative entropy functions are important in quantum information processing. DDS 2.0 accepts constraints involving these two functions. Let us start with the main definitions. Consider a function $f:\mathbb{R}\rightarrow\mathbb{R}\cup\{+\infty\}$ and let $X\in\mathbb{H}^{n}$ be a Hermitian matrix (with entries from $\mathbb{C}$) with a spectral decomposition $X=U\textup{Diag}(\lambda_{1},\ldots,\lambda_{n})U^{*}$, where Diag returns a diagonal matrix with the given entries on its diagonal and $U^{*}$ is the conjugate transpose of a unitary matrix $U$. We define the matrix extension $F$ of $f$ as $F(X):=U\textup{Diag}(f(\lambda_{1}),\ldots,f(\lambda_{n}))U^{*}$. Whenever we write $\phi(X):=\textup{Tr}(F(X))$, we mean a function $\phi:\mathbb{H}^{n}\rightarrow\mathbb{R}\cup\{+\infty\}$ defined as

[TABLE]

Study of such matrix functions go back to the work of Löwner as well as Von-Neumann (see [11], [21], and the references therein). A function $f:(a,b)\mapsto\mathbb{R}$ is said to be matrix monotone if for any two self-adjoint matrix $X$ and $Y$ with eigenvalues in $(a,b)$ that satisfy $X\succeq Y$, we have $F(X)\succeq F(Y)$. A function $f:(a,b)\mapsto\mathbb{R}$ is said to be matrix convex if for any pair of self-adjoint matrices $X$ and $Y$ with eigenvalues in $(a,b)$, we have

[TABLE]

Faybusovich and Tsuchiya [14] utilized the connection between the matrix monotone functions and self-concordant functions. Let $f$ be a continuously differentiable function whose derivative is matrix monotone on the positive semi-axis and let us define $\phi$ as (208). Then, the function

[TABLE]

is a $(n+1)$-s.c. barrier for the epigraph of $\phi(X)$. This convex set has many applications. Many optimization problems arising in quantum information theory and some other areas requires dealing with the so-called quantum or von Neumann entropy $qe:\mathbb{H}^{n}\rightarrow\mathbb{R}\cup\{+\infty\}$ which is defined as $qe(X):=\textup{Tr}(X\textup{ln}(X))$. If we consider $f(x):=x\textup{ln}(x)$, then $f^{\prime}(x)=1+\textup{ln}(x)$ is matrix monotone on $(0,\infty)$ (see, for instance [18]-Example 4.2) and so we have a s.c. barrier for the set

[TABLE]

We have to solve the optimization problem

[TABLE]

to calculate the LF conjugate of (210), which is done in Appendix C.

Another interesting function with applications in quantum information theory is quantum relative entropy $qre:\mathbb{H}^{n}\oplus\mathbb{H}^{n}\rightarrow\mathbb{R}\cup\{+\infty\}$ defined as

[TABLE]

This function is convex as proved in [34]. The epigraph of $qre$ is

[TABLE]

DDS 2.0 uses the following barrier (not yet known to be s.c.) for solving problems involving quantum relative entropy constraints:

[TABLE]

We can input constraints involving quantum entropy and quantum relative entropy into DDS as we explain in the following subsections. Appendix C contains some interesting theoretical results about quantum entropy functions and how to calculate the derivative and Hessian for the above s.c. barriers and also the LF conjugate given in (211).

9.1. Adding quantum entropy based constraints

Let $qe_{i}:\mathbb{S}^{n_{i}}\rightarrow\mathbb{R}\cup\{+\infty\}$ be quantum entropy functions and consider $\ell$ quantum entropy constraints of the form

[TABLE]

$A^{i}_{j}$’s are $n_{i}$-by-$n_{i}$ symmetric matrices. To input (212) to DDS as the $k$th block, we define:

[TABLE]

Example 9.1**.**

Assume that we want to find scalars $x_{1}$, $x_{2}$, and $x_{3}$ such that $2x_{1}+3x_{2}-x_{3}\leq 5$ and all the eigenvalues of $H:=x_{1}A_{1}+x_{2}A_{2}+x_{3}A_{3}$ are at least 3, for

[TABLE]

such that the quantum entropy $H$ is minimized. We can write this problem as

[TABLE]

For the objective function we have $\texttt{c}=(0,0,0,1)^{\top}$. Assume that the first and second blocks are LP and SDP as before. We define the third block of constraints as:

[TABLE]

If we run DDS, the answer we get is $(x_{1},x_{2},x_{3})=(4,-1,0)$ with $f(H)=14.63$.

9.2. Adding quantum relative entropy based constraints

The abbreviation we use for quantum relative entropy is QRE. Let $qre_{i}:\mathbb{S}^{n_{i}}\oplus\mathbb{S}^{n_{i}}\rightarrow\mathbb{R}\cup\{+\infty\}$ be quantum relative entropy functions and consider $\ell$ quantum entropy constraints of the form

[TABLE]

$A^{i}_{j}$’s and $B^{i}_{j}$’s are $n_{i}$-by-$n_{i}$ symmetric matrices. To input (212) to DDS as the $k$th block, we define:

[TABLE]

10. Hyperbolic polynomials

Hyperbolic programming problems form one of the largest classes of convex optimization problems which have the kind of special mathematical structure amenable to more robust and efficient computational approaches. DDS 2.0 accepts hyperbolic constraints. Let us first define a hyperbolic polynomial. A polynomial $p(x)\in\mathbb{R}[x_{1},\ldots,x_{m}]$ is said to be homogeneous if every term has the same degree $d$. A homogeneous polynomial $p:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is hyperbolic in direction $e\in\mathbb{R}^{m}$ if

• •

$p(e)>0$.

• •

for every $x\in\mathbb{R}^{m}$, the univariate polynomial $p(x-te)$ has only real roots.

Let $p(x)\in\mathbb{R}[x_{1},\ldots,x_{m}]$ be hyperbolic in direction $e$. We define the eigenvalue function $\lambda:\mathbb{R}^{m}\rightarrow\mathbb{R}^{d}$ with respect to $p$ and $e$ such that for every $x\in\mathbb{R}^{m}$, the elements of $\lambda(x)$ are the roots of the univariate polynomial $p(x-te)$. The underlying hyperbolicity cone is defined as

[TABLE]

Example 10.1**.**

The polynomial $p(x)=x_{1}^{2}-x_{2}^{2}-\cdots-x_{m}^{2}$ is hyperbolic in the direction $e:=(1,0,\ldots,0)^{\top}$ and the hyperbolicity cone with respect to $e$ is the second-order cone. The polynomial $p(X)=\det(X)$ defined on $\mathbb{S}^{n}$ is hyperbolic in the direction $I$, and the hyperbolicity cone with respect to $I$ is the positive-semidefinite cone.

By the above example, optimization over hyperbolicity cone is an extension of SOCP and SDP. The following theorem by Güler gives a s.c. barrier for the hyperbolicity cone.

Theorem 10.1** (Güler [17]).**

Let $p(x)$ be a homogeneous polynomial of degree $d$, which is hyperbolic in direction $e$. Then, the function $-\textup{ln}(p(x))$ is a $d$-LH s.c. barrier for $\Lambda_{+}(p,e)$.

DDS handles optimization problems involving hyperbolic polynomials using the above s.c. barrier. A computational problem is that, currently, we do not have a practical, efficient, algorithm to evaluate the LF conjugate of $-\textup{ln}(p(x))$. Therefore, DDS uses a primal-heavy version of the algorithm for these problems, as we explain in Section 12.

10.1. Different formats for inputting multivariate polynomials

To input constraints involving hyperbolic polynomials, we use a matrix named poly. In DDS, there are different options to input a multivariate polynomial:

Using monomials: In this representation, if $p(x)$ is a polynomial of $m$ variables, then poly is an $k$-by-$(m+1)$ matrix, where $k$ is the number of monomials. In the $j$th row, the first $m$ entries are the powers of the $m$ variables in the monomial, and the last entry is the coefficient of the monomial in $p(x)$. For example, if $p(x)=x_{1}^{2}-x_{2}^{2}-x_{3}^{2}$, then

[TABLE]

Note: In many applications, the above matrix is very sparse. DDS recommends that in the monomial format, poly should be defined as a sparse matrix.

Using straight-line program: Another way to represent a polynomial is by a straight-line program, which can be seen as a rooted acyclic directed graph. The leaves represent the variables or constants. Each node is a simple binary operation (such as addition or multiplication), and the root is the result of the polynomial. In this case, poly is a $k$-by-4 matrix, where each row represents a simple operation. Assume that $p(x)$ has $m$ variables, then we define

[TABLE]

The $\ell$th row of poly is of the form $[\alpha_{j}\ \ i\ \ j\ \ \Box]$, which means that

[TABLE]

Operations are indexed by 2-digit numbers as the following table:

[TABLE]

For example, if $p(x)=x_{1}^{2}-x_{2}^{2}-x_{3}^{2}$, we have the following representation:

[TABLE]

Note that straight-line program is not unique for a polynomial.

Determinantal representation: In this case, if possible, the polynomial $p(x)$ is written as

[TABLE]

where $H_{i}$, $i\in\{0,1,\ldots,m\}$ are in $\mathbb{S}^{m}$. We define

[TABLE]

For example, for $p(x)=x_{1}^{2}-x_{2}^{2}-x_{3}^{2}$, we can have

[TABLE]

Note: $H_{1},\ldots,H_{m}$ must be nonzero symmetric matrices.

10.2. Adding constraints involving hyperbolic polynomials

Consider a hyperbolic polynomial constraint of the form

[TABLE]

To input this constraint to DDS as the $k$th block, $A$ and $b$ are defined as before, and different parts of cons are defined as follows:

cons{k,1}=’HB’,

cons{k,2}= number of variables in $p(z)$.

cons{k,3} is the poly that can be given in one of the three formats of Subsection 10.1.

cons{k,4} is the format of polynomial that can be ’monomial’, ’straight_line’, or ’determinant’.

cons{k,5} is the direction of hyperbolicity or a vector in the interior of the hyperbolicity cone.

Example 10.2**.**

Assume that we want to give constraint (256) to DDS for $p(x)=x_{1}^{2}-x_{2}^{2}-x_{3}^{2}$, using the monomial format. Then, cons part is defined as

[TABLE]

11. Equality constraints

In the main formulation of Domain-Driven form (1), equality constraints are not included. However, the current version of DDS accepts equality constraints. In other form, users can input a problem of the form

[TABLE]

To input a block of equality constraints $Bx=d$, where $B$ is a $m_{EQ}$-by-$n$ matrix, as the $k$th block, we define

[TABLE]

Example 11.1**.**

If for a given problem with $x\in\mathbb{R}^{3}$, we have a constraint $x_{2}-x_{3}=2$, then we can add it as the $k$th block by the following definitions:

[TABLE]

In the main Domain-Driven formulation (1), for some theoretical and practical considerations, we did not consider the equality constraints. From a mathematical point of view, this is not a problem. For any matrix $Z$ whose columns form a basis for the null space of $B$, we may represent the solution set of $Bx=d$ as $x^{0}+Zw$, where $x^{0}$ is any solution of $Bx=d$. Then the feasible region in (259) is equivalent to:

[TABLE]

$D-Ax^{0}$ is a translation of $D$ with the s.c. barrier $\Phi(z-Ax^{0})$, where $\Phi$ is a s.c. barrier for $D$. Therefore, theoretically, we can reformulate (259) as (1), where $AZ$ replaces $A$. Even though this procedure is straightforward in theory, there might be numerical challenges in application. For example, if we have a nice structure for $A$, such as sparsity, multiplying $A$ with $Z$ may ruin the structure. In DDS, we use a combination of heuristics, LU and QR factorizations to implement this transition efficiently.

12. Primal-heavy version of the algorithm

For some class of problems, such as hyperbolic programming, a computable s.c. barrier $\Phi$ is available for the set $D$, while the LF of it is not available. For these classes, DDS uses a primal-heavy version of the algorithm. In the primal-heavy version, we approximate the primal-dual system of equations for computing the search directions by approximating the gradient and Hessian of $\Phi_{*}$. The approximations are based on the relations between the derivatives of $\Phi$ and $\Phi_{*}$: for every point $z\in\textup{int}D$ we have

[TABLE]

Instead of the primal-dual proximity measure defined in [19], we use the primal-heavy version:

[TABLE]

where $u:=Ax+\frac{1}{\tau}z^{0}$, $\tau$ is an artificial variable we use in the formulation of the central path (see [19] for details), and $\mu$ is the parameter of the central path. By [19]-Corollary 4.1, this proximity measures is “equivalent” to the primal-dual one, but (262) is less efficient computationally.

By using a primal-heavy version, we lose some of the desired properties of primal-dual setup, such as the ability to move in a wider neighbourhood of the central path. Moreover, in the primal-heavy version, we have to somehow make sure that the dual iterates $y$ are feasible (or at least the final dual iterate is). Another issue is with calculating the duality gap (4). For a general convex domain $D$, we need $\Phi_{*}^{\prime}$ to accurately calculate $\delta_{*}(y|D)$ as explained in [19]. Note that when $D$ is a shifted cone $D=K-b$, then we have

[TABLE]

To calculate the duality gap, we can write it as the summation of separate terms for the domains, and if a domain with only the primal barrier is a shifted cone, we can use (263). This is the case for the current version of DDS as all the set constraints using primal-heavy version are based on a shifted convex cone. To make sure that the dual iterates are feasible, we choose our neighborhoods to satisfy

[TABLE]

and by the Dikin ellipsoid property of s.c. functions, $y$ iterates stay feasible. We can specify in OPTIONS if we want to use a primal-heavy version of the algorithm in DDS by

OPTIONS.primal=1;

13. Numerical Results

In this section, we present some numerical examples of running DDS 2.0. We performed computational experiments using the software MATLAB R2018b, on a 4-core 3.2 GHz Intel Xeon X5672 machine with 96GB of memory. Many of the examples in this section are given as mat-files in the DDS package; in the folder Test_Examples.

DDS pre-processes the input data and returns the statistics before and after the pre-processing. This processing includes checking the structure of the input data removing obviously redundant constraints. For example, if the size of $A$ given for a block of constraints does not match the information in cons, DDS returns an error. The pre-processing also includes scaling the input data properly for each block of constraints.

As mentioned in Section 2, users have the option to give the starting point of the infeasible-start algorithm as input. Otherwise, DDS chooses an easy-to-find starting point for every input set constraint. The stopping criteria in DDS are based on the status determination result in [20], which studied how to efficiently certify different status of a given problem in Domain-Driven form, using duality theory. [20]-Section 9 discusses stopping criteria in a practical setting. We define the following parameters for the current point $(x,\tau,y)$:

[TABLE]

where $x$ and $y$ are the primal and dual points, and $\tau$ is the artificial variable we use (see [19] for details). DDS stops the algorithm and returns the status as “solved” when we have

[TABLE]

DDS returns that status as “infeasible” if the returned certificate $\bar{y}:=\frac{\tau}{\mu}y$ satisfies

[TABLE]

and returns the status as “unbounded” if

[TABLE]

DDS has different criteria for returning “ill-conditioned” status, which can be because of numerical inaccuracy issues.

In the following, we see the numerical performance of DDS for different combinations of set constraints.

13.1. LP-SOCP-SDP

In this subsection, we consider LP-SOCP-SDP instances mostly from the Dimacs library [28]. Note that the problems in the library are for the standard equality form and we solve the dual of the problems. Table 3 shows the results. DDS has a built-in function read_sedumi_DDS to read problems in SeDuMi input format. You can use either of the following commands based on the SeDuMi format file

[c,A,b,cons]=read_sedumi_DDS(At,b,c,K); [c,A,b,cons]=read_sedumi_DDS(A,b,c,K);

13.2. EXP cone optimization problems from CBLIB

Conic Benchmark Library (CBLIB) [15] is a collection of benchmark problems for conic mixed-integer and continuous optimization. It contains a collection of optimization problems involving exponential cone, they show as EXP. Exponential cone is a special case of vector relative entropy we discussed in Section 8 when $\ell=1$. Table 4 show the results of running DDS to solve problems with EXP cone from CBLIB. The files converted into DDS format can be found in the Test_Examples folder.

13.3. Minimizing Nuclear Norm

In this subsection, we are going to use DDS to solve problem (134); minimizing nuclear norm of a matrix subject to linear constraints. We are interested in the case that $X$ is an $n$-by-$m$ matrix with $n\gg m$, which makes the s.c. barrier for the epigraph of a matrix norm more efficient than the SDP representation. We compare DDS with convex modeling system CVX, which has a function norm_nuc for nuclear norm. For numerical results, we assume (134) has two linear constraints and $U_{1}$ and $U_{2}$ are 0-1 matrices generated randomly. Let us fix the number of columns $m$ and calculate the running time by changing the number of rows $n$. Figure 1 shows the plots of running time for both DDS and CVX. For epigraph of a matrix norm constraints, DDS does not form matrices of size $m+n$, and as can be seen in the figure, the running time is almost linear as a function of $n$. On the other hand, the CVX running time is clearly super-linear.

13.4. Quantum Entropy and Quantum Relative Entropy

In this subsection, we see optimization problems involving quantum entropy. As far as we know, there is no library for quantum entropy and quantum relative entropy functions. We have created examples of combinations of quantum entropy constraints and other types such as LP, SOCP, and $p$-norm. Problems for quantum entropy are of the form

[TABLE]

where $qe=\textup{Tr}{X\textup{ln}(X)}$ and $A_{i}$’s are sparse symmetric matrices. The problems for $qre$ are also in the same format. Our collection contains both feasible and infeasible cases. The results are shown in Tables 5 and 6 respectively for quantum entropy and quantum relative entropy. The MATLAB files in the DDS format for some of these problems are given in the Test_Examples folder. We compare our results to CVXQUAD, which is a collection of matrix functions to be used on top of CVX. The package is based on the paper [13], which approximates matrix logarithm with functions that an be described as feasible region of SDPs. Note that these approximations increase the dimension of the problem, which is why CVX fails for the larger problems in Table 5.

13.5. Hyperbolic Polynomials

We have created a library of hyperbolic cone programming problems for our experiments. Hyperbolic polynomials are defined in Section 10. We discuss three methods in this section for the creation of our library.

13.5.1. Hyperbolic polynomials from matriods

Let us first define a matroid.

Definition 13.1**.**

*A ground set $E$ and a collection of independent sets * $\mathcal{I}\subseteq 2^{E}$ form a matroid $M=(E,\mathcal{I})$ if:

• •

$\emptyset\in\mathcal{I}$,

• •

if $A\in\mathcal{I}$ and $B\subset A$, then $B\in\mathcal{I}$,

• •

if $A,B\in\mathcal{I}$, and $|B|>|A|$, then there exists $b\in B\setminus A$ such that $A\cup\{b\}\in\mathcal{I}$.

The independent sets with maximum cardinality are called the bases of the matroid, we denote the set of bases of the matroid by $\mathcal{B}$. We can also assign a rank function $r_{M}:2^{E}\rightarrow\mathbb{Z}_{+}$ as:

[TABLE]

Naturally, the rank of a matroid is the cardinality of any basis of it. The basis generating polynomial of a matroid is defined as

[TABLE]

A class of hyperbolic polynomials are the basis generating polynomials of certain matroids with the half-plane property [9]. A polynomials has the half-plane property if it is nonvanishing whenever all the variables lie in the open right half-plane. We state it as the following lemma:

Lemma 13.1**.**

Assume that $M$ is a matroid with half-plane property. Then its basis generating polynomial is hyperbolic in any direction in the positive orthant.

Several classes of matroids have been proven to have the half-plane property [9, 3, 1]. As the first example, we introduce the Vamos matroid. The ground set has size 8 (can be seen as 8 vertices of the cube) and the rank of matroid is 4. All subsets of size 4 are independent except 5 of them. Vamos matroid has the half-plane property [37]. The basis generating polynomial is:

[TABLE]

where $E_{d}(x_{1},\ldots,x_{m})$ is the elementary symmetric polynomial of degree $d$ with $m$ variables. Note that elementary symmetric polynomials are hyperbolic in direction of the vector of all ones. Some extensions of Vamos matroid also satisfy the half-plane property [6]. These matroids give the following Vamos-like basis generating polynomials:

[TABLE]

Vamos like polynomials in (272) have an interesting property of being counter examples to one generalization of Lax conjecture [4, 6]. Explicitly, there is no power $k$ and symmetric matrices $H_{2},\ldots,H_{2m}$ such that

[TABLE]

In the DDS package, we have created a function

poly = vamos(m)

that returns a Vamos like polynomial as in (272) in the ’monomial’ format, which can be given as cons{k,3} for inputting a hyperbolic constraint. vamos(4) returns the Vamos polynomial with 8 variables.

Graphic matroids also have the half-plane property [9]. However, the hyperbolicity cones arising from graphic matroids are isomorphic to positive semidefinite cones. This can be proved using a classical result that the characteristic polynomial of the bases of graphic matroids is the determinant of the Laplacian (with a row and a column removed) of the underlying graph [5]. Consider a graph $G=(V,E)$ and let $\mathcal{T}$ be the set of all spanning trees of $G$. Then the generating polynomial defined as

[TABLE]

is a hyperbolic polynomial. Several matroid operations preserve the half-plane property as proved in [9], including taking minors, duals, and direct sums.

We have created a set of problems with combinations of hyperbolic polynomial inequalities and other types of constraints such as those arising from LP, SOCP, and entropy function.

13.5.2. Derivatives of products of linear forms

Hyperbolicity of a polynomial is preserved under directional derivative operation (see [30]):

Theorem 13.1**.**

Let $p(x)\in\mathbb{R}[x_{1},\ldots,x_{m}]$ be hyperbolic in direction $e$ and of degree at least 2. Then the polynomial

[TABLE]

is also hyperbolic in direction $e$. Moreover,

[TABLE]

Assume that $p(x):=l_{1}(x)\cdots l_{\ell}(x)$, where $l_{1}(x),\ldots,l_{\ell}(x)$ are linear forms. Then, $p$ is hyperbolic in the direction of any vector $e$ such that $p(e)\neq 0$. As an example, consider $E_{m}(x_{1},\ldots,x_{m})=x_{1}\cdots x_{m}$. Recursively applying Theorem 13.1 to such polynomials $p$ leads to many hyperbolic polynomials, including all elementary symmetric polynomials. For some properties of hyperbolicity cones of elementary symmetric polynomials, see [38]. For some preliminary computational experiments on such hyperbolic programming problems, see [24].

For the polynomials constructed by the products of linear forms, it is more efficient to use the straight-line program. For a polynomial $p(x):=l_{1}(x)\cdots l_{\ell}(x)$, let us define a $\ell$-by-$m$ matrix $L$ where the $j$th row contains the coefficients of $l_{j}$. In DDS, we have created a function

[poly,poly_grad] = lin2stl(M,d)

which returns two polynomials in ”straight-line” form, one as the product of linear forms defined by $M$, and the other one its directional derivative in the direction of $d$. For example, if we want the polynomial $p(x):=(2x_{1}-x_{3})(x_{1}-3x_{2}+4x_{3})$, then our $M$ is defined as

[TABLE]

Table 8 shows some results of problem where the hyperbolic polynomial is either a product of linear forms or their derivatives.

13.5.3. Derivatives of determinant of matrix pencils

Let $H_{1},\ldots,H_{m}\in\mathbb{H}^{n}$ be Hermitian matrices and $e\in\mathbb{R}^{m}$ be such that $e_{1}H_{1}+\cdots+e_{m}H_{m}$ is positive definite. Then, $p(x):=\det(x_{1}H_{1}+\cdots+x_{m}H_{m})$ is hyperbolic in direction $e$. Now, by using Theorem 13.1, we can generate a set of interesting hyperbolic polynomials by taking the directional derivative of this polynomial up to $n-3$ times at direction $e$.

Appendix A Calculating the predictor and corrector steps

As discussed in [19], for both the predictor and corrector steps, the theoretical algorithm solves a linear system with the LHS matrix of the form

[TABLE]

where $U$ is a matrix that contains the linear transformations we need:

[TABLE]

$H$ and $G$ are positive definite matrices based on the Hessians of the s.c. functions, $F$ is a matrix whose rows form a basis for $\textup{null}(A)$, $c_{A}$ is any vector that satisfies $A^{\top}c_{A}=c$, and $\eta_{*}$, $\zeta$, $h$, and $h_{*}$ are scalars and vectors defined in [19]. A practical issue of this system is that calculating $F$ is not computationally efficient, and DDS uses a way around it. If we expand the system given in (286), the linear systems DDS needs to solve become of the form:

[TABLE]

At the end, we are interested in $F^{\top}d_{v}$ to calculate our search directions. If we consider the last equation, we can remove the matrix $F$ multiplied from the left to all the terms as

[TABLE]

where $\bar{G}:=G+\eta_{*}h_{*}h_{*}^{\top}$. Now, we multiply the last equation by $A^{\top}$ from the left and eliminate $d_{v}$ to get

[TABLE]

By using the equations in (293) and (294), we can get a linear system without $F$ as

[TABLE]

By Sherman–Morrison formula, we can write

[TABLE]

Let us define

[TABLE]

Then we can see that the LHS matrix in (304) can be written as

[TABLE]

This matrix is a $(2n+1)$-by-$(2n+1)$ matrix $\tilde{H}$ plus a rank one update. If we have the Cholesky or LU factorization of $A^{\top}HA$ and $A^{\top}G^{-1}A$ (in the case that $G:=\mu^{2}H$, we need just one such factorization), then we have such a factorization for the $2n$-by-$2n$ leading minor of $\tilde{H}$ and we can easily extend it to a factorization for the whole $\tilde{H}$. To solve the whole system, we can then use Sherman-Morrison formula.

Appendix B Implementation details for SDP and generalized epigraphs of matrix norms

DDS accepts many constraints that involve matrices. These matrices are given into DDS as vectors. Calculating the derivatives and implementing the required linear algebra operations efficiently is critical in the performance of DDS. In DDS 2.0, many ideas and tricks have been used for matrix functions. In this section, we mention some of these techniques for SDP and matrix norm constraints.

B.1. SDP

The primal and dual barrier functions being used for SDP constraints (53) are:

[TABLE]

For function $f=-\textup{ln}(\det(X))$, we have:

[TABLE]

To implement our algorithm, for each matrix $X$, we need to find the corresponding gradient $g_{X}$ and Hessian $H_{X}$, such that for any symmetric positive semidefinite matrix $X$ and symmetric matrix $H$ we have:

[TABLE]

It can be shown that $g_{X}=\texttt{sm2vec}(X^{-1})$ and $H_{X}=X^{-1}\otimes X^{-1}$, where $\otimes$ stands for the Kronecker product of two matrices. Although this representation is theoretically nice, there are two important practical issues with it. First, it is not efficient to calculate the inverse of a matrix explicitly. Second, forming and storing $H_{X}$ is not numerically efficient for large scale problems. DDS does not explicitly form the inverse of matrices. An important matrix in calculating the search direction is $A^{\top}\Phi^{\prime\prime}(u)A$, where $\Phi$ is the s.c. barrier for the whole input problem. In DDS, there exist an internal function hessian_A that directly returns $A^{\top}\Phi^{\prime\prime}(u)A$ in an efficient way, optimized for all the set constraints, including SDP. For transitions from matrices to vectors, we use the properties of Kronecker product that for matrices $A$, $B$, and $X$ of proper size, we have

[TABLE]

Similar to hessian_A, there are other internal functions in DDS for efficiently calculating matrix-vector products, such as hessian_V for evaluating the product of Hessian with a given vector of proper size.

B.2. Generalized epigraphs of matrix norms

Let us see how to calculate the first and second derivatives of functions in (113) and (114).

Proposition B.1**.**

(a) Consider $\Phi(X,U)$ defined in (113). Let, for simplicity, $\bar{X}:=X-UU^{\top}$, then, we have

[TABLE]

(b) Consider $\Phi_{*}(Y,V)$ defined in (114), we have

[TABLE]

Proof.

For the proof we use the fact that if $g=-\textup{ln}(\det(X))$, then $g^{\prime}(X)[H]=\textup{Tr}(X^{-1}H)$. Also note that if we define

[TABLE]

then

[TABLE]

and similarly for $\Phi_{*}(Y,V)$. We do not provide all the details, but we show how the proof works. For example, let us define

[TABLE]

and we want to calculate $f^{\prime}(0)$. We have

[TABLE]

∎

Note that all the above formulas for the derivatives are in matrix form. Let us explain briefly how to convert them to the vector form for the code. We explain it for the derivatives of $\Phi(X,U)$ and the rest are similar. From (B.1) we have

[TABLE]

Hence, if $g$ is the gradient of $\Phi(X,U)$ in the vector form, we have

[TABLE]

The second derivatives are trickier. Assume that for example we want the vector form $h$ for $\Phi^{\prime\prime}(X,U)[(d_{X},d_{U})]$. By using (B.1) we can easily get each entry of $h$; consider the identity matrix of size $m^{2}+mn$. If we choose $(\bar{d}_{X},\bar{d}_{U})$ to represent the $j$th column of this identity matrix, we get $h(j)$. Practically, this can be done by a $for$ loop, which is not efficient. What we did in the code is to implement this using matrix multiplication.

Appendix C Quantum entropy and quantum relative entropy

The s.c. barrier DDS uses for quantum entropy is given in (210). To derive the LF conjugate, we solve the optimization problem in (211); let us write the first order optimality conditions for (211). Section 3.3 of the book [18] is about the derivation of matrix-values functions. For the first derivative, we have the following theorem:

Theorem C.1**.**

Let $X$ and $H$ be self-adjoint matrices and $f:(a,b)\mapsto\mathbb{R}$ be a continuously differentiable function defined on an interval. Assume that the eigenvalues of $X+\alpha H$ are in $(a,b)$ for an interval around $\alpha_{0}\in\mathbb{R}$. Then,

[TABLE]

The first-order optimality conditions for (211) can be written as

[TABLE]

If we substitute the first equation in the second one, we get

[TABLE]

This equation implies that $Y$ and $X$ are simultaneously diagonalizable and if we have $Y=U\textup{Diag}(\lambda_{1}(Y),\ldots,\lambda_{n}(Y))$, then we have $X=U\textup{Diag}(\lambda_{1}(X),\ldots,\lambda_{n}(X))$ and so

[TABLE]

Here, we focus on the case that $f(x)=x\textup{ln}(x)$. This matrix function is related to quantum relative entropy and Von-Neumann entropy optimization problems (see [7] for a review of the applications). In this case, we can use results for type 3 univariate function in Table 2 and use the $\theta$ function we defined in (344). The LF conjugate of (211) is given in the following lemma:

Lemma C.1**.**

Assume that $f(x)=x\textup{ln}(x)$. For a given $\eta<0$ and a symmetric matrix $Y\in\mathbb{S}^{n}$, the function defined in (211) becomes

[TABLE]

where $\Theta:=\Theta(\frac{1}{\eta}Y+(1-\textup{ln}(-\eta))I)$. $\Theta$ is the matrix extension of $\theta$ defined in (344) as described in Section 9.

Proof.

Assume that for a given $(\eta,Y)$, $(t,X)$ is the optimal solution for (211). If we use theorem C.1, we have $f^{\prime}(X)=I+\textup{ln}(X)$. By substituting this in the first order optimality condition (332) we get

[TABLE]

By the first equation in (C) and using (335) we can write

[TABLE]

If we substitute $t$ and $X$ in (211), we get the result. ∎

To implement our primal-dual techniques, we need the gradient and Hessian of $\Phi(t,X)$ and $\Phi_{*}(\eta,Y)$. We already saw in Theorem C.1 how to calculate the gradient. The following theorem gives us a tool to calculate the Hessian.

Theorem C.2** ([18]-Theorem 3.25).**

Assume that $f:(a,b)\mapsto\mathbb{R}$ is a $\mathcal{C}^{1}$-function and $T=\textup{Diag}(t_{1},\ldots,t_{n})$ with $t_{i}\in(a,b),\ i\in\{1,\ldots,n\}$. Then, for a Hermitian matrix $H$, we have

[TABLE]

where $\odot$ is the Hadamard product and $T_{f}$ is the divided difference matrix:

[TABLE]

$T$ is diagonal in the statement of the theorem, which is without loss of generality. Note that by the definition of functional calculus in (208), for a Hermitian matrix $X$ and a unitary matrix $U$, we have

[TABLE]

Therefore, for a matrix $T=U\textup{Diag}(t_{1},\ldots,t_{n})U^{*}$, we can update (337)

[TABLE]

where we extend the definition of $T_{f}$ in (340) to non-diagonal matrices. Now we can use Theorems C.2 and C.1 to calculate the Hessian of a matrix function.

Corollary C.1**.**

Let $X$, $H$, and $\tilde{H}$ be self-adjoint matrices and $f:(a,b)\mapsto\mathbb{R}$ be a continuously differentiable function defined on an interval. Assume that the eigenvalues of $X+tH$ and $X+t\tilde{H}$ are in $(a,b)$ for an interval around $t=0$. Assume that $X=U\textup{Diag}(\lambda_{1},\ldots,\lambda_{n})U^{*}$. Then,

[TABLE]

Let us calculate the gradient and Hessian for our functions for $\phi(x)=x\textup{ln}(x)$. Let $X=U\textup{Diag}(\lambda_{1},\ldots,\lambda_{n})U^{*}$ in the following.

[TABLE]

For the second derivative, we can use the fact that

[TABLE]

Using this formula, we have ($\zeta:=\frac{1}{t-\textup{Tr}(X\textup{ln}X)}$)

[TABLE]

Now let us compute the gradient and Hessian for the conjugate function. Let $Y=U\lambda(Y)U^{*}$, by using Theorem C.1, the gradient of $\Phi_{*}(\eta,Y)$ is

[TABLE]

where $\bar{\Theta}$ is the matrix extension of $\left(\theta^{\prime}-\frac{\theta^{\prime}}{\theta^{2}}\right)$. For the second derivative, let us first define $\bar{Y}$ as in (340):

[TABLE]

In the expression of the Hessian, we use $\bar{\bar{\Theta}}$ as the matrix extension of $\left(\theta^{\prime\prime}-\frac{\theta^{\prime\prime}\theta-2(\theta^{\prime})^{2}}{\theta^{3}}\right)$. Then, we have

[TABLE]

Appendix D Univariate convex functions

DDS 2.0 uses the epigraph of six univariate convex function, given in Table 2, in the direct sum format to handle important constraints such as geometric and entropy programming. In these section, we show how to calculate the LF conjugate for 2-variable s.c. barriers. Some of these function are implicitly evaluated, but in a computationally efficient way. In the second part of this section, we show that gradient and Hessian of these functions can also be evaluated efficiently.

D.1. Legendre-Fenchel conjugates of univariate convex functions

The LF conjugates for the first three functions are given in Table 9.

Finding the LF conjugates for the first two functions can be handled with easy calculus. In the third row, $\theta(r)$, defined in [25], is the unique solution of

[TABLE]

It is easy to check by implicit differentiation that

[TABLE]

We can calculate $\theta(r)$ with accuracy $10^{-15}$ in few steps with the following Newton iterations:

[TABLE]

The s.c. barrier DDS uses for the set $|z|^{p}\leq t,\ p\geq 1$, is $\Phi(z,t)=-\textup{ln}(t^{\frac{2}{p}}-z^{2})-2\textup{ln}(t)$. To find the LF conjugate, we need to solve the following optimization problem:

[TABLE]

By writing the first order optimality conditions we have:

[TABLE]

By doing some algebra, we can see that $z$ and $t$ satisfy:

[TABLE]

Let us define $z(y,\eta)$ as the solution of the first equation in (D.1). For each pair $(y,\eta)$, we can calculate $z(y,\eta)$ by few iterations of Newton method. Then, the first and second derivative can be calculated in terms of $z(y,\eta)$. In DDS, we have two functions for these derivatives.

p1_TD(y,eta,p) % returns z p1_TD_der(y,eta,p) % returns [z_y z_eta z_y,y z_y,eta z_eta,eta]

For the set defined by $-z^{p}\leq t,\ 0\leq p\leq 1,z>0$, the corresponding s.c. barrier is $\Phi(z,t)=-\textup{ln}(z^{p}+t)-\textup{ln}(z)$. To calculate the LF conjugate, we need to solve the following optimization problem:

[TABLE]

By writing the first order optimality conditions we have:

[TABLE]

By doing some algebra, we can see that $z$ satisfies:

[TABLE]

Similar to the previous case, let us define $z(y,\eta)$ as the solution of the first equation in (351). For each pair $(y,\eta)$, we can calculate $z(y,\eta)$ by few iterations of Newton method. Then, the first and second derivatives can be calculated in terms of $z(y,\eta)$. The important point is that when we calculate $z(y,\eta)$, then the derivatives can be calculated by explicit formulas. In our code, we have two functions

p2_TD(y,eta,p) % returns z p2_TD_der(y,eta,p) % returns [z_y z_eta z_y,y z_y,eta z_eta,eta]

The inputs to the above functions can be vectors. Table 10 is the continuation of Table 9.

For the set defined by $\frac{1}{z}\leq t,z>0$, the corresponding s.c. barrier is $\Phi(z,t)=-\textup{ln}(zt-1)$. To calculate the LF conjugate, we need to solve the following optimization problem:

[TABLE]

At the optimal solution, we must have

[TABLE]

Since we have $z,t>0$, then we must have $y,\eta<0$ at a dual feasible point. By solving these systems we get

[TABLE]

D.2. Gradient and Hessian

In this subsection, we show how to calculate the gradient and Hessian for the s.c. barriers and their LF conjugates. Specifically for the implicit function, we show that the derivatives can also be calculated efficiently.

First consider the three pairs of functions in Table 9. Here are the explicit formulas for the first and second derivatives:

[TABLE]

For the primal function we have

[TABLE]

and for the dual function we have

[TABLE]

[TABLE]

For the primal function we have

[TABLE]

and for the dual function we have

[TABLE]

[TABLE]

For the primal function we have

[TABLE]

For the dual function, since the argument of the function $\theta(\cdot)$ is always $1+\frac{y}{\eta}-\textup{ln}(-\eta)$, we ignore that in the following formulas and use $\theta$, $\theta^{\prime}$, and $\theta^{\prime\prime}$ for the function and its derivative.

[TABLE]

where

[TABLE]

[TABLE]

where $z(y,\eta)$ is the solution of

[TABLE]

For simplicity, we drop the arguments of $z(y,\eta)$ and denote it as $z$. We denote the first derivatives with respect to $y$ and $\eta$ as $z^{\prime}_{y}$ and $z^{\prime}_{\eta}$, respectively. Similarly, we use $z^{\prime\prime}_{yy}$, $z^{\prime\prime}_{\eta y}$, and $z^{\prime\prime}_{\eta\eta}$ for the second derivatives. We have

[TABLE]

For calculating the second derivatives of $\Phi_{*}$, we need the derivatives of $B$:

[TABLE]

Then we have

[TABLE]

[TABLE]

[TABLE]

By the above definitions of $S$, $M$, and $T$, we have

[TABLE]

The first and second derivatives of $\Phi$ are calculated as follows:

[TABLE]

The first and second derivatives of $\Phi_{*}$ are messier. For the first derivative we have

[TABLE]

For calculating the second derivative, we also need the second derivatives of $B$:

[TABLE]

Using the second derivatives of $B$, we have

[TABLE]

where

[TABLE]

$f_{21}=f_{12}$ and $f_{22}$ have similar formulations.

[TABLE]

where $z$ is the solution of

[TABLE]

Similar to the previous case, for simplicity, we drop the arguments of $z(y,\eta)$ and denote it as $z$. We denote the first derivatives with respect to $y$ and $\eta$ as $z^{\prime}_{y}$ and $z^{\prime}_{\eta}$, respectively. By implicit differentiation, we have

[TABLE]

For the second derivatives of $z$, by using

[TABLE]

we have

[TABLE]

The first and second derivatives of $\Phi$ are calculated as follows:

[TABLE]

The first derivative of $\Phi_{*}$ is equal to

[TABLE]

and the second derivatives of $\Phi_{*}$ is equal to

[TABLE]

where

[TABLE]

[TABLE]

For the primal function we have

[TABLE]

and for the dual function we have

[TABLE]

Appendix E A s.c. barrier for vector relative entropy

In section, we prove that the function

[TABLE]

is a $(2\ell+1)$-LH-s.c. barrier for the epigraph of vector relative entropy function $f:\mathbb{R}_{++}^{\ell}\oplus\mathbb{R}_{++}^{\ell}\rightarrow\mathbb{R}$ defined as

[TABLE]

First note that $\Phi$ is $(2\ell+1)$-LH, so we just need to prove self-concordance.

The proof is based on the compatibility results of [27]. We first consider the case $\ell=1$ and then generalize it. Let us define the map $A:\mathbb{R}_{++}^{2}\rightarrow\mathbb{R}$ as

[TABLE]

For a vector $d:=(d_{u},d_{z})$, we can verify

[TABLE]

(E) implies that $A$ is concave with respect to $\mathbb{R}_{+}$. We claim that $A$ is $(\mathbb{R}_{+},1)$-compatible with the barrier $-\textup{ln}(u)-\textup{ln}(z)$ ([27]- Definition 5.1.2). For this, we need to show

[TABLE]

By using (E) and canceling out the common terms from both sides, (395) reduces to

[TABLE]

We can assume that the LHS is nonnegative, then by taking the square of both side and reordering, we get the obvious inequality

[TABLE]

Therefore $A(u,z)$ is $(\mathbb{R}_{+},1)$-compatible with the barrier $-\textup{ln}(u)-\textup{ln}(z)$. Also note that its summation with a linear term $t+A(u,z)$ is also $(\mathbb{R}_{+},1)$-compatible with the barrier $-\textup{ln}(z)-\textup{ln}(u)$. Hence, $-\textup{ln}(t+A)-\textup{ln}(u)-\textup{ln}(z)$ is a 3-s.c. barrier by [27]-Proposition 5.1.7.

For the general case, consider the map $\bar{A}:[\mathbb{R}_{++}^{2}]^{\ell}\rightarrow\mathbb{R}$ as

[TABLE]

For a vector $d$ of proper size, we have

[TABLE]

We claim that $\bar{A}$ is $(\mathbb{R}_{+},1)$-compatible with the barrier $-\sum_{i}\textup{ln}(u_{i})-\sum_{i}\textup{ln}(z_{i})$. First note that $\bar{A}$ is concave with respect to $\mathbb{R}_{+}$. We need to prove a similar inequality as (395) for $\bar{A}$. Clearly we have

[TABLE]

Using inequality (395) and (400) for all $j$ and adding them together yields the inequality we want for $\bar{A}$. Therefore, $\bar{A}$ is $(\mathbb{R}_{+},1)$-compatible with the barrier $-\sum_{i}\textup{ln}(u_{i})-\sum_{i}\textup{ln}(z_{i})$, and by [27]-Proposition 5.1.7, (392) is a $(2\ell+1)$-s.c. barrier.

Appendix F Comparison with some related solvers and modeling systems

In this section, we take a look at the input format for some other well-known solvers and modeling systems. [22] is a survey by Mittelmann about solvers for conic optimization, which gives an overview of the major codes available for the solution of linear semidefinite (SDP) and second-order cone (SOCP) programs. Many of these codes also solve linear programs (LP). We mention the leaders MOSEK, SDPT3, and SeDuMi from the list. We also look at CVX, a very user-friendly interface for convex optimization. CVX is not a solver, but is a modeling system that (by following some rules) detects if a given problem is convex and remodels it as a suitable input for solvers such as SeDuMi.

F.1. MOSEK [23]

MOSEK is a leading commercial solver for not just optimization over symmetric cones, but also many other convex optimization problems. The most recent version, MOSEK 9.0, for this state-of-the-art convex optimization software handles, in a primal-dual framework, many convex cone constraints which arise in applications [10]. There are different options for the using platform that can seen in MOSEK’s website [23].

F.2. SDPT3 [36, 33]

SDPT3 is a MATLAB package for optimization over symmetric cones, and it solves a conic optimization problem in the equality form as

[TABLE]

where our cone $K$ can be a direct sum of nonnegative rays (leading to LP problems), second-order cones or semidefinite cones. The Input for SDPT3 is given in the cell array structure of MATLAB. The command to solve SDPT3 is of he form

[obj,X,y,Z,info,runhist] = sqlp(blk,At,C,b,OPTIONS,X0,y0,Z0).

The input data is given in different blocks, where for the $k$th block, blk{k,1} specifies the type of the constraint. Letters ’l’, ’q’, and ’s’ are representing linear, quadratic, and semidefinite constraints. In the $k$th block, At{k}, C{k}, … contain the part of the input related to this block.

F.3. SeDuMi [32]

SeDuMi is also a MATLAB package for optimization over symmetric cones in the format of (404). For SeDuMi, we give as the input $A$, $b$ and $c$ and a structure array $K$. The vector of variables has a “direct sum” structure. In other words, the set of variables is the direct sum of free, linear, quadratic, or semidefinite variables. The fields of the structure array $K$ contain the number of constraints we have from each type and their sizes. SeDuMi can be called in MATLAB by the command

[x,y] = sedumi(A,b,c,K);

and the variables are distinguished by $K$ as follows:

1. (1)

$K.f$ is the number of free variables, i.e., in the variable vector $x$, x(1:K.f) are free variables. 2. (2)

$K.l$ is the number of nonnegative variables. 3. (3)

$K.q$ lists the dimension of Lorentz constraints. 4. (4)

$K.s$ lists the dimensions of positive semidefinite constraints.

For example, if K.l=10, K.q=[3 7] and K.s=[4 3], then x(1:10) are non-negative. Then we have x(11) >= norm(x(12:13)), x(14) >= norm(x(15:20)), and mat(x(21:36),4) and mat(x(37:45),3) are positive semidefinite. To insert our problem into SeDuMi, we have to write it in the format of (404) . We also have the choice to solve the dual problem because all of the above cones are self-dual.

F.4. CVX [16]

CVX is an interface that is more user-friendly than solvers like SeDuMi. It provides many options for giving the problem as an input, and then translates them to an eligible format for a solver such as SeDuMi. We can insert our problem constraint-by-constraint into CVX, but they must follow a protocol called Disciplined convex programming (DCP). DCP has a rule-set that the user has to follow, which allows CVX to verify that the problem is convex and convert it to a solvable form. For example, we can write a $<=$ constraint only when the left side is convex and the right side is concave, and to do that, we can use a large class of functions from the library of CVX.