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

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

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.

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);

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:

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

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

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.

Suppose we want to add the following constraints to DDS:

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

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

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:

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.

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

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:

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:

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:

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

Assume that we have constraints of the form

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:

with LF conjugate

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

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

DDS accepts constraints of the form

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

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

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

  $\displaystyle\texttt{cons\{k,2\}\{j,1\}}=\left[\begin{array}[]{ccc}2&l_{2}^{j}\\ 3&l_{3}^{j}\end{array}\right].$

  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

  $\displaystyle\texttt{cons\{k,2\}\{j,2\}}=[\alpha_{1}^{j,2},\alpha_{2}^{j,2},\alpha_{1}^{j,3},\cdots].$

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

Let us see an example: