HDSDP: Software for Semidefinite Programming

Dual-scaling method [4] works under three conditions. 1) the data $\left\{\mathbf{A}_{i}\right\}$ are linearly independent. 2) Both $(P)$ and $(D)$ admit an interior point solution. 3) an interior dual feasible solution $\left(\mathbf{y}^{0},\mathbf{S}^{0}\right)\in\mathcal{F}^{0}(D)$ is known. The first two conditions imply strong duality and the existence of an optimal primal-dual solution pair $\left(\mathbf{X}^{\ast},\mathbf{y}^{\ast},\mathbf{S}^{\ast}\right)$ satisfying complementarity $\left\langle\mathbf{X}^{\ast},\mathbf{S}^{\ast}\right\rangle=0$. Also, the central path $\mathcal{C}(\mu):=\left\{\left(\mathbf{X},\mathbf{y},\mathbf{S}\right)\in\mathcal{F}^{0}(P)\times\mathcal{F}^{0}(D):\mathbf{X}\mathbf{S}=\mu\mathbf{I}\right\}$, which serves as a foundation of path-following approaches, is well-defined. Given a centrality parameter $\mu$, dual method starts from $\left(\mathbf{y},\mathbf{S}\right)\in\mathcal{F}^{0}(D)$ and takes Newton’s step towards the solution of the perturbed KKT system $\mathcal{A}\mathbf{X}=\mathbf{b},\mathcal{A}^{\ast}\mathbf{y}+\mathbf{S}=\mathbf{C}$ and $\mathbf{X}\mathbf{S}=\mu\mathbf{I}$

where the last relation linearizes $\mathbf{X}=\mu\mathbf{S}^{-1}$ instead of $\mathbf{X}\mathbf{S}=\mu\mathbf{I}$ and $\mathbf{S}^{-1}$ is known as a scaling matrix. By geometrically driving $\mu$ to 0, dual-scaling eliminates (primal) infeasibility, approaches optimality, and finally solves the problem to some $\varepsilon$-optimal solution $(\mathbf{y_{\varepsilon}},\mathbf{S_{\varepsilon}})\in\mathcal{F}^{0}(D)$ such that $\langle\mathbf{b},\mathbf{y_{\varepsilon}}\rangle\geq\langle\mathbf{b},\mathbf{y^{*}}\rangle-\varepsilon$. Theory of dual potential reduction shows an $\varepsilon$-optimal solution can be obtained in $\mathcal{O}(\log(1/\varepsilon))$ iterations ignoring dimension dependent constants.

An important feature of dual-scaling is that $\mathbf{X}$ and $\Delta\mathbf{X}$ vanish in the Schur complement of (3.1)

and the dual algorithm thereby avoids explicit reference to the primal variable $\mathbf{X}$. For brevity we denote the left-hand side matrix in (6) by $\mathbf{M}$. If $\left\{\mathbf{A}_{i}\right\},\mathbf{C}$ are sparse, any feasible dual slack $\mathbf{S}=\mathbf{C}-\mathcal{A}^{\ast}\mathbf{y}$ inherits the aggregated sparsity pattern from the data and it is computationally cheaper to iterate in the dual space. Another feature of dual-scaling is the availability of primal solution by solving a projection subproblem at the end of the algorithm [2]. The aforementioned properties characterize the behavior of the dual-scaling method.

However, the desirable theoretical properties of the dual method is not free. First, an initial dual feasible solution is needed, while obtaining such a solution is often as difficult as solving the original problem. Second, due to a lack of information from primal space, dual-scaling method has to be properly guided to avoid deviating from the central path. To overcome the aforementioned difficulties, DSDP introduces artificial variables with big-$M$ penalty to ensure a nonempty interior and a trivial dual feasible solution. Moreover, a potential function is introduced to guide the dual iterations. This works well in practice and enhance performance of DSDP. We refer the interested readers to [3, 2] for more implementation details.

