2022

[1]Xin-Yee Lam These authors contributed equally to this work.

These authors contributed equally to this work.

[1]Department of Mathematics, National University of Singapore, 10 Lower Kent
Ridge Road, 119076, Singapore

2]Institute of Operations Research and Analytics, National University of
Singapore, 3 Research Link, 117602, Singapore

# On proximal augmented Lagrangian based decomposition methods for dual block-
angular convex composite programming problems

Kuang-Yu Ding<EMAIL_ADDRESS><EMAIL_ADDRESS>Kim-Chuan Toh
<EMAIL_ADDRESS>*  [

###### Abstract

We design inexact proximal augmented Lagrangian based decomposition methods
for convex composite programming problems with dual block-angular structures.
Our methods are particularly well suited for convex quadratic programming
problems arising from stochastic programming models. The algorithmic framework
is based on the application of the abstract inexact proximal ADMM framework
developed in [Chen, Sun, Toh, Math. Prog. 161:237–270] to the dual of the
target problem, as well as the application of the recently developed symmetric
Gauss-Seidel decomposition theorem for solving a proximal multi-block convex
composite quadratic programming problem. The key issues in our algorithmic
design are firstly in designing appropriate proximal terms to decompose the
computation of the dual variable blocks of the target problem to make the
subproblems in each iteration easier to solve, and secondly to develop novel
numerical schemes to solve the decomposed subproblems efficiently. Our inexact
augmented Lagrangian based decomposition methods have guaranteed convergence.
We present an application of the proposed algorithms to the doubly nonnegative
relaxations of uncapacitated facility location problems, as well as to two-
stage stochastic optimization problems. We conduct numerous numerical
experiments to evaluate the performance of our method against state-of-the-art
solvers such as Gurobi and MOSEK. Moreover, our proposed algorithms also
compare favourably to the well-known progressive hedging algorithm of
Rockafellar and Wets.

###### keywords:

dual block angular structure, progressive hedging algorithm, symmetric Gauss-
Seidel decomposition, ADMM, proximal augmented Lagrangian method, convex
composite programming

##  1 Introduction

Let $\mathcal{X}$ and ${\cal X}_{i}$ ($i=1,\ldots,N$) be $n_{0}$-dimensional
and $n_{i}$-dimensional inner product spaces each endowed with the inner
product $\langle\cdot,\,\cdot\rangle$ and its induced norm $\|\cdot\|$,
respectively. Similarly, let ${\cal Y}$ and ${\cal Y}_{i}$ ($i=1,\ldots,N$) be
$m_{0}$-dimensional and $m_{i}$-dimensional inner product spaces,
respectively. We consider the following general convex composite programming
problem with the so called dual block-angular structure (DBA) (according to
the terminology used by Dantzig Dantzig  in the case of linear programming):

$\displaystyle\begin{array}[]{rl}\min&\theta(x)+\langle
c,\,x\rangle+\sum_{i=1}^{N}\color[rgb]{0,0,0}\left(\bar{\theta}_{i}(\bar{x}_{i})+\langle\bar{c}_{i},\,\bar{x}_{i}\rangle\right)\\\\[5.0pt]
\mbox{s.t.}&\left[\begin{array}[]{ccccc}A&\\\\[5.0pt] B_{1}&\bar{B}_{1}&\\\
B_{2}&&\bar{B}_{2}&\\\ \vdots&&&\ddots\\\
B_{N}&&&&\bar{B}_{N}\end{array}\right]\left[\begin{array}[]{c}x\\\
\bar{x}_{1}\\\ \vdots\\\
\bar{x}_{N}\end{array}\right]\;=\;\left[\begin{array}[]{c}b\\\ \bar{b}_{1}\\\
\vdots\\\ \bar{b}_{N}\end{array}\right]\\\\[40.0pt] &x\in{\cal
K},\quad\bar{x}_{i}\in{\cal K}_{i}\subset{\cal
X}_{i},\;\forall\;i=1,\ldots,N,\end{array}$ (17)

where $b\in{\cal Y}$, $c\in{\cal X}$, $\bar{b}_{i}\in{\cal Y}_{i}$ and
$c_{i}\in{\cal X}_{i}$ ($i=1,\ldots,N$) are given data; $A:{\cal X}\to{\cal
Y}$, $B_{i}:{\cal X}\to{\cal Y}_{i}$ and $\bar{B}_{i}:{\cal X}_{i}\to{\cal
Y}_{i}$ are given linear maps; $\theta:{\cal X}\to(-\infty,\infty]$ and
$\bar{\theta}_{i}:{\cal X}_{i}\to(-\infty,\infty]$ are given proper closed
convex functions; ${\cal K}\subset{\cal X}$ and ${\cal K}_{i}\subset{\cal
X}_{i}$ ($i=1,\ldots,N$) are given closed convex sets, which are typically
cones but not necessarily so. In the above problem, we allow the case where
$A$ is absent. We will call the form in (1) as the standard DBA problem. In
the literature, there are other block-angular forms such as the standard
primal block-angular form and other variants which are not expressed in the
standard form. Following the terminology in the literature, we will sometimes
call the variable $x$ as the first-stage variable, and $\bar{x}_{i}$
($i=1,\ldots,N$) as the second-stage variables.

The above problem (which is reformulated as (P) in next page) that we consider
in this paper is very general. For example, it includes block-angular convex
quadratic programming problems as a special case when $\theta(\cdot)$ and
$\bar{\theta}_{i}(\cdot)$ ($i=1,\ldots,N$) are convex quadratic functions,
i.e, $\theta(x)=\frac{1}{2}\langle x,\,Qx\rangle$ and
$\bar{\theta}_{i}(\bar{x}_{i})=\frac{1}{2}\langle\bar{x}_{i},\,\bar{Q}_{i}\bar{x}_{i}\rangle$
where $Q:{\cal X}\to{\cal X}$ and $\bar{Q}_{i}:{\cal X}_{i}\to{\cal X}_{i}$
are self-adjoint positive semidefinite linear operators. As we shall see
later, our standard DBA problem (1) is much more general than other forms that
have been considered so far in the literature.

The dual block-angular problem generally has a separable convex objective
function and convex constraints but the linear constraints are coupled by the
linking variable ($x$) across different constraints. The dimension of these
problems can easily grow very large especially when the number of blocks, $N$,
increases. In particular, the number of linear equality constraints is equal
to $\sum_{i=0}^{N}m_{i}$. The main goal of this paper is to develop a
specialized algorithm that can fully exploit the block-angular structure to
solve the large-scale problem efficiently rather than treating it as a generic
convex programming problem for which even a powerful commercial solver may not
be able to handle efficiently.

In designing decomposition based methods, it is important to clearly specify
the block-angular structure one is dealing with, since the decomposition
strategies and the efficiency of the methods are invariably dictated by the
structure. In this paper, we will only be concerned with the DBA structure
presented in (1), but not touch on other block-angular structures as they will
lead to different decomposition strategies and methods. Consequently, our
literature review will also be confined to the DBA problems. Even then, the
literature on DBA problems is still too vast for us to give a comprehensive
coverage. Thus we shall only touch on papers which are most directly relevant
to our current paper when we described the related literature in Section 1.1.

In this paper, we will use the notation $[P\mathchar 24635\relax\;Q]$ to
denote the matrix obtained by appending the matrix $Q$ to the last row of the
matrix $P$, assuming that they have the same number of columns. We will also
use the same notation symbolically for linear maps $P$ and $Q$ as long as they
have the compatible domain and co-domain. Let $\bar{m}=\sum_{i=1}^{N}m_{i}$
and $\bar{n}=\sum_{i=1}^{N}n_{i}$ with $n_{i}$ being the dimension of ${\cal
X}_{i}$. Define

$\displaystyle{\bar{\cal X}}={\cal X}_{1}\times\cdots\times{\cal
X}_{N},\quad{\bar{\cal Y}}={\cal Y}_{1}\times\cdots\times{\cal
Y}_{N},\quad{\bar{\cal K}}={\cal K}_{1}\times\cdots\times{\cal K}_{N},$
$\displaystyle\bar{x}=[\bar{x}_{1}\mathchar 24635\relax\;\dots\mathchar
24635\relax\;\bar{x}_{N}]\in{\bar{\cal X}},\;\;\bar{c}=[\bar{c}_{1}\mathchar
24635\relax\;\dots\mathchar 24635\relax\;\bar{c}_{N}]\in{\bar{\cal
X}},\;\;\bar{b}=[\bar{b}_{1}\mathchar 24635\relax\;\dots\mathchar
24635\relax\;\bar{b}_{N}]\in{\bar{\cal Y}},$ $\displaystyle B=[B_{1}\mathchar
24635\relax\;\dots\mathchar
24635\relax\;B_{N}],\;\;\bar{B}=\mbox{diag}(\bar{B}_{1},\ldots,\bar{B}_{N}),\quad\bar{\theta}(\bar{x})=\sum_{i=1}^{N}\bar{\theta}_{i}(\bar{x}_{i}).$
(18)

Note that the dimension of ${\bar{\cal X}}$ is $\bar{n}$ and that of
${\bar{\cal Y}}$ is $\bar{m}.$ For convenience, we will sometimes refer to
$\bar{B}_{i}$ as an $m_{i}\times n_{i}$ matrix to mean its matrix
representation with respect to the standard bases in ${\cal X}_{i}$ and ${\cal
Y}_{i}$. A similar convention may also be used for $A$ and $B_{i}$.

We can rewrite (17) as an instance of the following general convex programming
problem with a block-angular structure:

$\displaystyle\begin{array}[]{rcll}\mbox{(P)}\quad\min&\theta(x)+\delta_{\cal
K}(x)+\langle c,\,x\rangle&+\;\bar{\theta}(\bar{x})+\delta_{{\bar{\cal
K}}}(\bar{x})+\langle\bar{c},\,\bar{x}\rangle\\\\[3.0pt]
\mbox{s.t.}&Ax&&=b,\quad x\in{\cal X}\\\\[3.0pt]
&Bx&+\;\bar{B}\bar{x}&=\bar{b},\quad\bar{x}\in{\bar{\cal X}}\end{array}$

where $A:{\cal X}\to{\cal Y}$, $B:{\cal X}\to{\bar{\cal Y}}$,
$\bar{B}:{\bar{\cal X}}\to{\bar{\cal Y}}$, and $\delta_{\cal K}(\cdot)$,
$\delta_{\bar{\cal K}}(\cdot)$ denote the indicator functions over ${\cal K}$,
${\bar{\cal K}}$, respectively. The above structure is often referred to in
the literature as a dual block-angular structure. But we shall drop the word
“dual” in this paper.

We shall show in Section 3 that the dual of (P) is given by

$\displaystyle\begin{array}[]{rrl}\mbox{(D)}\quad-\min&\theta^{*}(-v)+\delta_{\cal
K}^{*}(-z)-\langle
b,\,y\rangle\;+\bar{\theta}^{*}(-\bar{v})+\delta_{{\bar{\cal
K}}}^{*}(-\bar{z})-\langle\bar{b},\,\bar{y}\rangle\\\\[5.0pt]
\mbox{s.t.}&A^{*}y+B^{*}\bar{y}+z+v=c,\;\;z\in{\cal X},\;v\in{\cal
X}\\\\[5.0pt]
&\;\;\;\bar{B}^{*}\bar{y}+\bar{z}+\bar{v}=\bar{c},\;\;\bar{z}\in{\bar{\cal
X}},\;\bar{v}\in{\bar{\cal X}}\end{array}$

where $\bar{z}=[\bar{z}_{1}\mathchar 24635\relax\;\dots\mathchar
24635\relax\;\bar{z}_{N}]\in{\bar{\cal X}},$ $\bar{v}=[\bar{v}_{1}\mathchar
24635\relax\;\ldots\mathchar 24635\relax\;\bar{v}_{N}]\in{\bar{\cal X}}$,
$\bar{y}=[\bar{y}_{1}\mathchar 24635\relax\;\dots\mathchar
24635\relax\;\bar{y}_{N}]\in{\bar{\cal Y}}$, $\delta^{*}_{\cal K}$ is the
conjugate function of $\delta_{\cal K}$ (similarly for $\delta^{*}_{\bar{\cal
K}}$). Observe that the problem (D) has six blocks of variables
$(y,\bar{y},z,\bar{z},v,\bar{v})$ coupled together by the linear constraints
$A^{*}y+B^{*}\bar{y}+z+v=c$ and $\bar{B}^{*}\bar{y}+\bar{z}+\bar{v}=\bar{c}$.

We made the following assumptions on the problem data throughout this paper.

######  Assumption 1.

$A\in\mathbb{R}^{m\times n}$ has full row rank (if it is present), and
$[B,\bar{B}]$ has full row rank.

######  Assumption 2.

The projection $\Pi_{\cal K}(x)$ and $\Pi_{\bar{\cal K}}(\bar{z})$ can be
computed efficiently such as in the case when ${\cal K}$ is a nonnegative
orthant, a second-order cone or $\mathbb{S}^{n}_{+}$ (the cone of symmetric
positive semidefinite matrices). We also assume that the proximal mappings for
$\theta$ and $\bar{\theta}$ can be computed analytically or very efficiently.

In practical applications, optimization problems with dual block-angular
structures arise in many contexts ZFEM ; FLS ; Wollenberg . In this paper, we
focus on two classes of problems, namely, two-stage stochastic convex
programming problems and doubly nonnegative relaxations of uncapacitated
facility location problems, for the purpose of evaluating the numerical
performance of our algorithm against other existing state-of-the-art solvers.

In the case of two-stage stochastic convex programming problems, the
optimization model often has a linear or quadratic objective function.
Commonly used methods to solve two-stage stochastic linear-quadratic
programming problems can mainly be classified into three types. The first type
is based on the cutting plane method, including Bender decomposition methods
Bender  and L-shaped methods Lshape . However, this class of methods is known
to have very slow convergence TLAG . The second type is based on interior
point methods, where details can be found in BeDeOZ99 ; birge95 ; zhaogy991 ;
zhao2001  and Mehrotra’s works SM1 ; SM2 . The last type is decomposition
methods based on augmented Lagrangian method (ALM), see for example RTRW ; CDZ
; ruszcz95 ; ruszcz99  and the references therein. Others include distributed
simplex methods ParSplx , and Karmarkar algorithms Karm ; BirQi88 . Among the
decomposition methods, the Progressive Hedging Algorithm (PHA) in RTRW , which
is a scenario decomposition based method, is perhaps the best known for
solving 2-stage stochastic linear-quadratic programming problems. For the
class of uncapacitated facility location problems, they are originally binary
integer programming problems. We relax them to doubly nonnegative (DNN) convex
programming problems and solve the resulting dual block-angular DNN problems.

The problem (P) is generally large scale and developing efficient algorithms
that are capable of handling such a large problem is of great interest. In
this paper, our main objective is to design and implement an efficient and
robust (distributed) algorithm for solving the problems (P) and (D) by fully
exploiting the dual block-angular structure. Our method is based on the
abstract inexact symmetric Gauss-Seidel proximal ADMM (sGS-ADMM) framework
developed in CST . Specifically, we develop a concrete sGS-ADMM to solve the
block-angular problem (P) via its dual problem (D). While our algorithm is
based on the abstract algorithmic framework in CST , the key issue in the
current algorithmic design is to fully exploit the block-angular structure in
our target problem to decompose the computation of the variable blocks, with
the goal of making the subproblems in each iteration easier to solve.
Moreover, we also need to develop the corresponding numerical schemes to solve
the subproblems efficiently. As we shall see later, our method can fully
exploit the block-angular structure in the problem and hence the subproblems
in each iteration can be solved efficiently. In addition, it also offers us
the flexibility to solve the subproblems in parallel.

Now we summarize our contributions in this paper.

  1. 1.

We derived concrete DBA convex composite conic programming problems arising
from DNN relaxations of common mixed integer programming problems such as
uncapacitated facility location problems.

  2. 2.

We proposed and implemented a specialized algorithm to solve large scale DBA
problems via their duals. In particular, by leveraging on the recently
developed sGS decomposition theorem LST-sGS , we have designed proximal ADMM
and proximal ALM frameworks that are able to fully exploit the underlying
block angular structure. This is done by designing appropriate proximal terms
within a proximal ADMM or ALM framework to enable the decomposition of the
dual block variables’ computation. By further designing novel numerical
schemes to efficiently solve the decomposed subproblems, each iteration of the
overall algorithm is relatively easy and inexpensive to execute. In addition,
our algorithm can conveniently take advantage of parallel computation to
further improve its computational efficiency.

  3. 3.

We conduct comprehensive numerical experiments to compare our algorithm with
state-of-the-art solvers. The results show that our algorithm outperforms the
others especially when the constraint set ${\cal K}_{i}$ is the cone
$\mathbb{S}^{n_{i}}_{+}$ or the quadratic coefficient matrix in the objective
function is dense. For two stage stochastic linear programming, we also
compare our algorithm with the progressive hedging algorithm (PHA) in RTRW .
The numerical results show that our algorithm can perform much better when the
number of scenarios is large.

In the current literature, as far as we are aware of, all the popular
augmented Lagrangian based decomposition methods (including the PHA) for DBA
linear programming problems directly target the primal problem (P). Here our
methods are designed to deal with the structure in the dual problem (D) while
at the same time solving the primal problem (P) as a by-product. The key
advantage in this dual approach is that we do not need to introduce any extra
variables and constraints, unlike the case of directly dealing with the primal
problem where a large number of extra variables and constraints must be
introduced to split the first-stage variable $x$ so as to decompose the
computation of the second-stage variables $\bar{x}_{i}$ ($i=1,\ldots,N$). In
contrast, for our dual approach, the decomposition of the computation of the
dual variables can be obtained naturally with the sGS decomposition framework.

This paper is organized as follows. We will describe some example classes of
convex programming problems with dual block-angular structures in section 2.
In section 3, we will derive the dual (D) of the primal problem (P). We
provide some preliminaries in section 4 before providing the algorithmic
details of our inexact symmetric Gauss-Seidel proximal ADMM and ALM for
solving (D) in section 5. In section 6, we discuss an extension where we
incorporate a semismooth Newton-CG method in our proposed algorithm. We
conduct numerical experiments in section 7 to evaluate the performance of our
algorithms. Finally, we conclude the paper in the last section.

Notation.

  * •

$\mathbb{R}^{d}$ denotes the $d$-dimensional Euclidean space, $\mathbb{S}^{d}$
denotes the space of the $d$-dimensional symmetric matrices,
$\mathbb{S}^{d}_{+}$ denotes the $d$-dimensional positive semidefinite cone.

  * •

Let ${\cal T}:{\cal X}\rightarrow{\cal X}$ be a self-adjoint positive
semidefinite map. We denote $\|x\|_{\cal T}:=\sqrt{\langle x,\,{\cal
T}x\rangle}$ as the semi-norm induced by ${\cal T}$. For any linear map ${\cal
B}:{\cal X}\to{\cal Y}$, where ${\cal X}$ and ${\cal Y}$ are finite-
dimensional inner product spaces, its induced operator norm $\|{\cal
B}\|=\max_{x\not=0}\|{\cal B}x\|/\|x\|.$

  * •

Let $f:{\cal X}\rightarrow(-\infty,+\infty]$ be a closed proper convex
function. Then, we denote dom$f$ and $\partial f$ as its effective domain and
subdifferential, respectively.

  * •

We denote $f^{*}$ as the Fenchel conjugate function of a convex closed
function $f$, which is defined by $f^{*}(x):=\sup_{y}\\{\langle
x,\,y\rangle-f(y)\\}$.

  * •

The Moreau-Yosida envelope of $f$ is defined by
$M_{f}(x):=\min_{y}\\{f(y)+\frac{1}{2}\|x-y\|^{2}\\}$, and proximal mapping of
$f$ is defined by
$\text{Prox}_{f}(x):=\arg\min_{y}\\{f(y)+\frac{1}{2}\|x-y\|^{2}\\}$.

  * •

The following formulas are frequently used in this paper:

$\displaystyle\nabla M_{f^{*}/t}(x)=t^{-1}\text{Prox}_{tf}(tx)\quad
and\quad\nabla M_{tf}(x)=t\text{Prox}_{f^{*}/t}(x/t)$ $\displaystyle
x=\text{Prox}_{tf}(x)+t\text{Prox}_{f^{*}/t}(x/t)\quad
and\quad\frac{1}{2}\|x\|^{2}=M_{tf}(x)+t^{2}M_{f^{*}/t}(x/t)$

  * •

For given square matrices $P_{1},\ldots,P_{\ell}$, we use ${\rm
diag}(P_{1},\ldots,P_{\ell})$ to denote the block diagonal matrix with
diagonal blocks given by $P_{1},\ldots,P_{\ell}$.

###  1.1 Related literature

DBA linear-quadratic problems are well-known in the context of two-stage
stochastic programming. In the pioneering work of RTRW , a multistage
programming problem was considered where the variables are vectors in some
feasible set in the Euclidean spaces, and the objective function is linear-
quadratic. In that paper, the progressive hedging algorithm (PHA) was
proposed, which is an algorithmic framework based on a proximal augmented
Lagrangian method whereby at each iteration, each scenario requires the
computation of one linear-quadratic programming subproblem. The PHA can be
costly when the number of scenarios is large and/or when the subproblems are
difficult to solve. Subsequently, augmented Lagrangian based methods like
Diagonal Quadratic Approximation method (DQA) in ruszcz95  and Accelerated
Distributed Augmented Lagrangian method (ADAL) in CDZ  have been proposed.
These two methods are initially designed to solve large-scale convex
optimization problems with primal block angular (PBA) structure that have many
linking constraints. Both methods aim to make their subproblems easy to solve
by modifying the augmented Lagrangian function and try to make it separable.
However, most of these methods are still difficult to be applied to a large-
scale problem when the underlying feasible set is nontrivial. Just like the
PHA, the computational cost of solving those separable smaller linear-
quadratic programming subproblems can become expensive when the number of
subproblems is large and/or when the subproblems are difficult to solve.

Other than two-stage stochastic programming problem, there are less literature
devoted to designing specialized methods for directly solving dual block
angular problems. Most of them construct a general framework and solve DBA
problems as a by-product. For example, in Kontogiorgis , the authors developed
a framework based on alternating directions method for solving the
minimization of the sum of two convex functions subject to linear constraints.
It is then applied successfully to solve the Fermat-Weber problem recasted
into a dual-block angular problem.

While many papers have been written on augmented Lagrangian based
decomposition methods for solving DBA linear and quadratic programming
problems, there are far fewer papers on other classes of DBA problems such as
DBA convex conic programming problems where the cones are more general than
the nonnegative orthants. In SM1 , the authors proposed a log-barrier
decomposition interior-point method (IPM) for a 2-stage stochastic SDP that is
expressed as the dual of the primal block-angular form, wherein the algorithm
is based on the pioneering work of Zhao zhao2001  on a log-barrier Benders
decomposition IPM for a 2-stage stochastic LP. But note that such a
decomposition method is not directly designed for the DBA problem (1), so we
shall not touch on such a method any further.

In SM2 , a log-barrier Benders decomposition IPM is proposed for a 2-stage
stochastic convex QP expressed in the dual block-angular form, and the
polynomial iteration complexity of the method is analysed, but there are no
numerical implementation. The IPM is extended in Chen-Mehrotra  for a 2-stage
stochastic convex optimization problem where the nonnegative orthants are
replaced by compact convex sets whose barrier functions (and the corresponding
gradients and Hessians) are assumed to be readily computable.

Our algorithms are applicable to the very general standard DBA problem (1). In
Chen-Mehrotra , a theoretical log-barrier decomposition IPM has been proposed
to solve dual block-angular convex programming problems with linear objective
functions (i.e, with $\theta(\cdot)\equiv 0$ and $\bar{\theta}(\cdot)\equiv
0$). However, the implementation and practical evaluation of the methods so
far have not been performed.

In LST-PBA , the authors considered convex optimization problems with the
primal block angular (PBA) structure, which has a different structure from the
dual block-angular structure considered in this paper. They also consider the
dual problem to deal with the PBA structure by sGS-ADMM, but the design of the
proximal terms are different from ours in this paper. Moreover, numerical
experiments on more general cone like SDP cone are not performed in LST-PBA .

##  2 Example classes of (17)

###  2.1 Examples from 2-stage stochastic convex programming problems

Here we show that the dual block-angular structure shown in (17) naturally
arises from a 2-stage stochastic convex programming problem. Consider the
following 2-stage stochastic optimization problem:

$\displaystyle\min_{x\in{\cal X}}\Big{\\{}\theta(x)+\langle
c,\,x\rangle+\color[rgb]{0,0,0}E_{\xi}\left(Q(x\mathchar
24635\relax\;\xi)\right)\mid Ax\;=\;b,\;x\in{\cal K}\Big{\\}}.$ (21)

Here $\xi$ is a random vector, $E_{\xi}(\cdot)$ denotes the expectation with
respect to $\xi$, and

$\displaystyle Q(x\mathchar 24635\relax\;\xi)\;=\;\min_{\tilde{x}\in{\cal
X}_{1}}\left\\{\tilde{\theta}_{\xi}(\tilde{x})+\langle\tilde{c}_{\xi},\,\tilde{x}\rangle\,:\,\tilde{B}_{\xi}\tilde{x}=\tilde{b}_{\xi}-B_{\xi}x,\;\tilde{x}\in{\cal
K}_{1}\right\\},$

where for each given scenario $\xi$, $\tilde{c}_{\xi}\in{\cal X}_{1}$,
$\tilde{b}_{\xi}\in{\cal Y}_{1}$ are given data, $B_{\xi}:{\cal X}\to{\cal
Y}_{1}$, $\tilde{B}_{\xi}:{\cal X}_{1}\to{\cal Y}_{1}$ are given linear maps,
and ${\cal K}_{1}\subset{\cal X}_{1}$ is a given closed convex set,
$\tilde{\theta}_{\xi}:{\cal X}_{1}\to(-\infty,\infty]$ is a given proper
closed convex function. Note that there is no loss of generality in
considering only equality constraints in (21) since constraints of the form
$b-Ax\in{\cal Q}$, $x\in{\cal K}$ can always be reformulated as
$[A,I](x\mathchar 24635\relax\;s)=b$, $(x\mathchar 24635\relax\;s)\in{\cal
K}\times{\cal Q}$.

By sampling $N$ scenarios for $\xi$, one may approximate $E_{\xi}Q(x\mathchar
24635\relax\;\xi)$ by the sample mean $\sum_{i=1}^{N}p_{i}Q(x\mathchar
24635\relax\;\xi_{i})$, with $p_{i}$ being the probability of occurrence of
the $i$th scenario, and hence approximately solve (21) via the following
deterministic optimization problem which has exactly the same form as (17):

$\displaystyle\min\left\\{\theta(x)+\langle
c,\,x\rangle+\sum_{i=1}^{N}\color[rgb]{0,0,0}\left(\bar{\theta}_{i}(\bar{x}_{i})+\langle\bar{c}_{i},\,\bar{x}_{i}\rangle\right)\Bigm{|}\begin{array}[]{l}Ax=b,\;x\in{\cal
K}\\\\[5.0pt] B_{i}x+\bar{B}_{i}\bar{y}_{i}=\bar{b}_{i},\;\bar{y}_{i}\in{\cal
K}_{i},\;\forall\;i=1,...,N\end{array}\right\\}\;$ (24)

where $\bar{c}_{i}=p_{i}\tilde{c}_{\xi_{i}}$ and
$\bar{\theta}_{i}=p_{i}\tilde{\theta}_{\xi_{i}}$,
$B_{i}=B_{\xi_{i}},\bar{B}_{i}=\tilde{B}_{\xi_{i}},\bar{b}_{i}=\tilde{b}_{\xi_{i}}$
are the data, and $\bar{x}_{i}=\tilde{x}_{\xi_{i}}\in{\cal X}_{1}$ is the
second stage decision variable associated with the $i$th scenario.

###  2.2 Examples from doubly nonnegative relaxations of uncapacitated
facility location problems

The uncapacitated facility location (UFL) problem is one of the most studied
problems in operations research. Here we consider the general UFL with linear
and/or separable convex quadratic allocation costs that was introduced in
Gunluk :

$\displaystyle\begin{array}[]{rl}\mbox{(UFL)}\quad\min&\sum_{i=1}^{p}c_{i}u_{i}+\sum_{i=1}^{p}\sum_{j=1}^{q}\color[rgb]{0,0,0}\left(p_{ij}s_{ij}+\frac{1}{2}q_{ij}s_{ij}^{2}\right)\\\\[5.0pt]
\mbox{s.t}&\sum_{i=1}^{n}s_{ij}=1,\;\forall\;j=1,\ldots,q\\\\[5.0pt] &0\leq
s_{ij}\leq u_{i},\;\forall\;i=1,\ldots,p,\;j=1,\ldots,q\\\\[5.0pt]
&u_{i}\in\\{0,1\\}\;\forall\;i=1,\ldots,p,\end{array}$ (29)

where $c_{i}\geq 0$ for all $i$, and $\color[rgb]{0,0,0}p_{ij},q_{ij}\geq 0$
for all $i,j$ are given data. Observe that the allocation cost for customer
$j$ is a convex quadratic function of his demand $s_{ij}$ being served by the
opened facility $i$.

Let $U=uu^{T}$. One can see that the constraint $u_{i}\in\\{0,1\\}$ is
equivalent to the constraint that $u_{i}^{2}=U_{ii}=u_{i}$. Also, by
introducing the nonnegative slack variable $z_{ij}$, we can convert the
inequality constraint $s_{ij}\leq u_{i}$ to the equality constraint
$s_{ij}+z_{ij}=u_{i}$.

Let $S=(s_{ij})\in\mathbb{R}^{p\times q}$ and
$Z=(z_{ij})\in\mathbb{R}^{p\times q}$. Also let $e\in\mathbb{R}^{p}$ be the
vector of all ones and $S_{j}$ be the $j$th column of $S$. Then one can
express the problem (29) equivalently as follows:

$\displaystyle\begin{array}[]{rl}\min&\sum_{i=1}^{p}c_{i}u_{i}+\sum_{i=1}^{p}\sum_{j=1}^{q}\color[rgb]{0,0,0}\left(p_{ij}s_{ij}+\frac{1}{2}q_{ij}s_{ij}^{2}\right)\\\\[5.0pt]
\mbox{s.t}&e^{T}S_{j}=1,\;\;\forall\;j=1,\ldots,q\\\\[5.0pt]
&S_{j}+Z_{j}=u,\;\;S_{j},Z_{j}\geq 0,\;\;\forall\;j=1,\ldots,q\\\\[5.0pt]
&u-\mbox{diag}(U)=0,\quad U=uu^{T},\;u\geq 0,\;U\in\mathbb{S}^{p}.\end{array}$

Let $\color[rgb]{0,0,0}P=(p_{ij})\in\mathbb{R}^{p\times q}$ and
$Q=(q_{ij})\in\mathbb{R}^{p\times q}.$ By relaxing the rank-1 constraint
$U=uu^{T}$ to $U\succeq uu^{T}$ and using the equivalence that $U\succeq
uu^{T}$ if and only if $[1,\;u^{T}\mathchar 24635\relax\;u,\;U]\succeq 0$, we
get the following doubly nonnegative (DNN) relaxation of (29):

$\displaystyle\begin{array}[]{rl}\min&\langle[0\mathchar
24635\relax\;c],\,[\alpha\mathchar
24635\relax\;u]\rangle+\sum_{j=1}^{q}\color[rgb]{0,0,0}\left(\langle
P_{j},\,S_{j}\rangle+\frac{1}{2}\langle
S_{j},\,\mbox{diag}(Q_{j})S_{j}\rangle\right)\\\\[8.0pt]
\mbox{s.t}&\left[\begin{array}[]{c}0\\\
u\end{array}\right]+\left[\begin{array}[]{cc}e^{T}&0^{T}\\\
-I_{p}&-I_{p}\end{array}\right]\left[\begin{array}[]{c}S_{j}\\\
Z_{j}\end{array}\right]=\left[\begin{array}[]{c}1\\\ 0\end{array}\right],\quad
j=1,\ldots q\\\\[15.0pt] &u-\mbox{diag}(U)=0,\;\;\alpha=1,\\\\[5.0pt]
&\mbox{\boldmath$U$}:=\left[\begin{array}[]{cc}\alpha&u^{T}\\\\[3.0pt]
u&U\end{array}\right]\in\mathbb{S}^{1+p}_{+},\;\;\mbox{\boldmath$U$}\geq
0,\;\;S_{j},Z_{j}\geq 0,\;j=1,\ldots,q.\end{array}$ (45)

Now we can express the above DNN programming problem in the form (17) by
defining the following inner product spaces, functions, cones and linear maps:

$\displaystyle{\cal X}=\mathbb{S}^{1+p},\;\;{\cal
Y}=\mathbb{R}^{1+p},\;\;{\cal X}_{j}=\mathbb{R}^{2q},\;\;{\cal
Y}_{j}=\mathbb{R}^{1+p},\;j=1,\ldots,q,$
$\displaystyle\theta(\mbox{\boldmath$U$})=\delta_{{\cal
N}}(\mbox{\boldmath$U$}),\;{\cal K}=\mathbb{S}^{1+p}_{+},$
$\displaystyle\bar{\theta}_{j}([S_{j}\mathchar
24635\relax\;Z_{j}])=\frac{1}{2}\langle
S_{j},\,\mbox{diag}(Q_{j})S_{j}\rangle,\;{\cal
K}_{j}=\mathbb{R}^{2q}_{+},\;j=1,\ldots,q,$ $\displaystyle
b=\left[\begin{array}[]{c}1\\\\[5.0pt]
0\end{array}\right]\in\mathbb{R}^{1+p},\;\;A(\mbox{\boldmath$U$})=\left[\begin{array}[]{c}\alpha\\\\[5.0pt]
u-\mbox{diag}(U)\end{array}\right]=\left[\begin{array}[]{c}\langle
A_{0},\,\mbox{\boldmath$U$}\rangle\\\ \langle
A_{1},\,\mbox{\boldmath$U$}\rangle\\\ \vdots\\\ \langle
A_{p},\,\mbox{\boldmath$U$}\rangle\end{array}\right],$ $\displaystyle
B_{j}(\mbox{\boldmath$U$})=\left[\begin{array}[]{c}0\\\\[5.0pt]
u\end{array}\right]=\left[\begin{array}[]{c}\langle
O,\,\mbox{\boldmath$U$}\rangle\\\ \langle
E_{1},\,\mbox{\boldmath$U$}\rangle\\\ \vdots\\\ \langle
E_{p},\,\mbox{\boldmath$U$}\rangle\end{array}\right],\;\;\bar{B}_{j}=\left[\begin{array}[]{cc}e^{T}&0^{T}\\\
-I_{p}&-I_{p}\end{array}\right]\in\mathbb{R}^{(1+p)\times 2q},$
$\displaystyle\bar{b}_{j}=\left[\begin{array}[]{c}1\\\
0\end{array}\right]\in\mathbb{R}^{1+p},\;\;j=1,\ldots,q,$

where ${\cal N}=\\{X\in\mathbb{S}^{1+p}\mid X\geq 0\\}$, $Q_{j}$ is the $j$th
column of $Q$, and

$\displaystyle A_{0}=\left[\begin{array}[]{cc}1&0^{T}\\\\[3.0pt] 0&0_{p\times
p}\end{array}\right],\quad
A_{i}=\left[\begin{array}[]{cc}0&\frac{1}{2}e_{i}^{T}\\\\[3.0pt]
\frac{1}{2}e_{i}&-e_{i}e_{i}^{T}\end{array}\right],\quad
E_{i}=\left[\begin{array}[]{cc}0&\frac{1}{2}e_{i}^{T}\\\\[3.0pt]
\frac{1}{2}e_{i}&0_{p\times p}\end{array}\right],\;\;i=1,\ldots,p.$

In the above, $e_{i}$ is the $i$th unit vector in $\mathbb{R}^{p}.$ Observe
that for the problem (45), we have that

$B_{1}=\cdots=B_{q},\quad\bar{B}_{1}=\cdots=\bar{B}_{q}.$

Thus given any $\mbox{\boldmath$U$}\in{\cal X}$ and $\bar{y}\in{\cal
X}_{1}\times\cdots\times{\cal X}_{q}$, the evaluation of
$B(\mbox{\boldmath$U$})$ and $B^{*}\bar{y}$ can be done very efficiently since
$B(\mbox{\boldmath$U$})=[B_{1}(\mbox{\boldmath$U$})]_{j=1}^{q}$ and
$B^{*}\bar{y}=B_{1}^{*}(\bar{y}_{1}+\cdots+\bar{y}_{q})$. For later purpose,
we note that $\bar{B}_{j}\bar{B}_{j}^{*}$ has a very simple inverse given by

$\displaystyle(\bar{B}_{j}\bar{B}_{j}^{*})^{-1}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}0&0^{T}\\\
0&\frac{1}{2}I_{p}\end{array}\right]+\frac{1}{2p}\left[\begin{array}[]{c}2\\\
e\end{array}\right][2,\;e^{T}],\quad j=1,\ldots,q.$ (54)

