SDPNAL$+$: A Majorized Semismooth Newton-CG Augmented Lagrangian Method for Semidefinite Programming with Nonnegative Constraints

Introduction

Let ${\cal K}$ be a pointed closed convex cone whose interior $int({\cal K})\neq\emptyset$ and ${\cal P}$ be a polyhedral convex cone in a finite-dimensional Euclidean space ${\cal X}$ such that ${\cal K}\cap{\cal P}$ is non-empty. For any cone ${\cal C}\subseteq{\cal X}$, we denote the dual cone of ${\cal C}$ by ${\cal C}^{*}$. For any closed convex set ${\cal C}\subseteq{\cal X}$, we denote the metric projection of ${\cal X}$ onto ${\cal C}$ by $\Pi_{\cal C}(\cdot)$ and the tangent cone of ${\cal C}$ at $X\in{\cal C}$ by ${\cal T}_{{\cal C}}(X)$, respectively. We will make extensive use of the Moreau decomposition theorem in , which states that $X=\Pi_{\cal C}(X)-\Pi_{{\cal C}^{*}}(-X)$ for any $X\in{\cal X}$ and any closed convex cone ${\cal C}\subseteq{\cal X}$. Let ${\cal S}^{n}$ be the space of $n\times n$ real symmetric matrices and ${\cal S}^{n}_{+}$ be the cone of positive semidefinite matrices in ${\cal S}^{n}$. In this paper, we focus on the case where ${\cal X}={\cal S}^{n}$, ${\cal K}={\cal K}^{*}={\cal S}^{n}_{+}$. We are particularly interested in the case where ${\cal P}={\cal S}^{n}_{\geq 0}$, the cone of $n\times n$ real symmetric matrices whose elements are all nonnegative, though the algorithm which we will design later is also applicable to other cases. For any matrix $X\in{\cal S}^{n}$, we use $X\succ 0$ to indicate that $X$ is a real symmetric positive definite matrix.

Consider the semidefinite programming (SDP) with an additional polyhedral cone constraint, which we name as SDP$+$:

where $b\in\Re^{m}$ and $C\in{\cal X}$ are given data, ${\cal A}:{\cal X}\rightarrow\Re^{m}$ is a given linear map whose adjoint is denoted as ${\cal A}^{*}$. Note that ${\cal P}={\cal X}$ is allowed in (1), in which case there is no additional polyhedral cone constraint imposed on $X$. We assume that the matrix ${\cal A}{\cal A}^{*}$ is invertible, i.e., ${\cal A}$ is surjective. The dual of (P) is given by

The optimality conditions (KKT conditions) for (P) and (D) can be written as follows:

In order for the KKT conditions (5) to have solutions, throughout this paper we make the following blanket assumption.

(a) For problem (P), there exists a feasible solution $X_{0}\in{\cal S}^{n}_{+}$ such that

(b) For problem (D), there exists a feasible solution $(y_{0},S_{0},Z_{0})\in\Re^{m}\times{\cal S}^{n}_{+}\times{\cal S}^{n}$ such that

It is known from convex analysis (e.g, [2, Corollary 5.3.6]) that under Assumption 1, the strong duality for (P) and (D) holds and the KKT conditions (5) have solutions.

For a given $\sigma>0$, define the augmented Lagrangian function for the dual problem (D) as follows:

where $X\in{\cal X},y\in\Re^{m},S\in{\cal K}^{*},Z\in{\cal P}^{*}$. We can consider the following inexact augmented Lagrangian method to solve (D). Specifically, given $\sigma_{0}>0$, $(y^{0},S^{0},Z^{0})\in\Re^{m}\times{\cal K}^{*}\times{\cal P}^{*}$, perform the following steps at the $(k+1)$-th iteration:

where $\sigma_{k}\in(0,+\infty)$, $k=0,1,\ldots$. For a general discussion on the convergence of the augmented Lagrangian method for solving convex optimization problems and beyond, see .

Note that problem (P) can be reformulated as a standard SDP in the primal form by replacing the constraint $X\in{\cal P}$ with two constraints $X-Y=0$ and $Y\in{\cal P}$. In , SDPNAL introduced by Zhao, Sun and Toh is applied to solve such a reformulated problem. It works quite well for nondegenerate SDPs, especially those without the constraint $X\in{\cal P}$. However, many of the tested SDPs (with the constraint $X\in{\cal P}$) in are degenerate and SDPNAL is unable to solve those problems efficiently. Motivated by our desire to overcome the aforementioned difficulty in solving degenerate SDPs and to improve the performance of SDPNAL, we present here a majorized semismooth Newton-CG augmented Lagrangian method by directly working on (P) instead of its reformulated problem. We call this new method SDPNAL$+$ since it is a much enhanced version of SDPNAL and it is designed for SDP$+$ problems (P).