Although DSDP proves efficient in real practice, the big-$M$ method requires prior estimation of the optimal solution to avoid losing optimality. Also, a large penalty often leads to numerical difficulties and mis-classification of infeasible problems when the problem is ill-conditioned. Therefore it is natural to seek better alternatives to the big-$M$ method, and the self-dual embedding is an ideal candidate.

Given a standard form SDP, its homogeneous self-dual model is a skew symmetric system containing the original problem data, whose non-trivial interior point solution certificates primal-dual feasibility. HDSDP adopts the following simplified embedding [20].

where $\kappa,\tau$ are homogenizing variables for infeasibility detection. The central path of barrier parameter $\mu$ is given by

Here $\left(\mathbf{y},\mathbf{S},\tau\right)$ are jointly considered as dual variables. Given an (infeasible) dual point $\left(\mathbf{y},\mathbf{S},\tau\right)$ such that $-\mathcal{A}^{\ast}\mathbf{y}+\mathbf{C}\tau-\mathbf{S}=\mathbf{R}$, HDSDP selects a damping factor $\gamma\in[0,1]$ and takes Newton’s step towards

where, similar to DSDP, we modify (8) and linearize $\mathbf{X}=\mu\mathbf{S}^{-1},\kappa=\mu\tau^{-1}$. We use this damping factor $\gamma$ and the barrier parameter $\mu$ to manage a trade-off between dual infeasibility, centrality and optimality. If we set $\gamma=0$, then the Newton’s direction $\left(\Delta\mathbf{y},\Delta\tau\right)$ is computed through the following Schur complement.

In practice, HDSDP solves $\Delta\mathbf{y}_{1}:=\mathbf{M}^{-1}\mathbf{b},\Delta\mathbf{y}_{2}:=\mathbf{M}^{-1}\mathcal{A}\mathbf{S}^{-1},\Delta\mathbf{y}_{3}:=\mathbf{M}^{-1}\mathcal{A}\mathbf{S}^{-1}\mathbf{R}\mathbf{S}^{-1},\Delta\mathbf{y}_{4}=\mathbf{M}^{-1}\mathcal{A}\mathbf{S}^{-1}\mathbf{C}\mathbf{S}^{-1}$, plugs the solution into $\Delta\mathbf{y}=\frac{\tau+\Delta\tau}{\mu}\Delta\mathbf{y}_{1}-\Delta\mathbf{y}_{2}+\Delta\mathbf{y}_{3}+\Delta\mathbf{y}_{4}\Delta\tau$ to get $\Delta\tau$, and finally assembles $\Delta\mathbf{y}$. With the above relations, $\mathbf{X}(\mu):=\mu\mathbf{S}^{-1}(\mathbf{S}-\mathbf{R}+\mathcal{A}^{\ast}\Delta\mathbf{y}-\mathbf{C}\Delta\tau)\mathbf{S}^{-1}$ satisfies $\mathcal{A}\mathbf{X}(\mu)-\mathbf{b}(\tau+\Delta\tau)=\textbf{0}$ and $\mathbf{X}(\mu)\succeq\textbf{0}$ if and only if the backward Newton step $-\mathcal{A}^{\ast}(\mathbf{y}-\Delta\mathbf{y})+\mathbf{C}(\tau-\Delta\tau)-\mathbf{R}\succeq\textbf{0}$. When $-\mathcal{A}^{\ast}(\mathbf{y}-\Delta\mathbf{y})+\mathbf{C}(\tau-\Delta\tau)-\mathbf{R}$ is positive definite, a primal upper-bound follows from simple algebraic manipulation

and dividing both sides by $\tau$. Alternatively, HDSDP extracts a primal objective bound from the projection problem

whose optimal solution is $\mathbf{X}^{\prime}(\mu):=\mu\mathbf{S}^{-1}(\mathbf{C}\tau-\mathcal{A}^{\ast}(\mathbf{y}-\Delta^{\prime}\mathbf{y})-\mathbf{R})\mathbf{S}^{-1}$. Here $\Delta\mathbf{y}^{\prime}=\frac{\tau}{\mu}\Delta\mathbf{y}_{1}-\Delta\mathbf{y}_{2}$ and