######  Remark 1.

If the problem (29) has additional constraints in $u$, say $Hu=\beta$, then
the corresponding constraints in the DNN relaxation of the problem would
become

$\displaystyle\langle
H_{k},\,\mbox{\boldmath$U$}\rangle=\beta_{k},\quad\langle\color[rgb]{0,0,0}\widehat{H}_{k},\,\mbox{\boldmath$U$}\rangle=\beta_{k}^{2},\quad\quad
k=1,\ldots m$

where given the $k$th row $h_{k}^{T}$ of $H$,

$H_{k}=\left[\begin{array}[]{cc}0&h_{k}^{T}/2\\\\[3.0pt]
h_{k}/2&0\end{array}\right],\quad\widehat{H}_{k}=\left[\begin{array}[]{cc}0&0^{T}\\\\[3.0pt]
0&h_{k}h_{k}^{T}\end{array}\right],\quad k=1,\ldots,m.$

##  3 The derivation of (D) and KKT conditions

By introducing the auxiliary variables, $u=x$, $\bar{u}=\bar{x}$, $s=x$,
$\bar{s}=\bar{x}$ and replacing $\theta(x)$, $\bar{\theta}(\bar{x})$,
$\delta_{{\cal K}}(x)$, $\delta_{{\bar{\cal K}}}(\bar{x})$ by $\theta(u)$,
$\bar{\theta}(\bar{u})$, $\delta_{{\cal K}}(s)$, $\delta_{{\bar{\cal
K}}}(\bar{s})$ respectively, one can derive the dual of (P). Consider the
Lagrangian function for (P):

$\displaystyle{\cal L}(x,\bar{x},u,\bar{u},s,\bar{s}\mathchar
24635\relax\;y,\bar{y},v,\bar{v},z,\bar{z})$ $\displaystyle:=$
$\displaystyle\theta(u)+\langle
c,\,x\rangle+\bar{\theta}(\bar{u})+\langle\bar{c},\,\bar{x}\rangle+\delta_{{\cal
K}}(s)+\delta_{{\bar{\cal K}}}(\bar{s})-\langle
y,\,Ax-b\rangle-\langle\bar{y},\,Bx+\bar{B}\bar{x}-\bar{b}\rangle$
$\displaystyle-\langle
v,\,x-u\rangle-\langle\bar{v},\,\bar{x}-\bar{u}\rangle-\langle
z,\,x-s\rangle-\langle\bar{z},\,\bar{x}-\bar{s}\rangle$ $\displaystyle=$
$\displaystyle\langle
c-A^{*}y-B^{*}\bar{y}-v-z,\,x\rangle+\color[rgb]{0,0,0}\langle\bar{c}-\bar{B}^{*}\bar{y}-\bar{v}-\bar{z},\,\bar{x}\rangle+\theta(u)+\langle
v,\,u\rangle+\bar{\theta}(\bar{u})+\langle\bar{v},\,\bar{u}\rangle$
$\displaystyle+\delta_{{\cal K}}(s)+\langle z,\,s\rangle+\delta_{{\bar{\cal
K}}}(\bar{s})+\langle\bar{z},\,\bar{s}\rangle+\langle
b,\,y\rangle+\langle\bar{b},\,\bar{y}\rangle,$

where $x,u,s,v,z\in{\cal X}$,
$\bar{x},\bar{u},\bar{s},\bar{v},\bar{z}\in{\bar{\cal X}}$, $y\in{\cal Y}$ and
$\bar{y}\in{\bar{\cal Y}}$. We have the following:

$\displaystyle\inf_{x}{\cal L}(x,\bar{x},u,\bar{u},s,\bar{s}\mathchar
24635\relax\;y,\bar{y},v,\bar{v},z,\bar{z})$
$\displaystyle\Longleftrightarrow$ $\displaystyle\inf_{x}\\{\langle
c-A^{*}y-B^{*}\bar{y}-v-z,\,x\rangle\\}$ $\displaystyle=$
$\displaystyle\begin{cases}0,&\text{if }c-A^{*}y-B^{*}\bar{y}-v-z=0,\\\
-\infty,&\text{otherwise}.\end{cases}$ $\displaystyle\inf_{\bar{x}}{\cal
L}(x,\bar{x},u,\bar{u},s,\bar{s}\mathchar
24635\relax\;y,\bar{y},v,\bar{v},z,\bar{z})$
$\displaystyle\Longleftrightarrow$
$\displaystyle\inf_{\bar{x}}\\{\langle\bar{c}-\bar{B}^{*}\bar{y}-\bar{v}-\bar{z},\,\bar{x}\rangle\\}$
$\displaystyle=$ $\displaystyle\begin{cases}0,&\text{if
}\bar{B}^{*}\bar{y}+\bar{v}+\bar{z}=\bar{c},\\\
-\infty,&\text{otherwise}.\end{cases}$ $\displaystyle\inf_{u}{\cal
L}(x,\bar{x},u,\bar{u},s,\bar{s}\mathchar
24635\relax\;y,\bar{y},v,\bar{v},z,\bar{z})$
$\displaystyle\Longleftrightarrow$ $\displaystyle\inf_{u}\\{\theta(u)+\langle
v,\,u\rangle\\}=-\theta^{*}(-v)$ $\displaystyle\inf_{\bar{u}}{\cal
L}(x,\bar{x},u,\bar{u},s,\bar{s}\mathchar
24635\relax\;y,\bar{y},v,\bar{v},z,\bar{z})$
$\displaystyle\Longleftrightarrow$
$\displaystyle\inf_{\bar{u}}\\{\bar{\theta}(\bar{u})+\langle\bar{v},\,\bar{u}\rangle\\}=-\bar{\theta}^{*}(-\bar{v})$
$\displaystyle\inf_{s}{\cal L}(x,\bar{x},u,\bar{u},s,\bar{s}\mathchar
24635\relax\;y,\bar{y},v,\bar{v},z,\bar{z})$
$\displaystyle\Longleftrightarrow$ $\displaystyle\inf_{u}\\{\delta_{{\cal
K}}(s)+\langle z,\,s\rangle\\}=-\delta_{{\cal K}}^{*}(-z)$
$\displaystyle\inf_{\bar{s}}{\cal L}(x,\bar{x},u,\bar{u},s,\bar{s}\mathchar
24635\relax\;y,\bar{y},v,\bar{v},z,\bar{z})$
$\displaystyle\Longleftrightarrow$
$\displaystyle\inf_{\bar{s}}\\{\delta_{{\bar{\cal
K}}}(\bar{s})+\langle\bar{z},\,\bar{s}\rangle\\}=-\delta_{{\bar{\cal
K}}}^{*}(-\bar{z}),$

where "$\Longleftrightarrow$" means equivalence by removing the irrelevant
constants.

Hence the dual problem of (P) is given by