The remaining parts of this paper are organized as follows. In Section 2, we introduce a majorized semismooth Newton-CG method for solving the inner minimization problems of the augmented Lagrangian method and analyze the convergence for solving these inner problems. Section 3 presents the SDPNAL$+$ dual approach. Section 4 is on numerical issues. There we report numerical results for a variety of SDP$+$ and SDP problems. We make an extensive numerical comparison with two other competitive first order methods based codes: (1) an alternating direction method of multiplier (ADMM) based solver called SDPAD by Wen et al. and (2) a two-easy-block-decomposition hybrid proximal extragradient method called 2EBD-HPE by Monteiro et al. . Numerical results show that SDPNAL$+$ is both fast and robust in achieving accurate solutions.

For the first time, we are able to solve all the $95$ difficult SDP problems arising from the relaxations of quadratic assignment problems (QAPs) tested in SDPNAL to an accuracy of $10^{-6}$ efficiently, while SDPAD and 2EBD-HPE successfully solve 30 and 16 problems, respectively. In addition, SDPNAL$+$ appears to be the only viable method currently available to solve large scale SDPs arising from rank-1 tensor approximation problems constructed by Nie and Wang . The largest rank-1 tensor approximation problem solved is nonsym(21,4), in which its resulting SDP problem has matrix dimension $n=9,261$ and the number of equality constraints $m=12,326,390$. Finally, in order to demonstrate the power of the proposed majorized semismooth Newton-CG procedure, we list the numerical results by only running the convergent ADMM with $3$-block constraints (ADMM+ in short) introduced by Sun et al. . As one may observe, although ADMM+ outperforms both SDPAD and 2EBD-HPE, it can still encounter numerical difficulty in solving some hard problems such as those arising from QAPs to high accuracy. The superior numerical performance of SDPNAL$+$ over solvers based purely on first order methods such as SDPAD and 2EBD-HPE clearly shows the necessity of exploiting second order methods such as the semismooth Newton-CG method in order to solve hard SDP$+$ and SDP problems to high accuracy efficiently. While there has been a recent focus on using first order methods such as those based on ADMM or accelerated proximal gradient methods to solve structured convex optimization problems arising from machine learning and statistics, the extensive numerical results we obtained here for matrix conic programming problems serve to demonstrate that second order methods with good local convergence property are essential, if used wisely, for mitigating the inherent slow local convergence of first order methods, especially on difficult problems.

A Majorized Semismooth Newton-CG Method for Inner Problems

Let $\sigma>0$ and $\widetilde{X}\in{\cal S}^{n}$ be fixed. In this section we will present a majorized semismooth Newton-CG method for solving the following inner problems involved in the augmented Lagrangian method (9a):

Note that problem (10) is the dual of the following problem:

Since the objective function in (11) is strongly concave, (11) has a unique optimal solution. In order for its dual problem (10) to have a bounded solution set, we need the following generalized Slater condition.

There exists a positive definite matrix $X_{0}\in{\cal S}^{n}_{+}\cap relint({\cal P})$ such that

where $relint({\cal P})$ denotes the relative interior of ${\cal P}$.

Note that in (12), ${\cal T}_{{\cal P}}(X_{0})$ is actually a linear subspace of ${\cal S}^{n}$ as $X_{0}$ is assumed to be in the relative interior part of the polyhedral cone ${\cal P}$. When ${\cal P}={\cal S}^{n}$, Assumption 2 is equivalent to saying that

From [16, Theorems $17$ and $18$], we have the following useful lemma.

Suppose that Assumption 2 holds. Then for any $\alpha\in\Re$, the level set $\mathcal{L}_{\alpha}:=\{(y,S,Z)\in\Re^{m}\times{\cal K}^{*}\times{\cal P}^{*}\,\mid\,\phi(y,S,Z)\leq\alpha\}$ is a closed and bounded convex set.

In order to introduce our majorized semismooth Newton-CG method for solving (16), we need to majorize the second part of the objective function in (16) by a convex, but not necessarily strongly convex, quadratic function. Specifically, for given $(y^{l},S^{l})\in\Re^{m}\times{\cal K}^{*}$ and $l\geq 0$, since

where $Z^{l}:=\Pi_{{\cal P}^{*}}(\widehat{C}-{\cal A}^{*}y^{l}-S^{l})$, we know that for $(y,S)\in\Re^{m}\times{\cal S}^{n}$,

where $\widetilde{C}^{l}:=C-Z^{l}$. Thus $\varPsi_{l}$ is a majorization function of $\varPhi$ at $(y^{l},S^{l})$ because $\varPsi_{l}(y^{l},S^{l})=\varPhi(y^{l},S^{l})$ and $\varPsi_{l}(y,S)\geq\varPhi(y,S)$ $\forall(y,S)\in\Re^{m}\times{\cal S}^{n}$. In order to find an optimal solution for problem (16), for $l=0,1,\ldots$, we solve the following problem