Using the Newton’s direction $\Delta\mathbf{y}$, HDSDP applies a simple ratio test

through a Lanczos procedure [24]. But to determine the aforementioned damping factor $\gamma$, we resort to the following adaptive heuristic: assuming that $\mu\rightarrow\infty,\tau=1$ and fixing $\Delta\tau=0$, the dual update can be decomposed by

where the first term $\mathbf{S}+\alpha\mathcal{A}^{\ast}\Delta\mathbf{y}_{2}$ is independent of $\gamma$ and improves centrality of the iteration. HDSDP conducts a Lanczos line-search to find $\alpha_{c}$ such that the logarithmic barrier function is improved

Then given $\alpha=\alpha_{c}$, by a second Lanczos line-search we find the maximum possible $\gamma\leq 1$ such that $\mathbf{S}+\alpha_{c}\Delta\mathbf{S}=\mathbf{S}+\alpha_{c}\mathcal{A}^{\ast}\Delta\mathbf{y}_{2}+\alpha_{c}\gamma(\mathbf{R}-\mathcal{A}^{\ast}\Delta\mathbf{y}_{3})\succeq 0$. Finally, a full Newton’s direction is determined by $\gamma$ and a third Lanczos procedure finally carries out the ratio test (10). The intuition of the above heuristic is to eliminate infeasibility under centrality restrictions, so that the iterations stay away from the boundary of the cone. After the ratio test, HDSDP updates $\mathbf{y}\leftarrow\mathbf{y}+0.95\alpha\Delta\mathbf{y}$ and $\tau\leftarrow\tau+0.95\alpha\Delta\tau$ to ensure the feasibility of the next iteration. To make further use of the Schur complement $\mathbf{M}$, the above procedure is repeated several times in an iteration with $\mathbf{M}$ unchanged.

In HDSDP, the barrier parameter $\mu$ also significantly affects the algorithm. At the end of each iteration, HDSDP updates the barrier parameter $\mu$ by $(z-\mathbf{b}^{\top}\mathbf{y}+\theta\left\|\mathbf{R}\right\|_{F})/\rho n$, where $z$ is the best primal bound so far and $\rho,\theta$ are pre-defined parameters. By default $\rho=4$ and $\theta=10^{8}$. Heuristics also adjust $\mu$ based on $\alpha_{c}$ and $\gamma$. To get the best of both worlds, HDSDP implements the same dual-scaling algorithm as in DSDP5.8. If dual infeasibility satisfies $\left\|\mathcal{A}^{\ast}\mathbf{y}+\mathbf{S}-\mathbf{C}\right\|_{F}\leq\varepsilon\tau$ and $\mu$ is sufficiently small, HDSDP fixes $\tau=1$, re-starts with $(\mathbf{y}/\tau,\mathbf{S}/\tau)$ and instead applies dual potential function to guide convergence. To sum up, HDSDP implements strategies and computational tricks tailored for the embedding and can switch to DSDP5.8 once a dual feasible solution is available.

One new feature of HDSDP is a special pre-solving module designed for SDPs. It inherits techniques from DSDP5.8 and adds new strategies to work jointly with the optimization module. After the pre-solver is invoked, it first goes through $\{\mathbf{A}_{i}\}$ two rounds to detect the possible low-rank structure. The first round uses Gaussian elimination for rank-one structure and the second round applies eigenvalue decomposition from Lapack. Two exceptions are when the data matrix is too dense or sparse. Unlike in DSDP5.8 where eigen-decomposition is mandatory, matrices that are too dense to be eigen-decomposed efficiently will be skipped in HDSDP; if a matrix has very few entries, then a strategy from DSDP5.8 is applied: 1) a permutation gathers the non-zeros to a much smaller dense block. 2) dense eigen routines from Lapack applies. 3) the inverse permutation recovers the decomposition.