(D)
$\displaystyle\max_{y,\bar{y},v,\bar{v},z,\bar{z}}\;\inf_{x,\bar{x},u,\bar{u},s,\bar{s}}{\cal
L}(x,\bar{x},u,\bar{u},s,\bar{s}\mathchar
24635\relax\;y,\bar{y},v,\bar{v},z,\bar{z})$
$\displaystyle\begin{array}[]{rl}=\;\max_{y,\bar{y},v,\bar{v},z,\bar{z}}&-\theta^{*}(-v)-\bar{\theta}^{*}(-\bar{v})-\delta_{{\cal
K}}^{*}(-z)-\delta_{{\bar{\cal K}}}^{*}(-\bar{z})+\langle
b,\,y\rangle+\langle\bar{b},\,\bar{y}\rangle\\\\[5.0pt]
\mbox{s.t.}&c-A^{*}y-B^{*}\bar{y}-v-z=0,\\\\[5.0pt]
&\bar{c}-\bar{B}^{*}\bar{y}-\bar{v}-\bar{z}=0.\end{array}$

Consider the Lagrangian function for (D),

$\displaystyle{\cal L}(v,\bar{v},z,\bar{z},y,\bar{y}\mathchar
24635\relax\;x,\bar{x})$ $\displaystyle:=$
$\displaystyle\;\;\theta^{*}(-v)+\bar{\theta}^{*}(-\bar{v})+\delta^{*}_{{\cal
K}}(-z)+\delta_{\bar{\cal K}^{*}}(-\bar{z})-\langle
b,\,y\rangle-\langle\bar{b},\,\bar{y}\rangle$ $\displaystyle+\langle
x,\,c-A^{*}y-B^{*}\bar{y}-v-z\rangle+\langle\bar{x},\,\bar{c}-\bar{B}^{*}\bar{y}-\bar{v}-\bar{z}\rangle.$

Assume that the feasible regions of both the primal and dual problems have
nonempty interiors. Then the optimal solutions for both problems exist, and
they satisfy the following Karush-Kuhn-Tucker (KKT) conditions for the
problems (P) and (D):

$\displaystyle\left\\{\begin{array}[]{ccc}&Ax=b,\quad
Bx+\bar{B}\bar{x}=\bar{b},&\\\\[3.0pt]
&A^{*}y+B^{*}\bar{y}+z+v=c,\quad\bar{B}^{*}\bar{y}+\bar{z}+\bar{v}=\bar{c},&\\\\[5.0pt]
&0\in\partial\delta_{\cal K}^{*}(-z)-x\Leftrightarrow-z={\rm
Prox}_{\delta_{\cal K}^{*}}(x-z)\Leftrightarrow x=\Pi_{{\cal K}}(x-z),\\\
&0\in\partial\delta_{\bar{\cal
K}^{*}}(-\bar{z})-\bar{x}\Leftrightarrow-\bar{z}={\rm Prox}_{\delta_{\bar{\cal
K}^{*}}}(\bar{x}-\bar{z})\Leftrightarrow\bar{x}=\Pi_{\bar{\cal
K}}(\bar{x}-\bar{z}),\\\ &0\in\partial\theta^{*}(-v)-x\Leftrightarrow-v={\rm
Prox}_{\theta^{*}}(x-v)\Leftrightarrow x={\rm Prox}_{\theta}(x-v),\\\
&0\in\partial\bar{\theta}^{*}(-\bar{v})-\bar{x}\Leftrightarrow-\bar{v}={\rm
Prox}_{\bar{\theta}^{*}}(\bar{x}-\bar{v})\Leftrightarrow\bar{x}={\rm
Prox}_{\bar{\theta}}(\bar{x}-\bar{v}).\\\ \end{array}\right.$ (62)

where $\Pi_{{\cal K}}(\cdot)$ and $\Pi_{{\bar{\cal K}}}(\cdot)$ denote the
projection mappings onto ${\cal K}$ and ${\bar{\cal K}}$, respectively.

##  4 Preliminaries

In this section, we will first present the sGS decomposition theorem, which is
important for formulating the main algorithmic framework that is given in CST
. In short, this theorem helps to decompose a multi-block problem into smaller
subproblems with additional inexpensive steps but not sacrificing the
convergence guarantee. Then, we will present the abstract algorithmic template
by fitting our problem in the context of CST .

###  4.1 sGS decomposition theorem

Let ${\cal X}_{1},\ldots,{\cal X}_{s}$ ($s\geq 2$) be given finite-dimensional
inner product spaces, ${\cal Q}:{\cal X}\to{\cal X}$ be a given self-adjoint
positive semidefinite linear operator, and $c\in{\cal X}$, where ${\cal
X}={\cal X}_{1}\times\cdots\times{\cal X}_{s}.$ Suppose that ${\cal Q}$ is
partitioned according to ${\cal X}_{1}\times\cdots\times{\cal X}_{s}$ as

$\displaystyle{\cal Q}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccc}Q_{1,1}&\dots&Q_{1,s}\\\\[5.0pt]
\vdots&\vdots&\vdots\\\\[5.0pt]
Q_{1,s}^{*}&\dots&Q_{s,s}\end{array}\right]\;=\;{\cal U}+{\cal D}+{\cal
U}^{*},$ (66)

where ${\cal D}$ and ${\cal U}$ denote the block diagonal part and the
strictly upper block triangular part of ${\cal Q}$, respectively. That is,

$\displaystyle{\cal
U}=\left[\begin{array}[]{cccc}\mbox{\boldmath$0$}&Q_{1,2}&\dots&Q_{1,s}\\\
&\ddots&&\vdots\\\ &&\ddots&Q_{s-1,s}\\\\[5.0pt]
&&&\mbox{\boldmath$0$}\end{array}\right],\quad{\cal
D}=\left[\begin{array}[]{cccc}Q_{1,1}&&\\\\[3.0pt] &Q_{2,2}&&\\\ &&\ddots&\\\
&&&Q_{s,s}\end{array}\right].$ (75)

Assuming that $Q_{ii}$ is positive definite for all $i=1,\ldots,s$, then we
can define the following symmetric Gauss-Seidel (sGS) linear operator
associated with ${\cal Q}$:

$\displaystyle{\bf sGS}({\cal Q})={\cal U}{\cal D}^{-1}{\cal U}^{*}.$ (76)

For a given $x=(x_{1}\mathchar 24635\relax\;\ldots\mathchar
24635\relax\;x_{s})$, we define

$\displaystyle x_{\geq i}=(x_{i}\mathchar 24635\relax\;\ldots\mathchar
24635\relax\;x_{s}),\quad x_{\leq i}=(x_{1}\mathchar
24635\relax\;\ldots\mathchar 24635\relax\;x_{i}),\quad i=1,\ldots,s.$

We also define $x_{\geq s+1}=\emptyset$. We have the following theorem which
will be useful for us to design a decomposition method to solve (D).

######  Theorem 1.

(sGS decomposition theorem (LST-sGS, , Theorem 1)) Let $z\in{\cal X}$ be
given. Assume that ${\cal Q}\succeq 0$ and $Q_{i,i}$ are positive definite.
Suppose $\delta^{\prime},\delta\in{\cal X}$ are two given error vectors with
$\delta_{1}^{\prime}=\delta_{1}$. Define

$\displaystyle\Delta(\delta^{\prime},\delta):=\delta+{\cal U}{\cal
D}^{-1}(\delta-\delta^{\prime}).$ (77)

Let $p$ be a closed proper convex function, and $q(x)=\frac{1}{2}\langle
x,\,{\cal Q}x\rangle-\langle c,\,x\rangle\;\forall\;x\in{\cal X}.$ Consider
the following convex composite quadratic programming problem:

$\displaystyle x^{+}:=\mathop{\rm
argmin}\;\Big{\\{}p(x_{1})+q(x)+\frac{1}{2}\|x-z\|_{\cal S}^{2}-\langle
x,\,\Delta(\delta^{\prime},\delta)\rangle\;\mid\;x\in{\cal X}\Big{\\}},$ (78)

where ${\cal S}={\bf sGS}({\cal Q})$ is the sGS linear operator associated
with ${\cal Q}$. Then $x^{+}$ can be computed in the following symmetric
Gauss-Seidel fashion. For $i=s,\ldots,2,$ compute $x_{i}^{\prime}\in{\cal
X}_{i}$ defined by

$\displaystyle x^{\prime}_{i}:={}$ $\displaystyle\mathop{\rm
argmin}_{x_{i}\in{\cal X}_{i}}\;\big{\\{}p(z_{1})+q(z_{\leq i-1}\mathchar
24635\relax\;x_{i}\mathchar 24635\relax\;x^{\prime}_{\geq
i+1})-\langle\delta_{i}^{\prime},\,x_{i}\rangle\big{\\}}$ (79)
$\displaystyle={}$ $\displaystyle
Q_{i,i}^{-1}\big{(}c_{i}+\delta_{i}^{\prime}-\mbox{$\sum_{j=1}^{i-1}$}Q_{j,i}^{*}z_{j}-\mbox{$\sum_{j=i+1}^{s}$}Q_{i,j}x^{\prime}_{j}\big{)}.$

Then the optimal solution $x^{+}$ for (78) can be computed exactly via the
following steps:

$\left\\{\begin{aligned} x_{1}^{+}={}&\mathop{\rm argmin}_{x_{1}\in{\cal
X}_{1}}\;\big{\\{}p(x_{1})+q(x_{1}\mathchar 24635\relax\;x^{\prime}_{\geq
2})-\langle\delta_{1},\,x_{1}\rangle\big{\\}}\\\\[5.0pt]
x_{i}^{+}={}&\mathop{\rm argmin}_{x_{i}\in{\cal
X}_{i}}\;\big{\\{}p(x_{1}^{+})+q(x^{+}_{\leq i-1}\mathchar
24635\relax\;x_{i}\mathchar 24635\relax\;x^{\prime}_{\geq
i+1})-\langle\delta_{i},\,x_{i}\rangle\big{\\}}\\\\[5.0pt]
={}&Q_{i,i}^{-1}\big{(}c_{i}+\delta_{i}-\mbox{$\sum_{j=1}^{i-1}$}Q_{j,i}^{*}x_{j}^{+}-\mbox{$\sum_{j=i+1}^{s}$}Q_{i,j}x^{\prime}_{j}\big{)},\quad
i=2,\ldots,s.\end{aligned}\right.$ (80)

Note that the role of the error vectors $\delta^{\prime}$ and $\delta$ in the
above block sGS decomposition theorem is to allow for inexact computation in
(79) and (80). There is no need to know these error vectors in advance. We
should view the computed $x_{i}^{\prime}$ and $x_{i}^{+}$ from (79) and (80)
as approximate solutions to the minimization subproblems without the terms
involving $\delta_{i}^{\prime}$ and $\delta_{i}$. Once these approximate
solutions have been computed, the errors incurred would generate the vectors
$\delta_{i}^{\prime}$ and $\delta_{i}$ automatically. With these known error
vectors, we know that the computed approximate solutions are the exact
solutions to the slightly perturbed minimization problems in (79) and (80). In
particular, for $i=s,\ldots,2$, we know that $\delta_{i}^{\prime}$ is the
residual vector obtained when we solve the linear system of equations
$Q_{ii}x^{\prime}_{i}=c_{i}-\mbox{$\sum_{j=1}^{i-1}$}Q_{j,i}^{*}z_{j}-\mbox{$\sum_{j=i+1}^{s}$}Q_{i,j}x^{\prime}_{j}$
in (79) approximately. Similarly, $\delta_{i}$ ($i=2,\ldots,s$) is the
residual vector obtained when we solve the linear system of equations in (80).

###  4.2 An abstract inexact proximal ADMM template

To make use of the abstract inexact proximal ADMM template provided in CST ,
we first write (D) in the form that was presented in that paper:

$\displaystyle\min_{\mu\in\tilde{{\cal X}},\nu\in\tilde{{\cal
Y}}}\big{\\{}p(\mu_{1})+f(\mu_{1},\mu_{2})+q(\nu_{1})+g(\nu_{1},\nu_{2})\;\mid\;{\cal
F}^{*}\mu+{\cal G}^{*}\nu=h\big{\\}},$ (81)

where $\tilde{{\cal X}}=({\bar{\cal X}}\times{\cal X})\times{\cal Y}$,
$\tilde{{\cal Y}}=({\cal X}\times{\bar{\cal X}})\times{\bar{\cal Y}}$,

$\displaystyle\mu_{1}:=[\bar{z}\mathchar
24635\relax\;z],\quad\mu_{2}=y,\quad\nu_{1}=[v\mathchar
24635\relax\;\bar{v}],\quad\nu_{2}=\bar{y},$ $\displaystyle
p(\mu_{1}):=\delta_{{\bar{\cal K}}}^{*}(-\bar{z})+\delta_{{\cal
K}}^{*}(-z),\quad f(\mu_{1},\mu_{2}):=-\langle b,\,y\rangle,$ $\displaystyle
q(\nu_{1}):=\theta^{*}(-v)+\bar{\theta}^{*}(-\bar{v}),\quad
g(\nu_{1},\nu_{2}):=-\langle\bar{b},\,\bar{y}\rangle,$ $\displaystyle{\cal
F}^{*}:=\begin{bmatrix}0&I&{\cal A}^{*}\\\ I&0&0\end{bmatrix},\quad{\cal
G}^{*}:=\begin{bmatrix}I&0&B^{*}\\\ 0&I&\bar{B}^{*}\end{bmatrix},\quad
h:=\color[rgb]{0,0,0}\begin{bmatrix}c\\\ \bar{c}\end{bmatrix}.$

Given a penalty parameter $\sigma>0$, the augmented Lagrangian function
associated with (81) is given as follows. For
$\mu:=(\mu_{1},\mu_{2})\in({\bar{\cal X}}\times{\cal X})\times{\cal Y}$,
$\nu:=(\nu_{1},\nu_{2})\in({\cal X}\times{\bar{\cal X}})\times{\bar{\cal Y}}$,
$\omega=[x\mathchar 24635\relax\;\bar{x}]\in{\cal X}\times{\bar{\cal X}}$,

$\displaystyle{\cal L}_{\sigma}(\mu,\nu\mathchar 24635\relax\;\omega)$
$\displaystyle=$ $\displaystyle
p(\mu_{1})+f(\mu_{1},\mu_{2})+q(\nu_{1})+g(\nu_{1},\nu_{2})-\frac{1}{2\sigma}\|\omega\|^{2}$
$\displaystyle+\frac{\sigma}{2}\|{\cal F}^{*}\mu+{\cal
G}^{*}\nu-h+\sigma^{-1}\omega\|^{2}.$

Let $\partial_{\mu}{\cal L}_{\sigma}(\mu,\nu\mathchar 24635\relax\;\omega)$
and $\partial_{\nu}{\cal L}_{\sigma}(\mu,\nu\mathchar 24635\relax\;\omega)$
denote the $\mu$ and $\nu$ components of the subdifferential respectively. The
inexact proximal ADMM template for solving (81) is presented next.

Algorithm imsPADMM in CST : An abstract inexact proximal ADMM for solving
(81). Let $\tau\in(0,(1+\sqrt{5})/2)$ be the step-length and
$\\{\tilde{\epsilon}_{k}\\}_{k\geq 0}$ be a summable sequence of nonnegative
numbers. Choose an initial point $(\mu^{0},\nu^{0}\mathchar
24635\relax\;\omega^{0})\in({\rm dom}p\times{\cal Y})\times({\rm
dom}q\times{\bar{\cal Y}})\times({\cal X}\times{\bar{\cal X}})$. Let
$\tilde{{\cal S}},\tilde{{\cal T}}$ be positive semidefinite linear operators
that are to be chosen later but satisfy the property that ${\cal M}:={\cal
F}{\cal F}^{*}+\tilde{{\cal S}}\succ 0$ and ${\cal N}:={\cal G}{\cal
G}^{*}+\tilde{{\cal T}}\succ 0$. For $k=0,1,...$ perform the following steps:
Step 1. Compute $\displaystyle{\mu}^{k+1}\approx{\rm
argmin}_{\mu\in({\bar{\cal X}}\times{\cal X})\times{\cal Y}}\big{\\{}{\cal
L}_{\sigma}(\mu,\nu^{k}\mathchar
24635\relax\;\omega^{k})+\frac{\sigma}{2}\|\mu-\mu^{k}\|_{\tilde{{\cal
S}}}^{2}\big{\\}},$ such that there exists a vector $\delta^{k}$ satisfying
$\|{\cal M}^{-\frac{1}{2}}{\delta}^{k}\|\leq\tilde{\epsilon}_{k}$ and
${\delta}^{k}\in\partial_{\mu}{\cal L}_{\sigma}(\mu^{k+1},\nu^{k}\mathchar
24635\relax\;\omega^{k})+\sigma\tilde{{\cal S}}({\mu}^{k+1}-\mu^{k}).$ Step 2.
Compute $\displaystyle\nu^{k+1}\approx\mbox{argmin}_{\nu\in({\cal
X}\times{\bar{\cal X}})\times{\bar{\cal Y}}}\big{\\{}{\cal
L}_{\sigma}(\mu^{k+1},\nu\mathchar
24635\relax\;\omega^{k})+\frac{\sigma}{2}\|\nu-\nu^{k}\|_{\tilde{{\cal
T}}}^{2}\big{\\}},$ such that there exists a vector $\delta^{k}$ satisfying
$\|{\cal N}^{-\frac{1}{2}}{\gamma}^{k}\|\leq\tilde{\epsilon}_{k}$ and
$\gamma^{k}\in\partial_{\nu}{\cal L}_{\sigma}(\mu^{k+1},\nu^{k+1}\mathchar
24635\relax\;\omega^{k})+\sigma\tilde{{\cal T}}(\nu^{k+1}-\nu^{k}).$ Step 3.
Compute $\omega^{k+1}:=\omega^{k}+\tau\sigma({\cal F}^{*}\mu^{k+1}+{\cal
G}^{*}\nu^{k+1}-h)$.

The convergence of the above imsPADMM algorithm is proven in CST .

##  5 Inexact symmetric Gauss-Seidel proximal ADMM and ALM for solving (D)

We now apply the abstract inexact proximal ADMM (imsPADMM presented in the
last section) to the DBA problem. Explicitly, given a penalty parameter
$\sigma>0$, the augmented Lagrangian function associated with (D) is given as
follows: for $(\bar{z},z,y)\in{\bar{\cal X}}\times{\cal X}\times{\cal Y}$,
$(v,\bar{v},\bar{y})\in{\cal X}\times{\bar{\cal X}}\times{\bar{\cal Y}}$,
$(x,\bar{x})\in{\cal X}\times{\bar{\cal X}}$,

$\displaystyle{\cal L}_{\sigma}(\bar{z},z,y\mathchar
24635\relax\;v,\bar{v},\bar{y}\mathchar
24635\relax\;x,\bar{x})\;=\;\theta^{*}(-v)+\delta_{\cal K}^{*}(-z)-\langle
b,\,y\rangle+\bar{\theta}^{*}(-\bar{v})+\delta_{{\bar{\cal
K}}}^{*}(-\bar{z})-\langle\bar{b},\,\bar{y}\rangle$ $\displaystyle\hskip
85.35826pt-\frac{1}{2\sigma}\|x\|^{2}-\frac{1}{2\sigma}\|\bar{x}\|^{2}+\frac{\sigma}{2}\|A^{*}y+B^{*}\bar{y}+z+v-c+\sigma^{-1}x\|^{2}$
$\displaystyle\hskip
85.35826pt+\frac{\sigma}{2}\|\bar{B}^{*}\bar{y}+\bar{z}+\bar{v}-\bar{c}+\sigma^{-1}\bar{x}\|^{2}.$
(82)

Next we present the explicit algorithmic framework of our inexact sGS proximal
ADMM for solving the dual problem (D) of the block-angular problem (P). To
apply the sGS decomposition theorem within the inexact proximal ADMM
framework, we consider $(\bar{z},z,y)$ as the first group of block variables
and $(v,\bar{v},\bar{y})$ as the second group of block variables. Then we make
use of the sGS decomposition theorem in Step 1 of imsPADMM by considering the
quadratic term ${\cal F}{\cal F}^{*}$ corresponding to the block $\mu$ in
${\cal L}_{\sigma}(\mu,\nu^{k}\mathchar 24635\relax\;\omega^{k})$. Similarly,
we make use of the sGS decomposition in Step 2 of imsPADMM by considering
${\cal G}{\cal G}^{*}$ corresponding to the block $\nu$ in ${\cal
L}_{\sigma}(\mu^{k+1},\nu\mathchar 24635\relax\;\omega^{k})$.

The explicit template of the inexact sGS proximal ADMM for solving (D) is
given as follows.

Algorithm DBA-sGS-ADMM. Let $\\{\epsilon_{k}\\}$ be a given nonnegative
summable sequence, i.e. $\sum_{k=0}^{\infty}\epsilon_{k}<\infty$. Given the
initial points $(\bar{z}^{0},z^{0},y^{0})\in{\bar{\cal X}}\times{\cal
X}\times{\cal Y}$, $(v,\bar{v}^{0},\bar{y}^{0})\in{\cal X}\times{\bar{\cal
X}}\times{\bar{\cal Y}}$, $(x^{0},\bar{x}^{0})\in{\cal X}\times{\bar{\cal
X}}$, perform the following steps in each iteration.

Step 1.

Let $J:{\cal Y}\to{\cal Y}$ be a given self-adjoint positive semidefinite
linear operator, and ${\cal S}$ be a positive semidefinite linear operator
that is to be chosen later. Compute

$\displaystyle(\bar{z}^{k+1},z^{k+1},y^{k+1})\approx\mathop{\rm
argmin}\limits_{\bar{z},z,y}\left\\{\begin{array}[]{l}{\cal
L}_{\sigma}(\bar{z},z,y\mathchar
24635\relax\;v^{k},\bar{v}^{k},\bar{y}^{k}\mathchar
24635\relax\;x^{k},\bar{x}^{k})\;+\\\\[5.0pt]
\frac{\sigma}{2}\|y-y^{k}\|^{2}_{J}+\frac{\sigma}{2}\|(z\mathchar
24635\relax\;y)-(z^{k}\mathchar 24635\relax\;y^{k})\|_{\cal
S}^{2}\end{array}\Bigm{|}\begin{array}[]{l}y\in{\cal Y},\\\\[3.0pt] z\in{\cal
X},\bar{z}\in{\bar{\cal X}}\end{array}\right\\}.$

Next we explain how the above abstract subproblem can be solved in practice.
Let $R^{1,k}=A^{*}y^{k}+B^{*}\bar{y}^{k}+z^{k}+v^{k}-c^{k}$,
$R^{2,k}=\bar{B}^{*}\bar{y}^{k}+\bar{z}^{k}+\bar{v}^{k}-\bar{c}^{k}$, where
$c^{k}=c-\sigma^{-1}x^{k}$ and $\bar{c}^{k}=\bar{c}-\sigma^{-1}\bar{x}^{k}.$
Observe that $\bar{z}$ and $(z,y)$ are separable. Hence we can compute
$\bar{z}^{k+1}$ and $(z^{k+1},y^{k+1})$ separately by solving

$\displaystyle\bar{z}^{k+1}$ $\displaystyle=$ $\displaystyle\mathop{\rm
argmin}\limits_{\bar{z}}\Big{\\{}\frac{\sigma}{2}\|\bar{z}-\bar{z}^{k}+R^{2,k}\|^{2}+\delta_{{\bar{\cal
K}}}^{*}(-\bar{z})\bigm{|}\bar{z}\in{\bar{\cal X}}\Big{\\}}$ (84)
$\displaystyle(z^{k+1},y^{k+1})$ $\displaystyle\approx$
$\displaystyle\mathop{\rm
argmin}\limits_{z,y}\left\\{\begin{array}[]{l|}-\langle
b,\,y\rangle+\delta_{\cal K}^{*}(-z)\\\\[5.0pt]
+\frac{\sigma}{2}\|z+A^{*}y+B^{*}\bar{y}^{k}+v^{k}-c^{k}\|^{2}\\\\[5.0pt]
\frac{\sigma}{2}\|y-y^{k}\|^{2}_{J}+\frac{\sigma}{2}\|(z\mathchar
24635\relax\;y)-(z^{k}\mathchar 24635\relax\;y^{k})\|_{\cal
S}^{2}\end{array}\;\begin{array}[]{l}y\in{\cal Y},\\\\[3.0pt] z\in{\cal
X}\end{array}\right\\}.\qquad$ (90)

Notice that $\bar{z}^{k+1}$ can be computed by solving the following $N$
subproblems in parallel, i.e.,

$\displaystyle\bar{z}_{i}^{k+1}$ $\displaystyle=$ $\displaystyle\mathop{\rm
argmin}\limits_{\bar{z}_{i}}\Big{\\{}\frac{\sigma}{2}\|\bar{z}_{i}-\bar{z}^{k}_{i}+R^{2,k}_{i}\|^{2}+\delta_{{\cal
K}_{i}}^{*}(-\bar{z}_{i})\bigm{|}\bar{z}_{i}\in{\cal X}_{i}\Big{\\}}$
$\displaystyle=$ $\displaystyle-{\rm Prox}_{\sigma^{-1}\delta_{{\cal
K}_{i}}^{*}}\big{(}R^{2,k}_{i}-\bar{z}^{k}_{i}\big{)},\quad i=1,\ldots,N.$

In order to compute $(z^{k+1},y^{k+1})$ efficiently in a sequential manner, we
apply the sGS decomposition theorem to the problem (90). We regard $z$ and $y$
in (90) as $x_{1}$ and $x_{2}$ in the sGS decomposition theorem, respectively.
We then choose ${\cal S}={\bf sGS}({\cal Q}_{z,y})$ to be the sGS linear
operator (76) associated with the following positive semidefinite linear
operator to decompose the computation of $y^{k+1}$ and $z^{k+1}$:

$\displaystyle{\cal Q}_{z,y}:=\left[\begin{array}[]{cc}I_{n}&A^{*}\\\\[3.0pt]
A&AA^{*}+J\end{array}\right].$

In the context of imsPADMM, we have chosen $\tilde{{\cal S}}={\rm diag}(0,{\rm
diag}(0,J)+{\bf sGS}({\cal Q}_{z,y})).$ Now we can compute $(z^{k+1},y^{k+1})$
in a symmetric Gauss-Seidel fashion as follows:

$\displaystyle\begin{array}[]{lll}y^{k+1}_{\rm tmp}&=&\mathop{\rm
argmin}\limits_{y}\big{\\{}-\langle
b,\,y\rangle+\frac{\sigma}{2}\|A^{*}(y-y^{k})+R^{1,k}\|^{2}+\frac{\sigma}{2}\|y-y^{k}\|_{J}^{2}-\langle\delta^{k}_{\rm
tmp},\,y\rangle\big{\\}}\\\\[8.0pt]
&=&y^{k}+(AA^{*}+J)^{-1}\big{(}\sigma^{-1}(b+\delta^{k}_{\rm
tmp})-AR^{1,k}\big{)}\\\\[8.0pt] z^{k+1}&=&\mathop{\rm
argmin}\limits_{z}\big{\\{}\frac{\sigma}{2}\|z-z^{k}+A^{*}(y^{k+1}_{\rm
tmp}-y^{k})+R^{1,k}\|^{2}+\delta_{\cal K}^{*}(-z)\mid z\in{\cal
X}\big{\\}}\\\\[8.0pt] &=&-{\rm Prox}_{\sigma^{-1}\delta_{\cal
K}^{*}}\big{(}A^{*}(y^{k+1}_{\rm tmp}-y^{k})+R^{1,k}-z^{k}\big{)}\\\\[8.0pt]
y^{k+1}&=&\mathop{\rm argmin}\limits_{y}\big{\\{}-\langle
b,\,y\rangle+\frac{\sigma}{2}\|A^{*}(y-y^{k})+R^{1,k}+z^{k+1}-z^{k}\|^{2}+\frac{\sigma}{2}\|y-y^{k}\|_{J}^{2}-\langle\delta^{k},\,y\rangle\big{\\}}\\\\[5.0pt]
&=&y^{k}+(AA^{*}+J)^{-1}(\sigma^{-1}(b+\delta^{k})-A(R^{1,k}+z^{k+1}-z^{k}))\end{array}$

where the error vectors ${\delta}^{k}_{\rm tmp}$ and ${\delta}^{k}$ satisfy
the condition that $\max\\{\|{\delta}^{k}_{\rm
tmp}\|,\|\delta^{k}\|\\}\leq\epsilon_{k}$.

Step 2.

Let $\bar{J}$ be a given block diagonal linear operator with
$\bar{J}=\mbox{diag}(\bar{J}_{1},\ldots,\bar{J}_{N})$ such that each
$\bar{J}_{i}:{\cal Y}_{i}\to{\cal Y}_{i}$ is self-adjoint positive
semidefinite, and ${\cal T}:{\cal X}\times\bar{\cal X}\times\bar{\cal
Y}\to{\cal X}\times\bar{\cal X}\times\bar{\cal Y}$ be a positive semidefinite
linear operator that is to be chosen later. Compute

$\displaystyle(v^{k+1},\bar{v}^{k+1},\bar{y}^{k+1})\approx\mathop{\rm
argmin}\limits_{v,\bar{v},\bar{y}}\left\\{\begin{array}[]{l}{\cal
L}_{\sigma}(\bar{z}^{k+1},z^{k+1},y^{k+1}\mathchar
24635\relax\;v,\bar{v},\bar{y}\mathchar
24635\relax\;x^{k},\bar{x}^{k})\;+\\\\[5.0pt]
\frac{\sigma}{2}\|\bar{y}-\bar{y}^{k}\|_{\bar{J}}^{2}+\frac{\sigma}{2}\|(v-v^{k}\mathchar
24635\relax\;\bar{v}-\bar{v}^{k}\mathchar
24635\relax\;\bar{y}-\bar{y}^{k})\|_{\cal
T}^{2}\end{array}\Bigm{|}\begin{array}[]{l}\bar{y}\in\mathbb{R}^{\bar{m}},\\\
v\in{\cal X},\\\ \bar{v}\in{\bar{\cal X}}\end{array}\right\\}.$

Next we explain how the above abstract subproblem can be solved efficiently in
practice. Let

$\displaystyle M=BB^{*}+\bar{B}\bar{B}^{*}+\bar{J},$ (95)

and
$R^{3,k}=A^{*}y^{k+1}+B^{*}\bar{y}^{k}+z^{k+1}+v^{k}-c^{k}=R^{1,k}+A^{*}(y^{k+1}-y^{k})+(z^{k+1}-z^{k})$,
$R^{4,k}=\bar{B}^{*}\bar{y}^{k}+\bar{z}^{k+1}+\bar{v}^{k}-\bar{c}^{k}=R^{2,k}+(\bar{z}^{k+1}-\bar{z}^{k})$.
In order to compute $(v^{k+1},\bar{v}^{k+1},\bar{y}^{k+1})$ efficiently in a
sequential manner, we regard $(v,\bar{v})$ and $\bar{y}$ in the subproblem as
$x_{1}$ and $x_{2}$ in the sGS decomposition theorem, respectively. We then
choose ${\cal T}={\bf sGS}({\cal Q}_{v\bar{v},\bar{y}})$, which is the sGS
linear operator associated with the following positive semidefinite linear
operator:

$\displaystyle{\cal
Q}_{v\bar{v},\bar{y}}=\left[\begin{tabular}[]{cc|c}$I_{n}$&0&$B^{*}$\\\ {}
0&$I_{\bar{n}}$&\small{$\bar{B}^{*}$}\\\
\hline\cr$B$&\color[rgb]{0,0,0}\small{$\bar{B}$}&$M$\end{tabular}\right].$

In the context of imsPADMM, we have chosen $\tilde{{\cal T}}={\rm
diag}(0,0,\bar{J})+{\bf sGS}({\cal Q}_{v\bar{v},\bar{y}}).$ Then we can
decompose the computation of $(v^{k+1},\bar{v}^{k+1},\bar{y}^{k+1})$ in a
symmetric Gauss-Seidel fashion as follows:

$\displaystyle\begin{array}[]{l}\bar{y}^{k+1}_{\rm tmp}\;=\;\mathop{\rm
argmin}\limits_{\bar{y}}\left\\{\begin{array}[]{l}-\langle\bar{b},\,\bar{y}\rangle+\frac{\sigma}{2}\|B^{*}(\bar{y}-\bar{y}^{k})+R^{3,k}\|^{2}+\frac{\sigma}{2}\|\bar{B}^{*}(\bar{y}-\bar{y}^{k})+R^{4,k}\|^{2}\\\\[5.0pt]
+\frac{\sigma}{2}\|\bar{y}-\bar{y}^{k}\|_{\bar{J}}^{2}-\langle{\bar{\delta}^{k}_{\rm
tmp}},\,\bar{y}\rangle\end{array}\right\\}\\\\[20.0pt]
\qquad\;\;=\;\bar{y}^{k}+M^{-1}(\sigma^{-1}(\bar{b}+\bar{\delta}^{k}_{\rm
tmp})-BR^{3,k}-\bar{B}R^{4,k})\\\\[15.0pt]
\left\\{\begin{array}[]{rcl}v^{k+1}&=&\mathop{\rm
argmin}\limits_{v}\big{\\{}\theta^{*}(-v)+\frac{\sigma}{2}\|v-v^{k}+B^{*}(\bar{y}^{k+1}_{\rm
tmp}-\bar{y}^{k})+R^{3,k}\|^{2}\bigm{|}v\in{\cal X}\big{\\}}\\\\[8.0pt]
\bar{v}^{k+1}&=&\mathop{\rm
argmin}\limits_{\bar{v}}\big{\\{}\bar{\theta}^{*}(-\bar{v})+\frac{\sigma}{2}\|\bar{v}-\bar{v}^{k}+\bar{B}^{*}(\bar{y}^{k+1}_{\rm
tmp}-\bar{y}^{k})+R^{4,k}\|^{2}\bigm{|}\bar{v}\in{\bar{\cal
X}}\big{\\}}\end{array}\right.\\\\[25.0pt] \bar{y}^{k+1}\;=\;\mathop{\rm
argmin}\limits_{\bar{y}}\left\\{\begin{array}[]{l}-\langle\bar{b},\,\bar{y}\rangle+\frac{\sigma}{2}\|B^{*}(\bar{y}-\bar{y}^{k})+v^{k+1}-v^{k}+R^{3,k}\|^{2}\\\\[5.0pt]
+\frac{\sigma}{2}\|\bar{B}^{*}(\bar{y}-\bar{y}^{k})+\bar{v}^{k+1}-\bar{v}^{k}+R^{4,k}\|^{2}+\frac{\sigma}{2}\|\bar{y}-\bar{y}^{k}\|_{\bar{J}}^{2}-\langle{\bar{\delta}^{k}},\,\bar{y}\rangle\end{array}\right\\}\\\\[20.0pt]
\qquad\;\;=\;\bar{y}^{k}+M^{-1}\big{(}\sigma^{-1}(\bar{b}+\bar{\delta}^{k})-B(R^{3,k}+v^{k+1}-v^{k})-\bar{B}(R^{4,k}+\bar{v}^{k+1}-\bar{v}^{k})\big{)}\end{array}$
(108)

where the error vectors ${\bar{\delta}}^{k}_{\rm tmp}$ and
${\bar{\delta}}^{k}$ satisfy the condition that
$\max\\{\|{\bar{\delta}}^{k}_{\rm
tmp}\|,\|\bar{\delta}^{k}\|\\}\leq\epsilon_{k}$.

We should note that the computation of $\bar{v}^{k+1}$ in (108) can be done by
solving $N$ smaller subproblems in parallel, where

$\displaystyle\bar{v}^{k+1}_{i}$ $\displaystyle=$
$\displaystyle\mbox{argmin}\Big{\\{}\bar{\theta}_{i}^{*}(-\bar{v}_{i})+\frac{\sigma}{2}\|\bar{v}_{i}-\bar{v}_{i}^{k}+\bar{B}_{i}^{*}\big{(}(\bar{y}^{k+1}_{\rm
tmp})_{i}-\bar{y}_{i}^{k}\big{)}+{R^{4,k}_{i}}\|^{2}\bigm{|}\bar{v}_{i}\in{\cal
X}_{i}\Big{\\}}$ (109) $\displaystyle=$ $\displaystyle-{\rm
Prox}_{\sigma^{-1}\bar{\theta}^{*}_{i}}\Big{(}\bar{B}_{i}^{*}\big{(}(\bar{y}^{k+1}_{\rm
tmp})_{i}-\bar{y}_{i}^{k}\big{)}+R^{4,k}_{i}-\bar{v}_{i}^{k}\Big{)},\quad
i=1,\ldots,N.$

Step 3.

Compute

$\displaystyle x^{k+1}$ $\displaystyle=$ $\displaystyle
x^{k}+\tau\sigma(A^{*}y^{k+1}+B^{*}\bar{y}^{k+1}+z^{k+1}+v^{k+1}-c),$
$\displaystyle\bar{x}^{k+1}$ $\displaystyle=$
$\displaystyle\bar{x}^{k}+\tau\sigma(\bar{B}^{*}\bar{y}^{k+1}+\bar{z}^{k+1}+\bar{v}^{k+1}-\bar{c}),$

where $\tau\in(0,(1+\sqrt{5})/2)$ is the step-length that is usually set to be
$1.618$.

Let ${\cal U}_{\mu}$ and ${\cal D}_{\mu}$ be the upper block triangular part
and block diagonal part of ${\cal Q}_{z,y}$ respectively, ${\cal U}_{\nu}$ and
${\cal D}_{\nu}$ be the upper block triangular part and block diagonal part of
${\cal Q}_{v\bar{v},\bar{y}}$ respectively. Define

$\displaystyle\Delta^{k}_{\mu}=\left[\begin{array}[]{c}0\\\
\widetilde{\Delta}^{k}_{\mu}\end{array}\right],\;\widetilde{\Delta}^{k}_{\mu}=\left[\begin{array}[]{c}0\\\
\delta^{k}\end{array}\right]+{\cal U}_{\mu}{\cal
D}_{\mu}^{-1}\left(\left[\begin{array}[]{c}0\\\
\delta^{k}\end{array}\right]-\left[\begin{array}[]{c}0\\\ \delta_{\rm
tmp}^{k}\end{array}\right]\right),$ (110)
$\displaystyle\Delta_{\nu}^{k}=\left[\begin{array}[]{c}0\\\
\delta^{k}\end{array}\right]+{\cal U}_{\nu}{\cal
D}_{\nu}^{-1}\left(\left[\begin{array}[]{c}0\\\
\bar{\delta}^{k}\end{array}\right]-\left[\begin{array}[]{c}0\\\
\bar{\delta}_{\rm tmp}^{k}\end{array}\right]\right).$

Based on the above notations and the practical implementations of Step 1 and
Step 2 in DBA-sGS-ADMM, we have the following results, which is a direct
consequence of Theorem 80.

######  Lemma 1.

Using the notations in Section 4.2, let $\tilde{{\cal S}}={\rm diag}(0,{\rm
diag}(0,J)+{\bf sGS}({\cal Q}_{z,y}))$, $\tilde{{\cal T}}={\rm
diag}(0,0,\bar{J})+{\bf sGS}({\cal Q}_{v\bar{v},\bar{y}})$, suppose ${\cal
D}_{\mu}\succ 0$, ${\cal D}_{\nu}\succ 0$. The iterates
$\\{(\mu^{k},\nu^{k})\\}$ generated by DBA-sGS-ADMM is the unique solution to
the following perturbed proximal minimization problem:

$\displaystyle\mu^{k+1}$ $\displaystyle=\mathop{\rm argmin}_{\mu\in\bar{\cal
X}\times{\cal X}\times{\cal Y}}\left\\{{\cal L}_{\sigma}(\mu,\nu^{k}\mathchar
24635\relax\;\omega^{k})+\frac{\sigma}{2}\|\mu-\mu^{k}\|^{2}_{\tilde{\cal
S}}-\langle\Delta_{\mu}^{k},\,\mu\rangle\right\\}$ $\displaystyle\nu^{k+1}$
$\displaystyle=\mathop{\rm argmin}_{\nu\in{\cal X}\times\bar{\cal
X}\times\bar{\cal Y}}\left\\{{\cal L}_{\sigma}(\mu^{k+1},\nu\mathchar
24635\relax\;\omega^{k})+\frac{\sigma}{2}\|\nu-\nu^{k}\|^{2}_{\tilde{\cal
T}}-\langle\Delta_{\nu}^{k},\,\nu\rangle\right\\},$

which is equivalent to the following inclusion problem:

$\displaystyle\Delta_{\mu}^{k}$ $\displaystyle\in\partial_{\mu}{\cal
L}_{\sigma}(\mu,\nu^{k}\mathchar 24635\relax\;\omega^{k})+\sigma\tilde{\cal
S}(\mu-\mu^{k})$ $\displaystyle\Delta_{\nu}^{k}$
$\displaystyle\in\partial_{\nu}{\cal L}_{\sigma}(\mu^{k+1},\nu\mathchar
24635\relax\;\omega^{k})+\sigma\tilde{\cal T}(\nu-\nu^{k}).$

From the inexact conditions in DBA-sGS-ADMM and (110), we can obtain that
$\|\Delta^{k}_{\mu}\|$ and $\|\Delta_{\nu}^{k}\|$ are summable as long as
$\\{\epsilon_{k}\\}$ is summable.

###  5.1 An inexact symmetric Gauss-Seidel proximal ALM for solving (D)

When the functions $\theta(\cdot)$ and $\bar{\theta}_{i}(\cdot)$
($i=1,\ldots,N$) are absent in (17), we can design a symmetric Gauss-Seidel
proximal augmented Lagrangian method (sGS-ALM) to solve the problem. In this
case, the variables $v$ and $\bar{v}$ are absent in (D), and the associated
augmented Lagrangian function is given as follows: for
$\color[rgb]{0,0,0}(z,\bar{z},y,\bar{y})\in{\cal X}\times{\bar{\cal
X}}\times{\cal Y}\times{\bar{\cal Y}}$, $(x,\bar{x})\in{\cal
X}\times{\bar{\cal X}}$,

$\displaystyle{\cal L}_{\sigma}(z,\bar{z},y,\bar{y}\mathchar
24635\relax\;x,\bar{x})$ $\displaystyle=$ $\displaystyle-\langle
b,\,y\rangle-\langle\bar{b},\,\bar{y}\rangle+\delta_{\cal
K}^{*}(-z)+\delta_{\bar{\cal
K}}^{*}(-\bar{z})-\frac{1}{2\sigma}\|x\|^{2}-\frac{1}{2\sigma}\|\bar{x}\|^{2}$
(111)
$\displaystyle+\frac{\sigma}{2}\|A^{*}y+B^{*}\bar{y}+z-c+\sigma^{-1}x\|^{2}$
$\displaystyle+\frac{\sigma}{2}\|\bar{B}^{*}\bar{y}+\bar{z}-\bar{c}+\sigma^{-1}\bar{x}\|^{2}.$

Algorithm DBA-sGS-ALM. Let $\\{\epsilon_{k}\\}$ be a given nonnegative
summable sequence, i.e. $\sum_{k=0}^{\infty}\epsilon_{k}<\infty$. Given the
initial points $(y^{0},\bar{y}^{0},z^{0},\bar{z}^{0})\in{\cal
Y}\times{\bar{\cal Y}}\times{\cal X}\times{\bar{\cal X}}$,
$(x^{0},\bar{x}^{0})\in{\cal X}\times{\bar{\cal X}}$, perform the following
steps in each iteration.

Step 1.

Let $J:{\cal Y}\to{\cal Y}$ and $\bar{J}:{\bar{\cal Y}}\to{\bar{\cal Y}}$ be
given self-adjoint positive semidefinite linear operators, and ${\cal S}$ be a
positive semidefinite linear operator that is to be chosen later. Compute

$\displaystyle(\bar{z}^{k+1},z^{k+1},y^{k+1},\bar{y}^{k+1})\approx\mbox{argmin}\left\\{\begin{array}[]{l|}{\cal
L}_{\sigma}(z,\bar{z},y,\bar{y}\mathchar
24635\relax\;x^{k},\bar{x}^{k})\\\\[5.0pt]
+\frac{\sigma}{2}\|y-y^{k}\|^{2}_{J}+\frac{\sigma}{2}\|\bar{y}-\bar{y}^{k}\|^{2}_{\bar{J}}\\\\[5.0pt]
+\frac{\sigma}{2}\|(\bar{z}\mathchar 24635\relax\;z\mathchar
24635\relax\;y\mathchar 24635\relax\;\bar{y})-(\bar{z}^{k}\mathchar
24635\relax\;z^{k}\mathchar 24635\relax\;y^{k}\mathchar
24635\relax\;\bar{y}^{k})\|_{\cal
S}^{2}\end{array}\;\;\begin{array}[]{l}y\in{\cal Y},\\\ \bar{y}\in{\bar{\cal
Y}},\\\ z\in{\cal X},\\\ \bar{z}\in{\bar{\cal X}}\end{array}\right\\}.$

The inexact conditions will be clarified later. Note that the quadratic term
of the above objective function involves the following positive semidefinte
linear operator

$\displaystyle{\cal
Q}_{\bar{z}z,y,\bar{y}}=\left[\begin{tabular}[]{cc|c|c}$I_{\bar{n}}$&0&0&$\bar{B}^{*}$\\\
0&$I_{n}$&$A^{*}$&$B^{*}$\\\ \hline\cr 0&$A$&$AA^{*}+J$&$AB^{*}$\\\
\hline\cr\color[rgb]{0,0,0}\small{$\bar{B}$}&$B$&$BA^{*}$&\small{$BB^{*}+\color[rgb]{0,0,0}\bar{B}\bar{B}^{*}+\bar{J}$}\end{tabular}\right].$

In order to compute $(\bar{z}^{k+1},z^{k+1},y^{k+1},\bar{y}^{k+1})$
efficiently in a sequential manner, we regard $(\bar{z},z)$, $y$ and $\bar{y}$
in the subproblem as $x_{1}$, $x_{2}$ and $x_{3}$ in the sGS decomposition
theorem, respectively. By choosing ${\cal S}={\bf sGS}({\cal
Q}_{\bar{z}z,y,\bar{y}})$, the sGS linear operator associated with ${\cal
Q}_{\bar{z}z,y,\bar{y}}$, then we can decompose the computation of
$(\bar{z}^{k+1},z^{k+1},y^{k+1},\bar{y}^{k+1})$ in a symmetric Gauss-Seidel
fashion as follows. Let $c^{k}=c-\sigma^{-1}x^{k}$ and
$\bar{c}^{k}=\bar{c}-\sigma^{-1}\bar{x}^{k}.$ Compute

$\displaystyle\bar{y}^{k+1}_{\rm
tmp}\;=\;\mbox{argmin}\left\\{\begin{array}[]{l}-\langle\bar{b},\,\bar{y}\rangle+\frac{\sigma}{2}\|B^{*}\bar{y}+A^{*}y^{k}+z^{k}-c^{k}\|^{2}+\frac{\sigma}{2}\|\bar{B}^{*}\bar{y}+\bar{z}^{k}-\bar{c}^{k}\|^{2}\\\\[3.0pt]
+\frac{\sigma}{2}\|\bar{y}-\bar{y}^{k}\|_{\bar{J}}^{2}-\langle\bar{\delta}^{k}_{\rm
tmp},\,\bar{y}\rangle\end{array}\right\\}$ $\displaystyle y^{k+1}_{\rm
tmp}\;=\;\mbox{argmin}\big{\\{}-\langle
b,\,y\rangle+\frac{\sigma}{2}\|A^{*}y+B^{*}\bar{y}^{k+1}_{\rm
tmp}+z^{k}-c^{k}\|^{2}+\frac{\sigma}{2}\|y-y^{k}\|_{J}^{2}-\langle\delta^{k}_{\rm
tmp},\,y\rangle\big{\\}}$
$\displaystyle\left\\{\begin{array}[]{rcl}z^{k+1}&=&\mbox{argmin}\big{\\{}\frac{\sigma}{2}\|z+A^{*}y^{k+1}_{\rm
tmp}+B^{*}\bar{y}_{\rm tmp}^{k+1}-c^{k}\|^{2}+\delta_{\cal K}^{*}(-z)\mid
z\in{\cal X}\big{\\}}\\\\[8.0pt] &=&-{\rm Prox}_{\sigma^{-1}\delta_{\cal
K}^{*}}\big{(}A^{*}y^{k+1}_{\rm tmp}+B^{*}\bar{y}_{\rm
tmp}^{k+1}-c^{k}\big{)}\\\\[8.0pt]
\bar{z}^{k+1}_{i}&=&\mbox{argmin}\big{\\{}\frac{\sigma}{2}\|\bar{z}_{i}+\bar{B}_{i}^{*}(\bar{y}_{\rm
tmp}^{k+1})_{i}-\bar{c}^{k}_{i}\|^{2}+\delta_{{\cal
K}_{i}}^{*}(-\bar{z}_{i})\mid\bar{z}_{i}\in{\cal X}_{i}\big{\\}}\\\\[5.0pt]
&=&-{\rm Prox}_{\sigma^{-1}\delta_{{\cal
K}_{i}}^{*}}\big{(}\bar{B}_{i}^{*}(\bar{y}_{\rm
tmp}^{k+1})_{i}-\bar{c}^{k}_{i}\big{)}\quad i=1,\ldots,N\end{array}\right.$
$\displaystyle y^{k+1}\;=\;\mbox{argmin}\big{\\{}-\langle
b,\,y\rangle+\frac{\sigma}{2}\|A^{*}y+B^{*}\bar{y}^{k+1}_{\rm
tmp}+z^{k+1}-c^{k}\|^{2}+\frac{\sigma}{2}\|y-y^{k}\|_{J}^{2}-\langle\delta^{k},\,y\rangle\big{\\}}$
$\displaystyle\bar{y}^{k+1}\;=\;\mbox{argmin}\left\\{\begin{array}[]{l}-\langle\bar{b},\,\bar{y}\rangle+\frac{\sigma}{2}\|B^{*}\bar{y}+A^{*}y^{k+1}+z^{k+1}-c^{k}\|^{2}\\\\[5.0pt]
+\frac{\sigma}{2}\|\bar{B}^{*}\bar{y}+\bar{z}^{k+1}-\bar{c}^{k}\|^{2}+\frac{\sigma}{2}\|\bar{y}-\bar{y}^{k}\|_{\bar{J}}^{2}-\langle\bar{\delta}^{k},\,\bar{y}\rangle\end{array}\right\\}$

where the error vectors $\delta_{\rm tmp}^{k}$, $\bar{\delta}_{\rm tmp}^{k}$,
$\delta^{k}$, $\bar{\delta}^{k}$ satisfy $\max\\{\|\delta_{\rm
tmp}^{k}\|,\|\bar{\delta}_{\rm
tmp}^{k}\|,\|\delta^{k}\|,\|\bar{\delta}^{k}\|\\}\leq\epsilon_{k}$.

Step 2.

Compute

$\displaystyle x^{k+1}$ $\displaystyle=$ $\displaystyle
x^{k}+\tau\sigma(A^{*}y^{k+1}+B^{*}\bar{y}^{k+1}+z^{k+1}-c),$
$\displaystyle\bar{x}^{k+1}$ $\displaystyle=$
$\displaystyle\bar{x}^{k}+\tau\sigma(\bar{B}^{*}\bar{y}^{k+1}+\bar{z}^{k+1}-\bar{c}),$

where $\tau\in(0,2)$ is the step-length which is usually set to be $1.9$.

Let ${\cal U}$ and ${\cal D}$ be the upper block triangular part and the block
diagonal part of ${\cal Q}_{\bar{z}z,y,\bar{y}}$, respectively, and

$\Delta^{k}:=\left[\begin{array}[]{c}0\\\ \delta^{k}\\\
\bar{\delta}^{k}\end{array}\right]+{\cal U}{\cal
D}^{-1}\left(\left[\begin{array}[]{c}0\\\ \delta^{k}\\\
\bar{\delta}^{k}\end{array}\right]-\left[\begin{array}[]{c}0\\\ \delta_{\rm
tmp}^{k}\\\ \bar{\delta}_{\rm tmp}^{k}\end{array}\right]\right).$ (117)

The summable condition for $\epsilon_{k}$ and the inexact condition:
$\max\\{\|\delta_{\rm tmp}^{k}\|,\|\bar{\delta}_{\rm
tmp}^{k}\|,\|\delta^{k}\|,\|\bar{\delta}^{k}\|\\}\leq\epsilon_{k}$ in Step 1
can guarantee that the error $\|\Delta^{k}\|$ is also summable, since

$\|\Delta^{k}\|\leq\|\delta^{k}\|+\|\bar{\delta}^{k}\|+\|{\cal U}\|\|{\cal
D}^{-1}\|\left(\|\delta^{k}\|+\|\bar{\delta}^{k}\|+\|\delta^{k}_{\rm
tmp}\|+\|\bar{\delta}^{k}_{\rm tmp}\|\right),$

where the norms $\|{\cal U}\|$ and $\|{\cal D}^{-1}\|$ here are induced
operator norms defined in Section 1. Based on the above notations and the
practical implementation of Step 1, we have the following results, which is a
direct consequence of Theorem 80.

######  Lemma 2.

Suppose ${\cal D}\succ 0$, let
$w^{k+1}:=(z^{k+1},\bar{z}^{k+1},y^{k+1},\bar{y}^{k+1})$ and $\tilde{\cal
S}=\mbox{diag}(0,0,J,\bar{J})+{\bf sGS}({\cal Q}_{\bar{z}z,y,\bar{y}})$. The
iterates $w^{k+1}$ generated by Step 1 in DBA-sGS-ALM is the unique solution
to the following perturbed proximal minimization problem:

$w^{k+1}=\mathop{\rm argmin}_{w\in{\cal X}\times\bar{\cal X}\times{\cal
Y}\times\bar{\cal Y}}\left\\{{\cal L}_{\sigma}(w\mathchar
24635\relax\;x^{k},\bar{x}^{k})+\frac{\sigma}{2}\|w-w^{k}\|_{\tilde{\cal
S}}^{2}-\langle\Delta^{k},\,w\rangle\right\\},$

which is equivalent to the following inclusion problem

$\Delta^{k}\in\partial_{w}{\cal L}_{\sigma}(w\mathchar
24635\relax\;x^{k},\bar{x}^{k})+\sigma\tilde{\cal S}(w-w^{k}).$

###  5.2 Convergence of the inexact DBA-sGS-ADMM and DBA-sGS-ALM algorithms

The convergence theorem of the inexact DBA-sGS-ADMM and DBA-sGS-ALM algorithms
could be established directly using known results from CST , CLST  and LW-
imiPADMM . Here we would state the theorems again for the convenience of the
readers.

In order to understand the theorem, first we need some definitions.

######  Definition 1.

Let $\mathbb{F}:{\cal X}\rightrightarrows{\cal Y}$ be a multivalued mapping
and denote its inverse by $\mathbb{F}^{-1}$. The graph of the multivalued
function $\mathbb{F}$ is defined by ${\rm gph}\mathbb{F}:=\\{(x,y)\in{\cal
X}\times{\cal Y}\mid y\in\mathbb{F}(x)\\}$.

Denote $u:=(y,\bar{y},v,\bar{v},z,\bar{z},x,\bar{x})\in\mathcal{U}:={\cal
Y}\times{\bar{\cal Y}}\times{\cal X}\times{\bar{\cal X}}\times{\cal
X}\times{\bar{\cal X}}\times{\cal X}\times{\bar{\cal X}}$. The KKT mapping
$\mathcal{R}:\mathcal{U}\rightarrow\mathcal{U}$ of (17) is defined by

$\displaystyle\mathcal{R}(u):=\begin{pmatrix}Ax-b\\\
Bx+\bar{B}\bar{x}-\bar{b}\\\ A^{*}y+B^{*}\bar{y}+z+v-c\\\
\bar{B}^{*}\bar{y}+\bar{z}+\bar{v}-\bar{c}\\\ x-\Pi_{{\cal K}}(x-z)\\\
\bar{x}-\Pi_{{\bar{\cal K}}}(\bar{x}-\bar{z})\\\ x-{\rm Prox}_{\theta}(x-v)\\\
\bar{x}-{\rm Prox}_{\bar{\theta}}(\bar{x}-\bar{v})\end{pmatrix}.$ (118)

Denote the set of KKT points by $\bar{\Omega}$. The KKT mapping $\mathcal{R}$
is said to be metrically subregular at $(\bar{u},0)\in{\rm gph}\mathcal{R}$
with modulus $\eta>0$ if there exists a scalar $\rho>0$ such that

${\rm dist}(u,\bar{\Omega})\leq\eta\|\mathcal{R}(u)\|\;\;\forall
u\in\\{u\in\mathcal{U}:\|u-\bar{u}\|\leq\rho\\}.$

Now we are ready to present the convergence theorem of DBA-sGS-ADMM.

######  Theorem 2.

Let $\Delta_{\mu}^{k}$ and $\Delta_{\nu}^{k}$ be defined as (110), ${\cal
M}:={\cal F}{\cal F}^{*}+\tilde{{\cal S}}\succ 0$ and ${\cal N}:={\cal G}{\cal
G}^{*}+\tilde{{\cal T}}\succ 0$,
$\\{u^{k}:=(y^{k},\bar{y}^{k},v^{k},\bar{v}^{k},z^{k},\bar{z}^{k},x^{k},\bar{x}^{k})\\}$
be the sequence generated by DBA-sGS-ADMM, where $\tilde{{\cal S}}={\rm
diag}(0,{\rm diag}(0,J)+{\bf sGS}({\cal Q}_{z,y}))$, $\tilde{{\cal T}}={\rm
diag}(0,0,\bar{J})+{\bf sGS}({\cal Q}_{v\bar{v},\bar{y}})$. Assume the KKT
point set $\bar{\Omega}$ is nonempty. Then, we have the following results.

  1. (a)

The sequence $\\{(y^{k},\bar{y}^{k},v^{k},\bar{v}^{k},z^{k},\bar{z}^{k})\\}$
converges to an optimal solution of the dual problem (D) of (17) and the
sequence $\\{(x^{k},\bar{x}^{k})\\}$ converges to an optimal solution of the
primal problem (P) of (17).

  2. (b)

Suppose that the sequence $\\{u^{k}\\}$ converges to a KKT point
$\hat{u}:=(\hat{y},\hat{\bar{y}},\hat{v},\hat{\bar{v}},\hat{z},\hat{\bar{z}},\hat{x},\hat{\bar{x}})$
and the KKT mapping $\mathcal{R}$ is metrically subregular at
$(\hat{u},0)\in{\rm gph}\mathcal{R}$. If a relative error criterion is
employed instead:

$\displaystyle\|\Delta_{\mu}^{k}\|_{{\cal M}^{-1}}$
$\displaystyle\leq\eta_{k}\left(\|\mu^{k}-\mu^{k+1}\|^{2}_{\cal
M}+\|\mu^{k-1}-\mu^{k}\|_{\cal M}^{2}\right)^{\frac{1}{2}},$
$\displaystyle\|\Delta_{\nu}^{k}\|_{{\cal N}^{-1}}$
$\displaystyle\leq\eta_{k}\left(\|\nu^{k}-\nu^{k+1}\|^{2}_{\cal
N}+\|\nu^{k-1}-\nu^{k}\|_{\cal N}^{2}\right)^{\frac{1}{2}}.$

and $\eta_{k}\in(0,1)$, $\sum_{k\geq 0}\eta_{k}^{2}<\infty$, then the sequence
$\\{u^{k}\\}$ is linearly convergent to $\hat{u}$.

Proof.  The result of (a) follows directly by applying the convergence result
in (CST, , Theorem 5.1) to (D). The result of (b) follows directly from (LW-
imiPADMM, , Theorem 1). $\square$

######  Remark 2.

By Theorem 1 and Remark 1 in LSTc , we know that when (P) is a convex
programming problem where for each $i=1,\ldots,N$, $\theta$ and
$\bar{\theta}_{i}$ are piecewise linear-quadratic or strongly convex, and
${\cal K},{\cal K}_{i}$ are polyhedral, then $\cal R$ is metrically subregular
at $(\bar{u},0)\in{\rm gph}\mathcal{R}$ for any KKT point $\bar{u}.$

Implied by Lemma 2, DBA-sGS-ALM is equivalent to an inexact majorized
indefinite-proximal ALM. Therefore, the convergence of the DBA-sGS-ALM can be
established readily by using known results in (CLST, , Theorem 3.1 and Theorem
3.2). We state the convergence theorem here for the convenience of the
readers.

######  Theorem 3.

Let $\Delta^{k}$ be defined as (117), $\\{\eta^{k}\\}$ be a given sequence of
nonnegative numbers that converges to 0,
$\\{u^{k}:=(y^{k},\bar{y}^{k},z^{k},\bar{z}^{k},x^{k},\bar{x}^{k})\\}$ be the
sequence generated by DBA-sGS-ALM. Assume the KKT point set is nonempty. Then,
we have the following results.

  1. (a)

The sequence $\\{(y^{k},\bar{y}^{k},z^{k},\bar{z}^{k})\\}$ converges to an
optimal solution of the dual problem (D) of (17) and the sequence
$\\{(x^{k},\bar{x}^{k})\\}$ converges to an optimal solution of the primal
problem (P) of (17).

  2. (b)

Suppose that the sequence $\\{u^{k}\\}$ converges to a KKT point
$\hat{u}:=(\hat{y},\hat{\bar{y}},\hat{z},\hat{\bar{z}},\hat{x},\hat{\bar{x}})$
and the KKT mapping $\mathcal{R}$ is metrically subregular at
$(\hat{u},0)\in{\rm gph}\mathcal{R}$. If a relative error criterion is
employed instead:

$\|\Delta^{k}\|\leq\eta_{k}\|u^{k}-u^{k+1}\|,$

then the sequence $\\{u^{k}\\}$ is locally linearly convergent to $\hat{u}$.

###  5.3 The efficient computation of $\bar{y}^{k+1}_{\rm tmp}$ and
$\bar{y}^{k+1}$

As $A$ is a single block matrix, the computation of $y^{k+1}$ and $y_{\rm
tmp}^{k+1}$ involve solving the linear system of the form $(AA^{*}+J)y=h$,
which is usually not an issue. For example, if computing the Cholesky
factorization of $AA^{*}$ can be done cheaply, we can simply choose $J=0$. If
factorizing $AA^{*}$ is expensive, we can choose
$J=\lambda_{\max}(AA^{*})I-AA^{*}$ to efficiently compute $y$ at the expense
of perturbing the augmented Lagrangian function by the non-trivial proximal
term $\frac{\sigma}{2}\|y-y^{k}\|_{J}^{2}$.

On the other hand, since $B$ and $\bar{B}$ are multiblock matrices, to
implement DBA-sGS-ALM or DBA-sGS-ADMM efficiently, it is necessary to
carefully address the efficient computation of $\bar{y}^{k+1}_{\rm tmp}$ and
$\bar{y}^{k+1}$. To compute the solution $\bar{y}^{k+1}_{\rm tmp}$ in (108),
we need to solve a generally very large $\bar{m}\times\bar{m}$ system of
linear equations of the form:

$\displaystyle M\bar{y}=h,$ (119)

where the $\bar{m}\times\bar{m}$ matrix $M$ given in (95) has the following
structure:

$\displaystyle M=\left[\begin{array}[]{c}B_{1}\\\ \vdots\\\
B_{N}\end{array}\right]\left[\begin{array}[]{c}B_{1}\\\ \vdots\\\
B_{N}\end{array}\right]^{*}+\left[\begin{array}[]{ccc}\bar{B}_{1}\bar{B}_{1}^{*}&&\\\
&\ddots&\\\ &&\bar{B}_{N}\bar{B}_{N}^{*}\end{array}\right]+\bar{J}.$ (129)

Here we design several numerical strategies which are crucial for cutting down
the computational cost of solving (119).

(a)

Suppose we choose $\bar{J}=\mbox{diag}(\bar{J}_{1},\ldots,\bar{J}_{N})$ where
each $\bar{J}_{i}:{\cal X}_{i}\to{\cal X}_{i}$ is a self-adjoint positive
semidefinite linear operator. To compute the inverse of $M$ in (129), we may
make use of the Sherman-Morrison-Woodbury (SMW) formula which we shall
describe next. Let $\bar{D}=\mbox{diag}(\bar{D}_{1},\dots,\bar{D}_{N})$, where
$\bar{D}_{i}=\bar{B}_{i}\bar{B}_{i}^{*}+\bar{J}_{i}$, $i=1,\ldots,N.$ Assume
that $\bar{D}_{i}$ is invertible for all $i=1,\ldots,N$, then we have that

$\displaystyle
M^{-1}h=(\bar{D}+BB^{*})^{-1}h=\bar{D}^{-1}h-\bar{D}^{-1}BG^{-1}B^{*}\bar{D}^{-1}h,$

where $G:{\cal X}\to{\cal X}$ is given by

$\displaystyle G:=I_{n}+\sum_{i=1}^{N}B_{i}^{*}\bar{D}_{i}^{-1}B_{i}.$ (130)

Thus to solve (119), we only need to solve a linear system of the form
$Gx=g:=B^{*}\bar{D}^{-1}h.$

Observe that to compute $G$, the component matrices
$B_{i}^{*}\bar{D}_{i}^{-1}B_{i}$, $i=1,\ldots,N$, can be computed in parallel,
and operations such as $\bar{D}^{-1}_{i}h_{i}$, $i=1,\ldots,N$, can also be
done in parallel. Similarly, operations such as evaluating $Bx$ and
$B^{*}\bar{y}$ for given $x\in{\cal X}$ and $\bar{y}\in{\bar{\cal Y}}$ can
also be done in parallel.

(i) Generally, one would expect the $n\times n$ matrix
$B_{i}^{*}\bar{D}_{i}^{-1}B_{i}$ to be dense even if $B_{i}$ is sparse, unless
$\bar{D}_{i}$ has a nice structure such as being a diagonal matrix. If the
memory required to store the matrix $G$ and its Cholesky factor is
prohibitive, then one may choose $\bar{J}_{i}$ such that $\bar{D}_{i}$ is a
diagonal matrix. For example, we may choose

$\bar{J}_{i}=\lambda_{\max}(\bar{B}_{i}\bar{B}_{i}^{*})I_{m_{i}}-\bar{B}_{i}\bar{B}_{i}^{*},\quad
i=1,\ldots,N.$

Then $\bar{D}_{i}=\lambda_{\max}(\bar{B}_{i}\bar{B}_{i}^{*})I_{m_{i}}$. In
which case,
$G=I_{n}+\sum_{i=1}^{N}\lambda_{\max}(\bar{B}_{i}\bar{B}_{i}^{*})^{-1}B_{i}^{*}B_{i}$,
and it is more likely to be sparse. As an example, consider the case where all
$B_{i}$ are identical and $B_{i}^{*}B_{i}$ is sparse. In this scenario, $G$
will be sparse as well. However, if $\bar{J}_{i}=0$ and
$(\bar{B}_{i}\bar{B}^{*}_{i})^{-1}$ is dense, then $G$ may become dense even
if $B_{i}^{*}B_{i}$ is sparse. This technique of attempting to obtain a
sparser $G$ also explains the rationale for introducing the linear operator
$\bar{J}$ in Step 2 of DBA-sGS-ADMM and Step 1 of DBA-sGS-ALM.

(ii) In general, we can expect the symmetric positive definite matrix $G$ in
(130) to be well conditioned, especially if for each $i$, we choose
$\bar{J}_{i}=\alpha_{i}I_{m_{i}}$ for a sufficiently large positive scalar
$\alpha_{i}$. In this case, a preconditioned conjugate gradient (PCG) method
applied to solve $Gx=g$ is expected to converge in a small number of steps. In
addition, for a given $x$, the matrix-vector product $Gx$ can also be computed
very efficiently since the component vectors,
$\bar{B}_{i}^{*}\bar{D}_{i}^{-1}\bar{B}_{i}x$, $i=1,\ldots,N$, can be computed
in parallel.

(b)

It is also possible to solve (119) in a parallel fashion if we choose
$\bar{J}$ appropriately. For example, if we choose

$\displaystyle\bar{J}=\mbox{diag}(B_{i}B_{i}^{*}+\mbox{$\sum_{j=1,j\not=i}^{N}$}\|B_{i}B_{j}^{*}\|_{2}\,I_{m_{i}}\mid
i=1,\ldots,N)-BB^{*}\;\succeq\;0,$ (131)

then

$\displaystyle M=\mbox{diag}(\bar{E},\ldots,\bar{E}_{N}),$ (132)

where
$\bar{E}_{i}=\bar{B}_{i}\bar{B}_{i}^{*}+B_{i}B_{i}^{*}+\mbox{$\sum_{j=1,j\not=i}^{N}$}\|B_{i}B_{j}^{*}\|_{2}I_{m_{i}},$
$i=1,\ldots,N$. Then the solution $\bar{y}$ in (119) can be computed in
parallel by solving

$\displaystyle\bar{E}_{i}\bar{y}_{i}=h_{i},\quad i=1,\ldots,N.$

We should mention that in the literature, for example, LST-PBA  and references
therein, a standard way to decompose $M$ into a block diagonal form is to
choose $\bar{J}$ to be:

$\displaystyle\bar{J}^{\rm
std}=(N+1)\mbox{diag}(B_{1}B_{1}^{*},\ldots,B_{N}B_{N}^{*})-BB^{*}\succeq 0,$

where the positive semidefiniteness of $\bar{J}^{\rm std}$ can be shown by
using the inequality that $\|\sum_{i=1}^{N}\bar{x}_{i}\|^{2}\leq
N\sum_{i=1}^{N}\|\bar{x}_{i}\|^{2}$ for any $\bar{x}\in{\bar{\cal X}}$.
However, such a choice is usually more conservative than the one in (131),
especially when $\|B_{i}B_{j}^{*}\|_{2}\ll 1$ for many pairs in $\\{(i,j)\mid
1\leq i<j\leq N\\}$.

(c)

We should also take special consideration when the optimization problem (17)
has the property that $B_{1}=B_{2}=\cdots=B_{N}$. In this case

$G=I_{n}+B_{1}^{*}\big{(}\mbox{$\sum_{i=1}^{N}$}\bar{D}_{i}^{-1}\big{)}B_{1}.$

Thus if $m_{1}$ is not too large, one can precompute the $m_{1}\times m_{1}$
matrix $\sum_{i=1}^{N}\bar{D}_{i}^{-1}$ once at the beginning of the algorithm
so that the matrix-vector product $Gx$ can be evaluated very efficiently for
any given $x\in{\cal X}.$ If in addition, we also have that
$\bar{B}_{1}=\cdots=\bar{B}_{N}$ such as in the case of the DNN relaxation
(45) of the UFL problem, then we get

$G=I_{n}+NB_{1}^{*}\bar{D}_{1}^{-1}B_{1}.$

In this case, obviously there is tremendous saving in the computation of the
matrix-vector product $Gx$ for any given $x$.

##  6 Incorporating a semismooth Newton method in DBA-sSG-ADMM and DBA-sGS-ALM
for solving (D)

(a)

Recall that in computing $(z^{k+1},y^{k+1})$ in (90) for Algorithm DBA-sGS-
ADMM, we have added a proximal term based on the sGS linear operator ${\bf
sGS}({\cal Q}_{z,y})$ to decouple the computation of $z^{k+1}$ and $y^{k+1}$.
We should note that while the computation of the individual $z^{k+1}$ and
$y^{k+1}$ has been made easier, the negative effect of adding a possibly large
proximal term is that the overall algorithm may take more iterations to
converge to a required level of accuracy. Thus in practice, a judicious choice
must always be made to balance the two competing factors. In particular, if an
efficient solver is available to compute $(z^{k+1},y^{k+1})$ simultaneously in
(90) without the need to add a proximal term, i.e., $J=0$ and ${\cal S}=0$,
then one should always adopt this option. In this case, we have that

$\displaystyle(z^{k+1},y^{k+1})=\mathop{\rm
argmin}\limits_{z,y}\left\\{\begin{array}[]{l}-\langle
b,\,y\rangle+\delta_{\cal
K}^{*}(-z)-{\langle\delta^{k},\,(z,y)\rangle}\\\\[5.0pt]
+\frac{\sigma}{2}\|A^{*}y+z+B^{*}\bar{y}^{k}+v^{k}-c^{k}\|^{2}\end{array}\Bigm{|}\;\begin{array}[]{l}y\in{\cal
Y},\\\\[3.0pt] z\in{\cal X}\end{array}\right\\}.\quad$ (137)

In Step 1 of Algorithm DBA-sGS-ALM, to compute
$(\bar{z}^{k+1},z^{k+1},y^{k+1},\bar{y}^{k+1})$ efficiently in sequential
manner, we regard $(\bar{z},z,y)$ and $\bar{y}$ in (137) as $x_{1}$ and
$x_{2}$ in the sGS decomposition theorem, respectively. By choosing $J=0$, and
${\cal S}={\bf sGS}({\cal Q}_{\bar{z}zy,\bar{y}})$, the sGS linear operator
associated with the following positive semidefinite linear operator,

$\displaystyle{\cal
Q}_{\bar{z}zy,\bar{y}}=\left[\begin{tabular}[]{ccc|c}$I_{\bar{n}}$&0&0&$\bar{B}^{*}$\\\
0&$I_{n}$&$A^{*}$&$B^{*}$\\\ 0&$A$&$AA^{*}$&$AB^{*}$\\\
\hline\cr\color[rgb]{0,0,0}\small{$\bar{B}$}&$B$&$BA^{*}$&\small{$BB^{*}+\color[rgb]{0,0,0}\bar{B}\bar{B}^{*}+\bar{J}$}\end{tabular}\right],$

then we can compute $(\bar{z}^{k+1},z^{k+1},y^{k+1},\bar{y}^{k+1})$ in a
symmetric Gauss-Seidel fashion as follows. Let $c^{k}=c-\sigma^{-1}x^{k}$ and
$\bar{c}^{k}=\bar{c}-\sigma^{-1}\bar{x}^{k}.$ Compute

$\displaystyle\bar{y}^{k+1}_{\rm tmp}\;=\;\mathop{\rm
argmin}\limits_{\bar{y}}\left\\{\begin{array}[]{l}-\langle\bar{b},\,\bar{y}\rangle+\frac{\sigma}{2}\|B^{*}\bar{y}+A^{*}y^{k}+z^{k}-c^{k}\|^{2}+\frac{\sigma}{2}\|\bar{B}^{*}\bar{y}+\bar{z}^{k}-\bar{c}^{k}\|^{2}\\\\[5.0pt]
+\frac{\sigma}{2}\|\bar{y}-\bar{y}^{k}\|_{\bar{J}}^{2}-\langle\bar{\delta}^{k}_{\rm
tmp},\,\bar{y}\rangle\end{array}\right\\}$ (141)
$\displaystyle\left\\{\begin{array}[]{l}(z^{k+1},y^{k+1})\;=\;\mathop{\rm
argmin}\limits_{z,y}\left\\{\begin{array}[]{l}-\langle
b,\,y\rangle+\delta_{\cal K}^{*}(-z)+\\\\[5.0pt]
\frac{\sigma}{2}\|z+A^{*}y+B^{*}\bar{y}_{\rm
tmp}^{k+1}-c^{k}\|^{2}-{\langle\delta^{k},\,(z,y)\rangle}\end{array}\Bigm{|}\;\begin{array}[]{l}y\in{\cal
Y},\\\ z\in{\cal X}\end{array}\right\\}\\\\[18.0pt]
\bar{z}^{k+1}\;=\;\mathop{\rm
argmin}\limits_{\bar{z}}\big{\\{}\frac{\sigma}{2}\|\bar{z}+\bar{B}^{*}\bar{y}_{\rm
tmp}^{k+1}-\bar{c}^{k}\|^{2}+\delta_{{\bar{\cal
K}}}^{*}(-\bar{z})\mid\bar{z}\in{\bar{\cal
X}}\big{\\}}\end{array}\right.\qquad$ (148)
$\displaystyle\bar{y}^{k+1}\;=\;\mathop{\rm
argmin}\limits_{\bar{y}}\left\\{\begin{array}[]{l}-\langle\bar{b},\,\bar{y}\rangle+\frac{\sigma}{2}\|B^{*}\bar{y}+A^{*}y^{k+1}+z^{k+1}-c^{k}\|^{2}+\frac{\sigma}{2}\|\bar{B}^{*}\bar{y}+\bar{z}^{k+1}-\bar{c}^{k}\|^{2}\\\\[5.0pt]
+\frac{\sigma}{2}\|\bar{y}-\bar{y}^{k}\|_{\bar{J}}^{2}-\langle\bar{\delta}^{k},\,\bar{y}\rangle\end{array}\right\\}.$
(151)

In both (137) and (148), we need to solve a problem of the form:

$\displaystyle(z^{k+1},y^{k+1})\;=\;\mathop{\rm
argmin}\limits_{z,y}\Big{\\{}-\langle b,\,y\rangle+\delta_{\cal
K}^{*}(-z)+\frac{\sigma}{2}\|z+A^{*}y-\widehat{c}^{k}\|^{2}-{\langle\delta^{k},\,(z,y)\rangle}\Bigm{|}\begin{array}[]{l}y\in{\cal
Y},\\\ z\in{\cal X}\end{array}\Big{\\}}\quad$ (154)

for some given $\widehat{c}^{k}$. By making use of the projection onto ${\cal
K}$, we can first compute $y^{k+1}$ from the following minimization problem,
and obtain $z^{k+1}$ once $y^{k+1}$ is computed as follows:

$\displaystyle\left\\{\begin{array}[]{rcl}y^{k+1}&\approx&\mathop{\rm
argmin}\limits_{y}\left\\{\begin{array}[]{l}-\langle b,\,y\rangle+\sigma
M_{\delta_{{\cal
K}/\sigma}^{*}}(A^{*}y-\widehat{c}^{k})-\langle\delta_{y}^{k},\,y\rangle\end{array}\bigm{|}\;y\in{\cal
Y}\right\\}\\\\[5.0pt] z^{k+1}&=&\frac{1}{\sigma}\Pi_{\cal
K}(\sigma(A^{*}y^{k+1}-\hat{c}^{k}))-(A^{*}y^{k+1}-\hat{c}^{k}).\end{array}\right.$
(158)

The detail of deriving (158) can be found in Appendix 2. In the interest of
keeping the paper to be within a reasonable length, we do not delve into the
details on how one can compute $y^{k+1}$ via a semismooth Newton method,
interested readers can refer to, for example, LSTc  for a detailed explanation
on how to use a semismooth Newton method.

(b)

In the second step in DBA-sGS-ADMM, if we further remove the proximal term
with respect to $(v,\bar{v},\bar{y})$, i.e. $\bar{\cal J}=0,{\cal T}=0$, we
compute

$\displaystyle(v^{k+1},\bar{v}^{k+1},\bar{y}^{k+1})=\mathop{\rm
argmin}\limits_{v,\bar{v},\bar{y}}\left\\{\begin{array}[]{l}\theta^{*}(-v)+\bar{\theta}^{*}(-\bar{v})-\langle\bar{b},\,\bar{y}\rangle+\frac{\sigma}{2}\|B^{*}\bar{y}+z^{k+1}+\bar{v}-\bar{c}^{k}\|^{2}\\\\[3.0pt]
+\frac{\sigma}{2}\|A^{*}y^{k+1}+B^{*}\bar{y}+z^{k+1}+v-c^{k}\|^{2}-\langle\bar{\delta}^{k},\,\bar{y}\rangle\end{array}\right\\}.\quad$
(161)

Let $R^{5,k}=A^{*}y^{k+1}+z^{k+1}-\hat{c}^{k}$,
$R^{6,k}=\bar{z}^{k+1}-\hat{\bar{c}}^{k}$, then solution of (161) can be
obtained as follows:

$\displaystyle\left\\{\begin{array}[]{rcl}\bar{y}^{k+1}&\approx&\mathop{\rm
argmin}\limits_{\bar{y}\in\bar{\cal
Y}}\left\\{\begin{array}[]{l}-\langle\bar{b},\,\bar{y}\rangle+\sigma
M_{\theta^{*}/\sigma}(R^{5,k}+B^{*}\bar{y})+\sigma
M_{\bar{\theta}^{*}/\sigma}(R^{6,k}+\bar{B}^{*}\bar{y})-\langle\bar{\delta}^{k},\,\bar{y}\rangle\end{array}\right\\}\\\\[5.0pt]
\bar{v}^{k+1}&=&-\mathop{\rm
Prox}_{\bar{\theta}^{*}/\sigma}(R^{6,k}+\bar{B}^{*}\bar{y}^{k+1})\\\\[5.0pt]
v^{k+1}&=&-\mathop{\rm
Prox}_{\theta^{*}/\sigma}(R^{5,k}+B^{*}\bar{y}^{k+1}).\end{array}\right.$
(166)

The detail of deriving (166) can be found in Appendix 2. The scale of the
subproblem solving for $\bar{y}^{k+1}$ can be large when $N$ increases. On the
other hand, the update scheme of $\bar{v}^{k+1}$ shows that this subproblem
can be solved in a parallel manner. We leave the efficient computation of
these two subproblems as a future work.

##  7 Numerical experiments

In this section, we conduct numerical experiments to test the efficiency of
our algorithm against state-of-the-art solvers. We use data sets found in the
literature (such as two-stage stochastic programming problems and
uncapacitated facility location problems) as well as synthetic data. All the
experiments are run in MATLAB R2017b on a computer with 1.4 GHz Quad-Core
Intel Core i5 Processor (this processor has 4 cores and 8 threads) and 8 GB
memory, equipped with macOS.

###  7.1 Stopping conditions

Based on the optimality conditions in (LABEL:eq-KKT), we measure the accuracy
of a computed solution by the following relative KKT residue:

$\eta=\max\\{\eta_{P},\eta_{D},0.2\eta_{{\cal
K}},\color[rgb]{0,0,0}0.2\eta_{\theta},\eta_{\bar{P}},\eta_{\bar{D}},0.2\eta_{{\bar{\cal
K}}},\color[rgb]{0,0,0}0.2\eta_{\bar{\theta}}\\},$

where

$\displaystyle\begin{array}[]{cccc}\eta_{P}=\frac{\|Ax-b\|}{1+\|b\|},&\eta_{D}=\frac{\|A^{*}y+\bar{B}^{*}\bar{y}+z+v-c\|}{1+\|c\|},&\eta_{{\cal
K}}=\frac{\|x-\Pi_{{\cal
K}}(x-z)\|}{1+\|x\|+\|z\|},&\color[rgb]{0,0,0}{\eta_{\theta}=\frac{\|x-\mathop{\rm{Prox}}_{\theta}(x-v)\|}{1+\|x\|+\|v\|}},\\\\[5.0pt]
\eta_{\bar{P}}=\frac{\|Bx+\bar{B}\bar{x}-\bar{b}\|}{1+\|\bar{b}\|},&\eta_{\bar{D}}=\frac{\|\bar{B}^{*}\bar{y}+\bar{z}+\bar{v}-\bar{c}\|}{1+\|\bar{c}\|},&\eta_{{\bar{\cal
K}}}=\frac{\|\bar{x}-\Pi_{{\bar{\cal
K}}}(\bar{x}-\bar{z})\|}{1+\|\bar{x}\|+\|\bar{z}\|},&\color[rgb]{0,0,0}\eta_{\bar{\theta}}=\frac{\|\bar{x}-\mathop{\rm{Prox}}_{\bar{\theta}}(\bar{x}-\bar{v})\|}{1+\|\bar{x}\|+\|\bar{v}\|}.\end{array}$

In addition, we also compute the duality gap by:

$\displaystyle\eta_{\rm gap}=\frac{{\rm obj}_{\rm P}-{\rm obj}_{\rm
D}}{1+\mid{\rm obj}_{\rm P}\mid+\mid{\rm obj}_{\rm D}\mid},$

where

$\displaystyle{\rm obj}_{\rm P}$ $\displaystyle:=$
$\displaystyle\theta(x)+\langle
c,\,x\rangle+\bar{\theta}(\bar{x})+\langle\bar{c},\,\bar{x}\rangle,$
$\displaystyle{\rm obj}_{\rm D}$ $\displaystyle:=$
$\displaystyle-\theta^{*}(-v)-\delta_{\cal K}^{*}(-z)+\langle
b,\,y\rangle-\bar{\theta}^{*}(-\bar{v})-\delta_{{\bar{\cal
K}}}^{*}(-\bar{z})+\langle\bar{b},\,\bar{y}\rangle,$

are the primal and dual objective values, respectively. We terminate our
algorithm when $\eta\leq 10^{-5}$ and $\eta_{\rm gap}\leq 10^{-4}$.

###  7.2 2-stage stochastic programming problems

In this section, we aim to solve two-stage stochastic programming problems on
some real data sets. We solve linear problems in the first subsection and
quadratic problems in the second subsection. Since there are not much
quadratic stochastic programming test data available, we will follow the idea
mentioned in LauWomersley  to extend stochastic linear programming problems to
its quadratic counterparts by adding a small quadratic term in the objective.
Specifically, we choose to set
$\theta(x)=\frac{0.1}{2}\|x\|^{2},\quad\bar{\theta}(\bar{x})=\frac{0.1}{2}\|\bar{x}\|^{2}.$
${\cal K}$ and $\bar{\cal K}$ in these instances are box constraints.

The detailed description of each data set used in this section can be found in
the Appendix. In table 1, we only list down the name and source of the data
sets.

Name | Source
---|---
| assets
---
env
phone
AIRL
4node
http://www4.uwsp.edu/math/afelt/slptestset/download.html
| pltexp
---
storm
http://users.iems.northwestern.edu/~jrbirge/html/dholmes/post.html
gbd | http://pages.cs.wisc.edu/~swright/stochastic/sampling/
Table 1: Source of two-stage stochastic programming data sets.

We compare our proposed DBA-sGS-ADMM with Gurobi under its default setting
except for modifying the accuracy tolerance parameters: FeasibilityTol,
OptimalityTol,
BarConvTol111https://www.gurobi.com/documentation/9.5/refman/parameters.html.,
which are absolute primal feasibility, absolute dual feasibility and relative
gap between primal and dual values respectively. More specifically, we set
FeasibilityTol=$\left(1+\max\\{\|b\|,\|\bar{b}\|\\}\right)\times 10^{-5}$,
OptimalityTol=$\left(1+\max\\{\|c\|,\|\bar{c}\|\\}\right)\times 10^{-5}$, and
BarConvTol=$10^{-4}$. In this setting, Gurobi adopts concurrent
optimizers222https://www.gurobi.com/documentation/9.5/refman/concurrent_optimizer.html
for linear programming, which runs multiple solvers (the primal/dual simplex
method and the barrier method) simultaneously, and chooses the one that
finishes first. Concurrent strategy requires multiple (at least 3) threads,
and the first three threads are assigned to primal simplex, dual simplex and
barrier respectively. The extra threads are then assigned to the barrier
method. More threads will not make the simplex method any faster. However, the
matrix factorization operation in the barrier method is faster with more
threads. We use 8 threads for Gurobi in order to adopt the concurrent strategy
for linear programming. Therefore, when solving linear programming, our method
is almost compared with the best of the simplex method and the barrier method,
meanwhile the barrier method is accelerated by multiple threads. Even though
this comparison is a bit unfair to our method, the results later show that the
performance of our method is comparable to Gurobi, and even better than Gurobi
for some large-scale instances.

For quadratic programming, Gurobi uses the barrier method by default. In this
case, it does improve the efficiency with more threads, but when the number of
threads increases, the speed is not necessarily faster because the cost of
data transmission will also increase. Meanwhile, the computational efficiency
of our DBA-sGS-ADMM/ALM can be further improved by incorporating parallel
computation. Hence in quadratic programming, we only use a single thread in
Gurobi for a fair comparison.

We incorporate a semismooth Newton method in DBA-sGS-ADMM for some examples
for which the size of the matrix $A$ is small while the number of scenarios is
relatively large. More specifically, we adopt a semismooth Newton method
described in Section 6 when $m_{0}\leq 10,n_{0}\leq 20$ and $N\geq 100$ in our
numerical experiments.

We adopt different strategies to solve the linear system (119) efficiently for
different cases. For some simple cases, such as when the component blocks of
$B$ or $\bar{B}$ are identical, or the data size is not very large (more
specifically, when $\sum_{i=1}^{N}m_{i}<5000$), we directly solve the linear
system (119) by utilizing Cholesky factorization. When the data size is
relatively large (specifically $\sum_{i=1}^{N}m_{i}\geq 5000$), we make use of
SMW formula described in Section 5.4 to solve the linear system. In our
numerical experiment, when $\sum_{i=1}^{N}m_{i}<50000$, we use the SMW formula
without the proximal term $\bar{J}$, i.e. $\bar{J}=0$. When the data scale is
very large, more specifically $\sum_{i=1}^{N}m_{i}\geq 50000$, we use the SMW
formula by setting
$\bar{J}_{i}=\lambda_{\max}(\bar{B}_{i}\bar{B}_{i}^{*})I_{m_{i}}-\bar{B}_{i}\bar{B}_{i}^{*}$,
$i=1,...,N$ to make the matrix $G$ in (130) more likely to be sparse.

For stochastic linear programming problems, we also compare our algorithm with
the well-known progressive hedging algorithm (PHA), which is applied to the
primal problem (P). We should emphasize that just like the PHA, the
subproblems associated with the second stage variables
$\bar{y},\bar{z},\bar{v}$ in our algorithm can also be solved in a parallel
manner, as mentioned in section 5 and section 6(b). However, in our case, the
cost of each parallel computation is already very cheap. If we implement the
distributed version using MATLAB’s parallel for-loop, the overhead cost will
become dominant and affect the performance adversely. Hence, we choose to
implement the algorithm serially in this paper. For a fair comparison, it is
reasonable to compare our non-parallel version with the non-parallel version
of PHA. Recent work in SXZ  has shown that one can choose the step length
$\tau$ in PHA to be in the range $(0,(1+\sqrt{5})/2)$. In our numerical
experiments, we set $\tau=1.618$ for both DBA-sGS-ADMM and PHA.

####  7.2.1 2-stage stochastic linear programming problems

Table 2:  Comparison of computational results between DBA-sGS-ADMM, Gurobi and PHA for two-stage stochastic linear programming problems. The results of Gurobi are obtained using 8 threads. The objective function value is shown to help us verifying the correctness easily. “-” means that the solver runs out of time.  | DBA-sGS-ADMM | Gurobi | PHA
---|---|---|---
Data | $m_{0}\mid m_{i}$ | $n_{0}\mid n_{i}$ | $N$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$
phone_1 |  1 $\mid$ 23 |  9 $\mid$ 93 | 1 | 2301 | 0.33 | 3.69e+01 | 25 | 0.03 | 3.69e+01 | 153 | 0.37 | 3.69e+01
AIRL_first |  2 $\mid$ 6 |  6 $\mid$ 12 | 25 | 12651 | 1.73 | 2.49e+05 | 146 | 0.01 | 2.49e+05 | 312 | 11.99 | 2.49e+05
AIRL_second |  2 $\mid$ 6 |  6 $\mid$ 12 | 25 | 14751 | 2.02 | 2.70e+05 | 157 | 0.01 | 2.70e+05 | 242 | 9.28 | 2.70e+05
4node_256 |  14 $\mid$ 74 |  60 $\mid$ 198 | 256 | 637 | 2.4 | 4.25e+02 | 25078 | 2.72 | 4.25e+02 | 504 | 492.39 | 4.25e+02
4node_512 |  14 $\mid$ 74 |  60 $\mid$ 198 | 512 | 2681 | 7.54 | 4.30e+02 | 49908 | 8.37 | 4.30e+02 | 813 | 1617.74 | 4.30e+02
4node_1024 |  14 $\mid$ 74 |  60 $\mid$ 198 | 1024 | 1450 | 7.89 | 4.34e+02 | 119424 | 11.10 | 4.34e+02 | 1022 | 3945.46 | 4.34e+02
4node_2048 |  14 $\mid$ 74 |  60 $\mid$ 198 | 2048 | 1287 | 13.88 | 4.42e+02 | 259486 | 27.42 | 4.42e+02 | - | - | -
4node_4096 |  14 $\mid$ 74 |  60 $\mid$ 198 | 4096 | 1628 | 36.64 | 4.47e+02 | 464119 | 66.96 | 4.47e+02 | - | - | -
4node_8192 |  14 $\mid$ 74 |  60 $\mid$ 198 | 8192 | 767 | 41.14 | 4.47e+02 | 141 | 239.78 | 4.47e+02 | - | - | -
4node_16384 |  14 $\mid$ 74 |  60 $\mid$ 198 | 16384 | 782 | 89.88 | 4.47e+02 | 155 | 474.69 | 4.47e+02 | - | - | -
4node_32768 |  14 $\mid$ 74 |  60 $\mid$ 198 | 32768 | 748 | 186.84 | 4.47e+02 | 164 | 869.61 | 4.47e+02 | - | - | -
4node_base_1024 |  16 $\mid$ 74 |  60 $\mid$ 198 | 1024 | 1881 | 10.39 | 4.23e+02 | 121 | 20.36 | 4.23e+02 | - | - | -
4node_base_2048 |  16 $\mid$ 74 |  60 $\mid$ 198 | 2048 | 1511 | 19.01 | 4.24e+02 | 125 | 53.30 | 4.25e+02 | - | - | -
4node_base_4096 |  16 $\mid$ 74 |  60 $\mid$ 198 | 4096 | 1290 | 33.63 | 4.25e+02 | 127 | 112.69 | 4.25e+02 | - | - | -
4node_base_8192 |  16 $\mid$ 74 |  60 $\mid$ 198 | 8192 | 898 | 54.21 | 4.25e+02 | 158 | 326.67 | 4.25e+02 | - | - | -
4node_base_16384 |  16 $\mid$ 74 |  60 $\mid$ 198 | 16384 | 968 | 122.31 | 4.25e+02 | 117 | 522.26 | 4.25e+02 | - | - | -
4node_base_32768 |  16 $\mid$ 74 |  60 $\mid$ 198 | 32768 | 840 | 234.75 | 4.25e+02 | 152 | 1510.31 | 4.25e+02 | - | - | -

From Table LABEL:table:SLP, we can observe that DBA-sGS-ADMM outperforms
Gurobi especially when the number of scenarios $N$ is large. For example, for
the instance 4node_32768 with about $3.27\times 10^{4}$ scenarios, we are
faster by at least 4 times. We also see that PHA is the least efficient, and
it can only solve moderate scale problems, while DBA-sGS-ADMM can solve
problem with up to $3.27\times 10^{4}$ scenarios and medium-size dimensional
decision variables. The inefficiency of the PHA could be due to the fact that
at each iteration, it needs to solve $N$ smaller LP subproblems whose total
computational cost are rather high. We have to mention that under the
framework of PHA, the cost of computing the projections onto the feasible sets
of the scenarios subproblems to check the primal feasibility condition can be
quite expensive. Therefore, in practice, one only checks the feasibility of
the nonanticipativity constraint and the relative error of iteration points or
function values. Therefore, another advantage of our algorithm compared to PHA
is that the stopping criterion (KKT residue) of DBA-sGS-ADMM is theoretically
more rigorous than that of the PHA.

####  7.2.2 2-stage stochastic quadratic programming problems

Table 3:  Comparison of computational results between DBA-sGS-ADMM and Gurobi for two-stage quadratic stochastic programming problem (with added 0.1$I$ in the quadratic objective). All the results are obtained using a 8 thread. The objective function value is shown to help us verifying the correctness easily. | DBA-sGS-ADMM | Gurobi
---|---|---
Data | $m_{0}\mid m_{i}$ | $n_{0}\mid n_{i}$ | $N$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$
assets_small |  5 $\mid$ 5 |  13 $\mid$ 13 | 100 | 577 | 0.47 | 1.86e+05 | 21 | 0.43 | 1.86e+05
assets_large |  5 $\mid$ 5 |  13 $\mid$ 13 | 37500 | 6801 | 395.38 | 4.92e+06 | 27 | 14.12 | 4.92e+06
env_aggr |  32 $\mid$ 48 |  69 $\mid$ 85 | 5 | 501 | 0.49 | 6.38e+04 | 23 | 0.20 | 6.38e+04
env_loose |  32 $\mid$ 48 |  69 $\mid$ 85 | 5 | 255 | 0.16 | 8.07e+04 | 21 | 0.49 | 8.07e+04
env_first |  32 $\mid$ 48 |  69 $\mid$ 85 | 5 | 423 | 0.17 | 6.93e+04 | 21 | 0.41 | 6.93e+04
env_imp |  32 $\mid$ 48 |  69 $\mid$ 85 | 15 | 1022 | 0.50 | 1.22e+05 | 25 | 0.20 | 1.22e+05
env_1200 |  32 $\mid$ 48 |  69 $\mid$ 85 | 1200 | 2872 | 14.38 | 6.68e+06 | 36 | 3.89 | 6.68e+06
env_1875 |  32 $\mid$ 48 |  69 $\mid$ 85 | 1875 | 4800 | 21.62 | 1.05e+07 | 37 | 6.02 | 1.05e+07
env_3780 |  32 $\mid$ 48 |  69 $\mid$ 85 | 3780 | 5099 | 44.45 | 2.11e+07 | 43 | 8.98 | 2.11e+07
env_5292 |  32 $\mid$ 48 |  69 $\mid$ 85 | 5292 | 5361 | 76.59 | 2.95e+07 | 42 | 11.55 | 2.95e+07
env_lrge |  32 $\mid$ 48 |  69 $\mid$ 85 | 8232 | 4611 | 125.11 | 4.58e+07 | 39 | 20.17 | 4.58e+07
env_xlrge |  32 $\mid$ 48 |  69 $\mid$ 85 | 32928 | 4553 | 1470.55 | 1.85e+08 | 38 | 121.71 | 1.85e+08
pltexpA2_6 |  62 $\mid$ 104 |  188 $\mid$ 272 | 6 | 20000 | 5.20 | 6.50e+04 | 25 | 0.48 | 6.50e+04
pltexpA2_16 |  62 $\mid$ 104 |  188 $\mid$ 272 | 16 | 20000 | 7.91 | 1.54e+05 | 24 | 0.58 | 1.54e+05
stormG2_8 |  58 $\mid$ 528 |  179 $\mid$ 1377 | 8 | 14401 | 28.59 | 1.83e+07 | 63 | 0.96 | 1.83e+07
stormG2_27 |  58 $\mid$ 528 |  179 $\mid$ 1377 | 27 | 8101 | 23.88 | 2.11e+07 | 66 | 5.12 | 2.11e+07
stormG2_125 |  58 $\mid$ 528 |  179 $\mid$ 1377 | 125 | 15501 | 127.37 | 3.04e+07 | 62 | 49.37 | 3.04e+07
gbd |  4 $\mid$ 5 |  21 $\mid$ 10 | 646425 | 1851 | 462.80 | 7.94e+08 | 21 | 118.87 | 8.21e+08
phone |  1 $\mid$ 23 |  9 $\mid$ 93 | 32768 | 501 | 48.89 | 2.74e+05 | 34 | 56.63 | 2.74e+05
phone_1 |  1 $\mid$ 23 |  9 $\mid$ 93 | 1 | 157 | 0.08 | 4.77e+01 | 12 | 0.37 | 4.77e+01
AIRL_first |  2 $\mid$ 6 |  6 $\mid$ 12 | 25 | 1551 | 0.26 | 3.61e+06 | 21 | 0.33 | 3.61e+06
AIRL_second |  2 $\mid$ 6 |  6 $\mid$ 12 | 25 | 1651 | 0.25 | 3.55e+06 | 23 | 0.35 | 3.55e+06
4node_1 |  14 $\mid$ 74 |  60 $\mid$ 198 | 1 | 1371 | 0.25 | 8.95e+03 | 24 | 0.31 | 8.95e+03
4node_2 |  14 $\mid$ 74 |  60 $\mid$ 198 | 2 | 1311 | 0.25 | 9.02e+03 | 24 | 0.40 | 9.02e+03
4node_256 |  14 $\mid$ 74 |  60 $\mid$ 198 | 256 | 881 | 1.10 | 1.25e+04 | 24 | 1.15 | 1.25e+04
4node_512 |  14 $\mid$ 74 |  60 $\mid$ 198 | 512 | 764 | 1.61 | 1.39e+04 | 24 | 2.05 | 1.39e+04
4node_1024 |  14 $\mid$ 74 |  60 $\mid$ 198 | 1024 | 691 | 2.97 | 1.64e+04 | 31 | 4.41 | 1.64e+04
4node_2048 |  14 $\mid$ 74 |  60 $\mid$ 198 | 2048 | 763 | 6.26 | 2.15e+04 | 50 | 17.31 | 2.15e+04
4node_4096 |  14 $\mid$ 74 |  60 $\mid$ 198 | 4096 | 873 | 17.12 | 3.14e+04 | 36 | 25.57 | 3.14e+04
4node_8192 |  14 $\mid$ 74 |  60 $\mid$ 198 | 8192 | 912 | 42.18 | 5.05e+04 | 53 | 107.71 | 5.05e+04
4node_16384 |  14 $\mid$ 74 |  60 $\mid$ 198 | 16384 | 839 | 84.11 | 8.64e+04 | 33 | 100.94 | 8.65e+04
4node_32768 |  14 $\mid$ 74 |  60 $\mid$ 198 | 32768 | 646 | 142.01 | 1.55e+05 | 29 | 164.33 | 1.55e+05
4node_base_1024 |  16 $\mid$ 74 |  60 $\mid$ 198 | 1024 | 501 | 2.19 | 1.60e+04 | 65 | 6.69 | 1.60e+04
4node_base_2048 |  16 $\mid$ 74 |  60 $\mid$ 198 | 2048 | 501 | 4.65 | 2.02e+04 | 48 | 61.26 | 2.02e+04
4node_base_4096 |  16 $\mid$ 74 |  60 $\mid$ 198 | 4096 | 501 | 10.86 | 2.85e+04 | 30 | 21.29 | 2.85e+04
4node_base_8192 |  16 $\mid$ 74 |  60 $\mid$ 198 | 8192 | 811 | 46.01 | 4.51e+04 | 48 | 45.12 | 4.51e+04
4node_base_16384 |  16 $\mid$ 74 |  60 $\mid$ 198 | 16384 | 931 | 111.32 | 7.82e+04 | 188 | 701.29 | 7.82e+04
4node_base_32768 |  16 $\mid$ 74 |  60 $\mid$ 198 | 32768 | 1271 | 313.84 | 1.44e+05 | 63 | 3886.62 | 1.44e+05

In table LABEL:table:SP, we observe that the computational performance of DBA-
sGS-ADMM and Gurobi are comparable for many instances, in which the runtime of
DBA-sGS-ADMM is shorter among $21/38\approx 55.3\%$ of the tested instances.
In particular, for a few of the 4node_base_xxxx instances, DBA-sGS-ADMM can
outperform Gurobi by a large margin. On the other hand, Gurobi can solve some
of the instances such as env_xxxx very efficiently. This is reasonable, since
we only add a very sparse term of $0.1I$ in the quadratic objective, and
Gurobi can take full advantage of the sparsity in the data. Later in section
7.4.1, we will see that when the matrix in the quadratic objective is dense,
we can expect extreme inefficiency of the barrier method used by Gurobi, and
DBA-sGS-ADMM has substantial computational advantage.

###  7.3 DNN relaxations of uncapacitated facility location (UFL) problems

####  7.3.1 DNN relaxations of linear UFL problems

In this subsection, we aim to solve some DNN relaxations of linear UFL
problems. We have two benchmark datasets of different scales, namely the small
scale Bilde-Krarup dataset and the large scale Koerkel-Ghosh dataset. Both of
these data instances could be downloaded from the UflLib at
http://resources.mpi-
inf.mpg.de/departments/d1/projects/benchmarks/UflLib/index.html. These are all
linear UFL problem, i.e. $q_{ij}=0\;\forall i=1,\ldots,p,j=1,\ldots,q$ in
(29).

We should emphasize that our intention here is not to investigate the
tightness of the DNN relaxation (45) of the UFL problem (29) but to treat the
DNN relaxation as an example of the DBA problem (17) with constraint set
${\cal K}$ being a positive semdefinite matrix cone $\mathbb{S}^{n}_{+}$.

We first solve the DNN relaxations of linear UFL problems described in (45).
The numerical results obtained by DBA-sGS-ADMM will be compared against MOSEK
for the small scale dataset Bilde-Krarup and the results are reported in Table
LABEL:table:SUFL-SDP. Note that MOSEK is the current best solver for solving
general SDP problems. When it is employed to solve the dual block angular
problem (45), MOSEK treats it as a general SDP problem since it does not
exploit the block-angular structure in the constraint matrix except for
sparsity in the data matrices.

Note that for UFL problems, we can make use of the special structure of $B$
and $\bar{B}$ in the constraint to further cut down the cost of solving the
linear system (119). Specifically, the component blocks of $B$ and $\bar{B}$
are repeated, i.e. $B_{1}=B_{2}=...=B_{N}$,
$\bar{B}_{1}=\bar{B}_{2}=...=\bar{B}_{N}$. Furthermore,
$\bar{B}_{j}\bar{B}_{j}^{*}$ has a very simple inverse given by (54). Thus we
can solve the linear system (119) efficiently by utilizing the SMW formula
described in Section 5.4(a).

Table 4:  Comparison of computational results for the DNN relaxation of UFL problem in the small scale Bilde-Krarup dataset using DBA-sGS-ADMM and MOSEK. The results of MOSEK are obtained using 8 threads. Under the column “Gap”, the relaxation gap is also reported. | Original IP | SDP via DBA-sGS-ADMM | SDP via MOSEK
---|---|---|---
Data | $m_{0}\mid m_{i}$ | $n_{0}\mid n_{i}$ | $N$ | Obj | Obj | Gap | Iter | Time(s) | Obj | Gap | Iter | Time(s)
B1.1 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 23468 | 23368 | 0.43 | 1883 | 4.08 | 23363 | 0.45 | 24 | 36.33
B1.10 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 21864 | 21193 | 3.07 | 1911 | 4.28 | 21187 | 3.10 | 24 | 36.27
C1.1 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 16781 | 16034 | 4.45 | 2151 | 5.29 | 16033 | 4.46 | 33 | 63.74
C1.10 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 17994 | 17151 | 4.68 | 1973 | 4.40 | 17148 | 4.70 | 31 | 46.04
E1.1 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 15042 | 13876 | 7.75 | 1969 | 3.32 | 13871 | 7.78 | 28 | 44.90
E1.10 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 14630 | 13323 | 8.93 | 1872 | 3.15 | 13319 | 8.96 | 30 | 44.99
E5.1 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 32377 | 29290 | 9.53 | 1603 | 2.52 | 29279 | 9.57 | 25 | 37.82
E5.10 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 34086 | 30122 | 11.63 | 1854 | 2.84 | 30114 | 11.65 | 26 | 39.18
E10.1 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 46832 | 42496 | 9.26 | 1321 | 2.03 | 42492 | 9.27 | 26 | 39.62
E10.10 |  51 $\mid$ 51 |  51 $\mid$ 100 | 100 | 47449 | 43179 | 9.00 | 1869 | 2.96 | 43177 | 9.00 | 24 | 37.82
Table 5:  Computational results for the DNN relaxation of UFL problem in the large scale Koerkel-Ghosh dataset using DBA-sGS-ADMM. Under the column“Gap”, the relaxation gap is also reported. | Original IP | SDP via DBA-sGS-ADMM
---|---|---
Data | $m_{0}\mid m_{i}$ | $n_{0}\mid n_{i}$ | $N$ | Obj | Obj | Gap | Iter | Time(s)
ga250a-1 |  251 $\mid$ 251 |  251 $\mid$ 500 | 250 | 2.5796e+05 | 2.5764e+05 | 0.12 | 1751 | 16.58
ga250b-2 |  251 $\mid$ 251 |  251 $\mid$ 500 | 250 | 2.7514e+05 | 2.7272e+05 | 0.88 | 1646 | 15.61
ga250c-3 |  251 $\mid$ 251 |  251 $\mid$ 500 | 250 | 3.3366e+05 | 3.2288e+05 | 3.23 | 1504 | 13.55
ga500a-4 |  501 $\mid$ 501 |  501 $\mid$ 1000 | 500 | 5.1105e+05 | 5.1043e+05 | 0.12 | 2601 | 127.43
ga500b-5 |  501 $\mid$ 501 |  501 $\mid$ 1000 | 500 | 5.3748e+05 | 5.3300e+05 | 0.83 | 2901 | 137.53
ga500c-1 |  501 $\mid$ 501 |  501 $\mid$ 1000 | 500 | 6.2136e+05 | 6.0286e+05 | 2.98 | 1796 | 83.73
ga750a-2 |  751 $\mid$ 751 |  751 $\mid$ 1500 | 750 | 7.6367e+05 | 7.6254e+05 | 0.15 | 3001 | 368.26
ga750b-3 |  751 $\mid$ 751 |  751 $\mid$ 1500 | 750 | 7.9613e+05 | 7.8965e+05 | 0.81 | 5601 | 795.05
ga750c-4 |  751 $\mid$ 751 |  751 $\mid$ 1500 | 750 | 9.0004e+05 | 8.7556e+05 | 2.72 | 2197 | 246.48
gs250a-1 |  251 $\mid$ 251 |  251 $\mid$ 500 | 250 | 2.5796e+05 | 2.5769e+05 | 0.10 | 2091 | 23.28
gs250b-2 |  251 $\mid$ 251 |  251 $\mid$ 500 | 250 | 2.7568e+05 | 2.7305e+05 | 0.95 | 1578 | 16.92
gs250c-3 |  251 $\mid$ 251 |  251 $\mid$ 500 | 250 | 3.3300e+05 | 3.2199e+05 | 3.31 | 1707 | 17.28
gs500a-4 |  501 $\mid$ 501 |  501 $\mid$ 1000 | 500 | 5.1114e+05 | 5.1042e+05 | 0.14 | 3101 | 143.43
gs500b-5 |  501 $\mid$ 501 |  501 $\mid$ 1000 | 500 | 5.3827e+05 | 5.3312e+05 | 0.96 | 3401 | 149.40
gs500c-1 |  501 $\mid$ 501 |  501 $\mid$ 1000 | 500 | 6.2004e+05 | 6.0198e+05 | 2.91 | 1871 | 81.54
gs750a-2 |  751 $\mid$ 751 |  751 $\mid$ 1500 | 750 | 7.6355e+05 | 7.6257e+05 | 0.13 | 3251 | 382.80
gs750b-3 |  751 $\mid$ 751 |  751 $\mid$ 1500 | 750 | 7.9659e+05 | 7.8995e+05 | 0.83 | 4601 | 510.33
gs750c-4 |  751 $\mid$ 751 |  751 $\mid$ 1500 | 750 | 9.0134e+05 | 8.7541e+05 | 2.88 | 2608 | 292.61

From Table LABEL:table:SUFL-SDP, we could observe that DBA-sGS-ADMM always
outperform MOSEK in solving the DNN problems (45) with dual block angular
structure. The runtime is about 10-20 times faster than that required by
MOSEK.

In Table LABEL:table:LUFL-SDP, we report the performance of DBA-sGS-ADMM in
solving (45) for the large scale dataset Koerkel-Ghosh. We can see that our
DBA-sGS-ADMM can solve these large scale dual block angular DNN problems very
efficiently. In particular, for the last instance gs750c-4 which has $564001$
linear constraints, a first stage $751$ by $751$ DNN matrix variables, and
$750$ second stage $1500$-dimensional nonnegative variables, we can solve it
in about 5 minutes.

We also tried to use MOSEK to solve the instances in Table LABEL:table:LUFL-
SDP, but they often run out of memory on our machine. In addition, due to the
inefficiency of MOSEK as observed in Table LABEL:table:SUFL-SDP in solving
this type of problems, we do not report the numerical results of MOSEK for
these large scale instances.

####  7.3.2 DNN relaxations of quadratic UFL problems

In this subsection, we generate some random quadratic UFL problems following a
procedure mentioned in Gunluk . We compare the performance of our DBA-sGS-ADMM
and MOSEK for solving the DNN relaxation of the quadratic UFL in (45) in Table
LABEL:table:QUFL-SDP. The strategy we adopt to solve the linear system (119)
is the same as that in Section 7.3.1.

We should take note that MOSEK could not accept the model with both quadratic
objective function and semidefinite constraint. Thus we reformulate (45) as
follows. First, we have the Cholesky factorization for each ${\rm
diag}(Q_{j}):=H_{j}^{T}H_{j}$, where $H_{j}={\rm
diag}\big{(}\left(\sqrt{Q_{j1}},...,\sqrt{Q_{jp}}\right)^{T}\big{)}$, then we
can reformulate (45) to be a convex conic programming model with rotated
quadratic cone constraints:

$\displaystyle\begin{array}[]{rl}\min&\langle[0\mathchar
24635\relax\;c],\,[\alpha\mathchar
24635\relax\;u]\rangle+\sum_{j=1}^{q}\langle\color[rgb]{0,0,0}P_{j},\,S_{j}\rangle+\sum_{j=1}^{q}t_{j}\\\\[8.0pt]
\mbox{s.t}&\left[\begin{array}[]{c}0\\\
u\end{array}\right]+\left[\begin{array}[]{cc}e^{T}&0^{T}\\\
-I_{p}&-I_{p}\end{array}\right]\left[\begin{array}[]{c}S_{j}\\\
Z_{j}\end{array}\right]=\left[\begin{array}[]{c}1\\\ 0\end{array}\right],\quad
j=1,\ldots q\\\\[15.0pt] &u-\mbox{diag}(U)=0,\;\;\alpha=1,\\\\[5.0pt]
&\mbox{\boldmath$U$}:=\left[\begin{array}[]{cc}\alpha&u^{T}\\\\[3.0pt]
u&U\end{array}\right]\in\mathbb{S}^{1+p}_{+},\;\;\mbox{\boldmath$U$}\geq
0,\;\;S_{j},Z_{j}\geq 0,\;j=1,\ldots,q,\\\\[15.0pt] &y_{j}=H_{j}S_{j},\quad
j=1,\ldots q,\\\\[5.0pt]
&\beta_{j}=1,\quad\left[\begin{array}[]{c}\beta_{j}\\\ t_{j}\\\
y_{j}\end{array}\right]\in{\cal
Q}_{j}:=\\{x\in\mathbb{R}^{p+2}:2x_{1}x_{2}\geq\sum_{i=3}^{p+2}x_{i}^{2},\;x_{1}\geq
0,\;x_{2}\geq 0\\}.\end{array}$

We solve the problem (45) by MOSEK via the above reformulation.

Table 6:  Comparison of computational results for the DNN relaxation of random large scale QUFL problem using DBA-sGS-ADMM and MOSEK. The results of MOSEK are obtained using 8 threads. “-” means that the solver runs out of time. | SDP via DBA-sGS-ADMM | SDP via MOSEK
---|---|---
Data | $m_{0}\mid m_{i}$ | $n_{0}\mid n_{i}$ | $N$ | Obj | Iter | Time(s) | Obj | Iter | Time(s)
randQUFL-m50-n10 |  11 $\mid$ 11 |  11 $\mid$ 20 | 50 | 83.73 | 101 | 0.16 | 83.73 | 7 | 0.46
randQUFL-m50-n20 |  21 $\mid$ 21 |  21 $\mid$ 40 | 50 | 68.11 | 107 | 0.10 | 68.11 | 10 | 0.77
randQUFL-m50-n30 |  31 $\mid$ 31 |  31 $\mid$ 60 | 50 | 45.78 | 173 | 0.15 | 45.78 | 13 | 1.83
randQUFL-m200-n50 |  51 $\mid$ 51 |  51 $\mid$ 100 | 200 | 73.52 | 402 | 0.54 | 73.52 | 9 | 67.24
randQUFL-m200-n100 |  101 $\mid$ 101 |  101 $\mid$ 200 | 200 | 46.97 | 385 | 1.01 | 46.97 | 16 | 1325.55
randQUFL-m200-n150 |  151 $\mid$ 151 |  151 $\mid$ 300 | 200 | 42.00 | 469 | 1.66 | 42.00 | 16 | 41258.92
randQUFL-m200-n200 |  201 $\mid$ 201 |  201 $\mid$ 400 | 200 | 38.35 | 501 | 2.52 | - | - | -
randQUFL-m1000-n200 |  201 $\mid$ 201 |  201 $\mid$ 400 | 1000 | 79.85 | 611 | 7.17 | - | - | -
randQUFL-m1000-n500 |  501 $\mid$ 501 |  501 $\mid$ 1000 | 1000 | 47.94 | 1021 | 55.13 | - | - | -
randQUFL-m1000-n1000 |  1001 $\mid$ 1001 |  1001 $\mid$ 2000 | 1000 | 35.29 | 791 | 159.92 | - | - | -
randQUFL-m1000-n2000 |  2001 $\mid$ 2001 |  2001 $\mid$ 4000 | 1000 | 24.95 | 651 | 701.38 | - | - | -

The comparison of DBA-sGS-ADMM and MOSEK in Table LABEL:table:QUFL-SDP are
similar to that reported in section 7.3.1. One can also observe that DBA-sGS-
ADMM is much faster than MOSEK when solving the DNN problems. MOSEK can only
solve small or moderate scale problems, and become extremely inefficient when
the scale of problems is large. In particular, for the last instance solved by
MOSEK, our method is over 24,000 times faster than MOSEK.

###  7.4 Randomly generated problems

To demonstrate the generality and superiority of our algorithm, in this
section, we conduct numerical experiments on some randomly generated dual
block angular problems. We have two classes of random problems, namely, the
DBA quadratic programming (QP) and the DBA semidefinite programming (SDP).

####  7.4.1 QP

Here we generate some random quadratic programming problems with dual block
angular structure. In particular, we have the following:

$\displaystyle\theta(x):=\frac{1}{2}\langle x,\,Qx\rangle\;\text{ and
}\bar{\theta}_{i}(\bar{x}_{i}):=\frac{1}{2}\langle\bar{x}_{i},\,{\bar{Q}}_{i}\bar{x}_{i}\rangle\;\forall
i=1,...,N,$ $\displaystyle\text{with the bounds }{\cal
K}:=\mathbb{R}^{n}_{+},\;{\bar{\cal K}}_{i}:=\mathbb{R}^{n_{i}}_{+}\;\forall
i=1,...,N.$

We randomly generate sparse matrices $A,B$ (by MATLAB command sprand) with
density $10/{n_{0}}$. Similarly, $\bar{B}$ is generated with the same command
with density $10/{n_{i}}$. Finally, we randomly generate sparse symmetric
positive definite matrices $Q$ and ${\bar{Q}}_{i}$ with density $2/{n_{0}}$
and $2/{n_{i}}$ respectively (using MATLAB command sprandsym). The strategy we
adopt to solve the linear system (119) is the same as that in Section 7.2.

Table 7:  Comparison of computational results between DBA-sGS-ADMM and Gurobi for randomly generated QP problem. All the results are obtained using a 8 thread. The objective function value is shown to help us verifying the correctness easily. “-” means that the solver runs out of time. | DBA-sGS-ADMM | Gurobi
---|---|---
Data | $m_{0}\mid m_{i}$ | $n_{0}\mid n_{i}$ | $N$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$
randQP-m10-n20-N10 |  10 $\mid$ 10 |  20 $\mid$ 20 | 10 | 104 | 0.11 | 3.523e+02 | 11 | 0.28 | 3.523e+02
randQP-m50-n80-N10 |  50 $\mid$ 50 |  80 $\mid$ 80 | 10 | 151 | 0.15 | 1.806e+03 | 10 | 0.30 | 1.806e+03
randQP-m100-n200-N10 |  100 $\mid$ 100 |  200 $\mid$ 200 | 10 | 283 | 0.50 | 3.707e+03 | 11 | 0.51 | 3.707e+03
randQP-m200-n300-N10 |  200 $\mid$ 200 |  300 $\mid$ 300 | 10 | 374 | 2.77 | 7.048e+03 | 11 | 1.66 | 7.048e+03
randQP-m500-n800-N10 |  500 $\mid$ 500 |  800 $\mid$ 800 | 10 | 335 | 3.77 | 1.770e+04 | 10 | 17.34 | 1.770e+04
randQP-m100-n200-N50 |  100 $\mid$ 100 |  200 $\mid$ 200 | 50 | 810 | 2.50 | 1.713e+04 | 13 | 0.99 | 1.713e+04
randQP-m200-n300-N50 |  200 $\mid$ 200 |  300 $\mid$ 300 | 50 | 957 | 7.02 | 3.362e+04 | 10 | 3.66 | 3.362e+04
randQP-m500-n800-N50 |  500 $\mid$ 500 |  800 $\mid$ 800 | 50 | 791 | 32.86 | 8.494e+04 | 11 | 20.23 | 8.494e+04
randQP-m100-n200-N100 |  100 $\mid$ 100 |  200 $\mid$ 200 | 100 | 973 | 5.65 | 3.460e+04 | 13 | 7.45 | 3.460e+04
randQP-m200-n300-N100 |  200 $\mid$ 200 |  300 $\mid$ 300 | 100 | 1023 | 14.57 | 6.687e+04 | 12 | 11.53 | 6.687e+04
randQP-m500-n800-N100 |  500 $\mid$ 500 |  800 $\mid$ 800 | 100 | 601 | 7.30 | 1.688e+05 | 13 | 70.68 | 1.688e+05
randQP-m100-n200-N200 |  100 $\mid$ 100 |  200 $\mid$ 200 | 200 | 1169 | 13.35 | 6.897e+04 | 15 | 3.41 | 6.897e+04
randQP-m200-n300-N200 |  200 $\mid$ 200 |  300 $\mid$ 300 | 200 | 1187 | 33.57 | 1.324e+05 | 13 | 19.67 | 1.324e+05
randQP-m500-n800-N200 |  500 $\mid$ 500 |  800 $\mid$ 800 | 200 | 859 | 18.58 | 3.369e+05 | 13 | 118.72 | 3.369e+05
randQP-m1000-n2000-N200 |  1000 $\mid$ 1000 |  2000 $\mid$ 2000 | 200 | 1326 | 68.79 | 7.003e+05 | 22 | 869.92 | 7.003e+05
randQP-m2000-n3000-N200 |  2000 $\mid$ 2000 |  3000 $\mid$ 3000 | 200 | 916 | 114.02 | 1.351e+06 | 16 | 4237.40 | 1.351e+06
randQP-m5000-n8000-N200 |  5000 $\mid$ 5000 |  8000 $\mid$ 8000 | 200 | 901 | 466.62 | 3.403e+06 | - | - | -

From Table LABEL:table:randQP, we can observe that DBA-sGS-ADMM is much faster
than Gurobi when $m_{0},m_{i},n_{0},n_{i}$ are large. We note that
$n_{0},n_{i}$ represent the dimensions of first stage and second stage
decision variables respectively. Meanwhile, we generate denser $Q$ in the
quadratic objective than that in section 7.2.2 ($Q=0.1I$), making the
computational cost of Gurobi per iteration become more expensive. Therefore,
we could expect that our DBA-sGS-ADMM to be more efficient in solving two
stage stochastic quadratic programming problems when the matrix in quadratic
objective is dense and the dimensions of the decision variables are large.

####  7.4.2 SDP

We also generate a random data set of semidefinite programming problems with
dual block angular structure. In particular, we have

$\displaystyle\theta(x):=0,\;\bar{\theta}_{i}(\bar{x}_{i}):=0\;\forall
i=1,...,N,$ $\displaystyle\text{with the bounds }{\cal
K}:=S^{n}_{+},\;{\bar{\cal K}}_{i}:=S^{n_{i}}_{+}\;\forall i=1,...,N,$
$\displaystyle\text{and the mappings }A(x):=\begin{bmatrix}\langle
A_{1},\,x\rangle\\\ \langle A_{2},\,x\rangle\\\ \vdots\\\ \langle
A_{m},\,x\rangle\end{bmatrix},\;B_{i}(x):=\begin{bmatrix}\langle
B_{i,1},\,x\rangle\\\ \langle B_{i,2},\,x\rangle\\\ \vdots\\\ \langle
B_{i,m_{i}},\,x\rangle\end{bmatrix},$
$\displaystyle\bar{B}_{i}(\bar{x}_{i}):=\begin{bmatrix}\langle\bar{B}_{i,1},\,\bar{x}_{i}\rangle\\\
\langle\bar{B}_{i,2},\,\bar{x}_{i}\rangle\\\ \vdots\\\
\langle\bar{B}_{i,m_{i}},\,\bar{x}_{i}\rangle\end{bmatrix}.$

We generate the matrix representation of the linear mapping $A_{i}$ for
$i=1,\ldots,m$ (using MATLAB routine Ai = sprand(n0,n0,0.2); Ai = Ai*Ai’;).
Similarly, the matrix representation of $B_{i,k}$ and $\bar{B}_{i,k}$ for
$k=1,\ldots,{m}_{i}$, $i=1,\ldots,N$, is generated using the same routine
except with density $5/{n_{i}}$. The strategy we adopt to solve the linear
system (119) is the same as that in Section 7.2.

Table 8:  Comparison of computational results between DBA-sGS-ADMM and MOSEK for randomly generated SDP problem. The results of MOSEK are obtained using 8 threads.  | DBA-sGS-ADMM | MOSEK
---|---|---
Data | $m_{0}\mid m_{i}$ | $n_{0}\mid n_{i}$ | $N$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$ | Iter | Time(s) | $\mbox{obj}_{\rm P}$
randSDP-m50-n50-N5 |  50 $\mid$ 50 |  100 $\mid$ 100 | 5 | 544 | 8.79 | 2.529e+05 | 10 | 2.72 | 2.529e+05
randSDP-m100-n100-N5 |  100 $\mid$ 100 |  100 $\mid$ 100 | 5 | 202 | 5.56 | 4.682e+05 | 9 | 6.89 | 4.682e+05
randSDP-m100-n200-N5 |  100 $\mid$ 100 |  200 $\mid$ 200 | 5 | 501 | 36.22 | 2.769e+06 | 13 | 32.77 | 4.778e+06
randSDP-m200-n200-N5 |  200 $\mid$ 200 |  200 $\mid$ 200 | 5 | 201 | 30.35 | 4.778e+06 | 10 | 57.26 | 4.778e+06
randSDP-m200-n300-N5 |  200 $\mid$ 200 |  300 $\mid$ 300 | 5 | 256 | 70.95 | 1.521e+07 | 11 | 159.41 | 1.521e+07
randSDP-m200-n400-N5 |  200 $\mid$ 200 |  400 $\mid$ 400 | 5 | 277 | 125.52 | 3.710e+07 | 13 | 389.44 | 3.710e+07
randSDP-m400-n400-N5 |  400 $\mid$ 400 |  400 $\mid$ 400 | 5 | 265 | 248.69 | 6.399e+07 | 10 | 838.80 | 6.398e+07
randSDP-m10-n20-N10 |  10 $\mid$ 10 |  20 $\mid$ 20 | 10 | 2375 | 3.48 | 5.222e+03 | 9 | 0.49 | 5.222e+03
randSDP-m50-n80-N10 |  50 $\mid$ 50 |  80 $\mid$ 80 | 10 | 236 | 4.88 | 2.327e+05 | 10 | 3.37 | 2.327e+05
randSDP-m50-n100-N10 |  50 $\mid$ 50 |  100 $\mid$ 100 | 10 | 162 | 4.80 | 3.631e+05 | 11 | 5.58 | 3.631e+05
randSDP-m80-n100-N10 |  80 $\mid$ 80 |  100 $\mid$ 100 | 10 | 240 | 9.67 | 5.406e+05 | 10 | 10.69 | 5.406e+05
randSDP-m100-n200-N10 |  100 $\mid$ 100 |  200 $\mid$ 200 | 10 | 212 | 28.73 | 3.286e+06 | 11 | 48.13 | 3.286e+06
randSDP-m200-n300-N10 |  200 $\mid$ 200 |  300 $\mid$ 300 | 10 | 309 | 128.05 | 1.646e+07 | 13 | 379.36 | 1.646e+07
randSDP-m200-n400-N10 |  200 $\mid$ 200 |  400 $\mid$ 400 | 10 | 280 | 190.98 | 3.907e+07 | 12 | 773.04 | 3.907e+07
randSDP-m200-n400-N20 |  200 $\mid$ 200 |  400 $\mid$ 400 | 20 | 633 | 680.14 | 3.906e+07 | 12 | 2042.76 | 3.906e+07

In Table LABEL:table:randSDP, the number of scenarios $N$ is small. When the
number of scenarios $N$ is fixed and the dimensions of the decision variables
increase gradually, we can observe that the results are consistent with the
SDP relaxation problems in section 7.3. When $m_{0},m_{i},n_{0},n_{i}$
increase further, our DBA-sGS-ADMM is expected to outperform MOSEK by an even
larger margin.

##  8 Conclusion

We have proposed efficient augmented Lagrangian based decomposition algorithms
for solving convex composite programming problems with dual block-angular
structures. We show that various application problems such as 2-stage
stochastic convex programming problems, and DNN relaxations of uncapacitated
facility location problems can be recast into our framework. Numerical
experiments show that our method is especially competitive against MOSEK or
Gurobi for solving SDP problems and convex QP problems with dual block-angular
structures. Our algorithms are also demonstrated to be much more competitive
than the well-known PHA in solving 2-stage stochastic linear programming
problems. As mentioned in the paper, some of the subproblems can be solved in
a parallel manner. So as a natural future work, we may implement our algorithm
on a conducive parallel computing environment to further improve the
efficiency of our algorithms.

## Acknowledgements

We would like to thank the two anonymous referees for their careful reading of
this paper, and their insightful comments and suggestions which have helped to
improve the quality of this paper.

## Appendix 1: two-stage stochastic programming datasets

assets.

This network model represents the management of assets. Its nodes are asset
categories and its arcs are transactions. The problem is to maximize the
return of an investment from every stage with the balance of material at each
node.

env.

This model assists the Canton of Geneva in planning its energy supply
infrastructure and policies. The main objective is to minimize the
installation cost of various types of energy, while meeting all the supply-
demand at each node and satisfying several realistic constraints such as the
environmental constraints.

phone.

This is a problem which models the service of providing private lines to
telecommunication customers, often used by large corporations between business
locations for high speed, private data transmission. The goal is to minimize
the unserved requests, while staying within budget.

AIRL.

This dataset is used to schedule monthly airlift operations. The airlift
operation is scheduled so that the restriction on each aircraft such as the
number of flight hours available during the month can be satisfied. The cost
such as the available flight time to go unused, switching aircraft from one
route to another and buying commercial flights is the main objective to be
minimized.

4node.

This is a two stage network problem for scheduling cargo transportation. While
the flight schedule is completely determined in stage one, the amounts of
cargo to be shipped are uncertain and shall be determined in stage two.

pltexp.

This is a stochastic capacity expansion model that tries to allocate new
production capacity across a set of plants so as to maximize profit subject to
uncertain demand.

storm.

This is a two period freight scheduling problem described in Mulvey and
Ruszczynski. In this model, routes are scheduled to satisfy a set of demands
at stage 1, demands occur, and unmet demands are delivered at higher costs in
stage 2 to account for shortcomings.

gbd.

This is the aircraft allocation problem where aircraft of different types are
to be allocated to routes in a way that maximizes profit under uncertain
demand, and minimizes the cost of operating the aircraft as well as costs
associated with bumping passengers when the demand for seats outstrips the
capacity.

## Appendix 2: Derivation of solutions of (158) and (166)