Note that we can only solve problem (19a) inexactly by an iterative method. Here we will introduce a semismooth Newton-CG (SNCG) method for solving (19a). Specifically, for fixed $\widetilde{X},\widetilde{C}\in{\cal S}^{n}$, we need to consider the following problem of the form

The objective function in (20) is continuously differentiable and solving (20) is equivalent to solving the following nonsmooth equation:

Since $\Pi_{{\cal K}}(\cdot)$ is strongly semismooth , we can design a SNCG method as in to solve (21), and expect fast superlinear or even quadratic convergence.

where $``\circ"$ denotes the Hadamard product of two matrices and

where $\partial\Pi_{{\cal K}}(\widetilde{X}+\sigma(\mathcal{A}^{*}y-\widetilde{C}))$ is the Clarke subdifferential of $\Pi_{{\cal K}}(\cdot)$ at $\widetilde{X}+\sigma(\mathcal{A}^{*}y-\widetilde{C})$. Note that from , we know that

Now we will introduce the SNCG algorithm for solving (20). Choose $y^{0}\in\Re^{m}$. Then the algorithm can be stated as follows.

Algorithm SNCG: A Semismooth Newton-CG Algorithm (SNCG$(y^{0},\widetilde{X},\sigma)$). Given $\mu\in(0,1/2)$, $\bar{\eta}\in(0,1)$, $\tau\in(0,1]$, $\tau_{1},\tau_{2}\in(0,1)$, and $\delta\in(0,1)$. Perform the $j$th iteration as follows. Step 1. Given a maximum number of CG iterations $N_{j}>0$, compute $\eta_{j}:=\min(\bar{\eta},\|\nabla\varphi(y^{j})\|^{1+\tau}).$ Apply the conjugate gradient (CG) algorithm $(CG(\eta_{j},N_{j}))$, to find an approximation solution $d^{j}$ to $\displaystyle(V_{j}+\varepsilon_{j}I)\,d=-\nabla\varphi(y^{j}),$ (30) where $V_{j}\in\hat{\partial}^{2}\varphi(y^{j})$ is defined as in (27) and $\varepsilon_{j}:=\tau_{1}\min\{\tau_{2},\|\nabla\varphi(y^{j})\|\}$. Step 2. Set $\alpha_{j}=\delta^{m_{j}}$, where $m_{j}$ is the first nonnegative integer $m$ for which $\displaystyle\varphi(y^{j}+\delta^{m}d^{j})\leq\varphi(y^{j})+\mu\delta^{m}\langle\nabla\varphi(y^{j}),d^{j}\rangle.$ (31) Step 3. Set $y^{j+1}=y^{j}+\alpha_{j}\,d^{j}$.

The convergence results for the above SNCG algorithm are stated in Theorems 2.2 and 2.3 below. We shall omit the proofs as they can be proved in the same fashion as in [25, Theorems 3.4 and 3.5].

Suppose that Assumption 2 holds. Then Algorithm SNCG generates a bounded sequence $\{y^{j}\}$ and any accumulation point $\hat{y}$ of $\{y^{j}\}$ is an optimal solution to problem (20).

Suppose that Assumption 2 holds. Let $\hat{y}$ be an accumulation point of the infinite sequence $\{y^{j}\}$ generated by Algorithm SNCG for solving the problem (20). Suppose that at each step $j\geq 0$, when the CG algorithm terminates, the tolerance $\eta_{j}$ is achieved (e.g., when $N_{j}=m+1$), i.e.,

Assume that the constraint nondegenerate condition

holds at $\widehat{W}:=\Pi_{{\cal K}}(\widetilde{X}+\sigma(\mathcal{A}^{*}\hat{y}-\widetilde{C}))$, where ${\rm lin}(\mathcal{T}_{{\cal K}}(\widehat{W}))$ denotes the lineality space of $\mathcal{T}_{{\cal K}}(\widehat{W})$. Then the whole sequence $\{y^{j}\}$ converges to $\hat{y}$ and

Given $\xi_{1}\in(0,1)$ and $\xi_{2}\in(0,\infty)$, we will use the following stopping criteria for terminating Algorithm SNCG:

$\varphi_{l}(y^{l+1})\leq\varphi_{l}(y^{l})-\frac{\xi_{1}}{2}\left|\langle\nabla\varphi_{l}(y^{l}),\,y^{l+1}-y^{l}\rangle\right|$,

$\|\nabla\varphi_{l}(y^{l+1})\|\leq\xi_{2}(\varphi_{l}(y^{l})-\varphi_{l}(y^{l+1}))^{\frac{1}{2}}$.

We can now state our majorized semismooth Newton-CG method for solving (16) as follows:

Algorithm MSNCG: A Majorized Semismooth Newton-CG Algorithm (MSNCG$(y^{0},S^{0},Z^{0},\widetilde{X},\sigma)$). Given $\xi_{1}\in(0,1)$, $\xi_{2}\in(0,+\infty)$. Perform the $l$th iteration as follows. Step 1. Starting with $y^{l}$ as the initial point, apply Algorithm SNCG to minimize $\varphi_{l}(\cdot)$ to find $y^{l+1}={\rm SNCG}(y^{l},\widetilde{X},\sigma)$ satisfying (A1) and (A2). Step 2. Compute $S^{l+1}:=\Pi_{{\cal K}^{*}}(C-{\cal A}^{*}y^{l+1}-Z^{l}-\sigma^{-1}\widetilde{X})$ and $Z^{l+1}:=\Pi_{{\cal P}^{*}}(C-{\cal A}^{*}y^{l+1}-S^{l+1}-\sigma^{-1}\widetilde{X})$.

Next, we establish the convergence of Algorithm MSNCG. For notational convenience, for any $y\in\Re^{m}$ and $S\in{\cal S}^{n}$, let $Y:=(y,S)$. Define the linear map ${\cal M}:\Re^{m}\times{\cal S}^{n}\to{\cal S}^{n}$ by

Let $B:=(b,0)\in\Re^{m}\times{\cal S}^{n}$ and ${\cal C}:=\Re^{m}\times{\cal K}^{*}$. Then problem (16) is equivalent to

and the function $\varPsi_{l}(y,S)$ in (17) can be rewritten as

Furthermore, $\widehat{Y}=(\hat{y},\widehat{S})$ is an optimal solution of

if and only if $\hat{y}$ is an optimal solution of problem (19a) and $\widehat{S}=\Pi_{{\cal K}^{*}}(\widetilde{C}^{l}-{\cal A}^{*}\hat{y}-\sigma^{-1}\widetilde{X}).$ In addition, conditions (A1) and (A2) are equivalent to the following two conditions, respectively,

Suppose that Assumption 2 holds. Then for Algorithm MSNCG, (A1) and (A2) are achievable.

If $Y^{l}-\Pi_{\cal C}(Y^{l}-\nabla\varPsi_{l}(Y^{l}))=0$, then one can take $Y^{l+1}=Y^{l}$ to satisfy (A1) and (A2).

Next, we assume that $Y^{l}-\Pi_{\cal C}(Y^{l}-\nabla\varPsi_{l}(Y^{l}))\neq 0$. Then $Y^{l}$ is not an optimal solution of problem (37). Let $\widehat{Y}$ be an arbitrary optimal solution of problem (37). Then $\widehat{Y}=\Pi_{\cal C}(\widehat{Y}-\nabla\varPsi_{l}(\widehat{Y}))$. So $\langle Y^{l}-\widehat{Y},\,(\widehat{Y}-\nabla\varPsi_{l}(\widehat{Y}))-\widehat{Y}\rangle\leq 0,$ i.e.,

we obtain that $\langle\nabla\varPsi_{l}(Y^{l}),\,\widehat{Y}-Y^{l}\rangle+\frac{\sigma}{2}\|{\cal M}(\widehat{Y}-Y^{l})\|^{2}<0.$ This implies

Then by using (40), (41), (42) and the fact that $\widehat{Y}=\Pi_{\cal C}(\widehat{Y}-\nabla\varPsi_{l}(\widehat{Y}))$, we know that for given $\xi_{1}\in(0,1)$ and $\xi_{2}\in(0,\infty)$, there exists $\delta>0$ such that

Suppose that Assumption 2 holds. Let Algorithm MSNCG be executed with stopping criteria (A1) and (A2). Then it generates a bounded sequence $\{(y^{l},S^{l},Z^{l})\}$ and any accumulation point $(\hat{y},\widehat{S})$ of $\{(y^{l},S^{l})\}$ is an optimal solution to problem (16) and hence $(\hat{y},\widehat{S},\widehat{Z})$ is an optimal solution to problem (10), where $\widehat{Z}:=\Pi_{{\cal P}^{*}}(C-{\cal A}^{*}\hat{y}-\widehat{S}-\sigma^{-1}\widetilde{X})$. Furthermore, $\|Z^{l+1}-Z^{l}\|\to 0$ as $l\to+\infty$.

By (38), we have $\varPhi(Y^{l+1})\leq\varPsi_{l}(Y^{l+1})\leq\varPsi_{l}(Y^{l})-\frac{\xi_{1}}{2}\left|\langle\nabla\varPsi_{l}(Y^{l}),\,Y^{l+1}-Y^{l}\rangle\right|=\varPhi(Y^{l})-\frac{\xi_{1}}{2}\left|\langle\nabla\varPsi_{l}(Y^{l}),\,Y^{l+1}-Y^{l}\rangle\right|$. Hence, the sequence $\{\varPhi(Y^{l})\}$ is nonincreasing.