After detecting the hidden low-rank structures, the pre-solver moves on to collecting the sparsity and rank information of the matrices. The information will be kept for the KKT analysis to be described in Section 5.3.

Finally, the pre-solver scales down the large objective coefficients by their Frobenius norm and goes on to identify the following seven structures. 1). Implied trace: constraints imply $\operatorname{tr}(\mathbf{X})=\theta$. 2). Implied dual bound: constraints imply $\mathbf{l}\leq\mathbf{y}\leq\mathbf{u}$. 3). Empty primal interior: constraints imply $\operatorname{tr}(\mathbf{X}\mathbf{a}\mathbf{a}^{\top})\approx 0$. 4). Empty dual interior: constraints imply $\mathcal{A}^{*}\mathbf{y}=\mathbf{C}$. 5). Feasibility problem: $\mathbf{C}=\textbf{0}$. 6). Dense problem: whether all the constraint matrices are totally dense. 7). Multi-block problem: whether there are many SDP cones. For each of the cases the solver adjusts its internal parameters to enhance its numerical stability and convergence.

Underlie the procedure control of HDSDP is a two-phase method which integrates the embedding technique (Phase A) and dual-scaling (Phase B).  Phase A targets feasibility certificate, while Phase B aims to efficiently drive a dual-feasible solution to optimality.

The two phases share the same backend computation routines but are associated with different goals and strategies. HDSDP decides which strategy to use based on the behavior of algorithm iterations.

The most computationally intensive part in HDSDP is in the Schur complement matrix $\mathbf{M}$. In real-life applications, the SDP objective $\mathbf{C}$ often does not share the same sparsity or low-rank property as $\{\mathbf{A}_{i}\}$. Therefore, the additional computation from $\mathcal{A}\mathbf{S}^{-1}\mathbf{C}\mathbf{S}^{-1}$ and $\langle\mathbf{C},\mathbf{S}^{-1}\mathbf{C}\mathbf{S}^{-1}\rangle$ demands more efficient techniques to set up $\mathbf{M}$.

Efficiency of setting up $\mathbf{M}$ depends mainly on two aspects: 1). the structure of the cones participating in forming the Schur complement, and 2) the structure of the coefficients inside each cone. To exploit both structures, KKT solver in HDSDP implements two abstract interfaces, one from conic coefficient to conic operations, the other from cone to $\mathbf{M}$.

The first interface allows HDSDP to set up $\mathbf{M}$ faster, while the second allows the extension of dual-scaling to arbitrary cones. These two interfaces already appeared in DSDP5.8, but unifying them in a conic KKT solver makes the algorithm implementation more compact. Till the date when this manuscript is written, four cones and five types SDP coefficients have been implemented by HDSDP.

The interface from SDP coefficients to conic operations is directly relevant to the set up of $\mathbf{M}$. Using rank and sparsity information from the pre-solver, HDSDP analyzes the Schur system and generates a permutation $\sigma(m)$ of the rows of $\mathbf{M}$. The idea is to permute the most computationally-expensive matrices to the left-most position, and this is a direct extension of [9].

After deciding the permutation, for each row of the permuted system, HDSDP chooses one out of a pool of five candidate techniques. We denote them by M1 to M5, and they are efficient under different coefficient structures: Technique M1 and M2 are inherited from DSDP. They exploit the low-rank structure of the coefficient data and an eigen-decomposition of the data $\mathbf{A}_{i}=\sum_{r=1}^{\operatorname{rank}(\mathbf{A}_{i})}\lambda_{ir}\mathbf{a}_{i,r}\mathbf{a}^{\top}_{i,r}$ needs to be computed at the beginning of the algorithm; Technique M3, M4 and M5 exploit sparsity and need evaluation of $\mathbf{S}^{-1}$ [9] every iteration.