By Lemma 2.1, we know that the level set $\mathcal{L}:=\{Y\in{\cal C}\,\mid\,\varPhi(Y)\leq\varPhi(Y^{0})\}$ is a closed and bounded convex set. Then the sequence $\{Y^{l}\}$ is bounded and so is the sequence $\{Z^{l}\}$. Let $\widehat{Y}$ be an accumulation point of $\{Y^{l}\}$. Then $\varPhi(Y^{l})\to\varPhi(\widehat{Y})$ and $\langle\nabla\varPsi_{l}(Y^{l}),\,Y^{l+1}-Y^{l}\rangle\to 0$ as $l\to\infty$. Furthermore, $\varPsi_{l}(Y^{l})-\varPsi_{l}(Y^{l+1})\rightarrow 0$ as $l\to\infty$.

Since $\langle\nabla\varPsi_{l}(Y^{l}),\,Y^{l+1}-Y^{l}\rangle\to 0$ as $l\to\infty$, we obtain from (44) that

For any $l\geq 0$, denote $\Delta_{l}:=Y^{l}-\Pi_{\cal C}(Y^{l}-\nabla\varPhi(Y^{l}))$. Then we have

By direct computations, we have for $l\geq 1$,

Thus, by (45) and the fact that $(\varPsi_{l}(Y^{l})-\varPsi_{l}(Y^{l+1}))^{\frac{1}{2}}\rightarrow 0$ as $l\to\infty$, we derive that $\|\Delta_{l+1}\|\to 0$ as $l\to\infty$. Since $\widehat{Y}$ is an accumulation point of $\{Y^{l}\}$, we obtain that $\widehat{Y}-\Pi_{\cal C}(\widehat{Y}-\nabla\varPhi(\widehat{Y}))=0$. By the convexity of $\varPhi$, $\widehat{Y}$ is an optimal solution of problem (36).

Finally, by using (45) and (46), we know that $\|Z^{l+1}-Z^{l}\|\to 0$ as $l\to\infty$. ∎

A Majorized Semismooth Newton-CG Augmented Lagrangian Method

For any $k\geq 0$ and $(y,S,Z)\in\Re^{m}\times{\cal S}^{n}\times{\cal S}^{n}$, denote

Since the inner problems in (10) are solved inexactly, we will use the following standard stopping criteria considered in to terminate Algorithm MSNCG:

$\hat{\phi}_{k}(y^{k+1},S^{k+1},Z^{k+1})-\inf\hat{\phi}_{k}\leq\epsilon_{k}^{2}$, $\epsilon_{k}\geq 0$, $\sum^{\infty}_{k=0}\epsilon_{k}<\infty$.

$\hat{\phi}_{k}(y^{k+1},S^{k+1},Z^{k+1})-\inf\hat{\phi}_{k}\leq(\delta_{k}^{2}/2\sigma_{k})\|X^{k+1}-X^{k}\|$, $\delta_{k}\geq 0$, $\sum^{\infty}_{k=0}\delta_{k}<\infty$.

${\rm dist}(0,\partial\hat{\phi}_{k}(y^{k+1},S^{k+1},Z^{k+1}))\leq(\delta_{k}^{{}^{\prime}}/\sigma_{k})\|X^{k+1}-X^{k}\|$, $0\leq\delta^{{}^{\prime}}_{k}\rightarrow 0$.

Just like SDPNAL, each iteration of the MSNCG algorithm can be quite expensive. Thus it is crucial for us to find a reasonably good initial point to warm start Algorithm SDPNAL$+$. We can certainly do so by solving the inner problem (10) by using any gradient descent type method. However, for this purpose we find that ADMM+ introduced by Sun, Toh and Yang is usually more efficient than other choices. Now we can present our SDPNAL$+$ algorithm as follows.

Algorithm SDPNAL$+$: A Majorized Semismooth Newton-CG Augmented Lagrangian Algorithm (SDPNAL$+$$(y^{0},S^{0},Z^{0},X^{0},\sigma_{0})$) Stage 1. Use ADMM+ to generate an initial point $(y^{0},S^{0},Z^{0},X^{0},\sigma_{0})\leftarrow{\rm ADMM+}(y^{0},S^{0},Z^{0},X^{0},\sigma_{0}).$ Stage 2. For $k=0,\ldots,$ perform the $k$th iteration as follows: (a) Using $(y^{k},S^{k},Z^{k})$ as the initial point, apply Algorithm MSNCG to minimize $\hat{\phi}_{k}(\cdot)$ to find $(y^{k+1},S^{k+1},Z^{k+1})=$ MSNCG $(y^{k},S^{k},Z^{k},X^{k},\sigma_{k})$ and $X^{k+1}=X^{k}+\sigma_{k}({\cal A}^{*}y^{k+1}+S^{k+1}+Z^{k+1}-C)$ satisfying (B1), (B2) or (B3). (b) Update $\sigma_{k+1}=\rho\sigma_{k}$ for some $\rho>1$ or $\sigma_{k+1}=\sigma_{k}$.

As mentioned in the introduction, if (P) is reformulated as a standard SDP and Algorithm SDPNAL$+$ is applied to this reformulated form, then SDPNAL$+$ reduces to SDPNAL proposed in .

We can obtain similar theorems on the convergence of SDPNAL$+$ as SDPNAL ([25, Theorems 4.1 and 4.2]). The global convergence of Algorithm SDPNAL$+$ follows from Rockafellar [18, Theorem $1$] and [17, Theorem $4$] without much difficulty.

Suppose that Assumption 2 holds. Let Algorithm SDPNAL$+$ be executed with stopping criterion (B1). If there exists $(y_{0},S_{0},Z_{0})\in\Re^{m}\times{\cal S}^{n}_{+}\times{\cal S}^{n}$ such that

then the sequence $\{X^{k}\}\subset{\cal P}$ generated by Algorithm SDPNAL$+$ is bounded and $\{X^{k}\}$ converges to $\overline{X}$, where $\overline{X}$ is some optimal solution to (P), and $\{(y^{k},S^{k},Z^{k})\}$ is asymptotically minimizing for (D) with $\max({\rm P})=\inf({\rm D})$.

If $\{X^{k}\}$ is bounded, then the sequence $\{(y^{k},S^{k},Z^{k})\}$ is also bounded, and all of its accumulation points of the sequence $\{(y^{k},S^{k},Z^{k})\}$ are optimal solutions to (D).

Next we state the local linear convergence of Algorithm SDPNAL$+$.

Suppose that Assumption 2 holds. Let Algorithm SDPNAL$+$ be executed with stopping criteria (B1) and (B2). Assume that (D) satisfies condition (49). If the second order sufficient conditions (in the sense of the conditions in [1, Theorem $3.137$]) holds at $\overline{X}$, where $\overline{X}$ is an optimal solution to (P), then the generated sequence $\{X^{k}\}\subset{\cal P}$ is bounded and $\{X^{k}\}$ converges to the unique optimal solution $\overline{X}$ with $\max({\rm P})=\min({\rm D})$, and

for some $\theta_{\infty}\in[0,1)$ with the property that $\theta_{\infty}\ll 1$ if $\sigma_{k}\to\sigma_{\infty}$ for any sufficiently large $\sigma_{\infty}$. The conclusions of Theorem 3.1 about $\{(y^{k},S^{k},Z^{k})\}$ are also valid.

The conclusions of Theorem 3.2 follow from the results in [18, Theorem 2] and [17, Theorem 5 and Proposition 3] combined with [1, Theorem $3.137$]. ∎

Numerical Experiments

+ and SDP Problem Sets In our numerical experiments, we test the following SDP$+$ and SDP problem sets.

(i) SDP$+$ problems coming from the relaxation of a binary integer nonconvex quadratic (BIQ) programming:

This problem has been shown in that under some mild assumptions, it can equivalently be reformulated as the following completely positive programming (CPP) problem:

where ${\cal C}_{pp}^{n}$ denotes the $n$-dimensional completely positive cone. It is well known that even though ${\cal C}_{pp}^{n}$ is convex, it is computationally intractable. To solve the CPP problem, one would typically relax ${\cal C}_{pp}^{n}$ to ${\cal S}^{n}_{+}\cap{\cal S}^{n}_{\geq 0}$, and the relaxed problem has the form (P):

where the polyhedral cone ${\cal P}=\{X\in{\cal S}^{n}\mid X\geq 0\}$. In our numerical experiments, the test data for $Q$ and $c$ are taken from Biq Mac Library maintained by Wiegele, which is available at http://biqmac.uni-klu.ac.at/biqmaclib.html.

(ii) SDP and SDP$+$ problems arising from the relaxation of maximum stable set problems. Given a graph $G$ with edge set $\mathcal{E}$, the SDP and SDP$+$ relaxation $\theta(G)$ and $\theta_{+}(G)$ of the maximum stable set problem are given by

where $E_{ij}=e_{i}e_{j}^{\mbox{{\tiny{T}}}}+e_{j}e_{i}^{\mbox{{\tiny{T}}}}$ and $e_{i}$ denotes the $i$th column of the identity matrix, ${\cal P}=\{X\in{\cal S}^{n}\mid X\geq 0\}$. In our numerical experiments, we test the graph instances $G$ considered in , , and .

(iii) SDP$+$ relaxation for computing lower bounds for quadratic assignment problems (QAPs). Let $\Pi$ be the set of $n\times n$ permutation matrices. Given matrices $A,B\in{\cal S}^{n}$, the QAP is given by

For a matrix $X=[x_{1},\dots,x_{n}]\in\Re^{n\times n}$, we will identify it with the $n^{2}$-vector $x=[x_{1};\dots;x_{n}]$. For a matrix $Y\in R^{n^{2}\times n^{2}}$, we let $Y^{ij}$ be the $n\times n$ block corresponding to $x_{i}x_{j}^{T}$ in the matrix $xx^{T}$. It is shown in that $v^{*}_{\rm QAP}$ is bounded below by the following number generated from the SDP$+$ relaxation of (59):