This newly developed KKT solver improves the speed of HDSDP on a broad set of instances. And to our knowledge HDSDP is the first SDP software that simultaneously incorporates all the commonly-used Schur complement strategies.

Aside from the third-party routines, HDSDP itself implements a pre-conditioned conjugate gradient (CG) method with restart to solve the Schur complement system. The maximum number of iterations is chosen around $50/m$ and is heuristically adjusted. Either the diagonal of $\mathbf{M}$ or its Cholesky decomposition is chosen as pre-conditioner and after a Cholesky pre-conditioner is computed, HDSDP reuses it for the consecutive solves till a heuristic determines that the current pre-conditioner is outdated. When the algorithm approaches optimality, $\mathbf{M}$ might become ill-conditioned and HDSDP switches to LDL decomposition in case Cholesky fails.

We configure the experiment as follows

1. 1.

   (Testing platform). The tests are run on an intel i11700K PC with 128GB memory and 12 threads.We choose HDSDP, DSDP5.8 and COPT6.5 are selected as benchmark softwares.

2. 2.

   (Compiler and third-party dependency). Both HDSDP and DSDP5.8 are compiled using icc with -O2 optimization and linked with threaded version intel MKL. Executable of COPT is directly obtained from its official website.

3. 3.

   (Threading). We set the number of MKL threads to be 12 for DSDP5.8; the Threads parameter of COPT is also set to 12. To enhance reproducibility HDSDP uses 8 threads.

4. 4.

   (Dataset). Datasets are chosen from three sources: 1). SDP benchmark datasets from Hans Mittelmann’s website [18]; 2). SDP benchmark datasets from SDPLIB [6]. 3). Optimal diagonal pre-conditioner SDPs generated according to [22].

5. 5.

   (Algorithm). Default methods in HDSDP, DSDP5.8, and the primal-dual interior point method implemented in COPT are compared.

6. 6.

   (Tolerance). All the feasibility and optimality tolerances are set to $5\times 10^{-6}$.

7. 7.

   (Solution status). We adopt the broadly accepted DIMACs error to determine if a solution is qualified. According to [18], if any of the DIMACS errors exceeds $10^{-2}$, the solution is considered invalid.

8. 8.

   (Performance Metric). We use shifted geometric mean to compare the overall speed between different solvers.

The SDP relaxation of the max-cut problem is represented by

Let $\mathbf{e}_{i}$ be the $i$-th column of the identity matrix and the constraint $\operatorname{diag}\left(\mathbf{X}\right)=\mathbf{e}$ is decomposed into $\left\langle\mathbf{X},\mathbf{e}_{i}\mathbf{e}_{i}^{\top}\right\rangle=1,i=1,\ldots,n$. Note that $\mathbf{e}_{i}\mathbf{e}_{i}^{\top}$ is rank-one and has only one non-zero entry, M2 and M5 can greatly reduce the computation of the Schur matrix.

Computational experience suggests that on large-scale sparse max-cut instances, HDSDP is more than $5$ times faster than DSDP5.8. Two exceptions are G60mc and G60_mb, where the aggregated sparsity pattern of the dual matrix is lost due to the dense objective coefficient, and thereby the dense MKL routines take up most of the computation.

The SDP relaxation of the graph partitioning problem is given by

where 1 denotes the all-one vector and $k,\beta$ are the problem parameters. Although the dual $\mathbf{S}$ no longer enjoys sparsity, the low-rank structure is still available to accelerate convergence.

On graph partitioning instances, we see that HDSDP has comparable performance to DSDP but is more robust on some of the problems.

So far HDSDP is tested and tuned over a broad set of benchmarks including SDPLIB [6] and Hans Mittelmann’s sparse SDP benchmark [18]. Using geometric mean as the metric, compared to DSDP5.8, HDSDP achieves a speedup of 21% on SDPLIB and around 17% speedup on Hans Mittelmann’s benchmark dataset. Below we list some example benchmark datasets of nice structure for HDSDP to exploit.