where $E$ is the matrix of ones, and $\delta_{ij}=1$ if $i=j$, and otherwise, ${\cal P}=\{X\in{\cal S}^{n^{2}}\mid X\geq 0\}$. In our numerical experiments, the test instances $(A,B)$ are taken from the QAP Library .

(iv) SDP$+$ relaxations of clustering problems (RCPs) described in [14, eq. (13), up to a constant]:

where $W$ is the so-called affinity matrix whose entries represent the similarities of the objects in the dataset, $e$ is the vector of ones, and $K$ is the number of clusters, ${\cal P}=\{X\in{\cal S}^{n}\mid X\geq 0\}$. All the data sets we tested are from the UCI Machine Learning Repository (available at http://archive.ics.uci.edu/ml/datasets.html). For some large data instances, we only select the first $n$ rows. For example, the original data instance “spambase” has 4601 rows, we select the first 1500 rows to obtain the test problem “spambase-large.2” for which the number “2” means that there are $K=2$ clusters.

(v) SDP+ problems arising from the SDP relaxation of frequency assignment problems (FAPs) . The explicit description of the SDP in the form ($P$) is given in [4, eq. (5)]:

where $k>1$ is an integer, $L(G,W):={\rm Diag}(We)-W$ is the Laplacian matrix, $E^{ij}=e_{i}e_{j}^{T}+e_{j}e_{i}^{T}$ with $e_{i}\in\Re^{n}$ being the $i$th standard unit vector and $e\in\Re^{n}$ is the vector of all ones. Let

where ${\cal P}=\{X\in{\cal S}^{n}\mid X_{ij}=0,\forall(i,j)\in U;X_{ij}\geq 0,\forall(i,j)\in E\setminus U\}$.

We should mention that we can easily extend our algorithm to handle the following more general SDP$+$ problem:

where $M\in{\cal X}$ is a given matrix. Thus (73) can also be solved by our proposed algorithm.

(vi) SDP relaxations for rank-1 tensor approximations (R1TA) :

It is shown in that (76) can be transformed into a standard SDP (up to a constant):

where $C$ is a constant matrix and ${\cal A}$ is a linear map, which depend on $M,f,g$.

2 Numerical Results

In this subsection, we compare the performance of our SDPNAL$+$ algorithm with two other competitive publicly available first order methods based codes for solving large-scale SDP$+$ and SDP problems: an ADMM based solver, called SDPADhttp://www.bicmr.org/~wenzw/code/SDPAD-release-beta2.zip (release-beta2, released in December 2012) developed in and a two-easy-block-decomposition hybrid proximal extragradient method, which was called 2EBD-HPEwww2.isye.gatech.edu/~cod3/CamiloOrtiz/Software_files/2EBD-HPE_v0.2/2EBD-HPE_v0.2.zip (v0.2, released on May 31, 2013) and we call it 2EBD here, introduced in . Since we use the convergent ADMM with $3$-block constraints introduced by Sun et al. (which was called ADMM3c but we call it ADMM$+$ here to indicate that it is an enhanced version of ADMM with convergence guantantee) to warm start SDPNAL$+$, we also list the numerical results obtained by running ADMM$+$ alone for the purpose of demonstrating the power and the importance of the proposed majorized semismooth Newton-CG algorithm for solving difficult SDP$+$ and SDP problems.

All our computational results for the tested SDP$+$ and SDP problems are obtained by running Matlab on a Linux server having 6 cores with 12 Intel Xeon X5650 processors at 2.67GHz and 32G RAM.

Note that numerically it is difficult to compute ${\rm dist}(0,\partial\hat{\phi}_{k}(W^{k+1}))$ in the criterion (B3) for terminating Algorithm MSNCG directly, where $W^{k+1}=(y^{k+1},S^{k+1},Z^{k+1})$. Fortunately, by using the fact that for any closed convex cone ${\cal C}\subseteq{\cal S}^{n}$ and $X\in{\cal C}$, it holds that ${\cal C}\subseteq{\cal T}_{{\cal C}}(X)$, we have from that

where $R_{D}^{k+1}={\cal A}^{*}y^{k+1}+S^{k+1}+Z^{k+1}-C$. In order to avoid computing $\Pi_{{\cal K}^{*}}(-R_{D}^{k+1}-\sigma_{k}^{-1}X^{k})$, we majorize the second part of (78) by a simpler term. Specifically, by using the fact that for any $X\in{\cal K}$ and $Y\in{\cal S}^{n}$, $\|\Pi_{{\cal K}^{*}}(Y-X)\|=\|\Pi_{{\cal K}^{*}}(Y-X)-\Pi_{{\cal K}^{*}}(-X)\|\leq\|Y\|$, we have

where $\widetilde{Z}^{k}$ is the $Z$ at the penultimate iteration when we compute $(y^{k+1},S^{k+1},Z^{k+1})=$ MSNCG $(y^{k},S^{k},Z^{k},X^{k},\sigma_{k})$. Note that ${\cal A}^{*}y^{k+1}+S^{k+1}+\widetilde{Z}^{k}+\sigma_{k}^{-1}X^{k}-C=\Pi_{{\cal K}}({\cal A}^{*}y^{k+1}+\widetilde{Z}^{k}+\sigma_{k}^{-1}X^{k}-C)\in{\cal K}$. Thus for a given $\zeta\in(0,+\infty)$, we can replace (B3) by the following criterion for terminating Algorithm MSNCG:

\max\Big{\{}\|{\cal A}(X^{k}+\sigma_{k}R_{D}^{k+1})-b\|,\zeta\sigma_{k}\|Z^{k+1}-\widetilde{Z}^{k}\|\Big{\}}\leq(\delta^{{}^{\prime}}_{k}/\sigma_{k})\|X^{k+1}-X^{k}\|, $0\leq\delta^{{}^{\prime}}_{k}\downarrow 0$.

In our numerical experiments, we measure the accuracy of an approximate optimal solution $(X,y,S,Z)$ for (P) and (D) by using the following relative residual:

where $\eta_{P}=\frac{\|{\cal A}X-b\|}{1+\|b\|}$, $\eta_{D}=\frac{\|{\cal A}^{*}y+S+Z-C\|}{1+\|C\|}$, $\eta_{{\cal K}}=\frac{\|\Pi_{{\cal K}^{*}}(-X)\|}{1+\|X\|}$, $\eta_{{\cal P}}=\frac{\|\Pi_{{\cal P}^{*}}(-X)\|}{1+\|X\|}$, $\eta_{{\cal K}^{*}}=\frac{\|\Pi_{{\cal K}}(-S)\|}{1+\|S\|}$, $\eta_{{\cal P}^{*}}=\frac{\|\Pi_{{\cal P}}(-Z)\|}{1+\|Z\|}$, $\eta_{C1}=\frac{\left|\langle X,\,S\rangle\right|}{1+\|X\|+\|S\|}$, $\eta_{C2}=\frac{\left|\langle X,\,Z\rangle\right|}{1+\|X\|+\|Z\|}$. Additionally, we compute the relative gap by

Let $\varepsilon>0$ be a given accuracy tolerance. We terminate both SDPNAL$+$ and ADMM$+$ when

Note that SDPAD can be used to solve SDP$+$ problems of form (P) with ${\cal P}={\cal S}^{n}_{\geq 0}$ directly and we stop SDPAD when $\eta<\varepsilon,$ where $\eta$ is defined as in (79). However, it is shown recently that the direct extension of ADMM to the multi-block case is not necessarily convergent . Hence SDPAD, which is essentially an implementation of the direct extension of ADMM with the step length set at $1.618$ for solving the dual of SDP$+$ problems, does not have convergence guarantee in theory.

The implementation of 2EBD including its termination, along with ADMM$+$ and SDPAD, is done in the same way as in . For 2EBD, we reformulate QAP, RCP and R1TA problems as SDP problems in the standard form as these problems do not appear to have obvious two-easy blocks structures.

In our numerical experiments, we also use a restart strategy for SDPNAL$+$ if it is not able to achieve the required accuracy for the tested SDP$+$ problems. For some problems, even though $\eta_{P}$ and $\eta_{D}$ can reach the required accuracy tolerance, $\eta_{{\cal K}}$ or $\eta_{C1}$ may stay above the required tolerance or stagnate. This may happen, as in the case for SDPNAL, because many of these SDP$+$ problems are degenerate at the optimal solutions. One way to overcome this difficulty is to apply ADMM$+$ to (P) using the most recently computed $(y,S,Z,X,\sigma)$ to restart SDPNAL$+$ when its progress is not satisfactory. From this point of view, our proposed algorithm is quite flexible.

Table 4.2 shows the number of problems that have been successfully solved to the accuracy of $10^{-6}$ in $\eta$ by each of the four solvers SDPNAL$+$, ADMM$+$, SDPAD and 2EBD, with the maximum number of iterations set at $25000$ or the maximum computation time set at $99$ hours. As can be seen, only SDPNAL$+$ can solve all the problems to the accuracy of $10^{-6}$. In particular, for the first time, we are able to solve all the $95$ difficult SDP$+$ problems arising from QAP problems to an accuracy of $10^{-6}$ efficiently, while ADMM$+$, SDPAD and 2EBD can successfully solve 39, 30 and 16 problems, respectively.

The authors would like to thank Jiawang Nie and Li Wang for sharing their codes on semidefinite relaxations of rank-1 tensor approximation problems.

References