A Distributed Primal-Dual Interior-Point Method for Loosely Coupled Problems Using ADMM
Mariette Annergren$^{1
In this paper we propose an efficient distributed algorithm for solving loosely coupled convex optimization problems. The algorithm is based on a primal-dual interior-point method in which we use the alternating direction method of multipliers (ADMM) to compute the primal-dual directions at each iteration of the method. This enables us to join the exceptional convergence properties of primal-dual interior-point methods with the remarkable parallelizability of ADMM. The resulting algorithm has superior computational properties with respect to ADMM directly applied to our problem. The amount of computations that needs to be conducted by each computing agent is far less. In particular, the updates for all variables can be expressed in closed form, irrespective of the type of optimization problem. The most expensive computational burden of the algorithm occur in the updates of the primal variables and can be precomputed in each iteration of the interior-point method. We verify and compare our method to ADMM in numerical experiments.
main.tex
% gOMSguide.tex
% v4.4 released October 2008

\documentclass[10pt,a4paper]{article}



%\usepackage[caption=false]{subfig}
%\usepackage[dvips]{graphicx}
%\graphicspath{{Figures/}}
\usepackage{amsmath, amssymb}
%\usepackage{tikz} \usetikzlibrary{calc}
%\usepackage{psfrag}
\usepackage{pgfplots}
\usepackage{algorithm}     % <-- should be replaced with 'algorithms'
\usepackage{algpseudocode} % <-- should be replaced with 'algorithms'
\usepackage{arydshln}

\usepackage[square,numbers]{natbib}
\usepackage{color}
\usepackage{pstool}
\usepackage{graphicx}
\usepackage{psfrag}
\usepackage{subfigure}
\usepackage{tikz}
\usepackage{pgfplots}  

\usepackage{authblk}

\usepackage{amsthm}


\DeclareMathOperator*{\minimize}{minimize}
\DeclareMathOperator*{\maximize}{maximize}
\DeclareMathOperator*{\trace}{tr}
\DeclareMathOperator*{\supremum}{sup}
\DeclareMathOperator*{\lne}{ln}
\DeclareMathOperator*{\eigen}{eig}
\DeclareMathOperator*{\minimum}{minimum}
\DeclareMathOperator*{\subject}{subject\ to}
\DeclareMathOperator*{\determin}{det}
\DeclareMathOperator*{\maximum}{max}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator*{\argmax}{argmax}
\DeclareMathOperator*{\diag}{diag}
\DeclareMathOperator*{\relint}{relint}
\DeclareMathOperator*{\vectorize}{vec}
\DeclareMathOperator*{\infimum}{inf}
\DeclareMathOperator*{\fix}{fix}
\DeclareMathOperator*{\find}{find}
\DeclareMathOperator*{\adj}{adj}
\DeclareMathOperator*{\parent}{par}
\DeclareMathOperator*{\IQC}{IQC}
\DeclareMathOperator*{\children}{ch}
\DeclareMathOperator*{\degree}{deg}
\DeclareMathOperator*{\real}{Re}
\DeclareMathOperator*{\logdet}{log\ det}
\DeclareMathOperator*{\dom}{dom}
\DeclareMathOperator*{\ridom}{ri \ dom}
\DeclareMathOperator*{\rank}{rank}
\DeclareMathOperator*{\blkdiag}{blkdiag}
\DeclareMathOperator*{\leaves}{leaves}
\DeclareMathOperator*{\Ne}{Ne}



\setlength{\unitlength}{1.2mm}
\newcommand{\st}{\mbox{subj. to}}
\newcommand{\nt}{\text{nt}}
\newcommand{\Acal}{\mathcal{A}}
\newcommand{\Rcal}{\mathbb{R}}
\newcommand{\Bcal}{\mathcal{B}}
\newcommand{\Ccal}{\mathcal{C}}
\newcommand{\Dcal}{\mathcal{D}}
\newcommand{\Ecal}{\mathcal{E}}
\newcommand{\Gcal}{\mathcal{G}}
\newcommand{\Hcal}{\mathcal{H}}
%\newcommand{\qed}{\hfill \ensuremath{\Box}}
%\newcommand{\Rcal}{\mathcal{R}}
\newcommand{\Ucal}{\mathcal{U}}
\newcommand{\Ycal}{\mathcal{Y}}
\newcommand{\Wcal}{\mathcal{W}}
\newcommand{\PF}{\bf}
\newcommand{\FB}{\bf}
\newcommand{\new}{\color{Blue}}
\newcounter{ass}
\newcounter{lem}
\newcounter{ex}
\newcounter{cor}
%\newcounter{pro}
\newcounter{thm}
%\newcounter{rem}
\newtheorem{example}[ex]{Example}
\newtheorem{assumption}[ass]{Assumption}
\newtheorem{corollary}[cor]{Corollary}
\newtheorem{lemma}[lem]{Lemma}
\newtheorem{theorem}[thm]{Theorem}
\newtheorem{simplification}{Simplification}
%\newtheorem{proof}[pro]{Proof}
%\newtheorem{rem}[rem]{Remark}
\newcounter{remcount}
\newcounter{def}
\newtheorem{definition}[def]{Definition}
\definecolor{red}{rgb}{1,0,0}
%\newcommand{\etal}{\emph{et~al.}}
\newcommand{\ones}{\mathbf{1}}
\newcommand{\dist}{\mathrm{dist}}


%% TikZ stuff %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\figuretikz}[3]{\pgfplotsset{width=#1\columnwidth,height=#1\columnwidth,compat=newest,plot coordinates/math parser=false}\input{#1}}
\usetikzlibrary{shapes,arrows}
\tikzstyle{block} = [draw, fill=white, rectangle,
    minimum height=2.8em, minimum width=4.5em, node distance = 8em]
\tikzstyle{sum} = [draw, fill=white, circle, node distance=1cm, minimum size = 0.3cm]
\tikzstyle{input} = [coordinate]
\tikzstyle{output} = [coordinate]
\tikzstyle{pinstyle} = [pin edge={to-,thin,black}]
\tikzstyle{block2} = [draw,minimum width = 6em, minimum height=2.8em,text centered]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





\begin{document}

\title{A Distributed Primal-Dual Interior-Point Method for Loosely Coupled Problems Using ADMM} 

\author[*]{Mariette Annergren} 
\author[**]{Sina Khoshfetrat Pakazad} 
\author[**]{Anders Hansson}
\author[*]{Bo Wahlberg} 


\affil[*]{Department of Automatic Control, KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden. Email:\{mariette.annergren, bo.wahlberg\}@ee.kth.se}                                         
\affil[**]{Division of Automatic Control, Department of Electrical Engineering, Link\"oping University, SE-581 83 Link\"oping, Sweden. Email: \{sina.kh.pa, hansson\}@isy.liu.se.} 
%\author{Mariette Annergren$^{1}$\thanks{$^{1}$ Mariette Annergren and Bo Wahlberg are with the Department of Automatic Control, KTH Royal Institute of Technology, Stockholm, Sweden. Email:  \{mariette.annergren, bo.wahlberg\}@ee.kth.se.}, Sina Khoshfetrat Pakazad$^{2}$\thanks{$^{2}$ Sina Khoshfetrat Pakazad and Anders Hansson are with the Division of Automatic Control, Department of Electrical Engineering, Link\"oping University, Sweden. Email: \{sina.kh.pa, hansson\}@isy.liu.se.}, Anders Hansson$^{2}$ and Bo Wahlberg$^{1}$}\"{a}

\maketitle

\begin{abstract}
In this paper we propose an efficient distributed algorithm for solving loosely coupled convex optimization problems. The algorithm is based on a primal-dual interior-point method in which we use the alternating direction method of multipliers (ADMM) to compute the primal-dual directions at each iteration of the method. This enables us to join the exceptional convergence properties of primal-dual interior-point methods with the remarkable parallelizability of ADMM. The resulting algorithm has superior computational properties with respect to ADMM directly applied to our problem. The amount of computations that needs to be conducted by each computing agent is far less. In particular, the updates for all variables can be expressed in closed form, irrespective of the type of optimization problem. The most expensive computational burden of the algorithm occur in the updates of the primal variables and can be precomputed in each iteration of the interior-point method. We verify and compare our method to ADMM in numerical experiments.
\end{abstract}

%\begin{keywords}distributed optimization algorithm, loosely coupled problems, interior-point methods, alternating direction method of multipliers.\end{keywords}
%
%\begin{classcode}65K\end{classcode}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We are interested in solving convex optimization problems of the form
%
\begin{subequations}\label{eq:EOP}
\begin{align}
\minimize_x  \quad &f_1(x)+ \dots + f_N(x),\\
\subject \quad & G^i(x) \preceq 0, \quad i=1, \dots, N, \\
  & A^i x = b^i, \hspace{6mm} i=1, \dots, N,
\end{align}
\end{subequations}
%
where $f_i \colon \mathbb R^{n}\rightarrow \mathbb R$, $G^i \colon \mathbb R^{n}\rightarrow \mathbb R^{m_i}$ and $A^i \in \mathbb R^{p_i \times n}$ with $p_i < n$ and $\rank(A^i) = p_i$ for all $i = 1, \dots, N$. We assume that the function pairs $f_i$, $G^i$ and their corresponding $A^i$ for $i = 1, \dots, N$, depend only on a small subset of the elements of the variable $x$ and we denote the ordered set of the indices of these variables by $J_i$. We also denote  the ordered set of indices of triplets $f_i$, $G^i$, $A^i$ that depend on $x_i$ by $\mathcal I_i$, i.e., $\mathcal{I}_i = \{k \ | \ i \in J_k\}$. An optimization problem is called loosely coupled if $|\mathcal I_i|\ll N$ for all $i = 1, \dots, n$. We can explicitly express the coupling structure in \eqref{eq:EOP} using the so-called consistency or consensus constraints.

Centralized algorithms for solving optimization problems of the form \eqref{eq:EOP} can be unviable. This can be due to lack of powerful enough centralized computational units, or because the problem cannot be formed as a centralized optimization problem due to its structural constraints, such as privacy requirements. A sensible approach for circumventing such issues is to use distributed optimization algorithms, which rely on collaboration of multiple computing agents to solve the problem. In such a setting, each agent is assigned a local subproblem, and at every iteration it solves its subproblem and communicates or collaborates with certain other agents. This is done repeatedly until the network of agents arrives or agrees on a solution.

Distributed optimization methods have been studied for many years, and there are different approaches for devising such algorithms, see e.g. \cite{ber:97,eck:89,boyd:11,ned:09,ned:10}. One of the most common approaches for designing distributed algorithms is to apply first order or proximal point methods directly to the problem or some reformulation of it. In this class of distributed algorithms, the ones based on subgradient or gradient methods are perhaps among the simplest, see e.g. \cite{ned:09,ned:10}. The local computations that need to be performed by each agent are usually elementary. However, these algorithms are very sensitive to the scaling of the problem. They also generally require many iterations to converge to a solution with even medium accuracy, \citep{ber:97}. In order to alleviate these issues there has been a surge of interest to devise distributed algorithms based on proximal point methods, e.g. see \cite{ber:97,eck:89,boyd:11,com:11}. For certain classes of problems, for instance when the objective function of the equivalent unconstrained reformulation of the problem has two terms and/or is strongly convex, such algorithms commonly enjoy better convergence properties, \citep{gol:12,gols:12} and are less sensitive to the scaling of the problem. However, they are generally more complicated in that the local computational burden is higher, and the communication protocols are more sophisticated, see e.g. \cite{sum:12,ohl:13}. Moreover, extra care must be taken if one wishes to apply proximal point methods to more general  classes of problems, as these algorithms might even diverge, see e.g. \cite{che:13}. There have been suggestions on how to modify these methods to allow application to more general problems. However, the resulting algorithms can become overly complicated to implement, particularly in a distributed fashion, \citep{gol:12,han:12,hon:12}.

Another approach for designing distributed optimization algorithms is to use second order methods, e.g. see \cite{Chu:11,wei:13,nec:09}. For instance, in \cite{nec:09} the authors propose a distributed optimization method based on an interior-point method. The introduced algorithm is obtained by first performing a Lagrangian decomposition of the problem and then efficiently solving the subproblems using interior-point methods. However, in the proposed algorithm, the computational cost for solving the subproblems can still be considerable. The authors in \cite{Chu:11} propose a distributed Newton method for solving coupled unconstrained quadratic problems, which is used for anomaly detection in large populations. This distributed method is only applicable to unconstrained quadratic problems. In \cite{wei:13} a distributed Newton method for solving a network utility maximization problem is proposed. The cost function for such problems  is given by a summation of several terms where each term depends on a single scalar variable. This structure allows the authors to employ a matrix splitting method which in turn enables them to distribute the computations of the inexact Newton directions. However, this method relies on the special structure in the considered problem and hence can only be used in particular cases.

The approach presented in the latter paper falls in the class of inexact interior-point methods which have been studied thoroughly over the past two decades, e.g. see \cite{bel:98,fre:99,miz:99,han:00,kor:00,dur:03,bel:04,zho:04,bon:05,bon:07,toh:08,lu:09,al:09,cur:10,cur:12}. These methods combine primal or primal-dual interior point methods with iterative algorithms for solving linear systems of equations. This is motivated by the fact that we need to solve a linear system of equations in every iteration of a primal or primal-dual interior-point method, in order to compute primal or primal-dual directions. These methods provide bounds on the required accuracy of the computed directions at each iteration in order to guarantee convergence. The papers \cite{fre:99} and \cite{miz:99} consider Linear Programs (LPs) and focus on the design of these accuracy bounds. In particular, they provide bounds on primal and dual residuals and computed directions to assure convergence of their respective proposed inexact interior-point method. LPs are also considered in \cite{kor:00} where the author proposes an inexact interior-point method with Quasi-Minimal Residual (QMR) technique and Conjugate Gradient (CG) as inexact solvers of choice. Also in \cite{al:09}, the authors consider LPs and they focus on devising efficient pre-conditioners for CG algorithms for solving the underlying linear equations more efficiently, so-called Preconditioned Conjugate Gradient (PCG) algorithms. An inexact primal-dual method for solving robust optimal control problems is proposed in \cite{han:00} with QMR as the iterative solver of choice. The papers \cite{bel:04} and \cite{zho:04} consider semidefinite programs and propose inexact primal-dual interior point methods for solving the problem. The inexact solvers in these papers were PCG for which they both propose efficient pre-conditioners to improve the convergence properties. They also propose similar accuracy bounds on the computed directions that depend solely on the so-called complementarity gap. In \cite{toh:08} a quadratic semi-definite program is considered where the author uses a pre-conditioned QMR algorithm and proposes efficient pre-conditioners for further improvement of its convergence rate. Inexact interior-point methods have also been used for solving constrained nonlinear systems of equations, which can be considered as Karush-Kuhn-Tucker (KKT) optimality conditions for general optimization problems (not necessarily convex). For instance \cite{bel:98} proposes an inexact interior-point method for solving constrained nonlinear monotone systems of equations, under the assumption that the Jacobian of the system of equations is invertible at the solution. The authors in \cite{dur:03} put forth a similar framework for solving general constrained nonlinear systems of equations and they use the PCG algorithm for solving them with respect to primal-dual directions. In \cite{bon:05} an inexact interior-point method for solving constrained nonlinear system of equations is proposed which uses the so-called Hestenes' multipliers method for solving the underlying linear systems of equations. The authors further investigate the numerical properties of the proposed method and compare with the case when they use PCG as the iterative solver of choice in \cite{bon:07}.

Notice that design of distributed algorithms for solving optimization problems was not the focus of any of the works discussed in the previous paragraph. In this paper, we focus on devising a distributed optimization algorithm based on a primal-dual interior-point method for solving loosely coupled optimization problems. These constitute a more general class of problems than those considered by \cite{Chu:11,wei:13,nec:09}. To this end, we first exploit the coupling in the problem using consistency constraints and use proximal splitting methods, particularly Alternating Direction Method of Multipliers (ADMM), to compute the primal-dual directions in a distributed manner.

ADMM is a method for finding saddle points of an augmented Lagrangian and, as such, a method of finding a solution of an optimization problem \citep{gab:76}. In our approach, we use ADMM to solve the KKT conditions of a particular optimization problem that has the primal-dual directions as solution, see Section \ref{sec:PDIPM}. The benefits of using ADMM are several. The ADMM iterations
\begin{itemize}
\item converge to a solution under mild assumptions \citep{boyd:11}.
\item enable the solution to be calculated in a highly distributed way, see Section \ref{distributedmethod}.
\item consist of subproblems that are extremely cheap to solve, see Section \ref{distributedmethod}.
\end{itemize}

ADMM was first introduced in \cite{mar:75} for solving nonlinear Dirichlet problems. It was presented as a modified version of Uzawa's algorithm \citep{arr:64}. The method was developed further in \cite{gab:76}, where some convergence properties were stated.  In \cite{gab:83}, it was shown that ADMM is equivalent to Douglas-Rachford splitting for monotone operators \citep{dou:56} and similar to Peaceman-Rachford splitting \citep{pea:55}. ADMM is related to the method of multipliers, also known as Hestenes' multipliers method, \citep{ hes:69,pow:69}, and the proximal point algorithm \citep{eck:92}. For a detailed overview of ADMM and other related methods, see \cite{boyd:11}.

Our proposed distributed optimization algorithm has superior computational properties than other distributed solvers, and we believe that the key to achieving this has been the use of ADMM for computing the primal-dual directions. We are not aware of any other iterative solvers that would present the same characteristics as listed above. We illustrate the performance of the proposed algorithm using a numerical experiment.

\subsection*{Contribution}
We present a novel distributed optimization algorithm for solving loosely coupled problems of the form \eqref{eq:EOP}. The algorithm is a primal-dual interior-point method where ADMM is used to calculate the search direction in a distributed fashion. We also present an inexact version of the algorithm, where the search directions are calculated with an adaptive degree of accuracy. We formulate the conditions under which the inexact algorithm converges to a solution along with a formal proof thereof. In addition, we review how ADMM relates to other methods of solving linear system of equations in general and for our problem formulation specifically.

Our method exhibits several important qualities. Specifically, the method
\begin{itemize}
\item inherits the convergence properties of the primal-dual interior-point \\method.
\item inherits the ability of ADMM to distribute calculations.
\item has cheap search direction calculations, where the most expensive computational burden can be precomputed.
\end{itemize}

\subsection*{Outline}

First we define the notation, and in Section \ref{sec:LCP} we explain the problem formulation and the structure of loose coupling. In Section \ref{sec:PDIPM} we briefly describe a primal-dual interior-point method. We apply this method to loosely coupled problems in Section \ref{sec:DPDIPM} and describe the details of how we can devise a distributed algorithm for solving such problems by using ADMM. In order to increase the efficiency of the proposed algorithm we discuss the use of inexact primal-dual directions in the algorithm in sections \ref{sec:PDIIPM} and \ref{sec:IDP}. We then provide a connection to iterative saddle point solvers in Section \ref{sec:saddle}. Moreover, to further improve the convergence properties of the algorithm, over-relaxation and scaling for the ADMM iterations are briefly discussed in Section~\ref{sec:ADMM}. We illustrate the performance of the proposed algorithm using some numerical experiments in Section~\ref{sec:NE}. Conclusions and future work are stated in Section \ref{sec:C}. In Appendix \ref{app:global}, we provide a proof of global convergence of the inexact algorithm. In Appendix \ref{sec:Appendix2}, we derive the explicit relations between ADMM and Uzawa's method, and ADMM and fixed point iterations for the considered problem formulation.
\subsection*{Notation}
The set of real numbers is denoted by $\mathbb R$. The set of real $n$-dimensional vectors and $n\times m$ matrices are denoted by $\mathbb R^{n}$ and $\mathbb R^{n\times m}$, respectively, and the transpose of a matrix $A$ is denoted by $A^T$. Let $\mathbb{N}_p$ represent the ordered set of positive integers $\{1,2,\ldots,p\}$. Given a set $J \subset \{1,2,\ldots,n\}$, the matrix $E_J \in \mathbb{R}^{|J|\times n}$ is the matrix obtained by deleting the rows indexed by $\mathbb{N}_n \setminus J$ from an identity matrix of order $n$, where $|J|$ denotes the number of elements in set $J$. Consequently, $E_Jx$ is a $|J|$-dimensional vector with the components of $x$ that correspond to the elements in $J$, and we denote this vector $x_J$. We denote by $x^{i,(l)}_k$ the $k$th element of vector $x^i$ at the $l$th iteration. Given vectors $x^i$ and matrices $A^i$ for $i= 1, \dots, N$, the column vector $(x^1, \dots, x^N)$ is all of the given vectors stacked and $\blkdiag(A^1, \dots, A^N)$ represents a block-diagonal matrix with $A^i$ as its diagonal blocks. Similarly, given a vector $x \in \mathbb R^{n}$, $\diag(x_1, \dots, x_n)$ denotes a diagonal matrix with its diagonals  expressed by elements of $x$. The vector $\hat{e}$ is a vector of ones of appropriate dimensions given by the context. The minimum value of a set or of a function is denoted by ``$\min$'' and the minimizing argument of an optimization problem in denoted by ``$\argmin$''. The inequality $x\preceq y$, where $x,y\in \mathbb R^{n}$, means $x_i\leq y_i$ for $i=1,\ldots,n$. The standard uniform distribution over interval $[a,b]$ is denoted $U(a,b)$. To simplify notation we introduce
\begin{align*}
z^{(l)}&=(x^{(l)},s^{(l)},\lambda^{(l)},v^{(l)}),\\
\Delta z&=(\Delta x,\Delta s,\Delta \lambda,\Delta v),\\
z^{i,(l)}&=(x_{J_i}^{(l)},s^{i,(l)},\lambda^{i,(l)},v^{i,(l)}).
\end{align*}
Here $x$ denotes all primal variables, $s$ is the slack variable vector, $\lambda$ is the dual variable vector corresponding to inequality constraints and $v$ is the dual variable vector corresponding to equality constraints. When formulating a loosely coupled problem, we introduce an additional primal variable $w$ and an additional dual variable $v_c$, both correspond to the consistency constraint in the coupled problem. Given $z^{(l)}$ and $\Delta z$, we also define $f^{(l)}(\alpha)$ as $f(z^{(l)}+\alpha\Delta z)$ which yields $f^{(l)}(0)=f(z^{(l)})$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Loosely Coupled Problems}\label{sec:LCP}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 The problem in \eqref{eq:EOP} can be equivalently written as
%
\begin{subequations}\label{eq:DDEOP}
\begin{align}
\minimize_{W, x} \quad & \bar f_1(w^1) + \dots + \bar f_N(w^N),\\
\subject \quad & \bar G^i(w^i) \preceq 0,  \quad  i = 1, \dots, N,\\
  & \bar A^i w^i = b^i, \hspace{6mm} i=1, \dots, N, \label{eq:DEOPc}\\
  & \bar E x = W,\label{eq:DEOPb}
\end{align}
\end{subequations}
%
where $W = (w^1, \dots, w^N)$ and $\bar{E}  = \begin{bmatrix} E_{J_1}^T   &\cdots &  E_{J_N}^T \end{bmatrix}^T$ with $E_{J_i}$ as a $0$--$1$ matrix that is obtained from an identity matrix of order $n$ by deleting the rows indexed by $\mathbb{N}_n \setminus J$. We refer to the constraints in \eqref{eq:DEOPb} as consistency constraints. The functions $\bar f_i \colon \mathbb R^{|J_i|} \rightarrow \mathbb R$ are lower dimensional descriptions of the functions $f_i$ such that $f_i(x) = \bar f_i(E_{J_i}x)$ for all $x \in \mathbb R^n$ and $i = 1, \dots, N$. In this formulation, the functions $\bar G^i \colon \mathbb R^{|J_i|} \rightarrow \mathbb R^{m_i}$ are defined in the same manner as the functions $\bar f_i$, and the matrices $\bar A^i \in \mathbb R^{p_i \times |J_i|}$ are defined by removing unnecessary columns from $A^i$. We further assume that $p_i < |J_i|$ and that $\rank (\bar A^i) = p_i$ for all $i = 1, \dots, N$. In this paper, we intend to devise algorithms to solve problems of the form in \eqref{eq:EOP} or \eqref{eq:DDEOP} in a distributed manner, and we will investigate the possibility of using primal-dual interior-point methods, both exact and inexact. To ensure global convergence of our algorithm, when using an inexact interior point method, we further make the standard assumption that the functions $f_i$ and $G^i$, in addition to being convex, have Lipschitz continuous derivatives, see section \ref{distConRes}. Next, we briefly review primal-dual interior-point methods for solving convex problems.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Primal-Dual Interior-Point Methods}\label{sec:PDIPM}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%We layout a primal-dual interior-point method in Algorithm \ref{alg:alg1a}. The primal-dual step direction and step size calculations are described in the following.

Let us consider the convex optimization problem
%
\begin{equation}\label{eq:ConvexIneq}
\begin{split}
\minimize & \quad  F(x)\\
\subject & \quad  g_i(x) \leq 0, \quad i = 1, \dots, m, \\
&\quad  Ax = b,
\end{split}
\end{equation}
%
where $F \colon \mathbb R^n \rightarrow \mathbb R$, $g_i \colon \mathbb R^n \rightarrow \mathbb R$ and $A \in \mathbb R^{p \times n}$ with $p<n$ and $\rank (A) = p$. We introduce slack variables $s\in\mathbb R^m$ and reformulate \eqref{eq:ConvexIneq} as
%
\begin{equation}\label{eq:ConvexIneqSlack}
\begin{split}
\minimize & \quad  F(x)\\
\subject & \quad  g_i(x) +s_i= 0, \quad i = 1, \dots, m, \\
&\quad  Ax = b,\\
&\quad  s_i \geq 0,\quad i = 1, \dots, m.
\end{split}
\end{equation}
%
The problem in \eqref{eq:ConvexIneq} is equivalent to \eqref{eq:ConvexIneqSlack}. This means that $x$ is optimal for \eqref{eq:ConvexIneq} if and only if $(x,s)$ is optimal for \eqref{eq:ConvexIneqSlack} with $s_i=-g_i(x)$ for $i=1,\dots,m$, see \cite{boyd:04}. The KKT optimality conditions for Problem \eqref{eq:ConvexIneqSlack} can be written as
%
\begin{subequations}\label{eq:ConvexIneqKKTP1}
\begin{align}
\nabla F(x)  + \sum_{i = 1}^m \lambda_i\nabla g_i(x) + A^T v&= 0,\\ \lambda_i&\geq 0, \quad i= 1, \dots, m, \\ s_i&\geq 0, \quad i= 1, \dots, m, \\  \lambda_i s_i &= 0 , \quad i= 1, \dots, m,  \label{eq:ConvexIneqKKTP1-d}\\  g_i(x) +s_i&= 0, \quad i= 1, \dots, m,\\A x &= b.
\end{align}
\end{subequations}
%
The conditions are equivalent to those obtained for \eqref{eq:ConvexIneq} if $s_i$ is exchanged with $-g_i(x)$ for $i=1,\dots,m$. Primal-dual methods solve the problem in \eqref{eq:ConvexIneqSlack} by dealing with a sequence of modified versions of the optimality conditions in \eqref{eq:ConvexIneqKKTP1} where we perturb~\eqref{eq:ConvexIneqKKTP1-d}~as $\lambda_i s_i = \mu$ with $\mu>0$. Particularly, in a primal-dual framework and at each iteration, we get the primal and dual search directions by linearizing the perturbed KKT conditions and solving the resulting set of linear equations with respect to the search directions. The perturbed KKT conditions for \eqref{eq:ConvexIneqSlack} are
%
\begin{subequations}\label{eq:ConvexIneqKKTPD}
\begin{align}
\nabla F(x)  + \sum_{i = 1}^m \lambda_i\nabla g_i(x) + A^T v&= 0,\\ \lambda_i s_i &= \mu , \quad i= 1, \dots, m,  \label{eq:ConvexIneqKKTP-d}\\ g_i(x) +s_i&= 0, \quad i= 1, \dots, m,\\  A x &= b,
\end{align}
\end{subequations}
%
with $\lambda_i > 0$ and $s_i > 0$ for all $i =1, \dots, m$. Given the primal and dual iterates $x^{(l)}$, $s^{(l)}$, $\lambda^{(l)}$ and $v^{(l)}$ at iteration $l$ such that $\lambda_i^{(l)} > 0$ and $s_i^{(l)} > 0$ for all $i =1, \dots, m$, we linearize \eqref{eq:ConvexIneqKKTPD} which results in
%
\begin{subequations}\label{eq:QAppIneqKKTPD}
\begin{align}
\begin{split}
\!\left(\!\!\nabla^2 F(x^{(l)}) + \sum_{i = 1}^m \lambda_i^{(l)}\nabla^2 g_i(x^{(l)})\!\right)\!\Delta x+& \\\sum_{i = 1}^m \nabla g_i(x^{(l)})\Delta \lambda_i + A^T\Delta v &= -r^{(l)}_{\text{dual}},\end{split}\\
\label{eq:QAppIneqKKTPD1}
 -\lambda_i^{(l)} \Delta s_i +s_i^{(l)}\Delta \lambda_i &= -\left(r^{(l)}_{\text{cent}}\right)_i+\mu^{(l)},
  i= 1, \dots, m,\\
\label{eq:QAppIneqKKTPD2}
 \nabla g_i(x^{(l)})^T \Delta x +\Delta s_i&= -\left(r^{(l)}_{\text{primal,1}}\right)_i,
  i= 1, \dots, m, \\
 A \Delta x &= -r^{(l)}_{\text{primal,2}},
\end{align}
\end{subequations}
%
where
%
\begin{subequations}\label{eq:res}
\begin{align}
r_{\text{dual}}^{(l)} &= \nabla F(x^{(l)}) +  \sum_{i = 1}^m \lambda^{(l)}_i\nabla g_i(x^{(l)}) + A^T v^{(l)},\\
\left(r_{\text{cent}}^{(l)} \right)_i &= \lambda^{(l)}_i s_i^{(l)}, \quad i=1, \dots, m,\\
\left(r^{(l)}_{\text{primal,1}}\right)_i& = g_i(x^{(l)}) + s_i^{(l)}, \quad i=1, \dots, m,\\
r_{\text{primal,2}}^{(l)}& = Ax^{(l)} - b.
\end{align}
\end{subequations}
%
The linearized KKT conditions in~\eqref{eq:QAppIneqKKTPD} can be written in a compact form as
%
\begin{align}\label{eq:QAppIneqKKTPDCom}
H'(z^{(l)})\Delta z=- H(z^{(l)})+ \mu^{(l)}\begin{bmatrix}0 \\ 0\\ 0  \\\hat{e} \end{bmatrix},
\end{align}
%
where
%
\begin{align*}
H'(z^{(l)})=&\begin{bmatrix}\nabla^2 F(x^{(l)}) + \sum_{i = 1}^m \lambda_i^{(l)}\nabla^2 g_i(x^{(l)}) & 0&Dg(x^{(l)})^T & A^T \\  Dg(x^{(l)}) & I&0& 0 \\ A & 0 & 0&0\\0&\Lambda^{(l)}&S^{(l)}&0\end{bmatrix},\end{align*}
and
\begin{align*}
H(z^{(l)})=& \begin{bmatrix} r^{(l)}_{\text{dual}} \\ r^{(l)}_{\text{primal,1}}  \\ r^{(l)}_{\text{primal,2}}\\ r^{(l)}_{\text{cent}} \end{bmatrix},
\end{align*}
%
with
\begin{subequations}
\begin{align*}
Dg(x) &= \begin{bmatrix} \nabla g_1(x) & \dots & \nabla g_m(x) \end{bmatrix}^T\in\mathbb R^{m\times n},\\
\Lambda^{(l)}&=\diag(\lambda_1^{(l)}, \dots, \lambda_m^{(l)}) \in\mathbb R^{m\times m},\\
S^{(l)}&=\diag(s_1^{(l)}, \dots, s_m^{(l)})\in\mathbb R^{m\times m},\\
\hat{e}&=(1 \dots 1)^T \in\mathbb R^{m}.
\end{align*}
\end{subequations}
We assume that $H'(z^{(l)})$ is nonsingular, which is a standard assumption in an interior-point method. One way to solve \eqref{eq:QAppIneqKKTPDCom} is by first eliminating $\Delta s$ and $\Delta \lambda$ as
%
\begin{subequations}
\label{eq:PDsolution}
\begin{align}
\label{eq:PDS}
\Delta s &= -Dg(x^{(l)})\Delta x -   r^{(l)}_{\text{primal,1}},\\
\label{eq:PDLambda}
\Delta \lambda &= -(S^{(l)})^{-1}\left( \Lambda^{(l)}\Delta s -   r^{(l)}_{\text{cent}}+\mu^{(l)}\hat{e} \right).
\end{align}
\end{subequations}
%
We can then rewrite \eqref{eq:QAppIneqKKTPDCom} as
%
\begin{align}\label{eq:PD}
\begin{bmatrix} H^{(l)}_{\text{pd}} & A^T \\ A & 0 \end{bmatrix}\begin{bmatrix} \Delta x \\ \Delta v \end{bmatrix}  = - \begin{bmatrix}r^{(l)}\\ r^{(l)}_{\text{primal,2}}     \end{bmatrix}
\end{align}
%
where
%
\begin{align*}
H_{\text{pd}}^{(l)} = \nabla^2 F(x^{(l)}) + \sum_{i = 1}^m \lambda_i^{(l)}\nabla^2 g_i(x^{(l)})+ \sum_{i=1}^m \frac{\lambda_i^{(l)}}{s_i^{(l)}} \nabla g_i(x^{(l)})\nabla g_i(x^{(l)})^T,
\end{align*}
%
and
%
\begin{multline*}
r^{(l)} = r^{(l)}_{\text{dual}} + Dg(x^{(l)}) ^T(S^{(l)})^{-1} r^{(l)}_{\text{cent}}-Dg(x^{(l)}) ^T(S^{(l)})^{-1} \mu^{(l)}\hat{e}+\\Dg(x^{(l)}) ^T(S^{(l)})^{-1} \Lambda^{(l)}r^{(l)}_{\text{primal,1}}.
\end{multline*}
%
The key observation for our proposed algorithm is that the set of equations in \eqref{eq:PD} express the optimality conditions for the quadratic program

\begin{align}\label{eq:PDQP}
\begin{split}
\minimize & \quad \frac{1}{2}\Delta x^T H^{(l)}_{pd} \Delta x + (r^{(l)})^T \Delta x,\\
\subject & \quad A \Delta x = - r_{\text{primal,2}}^{(l)}.
\end{split}
\end{align}
%
Hence $\Delta x$ and $\Delta v$ can be computed through solving \eqref{eq:PDQP}. Based on the solution obtained, $\Delta s$ and $\Delta \lambda$ can be calculated using \eqref{eq:PDsolution}. With this, we lay out a primal-dual interior-point method in Algorithm \ref{alg:alg1a}.
%
\begin{algorithm}[tb]
\caption{Primal-Dual Interior-Point Method, \cite[]{boyd:04}.}\label{alg:alg1a}
\begin{algorithmic}[1]
\State{Given $l = 0$, $\sigma \in (0\ 1)$, $\epsilon>0$, $\epsilon_{\text{feas}}>0$, $\lambda_i^{(0)} > 0$, $s_i^{(0)} > 0$ for all $i = 1, \dots, m$ and $\hat \eta^{(0)} = \sum_{i=1}^m \lambda_i^{(0)}s_i^{(0)}$.}
\Repeat
\State{Set $\mu = \sigma \hat \eta^{(l)}/m$.}
\State{Given $\mu$, $x^{(l)}$, $s^{(l)}$,  $v^{(l)}$ and $\lambda^{(l)}$  compute $\Delta x^{(l+1)}$, $\Delta s^{(l+1)}$, $\Delta v^{(l+1)}$ and $\Delta \lambda^{(l+1)}$ by solving \eqref{eq:PD} and~\eqref{eq:PDsolution}.}
\State{Compute $\alpha^{(l+1)}$ using line search.}
\State {Update:
\begin{align*}
x^{(l+1)} &=  x^{(l)} + \alpha^{(l+1)}\Delta x^{(l+1)},\\
s^{(l+1)} &=  s^{(l)} + \alpha^{(l+1)}\Delta s^{(l+1)},\\
\lambda^{(l+1)}& =  \lambda^{(l)} + \alpha^{(l+1)}\Delta \lambda^{(l+1)},\\
v^{(l+1)} &=  v^{(l)} + \alpha^{(l+1)}\Delta v^{(l+1)},\\
l &= l + 1.
\end{align*}}
\State{Set $\hat \eta^{(l)} =  \sum_{i=1}^m \lambda_i^{(l)}s_i^{(l)}$.}
\Until{$\| (r^{(l)}_{\text{primal,1}}, r^{(l)}_{\text{primal,2}}) \|\leq \epsilon_{\text{feas}}$, $\| r^{(l)}_{\text{dual}} \| \leq \epsilon_{\text{feas}}$ and $\hat \eta^{(l)} \leq \epsilon$.}
\end{algorithmic}
\end{algorithm}
%
%
\subsubsection*{Remark 1}
\emph{We do not use \eqref{eq:QAppIneqKKTPDCom} for computing the primal-dual directions. This is because the coefficient matrix in \eqref{eq:QAppIneqKKTPDCom} is not symmetric, which limits our capability to solve \eqref{eq:QAppIneqKKTPDCom} efficiently. Instead we focus on the linear system of equations in \eqref{eq:PD}, which is sometimes referred to as the \emph{augmented system}. The structure in \eqref{eq:PD}, or equivalently in \eqref{eq:PDQP}, enables us to distribute the computations of primal-dual directions. Another approach to computing the primal-dual directions eliminates $\Delta x$ and $\Delta s$ and then solves a linear set of equations, referred to as the \emph{normal equations}, for computing $\Delta v$. This, however, generally destroys the inherent structure of the problem and inhibits us from devising distributed solutions.}
\newpage
%
%========================================================================
\subsection{Step Size Computations}\label{sec:PDIPMstep}
%========================================================================
We briefly review one of the ways to compute suitable step sizes to ensure convergence of the interior-point method. At each iteration, $l$, in order to have $s^{(l+1)}\succ 0$ and $\lambda^{(l+1)}\succ 0$, we first compute
%
\begin{align*}
\alpha_{\textrm{max}} = \textrm{min} \left\{ 1, \underset{i}{\textrm{min}} \left\{ -\lambda_i^{(l)}/\Delta \lambda_i^{(l+1)}  \ \big | \ \Delta \lambda_i^{(l+1)} < 0 \right\} \right\},
\end{align*}
%
and perform a backtracking line search as\\
%
\begin{algorithmic}
  \While{$\exists \ i  \colon s_i^{(l)} + \alpha^{(l+1)} \Delta s_i^{(l+1)} \leq 0 $}
    \State $\alpha^{(l+1)} = \beta \alpha^{(l+1)}$
  \EndWhile
\end{algorithmic}
%
with $\beta \in (0,1)$ and $\alpha^{(l+1)}$ initialized as $0.99 \alpha_{\textrm{max}}$. In order to guarantee convergence of primal and dual residuals to zero we continue the back tracking as
%
\begin{flushleft}
\begin{algorithmic}
 \While{$\left\| H^{(l)}(\alpha^{(l+1)}) \right\| > (1 - \gamma \alpha^{(l+1)})  \left\| H^{(l)}(0) \right\|$}
    \State  $\alpha^{(l+1)} = \beta \alpha^{(l+1)}$
  \EndWhile
\end{algorithmic}
\end{flushleft}
%
\noindent where $\gamma \in [0.01, 0.1]$. The resulting $\alpha^{(l+1)}$ ensures that the iterates remain feasible and that the norm of the KKT conditions, $\left\| H(z^{(l)}) \right\|$, is decreased consistently after each iteration, \cite{boyd:04}.
%
\subsubsection*{Remark 2}
\emph{The primal-dual method presented in this section is an implementation of a so-called infeasible long-step interior-point method. There are other variants of primal-dual methods, such as short-step, predictor-corrector and Mehrotra's predictor-corrector, that differ in their choice of primal-dual directions. One of the major differences among these variants is in the way they perturb the KKT conditions, i.e., the choice of $\mu$ in~\eqref{eq:ConvexIneqKKTPD}. This means that regardless of the choice of primal-dual interior point method the structure of the coefficient matrix in the resulting linear system of equations remains the same, and hence the discussions that follow can be extended to other variants of primal-dual methods.}
%
\\

\noindent Next, we apply the described primal-dual interior-point method to the loosely coupled problem in \eqref{eq:DDEOP} and propose a distributed algorithm for solving the problem.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{A Distributed Primal-Dual Interior-Point\\ Method For solving Loosely Coupled Problems}\label{sec:DPDIPM}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Let us now apply the primal-dual interior-point method in Algorithm \ref{alg:alg1a} to the problem in \eqref{eq:DDEOP}. As can be seen in Section \ref{sec:PDIPM}, the primal-dual directions computation is at the heart of a primal-dual interior-point method. Hence, the first step in devising a distributed primal-dual interior-point method for solving \eqref{eq:DDEOP} is to distribute the computations of these directions. To this end, we focus on the structure of \eqref{eq:PD} for the problem in~\eqref{eq:DDEOP}, which is given by
\vspace{5pt}
%
{
\renewcommand{\arraystretch}{1.8}
\begin{align}\label{sysOfEq}
\begin{bmatrix} \begin{array} {c:c:c:c}   \bar H^{(l)}_{\textrm{pd}} &0 & \bar A^T & I \\ \hdashline 0 & 0 & 0 & - \bar E^T \\\hdashline \bar A & 0 & 0 & 0 \\ \hdashline I & - \bar E  & 0 & 0 \end{array}\end{bmatrix} \begin{bmatrix} \begin{array}{c} \Delta w^1\\\vdots\\ \Delta w^N \\ \hdashline \Delta x \\ \hdashline \Delta v^1 \\ \vdots \\ \Delta v^N \\ \hdashline \Delta v_c  \end{array}  \end{bmatrix} = -\begin{bmatrix} \begin{array} {c} r^{1,(l)} \\ \vdots \\ r^{N,(l)} \\ \hdashline -\bar E^Tv_c^{(l)} \\ \hdashline r_{\text{primal,2}}^{1,(l)} \\ \vdots \\ r_{\text{primal,2}}^{N,(l)} \\ \hdashline \\[-17pt] r_c^{(l)} \end{array}  \end{bmatrix},
\end{align}
}%
\vspace{5pt}
where $\Delta v$ and $\Delta v_c$ are the dual variable directions for the constraints in \eqref{eq:DEOPc} and \eqref{eq:DEOPb}, respectively; $\bar H^{(l)}_{\textrm{pd}} = \blkdiag\left(H_{\text{pd}}^{1,(l)}, \dots, H_{\text{pd}}^{N,(l)}\right)$ with
%
\vspace{5pt}
\begin{align*}
H_{\text{pd}}^{i,(l)} =& \nabla^2 \bar f_i(w^{i,(l)}) +\sum_{j = 1}^{m_i}\! \lambda^{i,(l)}_j \nabla^2 \bar G_j^i(w^{i,(l)})+ \\ &\sum_{j=1}^{m_i} \frac{\lambda^{i,(l)}_j}{s^{i,(l)}}\nabla \bar G_j^i(w^{i,(l)})\left(\nabla \bar G_j^i(w^{i,(l)})\right)^T,
\end{align*}
\vspace{5pt}
%
$\bar A = \blkdiag\left(\bar A^1, \dots, \bar A^N\right)$ and
%
\vspace{5pt}
\begin{align*}
r^{i,(l)} =& \nabla \bar f_i (w^{i,(l)}) + \sum_{j=1}^{m_i} \lambda^{i,(l)}_j \nabla \bar G_j^i(w^{i,(l)})+ (\bar A^i)^T v^{i,(l)} +  v_c^{i,(l)}+\\&D\bar G^i(w^{i,(l)}) (S^{i,(l)})^{-1} r_{\text{cent}}^{i,(l)}-D\bar G^i(w^{i,(l)}) (S^{i,(l)})^{-1} \mu^{(l)}\hat{e}+\\&D\bar G^i(w^{i,(l)}) (S^{i,(l)})^{-1}\Lambda^{i,(l)} r_{\text{primal,1}}^{i,(l)}, \\
r_{\text{primal,2}}^{i,(l)} = &A^iw^{i,(l)} - b^i,\\
r_c^{(l)} =&  W^{(l)} - \bar E x^{(l)},
\end{align*}
%
\vspace{5pt}
with \begin{align*}r^{i,(l)}_{\text{cent}}& =  \Lambda^{i,(l)}s^{i,(l)},\\ r^{i,(l)}_{\text{primal,1}} &= \bar{G}^i(w^{i,(l)})+s^{i,(l)}.\end{align*}
%
\newpage
The system of equations in \eqref{sysOfEq} coincides with the necessary and sufficient optimality conditions for the optimization problem
%
\begin{subequations}\label{eq:QAppEqLogPD}
\begin{align}
\minimize_{\Delta W, \Delta x} & \quad \sum_{i = 1}^{N} \frac{1}{2}(\Delta w^i)^T H_{\text{pd}}^{i,(l)} \Delta w^i +  (r^{i,(l)})^T\Delta w^i - (v_c^{(l)})^T \bar E \Delta x,\\
\subject & \quad \bar A^i (\Delta w^i + w^{i,(l)}) = b^i, \quad i = 1, \dots, N,\label{eq:QAppEqLogb}\\
& \quad \Delta W - \bar E \Delta x = \bar E x^{(l)}-W^{(l)}. \label{eq:QAppEqLogc}
\end{align}
\end{subequations}
%

\noindent Note that \eqref{eq:QAppEqLogPD} has the same coupling structure as in~\eqref{eq:DDEOP} and can be solved in a distributed way. This enables us to compute the primal-dual directions in a distributed manner. In the following sections, we describe how to distribute the calculation of the search directions, perturbation parameter, step sizes and stopping criteria for the over-all method.
%=================================================================
\subsection{Distributed Primal-Dual Direction Computations}\label{distributedmethod}
%=================================================================
The problem in \eqref{eq:QAppEqLogPD} is of the form
%
\begin{equation}\label{eq:ADMM}
\begin{split}
\minimize_{\Delta W, \Delta x}  & \quad  F_1 (\Delta W) + F_2(\Delta x), \\
\subject & \quad  A\Delta  W + B\Delta x = c,
\end{split}
\end{equation}
%
which can be solved in a distributed fashion using proximal splitting methods, for example ADMM as described in Algorithm \ref{alg:alg2}, \cite{boyd:11}, \cite{com:11}.
%
\begin{algorithm}[tb]
\caption{ADMM, \cite{boyd:11}.}\label{alg:alg2}
Let $\bar{v}$ denote the scaled dual variable, that is, $\bar{v}=(1/\rho)v$.
\begin{algorithmic}[1]
\State{Given $k= 0$, $\rho > 0$, $\epsilon_{\textrm{pri}}>0$, $\epsilon_{\textrm{dual}}>0$, $x^{(0)}$ and $\bar{v}^{(0)}$.}
\Repeat
\State  $\Delta W^{(k+1)} = \argmin_{\Delta W} \left\{ F_1(\Delta W) + \frac{\rho}{2} \| A\Delta  W - B \Delta x^{(k)} - c + \bar{v}^{(k)} \|^2  \right\}$.
\State  $\Delta x^{(k+1)} = \argmin_{\Delta x} \left\{ F_2(\Delta x) + \frac{\rho}{2} \| A\Delta  W^{(k+1)} - B\Delta  x - c + \bar{v}^{(k)} \|^2  \right\}$.
\State  $\bar{v}^{(k+1)} =  \bar{v}^{(k)} + \left(A\Delta  W^{(k+1)} + B\Delta  x^{(k+1)} - c \right)$.
\State{$k = k+1$.}
\If  {$\| A\Delta W^{(k+1)} + B\Delta x^{(k+1)} - c  \|^2 < \epsilon_{\text{pri}}$ and $\| \rho A^TB(\Delta x^{(k+1)} - \Delta x^{(k)}) \|^2 < \epsilon_{\text{dual}}$.}
\State Terminate the algorithm.
\EndIf
\State{$k = k+1$.}
\Until{Algorithm is terminated.}
%
\end{algorithmic}
\end{algorithm}
%
In particular, \eqref{eq:QAppEqLogPD} can be written as
%
\begin{equation}\label{eq:ADMMp}
\begin{split}
\minimize_{\Delta W, \Delta x} &  \hspace{2mm} \underbrace{\sum_{i = 1}^{N}\left( \frac{1}{2}(\Delta w^i)^T H_{\text{pd}}^{i,(l)} \Delta w^i + (r^{i,(l)})^T\Delta w^i\right)}_{F_1(\Delta W)} + \underbrace{(-v_c^{(l)})^T \bar E \Delta x}_{F_2(\Delta x)},\\
\subject &  \\
& \hspace{-3mm}  \underbrace{\begin{bmatrix} \bar A^1 & 0 & \dots & 0\\ 0 & \bar A^2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \bar A^N \\ I & 0 & \dots & 0\\ 0 & I & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & I \end{bmatrix}}_{A} \underbrace{\begin{bmatrix} \Delta w^1 \\ \Delta w^2 \\ \vdots \\ \Delta w^N\end{bmatrix}}_{\Delta W} + \underbrace{\begin{bmatrix} 0 \\0 \\ \vdots \\ 0 \\ -\bar E \end{bmatrix}}_{B} \Delta x = \underbrace{-\begin{bmatrix} r_{\text{primal,2}}^{1,(l)} \\ r_{\text{primal,2}}^{2,(l)} \\ \vdots \\ r_{\text{primal,2}}^{N,(l)} \\ r_c^{(l)} \end{bmatrix}}_{c}.
\end{split}
\end{equation}
%
Applying ADMM to \eqref{eq:ADMMp} results in the following update rules for the primal variable directions:
%
\begin{multline*}
\Delta W^{(k+1)} = \argmin_{\Delta W} \left\{ \sum_{i=1}^{N}  \left( \frac{1}{2}(\Delta w^i)^T H_{\text{pd}}^{i,(l)} \Delta w^i +  (r^{i,(l)})^T\Delta w^i+ \right. \right.\\ \left. \left.\frac{\rho}{2} \| \Delta w^{i} - \Delta x_{J_i}^{(k)} + r_c^{i,(l)} + \Delta \bar v_c^{i,(k)} \|^2+ \frac{\rho}{2} \|\bar A^i\Delta w^{i} + r_{\text{primal,2}}^{i,(l)} + \Delta \bar v^{i,(k)} \|^2 \right)  \right\},
\end{multline*}
%
with $r_c^{i,(l)} = w^{i,(l)} - x_{J_i}^{(l)}$, and
%
\begin{align*}
\Delta x^{(k+1)} =  \argmin_{\Delta x} \left\{ (-v_c^{(l)})^T \bar E \Delta x+   \frac{\rho}{2} \| \Delta W^{(k+1)} - \bar E \Delta x +r_c^{(l)} + \Delta \bar v_c^{(k)} \|^2 \right\},
\end{align*}
%
which results in
%
\begin{align}\label{eq:xPD1}
\Delta x^{(k+1)} =  \left( \bar E^T \bar E  \right)^{-1}\bar E^T \left(  v_c^{(l)} + \Delta W^{(k+1)} +r_c^{(l)} + \Delta \bar v_c^{(k)} \right).
\end{align}
%
Note that the update for $\Delta W$ is highly parallelizable and can be rewritten as
%
\begin{multline*}
\Delta w^{i,(k+1)} = \argmin_{\Delta w^i} \left\{  \frac{1}{2}(\Delta w^i)^T H_{\text{pd}}^{i,(l)} \Delta w^i + \right.\\ \left. (r^{i,(l)})^T\Delta w^i +  \frac{\rho}{2} \| \Delta w^{i} - \Delta x_{J_i}^{(k)} + r_c^{i,(l)} + \Delta \bar v_c^{i,(k)} \|^2+ \right.\\ \left. \frac{\rho}{2} \| \bar A^i\Delta w^{i} + r_{\text{primal,2}}^{i,(l)} + \Delta \bar v^{i,(k)} \|^2  \right\},
\end{multline*}
%
which in turn results in
%
\begin{multline}\label{eq:SPD1}
\Delta w^{i,(k+1)} = -\left[  H^{i,(l)}_{pd} + \rho I + \rho (\bar A^i)^T\bar A^i  \right]^{-1}\times\\  \left[ r^{i,(l)} + \rho \left( r_c^{i,(l)} + \Delta \bar v_c^{i,(k)} -   \Delta x_{J_i}^{(k)} \right)+  \rho (\bar A^i)^T \left( r_{\text{primal,2}}^{i,(l)} + \Delta \bar v^{i,(k)}  \right)  \right],
\end{multline}
%
for $i = 1, \dots, N$. By considering the update in \eqref{eq:xPD1} and the structure in matrix $\bar E$, we see that each Agent $i$ can update their corresponding elements of $\Delta x$ (i.e. $\Delta x_{J_i}$) in a distributed manner, through communication with its neighbors defined by $\text{Ne}(i) = \left\{ j \ | \ J_i\cap J_j \neq \emptyset  \right\}.$ The updates for the dual variable directions are given by
%
\begin{align}\label{dualDir}
\begin{split}
\Delta \bar v^{i,(k+1)} &= \Delta \bar v^{i,(k)} + \left( \bar A^i\Delta w^{i,(k+1)} + r_{\text{primal,2}}^{i,(l)}\right) ,\\
\Delta \bar v_c^{i,(k+1)} &= \Delta \bar v_c^{i,(k)} + \left( \Delta w^{i,(k+1)} - \Delta x_{J_i}^{(k+1)} + r_c^{i,(l)}\right),
\end{split}
\end{align}
%
for $i = 1, \dots, N$. The dual variable directions \eqref{dualDir} are scaled, \cite{boyd:11}, and they have to be rescaled to give the actual dual variable directions, that is
%
\begin{align*}
\Delta \bar v= (1/\rho)\Delta v,\quad \Delta \bar v_c= (1/\rho)\Delta v_c.
\end{align*}
%
Having computed the directions $\Delta W$, $\Delta x$, $\Delta v_c$ and $\Delta v$, we can now compute $\Delta s$ and $\Delta \lambda$ as
\begin{align}
\label{extraDir}
\begin{split}
\Delta s^i\! &=\! -D\bar{G}^i(w^{(l)})\Delta w^i\!- \!  r^{i,(l)}_{\text{primal,1}},\\
\Delta \lambda^i\! &=\! -(S^{i,(l)})^{-1}\!\left(\! \Lambda^{i,(l)}\Delta s^i\! - \!  r^{i,(l)}_{\text{cent}}+\mu^{(l)}\hat{e} \right)\!,
\end{split}
\end{align}
%
for $i = 1, \dots, N$.
%
The distributed algorithm for computing the primal-dual directions is expressed in Algorithm \ref{alg:alg4PD}.
%
\subsubsection*{Remark 3}
\emph{The computational effort for each iteration of Algorithm \ref{alg:alg4PD} is dominated by the cost of updating the iterates $\Delta w^{i,(k+1)}$, which requires factorizing the matrices $H^{i,(l)}_{\text{pd}} + \rho\left( I + (\bar A^i)^T\bar A^i\right)$ for $i = 1, \dots, N.$ In case $\rho$ is chosen to be a constant, these matrices remain the same within each iteration of the algorithm, and hence the computational burden of each instance of Algorithm \ref{alg:alg4PD} can be significantly reduced by pre-caching the factorizations and reusing them in the subsequent iterations. In fact, even if $\rho$ is nonconstant we can adopt a procedure that would allow us to update the factorizations of these matrices without having to recompute them entirely \emph{\cite[Sec. 4.2]{liu:13}}.}
%
%
\subsubsection*{Remark 4}
\emph{We can use other proximal splitting methods than ADMM for solving \eqref{eq:QAppEqLogPD} in a distributed way and possibly get better convergence properties. However, other proximal splitting methods generally require a reformulation of \eqref{eq:QAppEqLogPD}, which can in turn complicate the recovery of the dual variable directions, $\Delta v$ and $\Delta v_c$. In order to keep the presentation simple, we have restricted ourselves to using ADMM in this paper.}

\begin{algorithm}[H]
\caption{ADMM-Based Primal-Dual Direction Computation}\label{alg:alg4PD}
\begin{algorithmic}[1]
\State{Given $k= 0$, $\rho > 0$, $\hat{\eta}>0$, $\epsilon_{\textrm{pri}}>0$, $\epsilon_{\textrm{dual}}>0$, $s^{i}$, $\lambda^{i}$, $m_i$ for $i=1,\dots,N$, $\Delta W^{(0)}$, $\Delta \bar v^{(0)}$ and $\Delta \bar v_c^{(0)}$.}
\For{$i = 1, 2, \dots,N$}
\State {Communicate with all agents $r$ belonging to $\text{Ne}(i)$.}
\For {all $j \in J_i$}
 \State $\Delta x_j^{(0)} = \frac{1}{|\mathcal{I}_j|}\sum_{q \in \mathcal{I}_j}^{} \left( E_{J_q}^T \Delta w^{q,(0)} \right)_j$.
 \EndFor
 \EndFor
 \Repeat
\For{$i = 1, 2, \dots,N$}
\State {\vspace*{-10pt}
\begin{multline*}
\hspace{6mm}\Delta w^{i,(k+1)} = -\left[  H^{i,(l)}_{\text{pd}} + \rho\left( I + (\bar A^i)^T\bar A^i\right)  \right]^{-1}\times\\  \left[ r^{i,(l)} + \rho \left( r_c^{i,(l)} + \Delta \bar v_c^{i,(k)} -   \Delta x_{J_i}^{(k)} \right)+ \right.  \\ \left.  \rho (\bar A^i)^T \left( r_{\text{primal,2}}^{i,(l)} + \Delta \bar v^{i,(k)}  \right)  \right].
\end{multline*}}
\State {Communicate with all agents $r$ belonging to $\text{Ne}(i)$.}
\For {all $j \in J_i$}
 \State $\Delta x_j^{(k+1)} = \frac{1}{|\mathcal{I}_j|}\sum_{q \in \mathcal{I}_j}^{} \left[ E_{J_q}^T\left( \Delta w^{q,(k+1)}+  \Delta \bar v_c^{q,(k)} + v_c^{q,(l)} + r_c^{q,(l)}   \right) \right]_j. $
 \EndFor
 \State $\Delta \bar v^{i,(k+1)} = \Delta \bar v^{i,(k)} + \left( \bar A^i\Delta w^{i,(k+1)} + r_p^{i,(l)}\right)$.
  \State $\Delta \bar v_c^{i,(k+1)} = \Delta \bar v_c^{i,(k)} + \left( \Delta w^{i,(k+1)} - \Delta x_{J_i}^{(k+1)} + r_c^{i,(l)}\right)$.
\State  Check whether $\| \Delta x_{J_i}^{(k+1)}-\Delta x_{J_i}^{(k)}\|^2 \leq \epsilon_{\text{dual}}/N$, $\| \Delta w^{i,(k+1)}-\Delta x_{J_i}^{(k+1)}+r_c^{i,(l)}\|^2 \leq \epsilon_{\text{pri}}/(2N)$ and $\| \bar A^i \Delta w^{i,(k+1)} + r_{\text{primal,2}}^{i,(l)}\|^2 \leq \epsilon_{\text{pri}}/(2N)$.
\EndFor
\If  {Condition in Step (17) satisfied for all $i = 1, \dots, N$.}
\State {Terminate the algorithm.}
\EndIf
\State{$k = k+1$.}
\Until{Algorithm is terminated.}
\For{$i = 1, 2, \dots,N$}
\State{\vspace*{-10pt}\begin{align*}
\Delta s^i &= -D\bar{G}^i(w^{(l)})\Delta w^i-   r^{i,(l)}_{\text{primal,1}},\\
\Delta \lambda^i &= -(S^{i,(l)})^{-1}\left( \Lambda^{i,(l)}\Delta s^i -   r^{i,(l)}_{\text{cent}}+\mu^{(l)}\hat{e} \right).
\end{align*}}
\EndFor
 \end{algorithmic}
 \end{algorithm}
%

%
%In order to provide a fully distributed algorithm, we next discuss distributed ways for computing the perturbation parameter and step size, and checking the algorithm.
%
%=================================================================
\subsection{Distributed Computations of Perturbation Parameter, Step Size and Stopping Criterion}\label{distStep}
%=================================================================
In the distributed case, we set $\mu$ at iteration $l$ to
\begin{align}\label{eq:Distmu}
\mu^{(l)} = \sigma \frac{\underset{i}{\textrm{min}}\big(\hat{\eta}^{i,(l)}\big)}{\sum_{i=1}^N m_i},
\end{align}
where $\hat{\eta}^{i,(l)}=(s^{i,(l)})^T\lambda^{i,(l)}$ and $\sigma \in (0,1)$ is a user-defined constant. We can distribute the required minimum value computations using algorithms such as \emph{min-consensus}, \cite{Iut:12}. To calculate the step size in a distributed way, we need to provide an alternative representation of $H(z^{(l)})$. For this purpose, we express $\|r_{\textrm{dual}}(z^{(l)})\|^2=\sum_{i=1}^N\|r_{\textrm{dual}}^{i}(z^{i,(l)})\|^2$ with
\begin{align}\label{rdual}
%\begin{split}
%r_{\textrm{dual}}^{i}(z^{i,(l)})=
%\begin{bmatrix}\begin{array}{c}0\\\vdots\\0\\
%\nabla \bar{f}_i(w^{i,(l)})+\sum_{j=1}^{m_i}\nabla \bar{G}_j^i(w^{i,(l)})\lambda_j^i+(\bar{A}^i)^Tv^{i,(l)}+v_c^{i,(l)}\\
%0\\\vdots\\0\\
%0^{(*)}\\
%\end{array}\end{bmatrix},
r_{\textrm{dual}}^{i}(z^{i,(l)})=\nabla \bar{f}_i(w^{i,(l)})+\sum_{j=1}^{m_i}\nabla \bar{G}_j^i(w^{i,(l)})\lambda_j^i+(\bar{A}^i)^Tv^{i,(l)}+v_c^{i,(l)}.
\end{align}
%
 Note that the sum $-\sum_{i=1}^NE_{J_i}^Tv_c^{i,(l)}$ is not included in the expression for the residual. This is because it is always zero for each iteration of the algorithm provided that we initialize $v_c^{i,(0)}$ such that $-\sum_{i=1}^NE_{J_i}^Tv_c^{i,(0)}=0$ and choose the same step size for all subsystems. To see this, consult the derivation of why the second block in residuals vector \eqref{residual2} is always equal to zero in the proof of Theorem \ref{connectStop}. We express the primal residuals as before,
\begin{align}\label{rprimal}
\begin{split}
(r_{\textrm{primal,1}}^{i}(z^{i,(l)}),r_{\textrm{primal,2}}^{i}(z^{i,(l)}))=
\begin{bmatrix}\begin{array}{c}
\bar{G}^i(w^{i,(l)})-s^i\\
\hdashline\\
\bar{A}^iw^{i,(l)}-b^i\\
w^{i,(l)}-x_{J_i}^{(l)}\\
\end{array}\end{bmatrix}.
\end{split}
\end{align}
We are now ready to express $\|H(z^{(l)})\|^2=\sum_{i=1}^N\|H^{i}(z^{i,(l)})\|^2$ with
\begin{align}\label{H}
\begin{split}
H^{i}(z^{i,(l)})=
\begin{bmatrix}\begin{array}{c}
r_{\textrm{dual}}^{i}(z^{i,(l)})\\
(r_{\textrm{primal,1}}^{i}(z^{i,(l)}),r_{\textrm{primal,2}}^{i}(z^{i,(l)}))\\
S^{i,(l)}\Lambda^{i,(l)}\hat{e}\\
\end{array}\end{bmatrix},
\end{split}
\end{align}
and are thus able to distribute the evaluation of $\|H(z^{(l)})\|^2$. At this point, let each agent $i$ compute its local step size $\alpha^{i,(l+1)}$ using the approach described in Section \ref{sec:PDIPMstep}, that is the agent first sets
%
\begin{align*}
\alpha_{\textrm{max}}^i = \textrm{min} \left\{ 1, \underset{j}{\textrm{min}} \left\{ -\lambda_j^{i,(l)}/\Delta \lambda_j^{i,(l+1)} \ \big | \ \Delta \lambda_j^{i,(l+1)}<0  \right\} \right\},
\end{align*}
%
for $j=1, \dots, m_i$, and perform a backtracking line search as\\
%
\begin{algorithmic}
  \While{$\exists \ j  \colon s_j^{i,(l)} + \alpha^{i,(l+1)} \Delta s_j^{i,(l+1)} \leq 0 $}
    \State $\alpha^{i,(l+1)} = \beta \alpha^{i,(l+1)}$
  \EndWhile
\end{algorithmic}
%
with $\beta \in (0,1)$ and $\alpha^{i,(l+1)}$ initialized as $0.99 \alpha_{\textrm{max}}^i$. In order to guarantee convergence of the KKT conditions to zero we continue the back tracking as
%
\begin{flushleft}
\begin{algorithmic}
 \While{$\left\| H^{i,(l)}(\alpha^{i,(l+1)}) \right\|^2> (1 - \gamma \alpha^{i,(l+1)})^2  \left\| H^{i,(l)}(0)\right\|^2   $}
    \State  $\alpha^{i,(l+1)} = \beta \alpha^{i,(l+1)}$
  \EndWhile
\end{algorithmic}
\end{flushleft}
%
\noindent where $\gamma \in [0.01, 0.1]$. The resulting $\alpha^{i,(l+1)}$ ensures that the local iterates remain feasible with respect to local inequality constraints and that the norm of the local KKT residuals is decreased consistently after each iteration. Once all agents have computed their local step sizes, we then choose the global step size as the smallest one, that is
\begin{align*}
\alpha^{(l+1)} = \underset{i}{\textrm{min}}\left\{\alpha^{i,(l+1)}\right\}, \textrm{ for } i = 1,\dots,N.
\end{align*}
This then allows us to guarantee the aforementioned properties consistently for all agents. As for the perturbation parameter, the minimum value computations can be performed in a distributed fashion using algorithms such as \emph{min-consensus}, \cite{Iut:12}.

It is also possible to check the stopping criterion at each iteration in a distributed way. For this purpose we need to distribute the check of whether the primal and dual residual norms together with the centrality residual are small enough. Note that, due to~\eqref{rdual} and~\eqref{rprimal}, we have the following implications:
\begin{align*}
\|(r_{\textrm{primal,1}}^{i,(l)},r_{\textrm{primal,2}}^{i,(l)})\|^2\leq\frac{\epsilon_{\textrm{feas}}^2}{N} \textrm{ for all } i=1,\dots,N \Rightarrow &\|(r_{\textrm{primal,1}}^{(l)}r_{\textrm{primal,2}}^{(l)})\|\!\leq\! \epsilon_{\textrm{feas}},\\
\|r_{\textrm{dual}}^{i,(l)}\|^2\leq\frac{\epsilon_{\textrm{feas}}^2}{N} \textrm{ for all } i=1,\dots,N \Rightarrow& \|r_{\textrm{dual}}^{(l)}\|\leq \epsilon_{\textrm{feas}}
\end{align*}
and
\begin{align*}
\hat{\eta}^{i,(l)}=(s^{i,(l)})^T\lambda^{i,(l)}\leq\frac{\epsilon}{N}\textrm{ for all } i=1,\dots,N \Rightarrow \hat{\eta}^{(l)}\leq\epsilon.
\end{align*}
Hence, in case for all agents $i = 1, \dots, N$ we have
\begin{align*}
\|(r_{\textrm{primal,1}}^{i,(l)},r_{\textrm{primal,2}}^{i,(l)})\|^2&\leq\frac{\epsilon_{\textrm{feas}}^2}{N},\\
\|r_{\textrm{dual}}^{i,(l)}\|^2&\leq\frac{\epsilon_{\textrm{feas}}^2}{N},\\
\hat{\eta}^{i,(l)}=(s^{i,(l)})^T\lambda^{i,(l)}&\leq\frac{\epsilon}{N},
\end{align*}
then the termination condition of the primal-dual algorithm is satisfied and the algorithm is terminated. Notice that this can be performed in a distributed way and to this end, each agent would merely need to declare whether its local termination conditions are satisfied or not.

At this point we can distribute the computations of the primal-dual directions $\Delta W$, $\Delta x$, $\Delta v$, $\Delta v_c$, $\Delta s$ and $\Delta \lambda$ using Algorithm \ref{alg:alg4PD}. The algorithm can be used in Step 4 of Algorithm~\ref{alg:alg1a}. Combining algorithms~\ref{alg:alg1a} and~\ref{alg:alg4PD} with the modifications discussed above result in Algorithm~\ref{alg:algNPDexact}, which is a distributed primal-dual interior-point method for solving \eqref{eq:DDEOP}.

Notice that for Algorithm \ref{alg:algNPDexact} to function consistently, it is important that the primal-dual directions are computed accurately. This in turn can require many ADMM iterations, particularly for the first iterations of the primal-dual method. It is, however, expected that the number of ADMM iterations would decrease as we progress through the primal-dual iterations, thanks to warm-starting of the ADMM iterations as outlined in Step 10 of Algorithm \ref{alg:algNPDexact}. One way of improving the efficiency of our proposed algorithm is to incorporate the use of inexact directions. This can potentially reduce the number of required ADMM iterations for computing the primal-dual directions. Analysis of interior-point methods under inexact directions has been investigated, and such methods are referred to as inexact interior-point methods. Studying the convergence results for such methods suggests that during the first iterations of an interior-point method, when we are far away from the optimal solution, it is not necessary to compute the search directions accurately, and the accuracy requirements become more stringent as we progress through the interior-point method iterations, e.g., see \cite{han:00}, \cite{al:09}, \cite{bel:98}. This means that such methods utilize an adaptive stopping criterion for the search direction calculations, which yields higher and higher accuracy as the interior-point iterates get closer and closer to a solution.
Next we investigate the possibility of devising a distributed version of one such method.

%Note that the required minimum value computations in steps 9 and~14 of Algorithm \ref{alg:algNPDexact} can be performed distributedly using algorithms such as \emph{min-consensus}, \cite{Iut:12}.
\begin{algorithm}[tb]
\caption{Distributed Primal-Dual Interior-Point Method}\label{alg:algNPDexact}
\begin{algorithmic}[1]
\State
%
Given $l=0$, $\sigma \in (0,1)$, $v_c^{(0)}$ such that $\bar{E}^Tv_c^{(0)}=0$, $W^{0}$, $(s^{(0)},\lambda^{(0)})\succ 0$, $\epsilon_{\textrm{feas}}>0$, $\epsilon>0$, $\epsilon_{\textrm{pri}}>0$, $\epsilon_{\textrm{dual}}>0$, $m_i$ for $i=1,\dots,N$, $\Delta W^{(0)}$, $\Delta \bar v^{(0)}$ and $\Delta \bar v_c^{(0)}$.
\For{$i = 1, 2, \dots,N$}
\State {Communicate with all agents $r$ belonging to $\text{Ne}(i)$.}
\For {all $j \in J_i$}
 \State $x_j^{(0)} = \frac{1}{|\mathcal{I}_j|}\sum_{q \in \mathcal{I}_j}^{} \left( E_{J_q}^T  w^{q,(0)} \right)_j.$
 \EndFor
 \EndFor
%
\Repeat
%
%
\State{Compute $\mu^{(l)}$ as in \eqref{eq:Distmu}.}
\State{Given $\mu^{(l)}$ and $z^{(l)}$ compute $\Delta z^{(l+1)}$ using Alg. \ref{alg:alg4PD} with the initial iterates $\Delta z^{(l)}$. }
%
\For{$i = 1, 2, \dots,N$}
\State{Compute local step size $\alpha^{i,(l+1)}$ using the approach presented in Section \ref{distStep}.}
\EndFor
\State{Let $\alpha^{l+1}=\underset{i}{\textrm{min}}\{\alpha^{i,(l+1)}\}$.}
\For{$i = 1, 2, \dots,N$}
\State{Set $z^{i,(l+1)}=z^{i,(l)}+\alpha^{(l+1)}\Delta_z^{i,(l+1)}$.}
\EndFor
\State{Set $l=l+1$.}
%
\Until{$\|(r_{\textrm{primal,1}}^{i}(z^{i,(l)}),r_{\textrm{primal,2}}^{i}(z^{i,(l)}))\|^2\leq \epsilon_{\textrm{feas}}^2/N$, $\|(r_{\textrm{dual}}^{i}(z^{i,(l)})))\|^2\leq \epsilon_{\textrm{feas}}^2/N$ and $\hat{\eta}^{i,(l)}\leq \epsilon/N$ for all $i=1,\dots,N$.}
\end{algorithmic}
\end{algorithm}
%
\vspace{32pt}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Primal-Dual Inexact Interior-Point Methods}\label{sec:PDIIPM}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%A common claim [reference] is that it is unnecessary to require high accuracy of the step directions when the interior-point iterates have low accuracy. Meaning, we want an adaptive stopping criteria for the step direction calculations, which yields higher and higher accuracy as the interior-point iterates get closer and closer to a solution. Such an implementation of a primal-dual interior point method is called a primal-dual \emph{inexact} interior-point method. The framework presented here is laid out in Algorithm \ref{alg:alg1b} and was introduced in \cite{bel:98}.
%------------------------------------------------------------------------------------------------------------
\begin{algorithm}[tb]
\caption{Primal-dual Inexact Interior-point Method, \cite[]{bel:98}.}\label{alg:alg1b}
\begin{algorithmic}[1]
\State
%
Given $(s^{(0)},\lambda^{(0)})\succ 0$, $\bar{\tau}_1=\textrm{min}(\Lambda^{(0)}S^{(0)}\hat{e})/((\lambda^{(0)})^Ts^{(0)}/m),$ $\bar{\tau}_2=(\lambda^{(0)})^Ts^{(0)}/\|R(z^{(0)})\|,$
$\gamma^{(l-1)}\in[1/2,1)$, $\eta_{\textrm{max}}\in(0,1)$, $\beta\in(0,1)$, $\theta\in(0,1)$ and $\epsilon>0$.
%
\Repeat
%
		\State{
%
Choose $\sigma^{(l)}$, $\hat{\eta}^{(l)}$ and $\gamma^{(l)}\in[1/2,\gamma^{(l-1)}]$ such that $(\sigma^{(l)}+\hat{\eta}^{(l)})\in(0,\eta_{\textrm{max}})$ and
\begin{align*}
\sigma^{(l)}>\textrm{max}\left( \frac{\sqrt{m}+\bar{\tau}_1\gamma^{(l)}}{\sqrt{m}(1-\bar{\tau}_1\gamma^{(l)})},\frac{\sqrt{m}+\bar{\tau}_2\gamma^{(l)}}{m} \right)\hat{\eta}^{(l)}.
\end{align*}
Put $\mu^{(l)}=\sigma^{(l)}(s^{(l)})^T\lambda^{(l)}/m$ and $\bar{\eta}^{(l)}=\sigma^{(l)}+\hat{\eta}^{(l)}$.
%
		\State
%
Compute $\Delta z^{(l)}$ by solving \eqref{eq:QAppIneqKKTPDComInexact} with $\|\hat{r}^{(l)}\|\leq \hat{\eta}^{(l)}(s^{(l)})^T\lambda^{(l)}/m.$
%
		\State
%
Choose $\bar{\alpha}_1^{(l)}$ such that
$$\textrm{min}\left(S^{(l)}(\bar{\alpha}_1^{(l)})\Lambda^{(l)}(\bar{\alpha}_2^{(l)})\hat{e}\right)\geq\bar{\tau}_1\gamma^{(l)}\left(s^{(l)}(\bar{\alpha}_1^{(l)})\right)^T\lambda^{(l)}(\bar{\alpha}_1^{(l)})/m.$$}
\State{Choose $\bar{\alpha}_2^{(l)}$ such that
$$\left(s^{(l)}(\bar{\alpha}_2^{(l)})\right)^T\lambda^{(l)}(\bar{\alpha}_2^{(l)})\geq\bar{\tau}_2\gamma^{(l)}\|R(z^{(l)}(\bar{\alpha}_2^{(l)}))\|.$$
%
		\State
%
Set $\alpha^{(l)}=\textrm{min}(\bar{\alpha}_1^{(l)},\bar{\alpha}_2^{(l)}).$
%
		\State
%
Set $\eta^{(l)}=1-\alpha^{(l)}(1-\bar{\eta}^{(l)})$.
%
			\While{
%
$\| H^{(l)}(\alpha^{(l)}) \|> (1-\beta(1-\eta^{(l)}))\| H^{(l)}(0) \|$}
%
				\State{
%
$\alpha^{(l)}=\theta \alpha^{(l)}$ and $\eta^{(l)}=1-\theta(1-\eta^{(l)})$.}
%
			\EndWhile
%
		\State
%
Set $z^{(l+1)}=z^{(l)}+\alpha^{(l)}\Delta z^{(l)}$ and $l=l+1$.
%
	}
%
\Until{$\|H(z^{(l)})\|\leq \epsilon$.}
\end{algorithmic}
\end{algorithm}
%------------------------------------------------------------------------------------------------------------
In an inexact interior-point method we only need to solve \eqref{eq:QAppIneqKKTPDCom} approximately for the primal-dual directions. That is, we solve
\vspace{-5pt}
\begin{align}\label{eq:QAppIneqKKTPDComInexact}
H'(z^{(l)})\Delta z =- H(z^{(l)})+\mu^{(l)}\begin{bmatrix}0 \\ 0\\ 0  \\\hat{e} \end{bmatrix}+\hat{r}^{(l)},
\end{align}
where $\hat{r}^{(l)}$ is the residual. However, in order to assure the convergence of the algorithm, it is necessary to make modifications to the framework. Specifically, consider the framework laid out in Algorithm \ref{alg:alg1b} as introduced in \cite{bel:98}. In order to ensure a global convergent inexact interior-point method, \cite{bel:98}, \cite{bakry:96}, it is necessary to
\begin{itemize}
\item impose restrictions on the problem formulation;
\item add requirements on the residual $\hat{r}$ in \eqref{eq:QAppIneqKKTPDComInexact};
\item change the choice of $\hat{\eta}$, $\sigma$ and $\alpha$.
\end{itemize}

For example, the total residual $\tilde{r}^{(l)}$ of the KKT system must fulfill
\begin{align}
\|\tilde{r}^{(l)}\|=\left\|\mu^{(l)}\begin{bmatrix}0 \\ 0\\ 0  \\\hat{e} \end{bmatrix}+\hat{r}^{(l)}\right\|\leq(\sigma^{(l)}+\hat{\eta}^{(l)})\|H(z^{(l)})\|,
\end{align}
see \cite{bel:98}.

We define the set $\Omega(\epsilon)$ for a given $\epsilon>0$ as
\begin{align}\label{omega:set}
\begin{split}
\Omega(\epsilon)=\bigg\{&z\in \mathbb{R}^{n+2m+p}|\epsilon\leq\|H(z)\|\leq\|H(z^{(0)})\|,\\ &\textrm{min}\left(S\Lambda\hat{e}\right)\geq \frac{\bar{\tau}_1}{m}\frac{1}{2}s^T\lambda,\ s^T\lambda\geq\bar{\tau}_2\frac{1}{2}\|R(z)\|\bigg\},
\end{split}
\end{align}
with $\bar{\tau}_1=\textrm{min}(\Lambda^{(0)}S^{(0)}\hat{e})/((\lambda^{(0)})^Ts^{(0)}/m)$ and $\bar{\tau}_2=(\lambda^{(0)})^Ts^{(0)}/\|R(z^{(0)})\|,$
and the following assumptions
\begin{enumerate}
\item [A1] $H$ is continuously differentiable in $\Omega(0)$.
\item [A2] $\{z^{(l)}\}$ is bounded.
\item [A3] $H'(z)$ is nonsingular in $\Omega({\epsilon})$ with $\epsilon>0$.
\item [A4] $R'$ is Lipschitz continuous in $\Omega(0)$ with constant $L$, where $$R=(r_{\textrm{dual}},r_{\textrm{primal,1}},r_{\textrm{primal,2}}).$$
\end{enumerate}
Then if $\{z^{(l)}\}$ is generated by Algorithm \ref{alg:alg1b} and assumptions A1--A4 are fulfilled, the sequence $\{\|H(z^{(l)})\|\}$ will converge to zero and $z^{(l)}$ will converge to the limit point of $\{z^{(l)}\}$, see Theorem 3.3 in \cite{bel:98} (with the additional assumption that $\sigma^{(l)}$ is bounded away from zero). Similar to Algorithm~\ref{alg:algNPDexact}, it is also possible to distribute the computations in Algorithm \ref{alg:alg1b}, and that is discussed in the next section.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{A Distributed Primal-Dual Inexact Interior-Point Method for Solving Loosely Coupled Problems}\label{sec:IDP}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Let us apply the primal-dual inexact interior-point method described in Section \ref{sec:PDIIPM} to the problem in \eqref{eq:DDEOP}. Similar to Algorithm \ref{alg:algNPDexact}, we can use Algorithm \ref{alg:alg4PD} for computing the inexact directions in a distributed fashion. Particularly, this algorithm can be used in Step~4 of Algorithm~\ref{alg:alg1b}. However, in order for the computed directions to satisfy the required accuracy in Step~4 of Algorithm \ref{alg:alg1b}, we need to establish a connection between the ADMM stopping criteria and the norm of the residuals in \eqref{eq:QAppIneqKKTPDComInexact}. This connection is established in the following theorem.
%
\begin{theorem}\label{connectStop}
It is possible to choose the thresholds $\epsilon_{\text{pri}}, \epsilon_{\text{dual}}>0$ such that the stopping criteria in Algorithm \ref{alg:alg4PD} and the residual conditions in Step 4 of Algorithm \ref{alg:alg1b} are equivalent.
\end{theorem}
%
\begin{proof}
Note that in our approach for solving \eqref{eq:QAppIneqKKTPDCom}, we in fact solve \eqref{eq:QAppIneqKKTPD1} and~\eqref{eq:QAppIneqKKTPD2} exactly, since $\Delta s$ and $\Delta \lambda$ are eliminated, see \eqref{eq:PDsolution}. The residuals in~\eqref{eq:QAppIneqKKTPDComInexact} for our approach and for the problem in \eqref{eq:DDEOP} are therefore given as
\begin{align}\label{residual2}\hat{r}^{(l)}\!\!=\!\!
\begin{bmatrix}\begin{array}{c}H^{1,(l)}_{pd}\Delta w^1+(\bar{A}^1)^T\Delta v^1+\Delta v_c^1+r^{1(l)} \\ \vdots\\ H^{N,(l)}_{pd}\Delta w^N+(\bar{A}^N)^T\Delta v^N+\Delta v_c^N +r^{N,(l)} \\\hdashline
\sum_{i=1}^N(-E_{J_i}^T\Delta v_c-E_{J_i}^Tv_c^{(l)})\\
\hdashline \\ 0\\\vdots\\0\\\hdashline \bar{A}^1\Delta w^1+r_{\textrm{primal,2}}^{1,(l)}\\
 \vdots\\
 \bar{A}^N\Delta w^N+r_{\textrm{primal,2}}^{N,(l)}\\
\\ \Delta w^1-\Delta x_{J_1}+r_c^{1,(l)}\\ \vdots\\
\Delta w^N-\Delta x_{J_N}+r_c^{N,(l)}\\\hdashline \\ 0\\\vdots\\0\end{array}\end{bmatrix}\!\!.
\end{align}
%
The norm of the fourth block of the right hand side of \eqref{residual2} is already included in the stopping criteria in Algorithm \ref{alg:alg4PD}. Furthermore, in the ADMM iterations, $\Delta w^{(k+1)}$ and $\Delta v^{(k+1)}$ are computed such that
\begin{align*}
0=&\left( H^{i,(l)}_{pd} + \rho I + \rho (\bar A^i)^T\bar A^i\right)\Delta w^{i,(k+1)}+ \rho \left( r_c^{i,(l)} + \frac{1}{\rho}\Delta v_c^{i,(k)} - \Delta x_{J_i}^{(k)} \right)+\\&r^{i,(l)}+    \rho (\bar A^i)^T \left( r_{\text{primal,2}}^{i,(l)} + \frac{1}{\rho}\Delta v^{i,(k)}  \right)\\
=&\left( H^{i,(l)}_{pd}\! + \rho I + \rho (\bar A^i)^T\bar A^i\right)\Delta w^{i,(k+1)}\!+\!\rho \left( -\Delta w^{i,(k+1)}\!+\Delta x_{J_i}^{(k+1)}\!-\Delta x_{J_i}^{(k)} \right)
\!+\\& r^{i,(l)} +
\rho (\bar A^i)^T \left( \frac{1}{\rho}\Delta v^{i,(k+1)} -\bar{A}^i\Delta w^{i,(k+1)} \right),
\end{align*}
which gives
\begin{align*}
H^{i,(l)}_{pd}\Delta w^{i,(k+1)}+r^{i,(l)}+(\bar A^i)^T \Delta v^{i,(k+1)}+\Delta v_c^{i,(k+1)} =
\rho(\Delta x_{J_i}^{(k)}-\Delta x_{J_i}^{(k+1)}).
\end{align*}
Consequently, the norm of the first block of the right hand side of \eqref{residual2} would be small if and only if $\|\rho(\Delta x_{J_i}^{(k)}-\Delta x_{J_i}^{(k+1)})\|^2$  would be small, which is also included in the stopping criteria of Algorithm \ref{alg:alg4PD}. The only remaining part of the residuals vector in \eqref{residual2} is the second block which is always equal to zero. This is because, in the ADMM iterations, $\Delta x^{(k+1)}$ and $\Delta v_c^{(k+1)}$ are chosen such that
\begin{align*}
0\!=\!-\bar{E}^T\!v_c^{(l)}\!-\!\rho\bar{E}^T\!\!\left(\!\!\Delta w^{(k+1)} \!-\!\bar{E}\Delta x^{(k+1)}\!+r_c^{(l)}\!+\frac{1}{\rho}\Delta v_c^{(k)}\!\!\right)\!=\!-\bar{E}^T\!v_c^{(l)}\!-\!\bar{E}^T\!\Delta v_c^{(k+1)}.
\end{align*}
As a result, we have
%
\begin{multline}
\label{residualhat}
\| \hat r^{(l)}  \|^2 = \sum_{i= 1}^N \left(\|\rho(\Delta x_{J_i}^{(k)}-\Delta x_{J_i}^{(k+1)})\|^2 +\right. \\ \left. \| \Delta w^{i,(k+1)}-\Delta x_{J_i}^{(k+1)}+r_c^{i,(l)}\|^2  +\| \bar A^i \Delta w^{i,(k+1)} + r_{\text{primal,2}}^{i,(l)}\|^2\right)
\end{multline}
and hence if the thresholds $\epsilon_{\text{pri}}, \epsilon_{\text{dual}}>0$ are chosen appropriately, then the stopping criteria in Algorithm \ref{alg:alg4PD} and the residual conditions in Step 4 of Algorithm \ref{alg:alg1b} will be equivalent.
%\begin{subequations}
%\label{conditionADMM}
%\begin{align}
%\label{conditionADMM1}
%\bar{A}_i\Delta w^i+r_{\textrm{primal,2}}^{i,(l)}&=0,\ i=1,\dots,N,\\
%\label{conditionADMM2}
%\Delta w^i-\Delta x_{J_i}+r_c^{i,(l)} &=0,\ i=1,\dots,N,\\
%\label{conditionADMM3}
%H_{\textrm{pd}}^{i,(l)}\Delta w^i+r^{i,(l)}+(\bar{A}^i)^T\Delta v^i+\Delta v_c^i&=0,\ i=1,\dots,N,\\
%\label{conditionADMM4}
%\sum_{i=1}^N (-E^T_{J_i}v_c^{i,(l)}-E^T_{J_i}\Delta v_c^i)&=0.
%\end{align}
%\end{subequations}
%Conditions \eqref{conditionADMM1} and \eqref{conditionADMM2} are considered fulfilled if the squared norm of their left hand side is smaller or equal to $\epsilon_{\text{pri}}/2N$ for some $\epsilon_{\text{pri}}>0$. See the second and third check of Step 17 in Algorithm \eqref{alg:alg4PD}. The left hand side of \eqref{conditionADMM3} is equivalent to $\rho(\Delta x_{J_i}^{(k)}-\Delta x_{J_i}^{(k+1)})$ since $\Delta w^{(k+1)}$ and $\Delta v^{k+1}$ are chosen in ADMM such that

%Consequently, \eqref{conditionADMM3} is considered fulfilled if $\|\Delta x_{J_i}^{(k+1)}-\Delta x_{J_i}^{(k)}\|^2\leq\epsilon_{\text{dual}}/N$ for some $\epsilon_{\text{dual}}>0$. See the first check of Step 17 in Algorithm \eqref{alg:alg4PD}. Condition \eqref{conditionADMM4} is always fulfilled since $\Delta x^{(k+1)}$ and $\Delta v_c^{k+1}$ are chosen in ADMM such that
%\begin{align*}
%0&=\!-\bar{E}^Tv_c^{(l)}\!-\rho\bar{E}^T\!\!\left(\!\!\Delta w^{(k+1)}-\bar{E}\Delta x^{(k+1)}+r_c^{(l)}+\frac{1}{\rho}\Delta v_c^{(k)}\!\!\right)\\
%&=-\bar{E}^Tv_c^{(l)}-\bar{E}^T\Delta v_c^{(k+1)}.
%\end{align*}
%Consequently, we do not have to check if this condition is fulfilled.
\end{proof}
Up to this point, we have illustrated that it is possible to use Algorithm \ref{alg:alg4PD} to distribute the computations of the inexact directions. Next we show how to choose the thresholds in the stopping criteria of Algorithm \ref{alg:alg4PD} so that the computed directions satisfy the necessary accuracy requirements. Moreover,  we describe how to distribute the update of the remaining iteration-dependent parameters and the remaining steps of Algorithm~\ref{alg:alg1b}.

\subsection{Distributed Computations of Perturbation Parameter, Step Size and Stop Criterion}\label{sec:diststep}
%
%Before we demonstrate how to calculate the perturbation parameter $\mu^{i,(l)}$ in the distributed case we

Let us first define $R^{i}(z^{i,(l)})$, $\sigma^{i,(l)}$ and $\hat{\eta}^{i,(l)}$ for each agent $i=1,\dots,N$. Particularly, we define $R^{i}(z^{i,(l)})$ such that $\|R(z^{(l)})\|^2=\sum_{i=1}^N\|R^{i}(z^{i,(l)})\|^2$ with
\begin{align}\label{R}
R^{i}(z^{i,(l)})=(r_{\textrm{dual}}^i(z^{i,(l)}),r_{\textrm{primal,1}}^i(z^{i,(l)}),r_{\textrm{primal,2}}^i(z^{i,(l)})),
\end{align}
%
and we choose $\sigma^{i,(l)}$ such that
\begin{align*}
\sigma^{i,(l)}>\bar{\tau}_2^i\gamma^{i,(l)}\hat{\eta}^{i,(l)}(s^{i,(l)})^T\lambda^{i,(l)}/\underset{i}{\textrm{min}}\left((s^{i,(l)})^T\lambda^{i,(l)}\right) + \epsilon_{\sigma},
\end{align*}
%
with $\bar{\tau}_2^i=(\lambda^{i,(0)})^Ts^{i,(0)}/\|R^i(z^{i,(0)})\|$, $\epsilon_{\sigma}\in (0,1)$ and $\hat{\eta}^{i,(l)}$ such that $(\sigma^{i,(l)}+\hat{\eta}^{i,(l)})\in (0,\eta_{\textrm{max}})$, $\gamma^{i,(l)}\in[1/2,\gamma^{i,(l-1)}]$ and $\gamma^{i,(0)}\in[1/2,1)$. Then choose $\hat{\eta}^{(l)} = \min_i\left\{ \hat{\eta}^{i,(l)} \right\}$, $\sigma^{(l)} = \max_i\left\{ \sigma^{i,(l)} \right\}$ and set $\bar{\eta}^{l}=\sigma^{(l)}+\hat{\eta}^{(l)}$. At this point, we can describe how to compute proper thresholds for the ADMM iterations termination criteria. In the distributed setting we set the residual norm condition in Step 4 in Algorithm~\ref{alg:alg1b} for the problem in~\eqref{eq:DDEOP} to
%
\begin{align}\label{stopIIP}
\|\hat{r}^{(l)}\|\leq \hat{\eta}^{(l)}\sum_{i=1}^N(s^{i,(l)})^T\lambda^{i,(l)}/m,
\end{align}
%
where $m = \sum_{i=1}^N m_i$ and the residual $\hat{r}^{(l)}$ is defined in \eqref{residual2}. There are several choices of $\epsilon_{\text{pri}}$ and $\epsilon_{\text{dual}}$ that ensure that criterion \eqref{stopIIP} is fulfilled when the stopping criteria of ADMM are satisfied. In this paper we set
\begin{align}\label{epsADMMInexact1}
\epsilon_{\text{pri}}^{i,(l)}=\frac{N}{2}\left(\hat{\eta}^{(l)}(s^{i,(l)})^T\lambda^{i,(l)}/m\right)^2
\end{align}
and
\begin{align}\label{epsADMMInexact2}
\epsilon_{\text{dual}}^{i,(l)}=\frac{N}{2}\left(\hat{\eta}^{(l)}(s^{i,(l)})^T\lambda^{i,(l)}/(\rho m)\right)^2.
\end{align}
Note that $\epsilon_{\text{pri}}$ and $\epsilon_{\text{dual}}$ are now subproblem-specific and they change with each iteration $l$ of the primal-dual inexact interior-point method.

Next we focus on computation of the perturbation parameter. In order to compute $\mu$ at iteration $l$, each agent first needs to compute
\begin{align*}
\mu^{i,(l)} = \sigma^{(l)}(s^{i,(l)})^T\lambda^{i,(l)}/m.
\end{align*}
Then the perturbation parameter is chosen as $\mu^{(l)} = \underset{i}{\textrm{min}} \left\{ \mu^{i,(l)} \right\}$.

It now remains to compute a proper step size for updating the iterates. As it is laid out in Algorithm \ref{alg:alg1b}, the process of computing the step size $\alpha^{(l)}$ at each iteration consists of two stages, namely, computation of an upper-bound on the step size and the line search. Similar to the approach we undertook for computing a step size in Section \ref{distStep}, each agent first needs to compute their local step size $\alpha^{i,(l)}$. To this end, each agent $i$ initially sets
%%
\begin{align*}
\alpha^{i,(l)}=\textrm{min}\left\{\bar{\alpha}_1^{i,(l)},\bar{\alpha}_2^{i,(l)}\right\},
\end{align*}
where
\begin{align*}
\bar \alpha_j^{i,(l)}=\underset{\alpha \in [0,1]}{\textrm{max}}\left\{\alpha|f_j(\alpha')\geq 0, \textrm{ for all }\alpha'\leq\alpha\right\},
\end{align*}
with $j \in \{1, 2 \}$, and
\begin{align*}f_1(\alpha)=\textrm{min}\left(S^{i,(l)}(\alpha)\Lambda^{i,(l)}(\alpha)\hat{e}\right)-\bar{\tau}^i_1\gamma^{i,(l)}\left(s^{i,(l)}(\alpha)\right)^T\lambda^{i,(l)}(\alpha)/m_i,\end{align*}
and
$$f_2(\alpha)=\left(s^{i,(l)}(\alpha)\right)^T\lambda^{i,(l)}(\alpha)-\bar{\tau}^i_2\gamma^{i,(l)}\|R^i(z^{i,(l)}(\alpha))\|.$$
Then the agent sets $\eta^{i,(l)}=1-\alpha^{i,(l)}(1-\bar{\eta}^{(l)})$ and performs a line search as
%
\begin{flushleft}
\begin{algorithmic}
 \While{$\left\| H^{i,(l)}(\alpha^{i,(l)}) \right\|> (1 - \beta ( 1 -  \eta^{i,(l+1)}))  \left\| H^{i,(l)}(0)\right\|   $}
    \State  $\alpha^{i,(l)} = \theta \alpha^{i,(l)}$
    \State  $\eta^{i,(l)}=1-\theta(1-\eta^{i,(l)})$
  \EndWhile
\end{algorithmic}
\end{flushleft}
%
When all agents are done computing their local step sizes, we then set $\alpha^{(l)}=\underset{i}{\textrm{min}}\{\alpha^{i,(l)}\}$. Finally, similar to the approach in Section \ref{distStep}, we check that $\| H^i(z^{i,(l)}) \|^2\leq \epsilon^2/N$ for $i=1,\dots,N$, which implies that $\| H(z^{(l)}) \|\leq \epsilon$, to decide whether to terminate the primal-dual iterations or not.

%=================================================================
\subsection{Distributed Primal-Dual Inexact Interior-Point\\ Method}
%=================================================================

Using Algorithm \ref{alg:alg4PD} with $\epsilon_{\text{pri}}^{i,(l)}$ and $\epsilon_{\text{dual}}^{i,(l)}$ set in accordance to \eqref{epsADMMInexact1} and \eqref{epsADMMInexact2}, respectively, we can distribute the computations of the primal-dual directions $\Delta W$, $\Delta x$, $\Delta v$, $\Delta v_c$, $\Delta s$ and $\Delta \lambda$. This algorithm can then be used in Step 4 of Algorithm~\ref{alg:alg1b} which distributes the major computations in the primal-dual inexact interior-point method. Combining algorithms \ref{alg:alg4PD} and~\ref{alg:alg1b} with the modifications discussed above result in Algorithm~\ref{alg:algNPD}, which is a distributed primal-dual inexact interior-point method for solving \eqref{eq:DDEOP}.

\begin{algorithm}[tb]
\caption{Distributed Primal-dual Inexact Interior-point Method}\label{alg:algNPD}
\begin{algorithmic}[1]
\State
%
Given $l=0$, $\rho>0$, $v_c^{(0)}$ such that $\bar{E}^Tv_c^{(0)}=0$, $W^{(0)}$, $(s^{(0)},\lambda^{(0)})\succ 0$,
 $\bar{\tau}_1^i=\textrm{min}(\Lambda^{i,(0)}S^{i,(0)}\hat{e})/((\lambda^{i,(0)})^Ts^{i,(0)}/m_i)$, $\bar{\tau}_2^i=(\lambda^{i,(0)})^Ts^{i,(0)}/\|R^i(z^{i,(0)})\|,$ $\gamma^{i,(0)}\in[1/2,1)$, $\eta_{\textrm{max}}^i\in(0,1)$, for $i=1,\dots,N$, $\epsilon_{\textrm{feas}}>0$, $\epsilon>0$, $\beta\in(0,1)$, $\theta\in(0,1)$, $\Delta W^{(0)}$,
$\Delta \bar v^{(0)}$ and $\Delta \bar v_c^{(0)}$.
%
\For{$i = 1, 2, \dots,N$}
\State {Communicate with all agents $r$ belonging to $\text{Ne}(i)$.}
\For {all $j \in J_i$}
 \State $x_j^{(0)} = \frac{1}{|\mathcal{I}_j|}\sum_{q \in \mathcal{I}_j}^{} \left( E_{J_q}^T  w^{q,(0)} \right)_j$.
 \EndFor
 \EndFor
\Repeat
%
%
\State Compute $\Delta z^{(l+1)}$ using Alg. \ref{alg:alg4PD} with the initial iterates $\Delta z^{(l)}$ with the stopping criteria thresholds given in \eqref{epsADMMInexact1} and \eqref{epsADMMInexact2}.
%\For{$i=1,\dots,N$}
%\State{Choose $\sigma^{i,(l)}$, $\hat{\eta}^{i,(l)}$ and $\gamma^{i,(l)}\in[1/2,\gamma^{i,(l-1)}]$ such that $(\sigma^{i,(l)}+\hat{\eta}^{i,(l)})\in(0,\eta_{\textrm{max}}^i)$ and $\sigma^{i,(l)}>\bar{\tau}_2^i\gamma^{i,(l)}\hat{\eta}^{i,(l)}/m_i.$ Put $\bar{\eta}^{i,(l)}=\sigma^{i,(l)}+\hat{\eta}^{i,(l)}$ and $\mu^{i,(l)}=\sigma^{i,(l)}(s^{i,(l)})^T\lambda^{i,(l)}/m_i$. Set $\epsilon_{\text{pri}}^{(l)}=N\left(\hat{\eta}^{i,(l)}(s^{i(l)})^T\lambda^{i(l)}/(\rho m)\right)^2$ and $\epsilon_{\text{dual}}^{(l)}=N\left(\hat{\eta}^{i,(l)}(s^{i(l)})^T\lambda^{i(l)}/m\right)^2$.}
%\State{Compute $\Delta z^{(l+1)}$ using Alg. \ref{alg:alg4PD} with the initial iterates $\Delta z^{(l)}$.}
%\EndFor
%
\For{$i=1,\dots,N$}
\State{Compute local step size $\alpha^{i,(l+1)}$ using the approach presented in Section \ref{distStep}.}
\EndFor
\State{Set $\alpha^{(l+1)}=\underset{i}{\textrm{min}}\{\alpha^{i,(l+1)}\}$.}
\State{Set $z^{(l+1)}=z^{(l)}+\alpha^{(l+1)}\Delta z^{(l+1)}$.}
\State{Set $l=l+1$.}
%
%\For{$i = 1, 2, \dots,N$}
%\State {Communicate with all agents $r$ belonging to $\text{Ne}(i)$}
%\For{all $j \in J_i$}
% \State $x_j^{(l)} = \frac{1}{|\mathcal{I}_j|}\sum_{q \in \mathcal{I}_j}^{} \left( E_{J_q}^T  s^{q,(l)} \right)_j$
% \EndFor
% \EndFor
\Until{$\| H^i(z^{i,(l)}) \|^2\leq \epsilon^2/N$ for all $i=1,\dots,N$.}
\end{algorithmic}
\end{algorithm}
%

%=================================================================
\subsection{Distributed Convergence Result}\label{distConRes}
%=================================================================
To ensure a global convergent inexact interior-point method, we modify the conditions stated in Section \ref{sec:PDIIPM}. In particular, we define the set $\Omega^i(\epsilon)$ for a given $\epsilon\geq0$ as
\begin{align}\label{omega:setinexact}
\begin{split}
\Omega^i(\epsilon)=\bigg\{&z^i\in \mathbb{R}^{|J_i|+2m_i+p_i}\big|\frac{\epsilon}{N}\leq\|H^i(z^i)\|\leq\|H^i(z^{i,(0)})\|,\\
&\textrm{min}\left(S^{i}\Lambda^{i}\hat{e}\right)\geq \frac{\bar{\tau}_1^i}{m_i}\frac{1}{2}(s^{i})^T\lambda^{i},\
(s^{i})^T\lambda^{i}\geq\bar{\tau}_2^i\frac{1}{2}\|R(z^{i})\|\bigg\},
\end{split}
\end{align}
with $$\bar{\tau}_1^i=\textrm{min}(\Lambda^{i,(0)}S^{i,(0)}\hat{e})/((\lambda^{i,(0)})^Ts^{i,(0)}/m_i),\ \bar{\tau}_2^i=(\lambda^{i,(0)})^Ts^{i,(0)}/\|R^i(z^{i,(0)})\|,$$ and the following assumptions
\begin{enumerate}
\item [B1] $H^i$ is continuously differentiable in $\Omega^i(0)$.
\item [B2] $\{z^{i,(l)}\}$ is bounded.
\item [B3] $H'$ is nonsingular in $\underset{i=1,\dots,N}{\prod}\Omega^i(\epsilon)$ with $\epsilon>0$.
\item [B4] $(R^i)'$ is Lipschitz continuous in $\Omega^i(0)$ with constant $L^i$, where $$R^i=(r_{\textrm{dual}}^i,r_{\textrm{primal,1}}^i,r_{\textrm{primal,2}}^i).$$
\end{enumerate}

If $\{z^{(l)}\}$ is generated by Algorithm \ref{alg:algNPD} and assumptions B1--B4 are fulfilled, then $\{\|H(z^{(l)})\|\}$ converges to zero and $z^{(l)}$ will converge to the limit point of $\{z^{(l)}\}$. For a proof of this see Appendix \ref{app:global}. Note that $\{z^{i,(l)}\}$ generated by Algorithm \ref{alg:algNPD} lies in $\Omega^i(0)$ for $i=1,\dots,N$ and all $l$.  Next we establish the connection of our proposed approach to that of iterative solvers and put forth suggestions on how to improve the convergence properties of the algorithm.
%=================================================================
\section{Iterative Solvers for Saddle Point Systems}\label{sec:saddle}
%=================================================================
The optimality conditions \eqref{sysOfEq} is a saddle point system where the solution, since strong duality holds, is a saddle point of the Lagrangian function of optimization problem \eqref{eq:ADMM}, see \cite{boyd:04}. In addition, a saddle point of the Lagrangian function is a saddle point of the augmented Lagrangian, and vice versa, \citep{gab:76}. Consequently, to find a solution of \eqref{sysOfEq}, we can instead consider the saddle point system corresponding to the augmented Lagrangian. The benefit of using the augmented Lagrangian is improved convergence properties when using dual methods for solving the saddle point system, \citep{arr:64,hes:69,for:83}.

%=================================================================
\subsection{Uzawa's Method and Fixed Point Iterations}\label{sec:UMAFPI}
%=================================================================
A well known algorithm for solving saddle point systems such as \eqref{sysOfEq} is Uzawa's method \citep{arr:64}. We solve the system of equations \eqref{sysOfEq} using ADMM (Algorithm \ref{alg:alg2}), which was originally derived as a modified version of Uzawa's method, see \citep{mar:75,gab:76}. ADMM applied to \eqref{eq:ADMM} is equivalent to Uzawa's method applied to the problem corresponding to the augmented Lagrangian of \eqref{eq:ADMM}, \citep{arr:64,hes:69,pow:69}, with one Gauss-Seidel iteration, \citep{saa:03}, in the update of the primal variables, \citep{mar:75,gab:76}. In addition, ADMM can be viewed as fixed point iterations of a pre-conditioned version of \eqref{sysOfEq}, see e.g. \cite{boy:14}. The similarities between ADMM, Uzawa's method and fixed point iterations are explored explicitly for our problem in Appendix \eqref{admom}.

Uzawa's method is also equivalent to the method of multipliers when applied to the augmented Lagrangian and the relaxation parameter in Uzawa's method is set to be equal to the penalty parameter in the method of multipliers, see \cite{ben:05}.
%
%=================================================================
\subsection{Other Iterative Methods}
%=================================================================
For a rigorous overview of the iterative methods available, we refer to \cite{saa:03,ben:05}. The problem in \eqref{sysOfEq} has an indefinite system matrix with an upper left block matrix that is singular, which limits the number of applicable methods or at least requires some pre-conditioning beforehand. For example, we could use the CG method applied to the normal equations of \eqref{sysOfEq}, see \cite{saa:03}. However, to the best of our knowledge, one has then destroyed the inherent structure of the problem which prevent us from distributing the calculations.




%=================================================================
\section{Improving Convergence Rate of ADMM}\label{sec:ADMM}
%=================================================================
We can improve the convergence rate of ADMM by using over-relaxation, warm starting the ADMM iterations, choosing the penalty parameter $\rho$ carefully and scaling the problem formulation appropriately.

In over-relaxation, we replace the primal quantity $A\Delta W^{k+1}$ with
%
\[\alpha_{\text{OR}} A\Delta W^{k+1}-(1-\alpha_{\text{OR}})(B \Delta x^{k}-c),\  \alpha_{\text{OR}} \in (1,2),\]
%
 in the update of $\Delta x^k$ and scaled dual variables $\Delta \bar{v}$ and $\Delta \bar{v}_c$. Empirical studies have shown that an $\alpha_{\text{OR}} \in [1.5,1.8]$ may improve the rate of convergence, see \citep{boyd:11}.

We can also improve the convergence rate by warm starting ADMM. That is, using the solution of the previous ADMM iteration as initial condition in the current iteration, see \citep{boyd:11}. The improvement in convergence rate is due to that the primal-dual directions of the interior point method do not change much as the iterates are approaching the solution.

In general it is an open problem how to choose the penalty parameter $\rho$ optimally. Certain heuristics suggest that $\rho$ should be chosen such that the primal and dual residuals converge at the same rate, \citep{boyd:11}. In \cite{2013arXiv1303.6680T}, the authors rescale the optimization problem using a block-diagonal matrix. For the scaled problem, they derive the $\rho$ and $\alpha_{\text{OR}}$ which guarantee the lowest worst-case amount of iterations in ADMM. The scaling matrix can be constructed in a distributed way and since it is block-diagonal it maintains the structure of our problem. However, we would have to recalculate the scaling matrix for each iteration of the interior-point method.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical Experiment}\label{sec:NE}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
To illustrate the proposed method, we apply it to a randomly generated optimization problem of the form  \eqref{eq:DDEOP}. 

In most of the examples, we warm start the ADMM algorithm. That is, we use the previous step direction as an initial point. We have tuned the penalty parameter $\rho$ of ADMM slightly, to provide a better balance of the convergence rate between the primal residual and the dual residual of the ADMM formulation. We do not use the scaling suggested by  \cite{2013arXiv1303.6680T}, although it is believed that this will improve the convergence rate of the search direction calculations.

\subsection{Simulation Set-Up 1}\label{setup1}
We consider fifty subproblems ($N=50$). The total number of variables ($n$), equality constraints ($p$) and inequality constraints ($m$) are 5091, 5089, and 1524, respectively. The number of local variables, local equality constraints and local inequality constraints  are drawn from the standard uniform distributions $U(55,65)$, $U(7,13)$ and $U(27,33)$, respectively. The indices defining the consistency constraints ($J_i$) are drawn from $U(0,900)$ and there are 3017 such constraints.

To ensure that the problem formulation is feasible, we first draw a global variable $x$ and slack variable $S$ from $U(-10,10)$ and $U(1,10)$, respectively. We then generate equality and inequality constraints from these values. The inequality constraints are affine, $A_{in}^ix+b_{in}^i\preceq0$ for $i=1,\dots,N$, where the elements of the matrices $A_{in}^i$ for $i=1,\dots,N$ are drawn from $U(0,1)$. The vectors $b_{in}^i$ for $i=1,\dots,N$ are then calculated. The equality constraints are constructed in the same way; the matrices $\bar{A}^i$ for $i=1,\dots,N$ are drawn from $U(0,1)$. The vectors $b^i$ for $i=1,\dots,N$ are then calculated.

The objective function is quadratic in the global variable with $f_i(x)=x^TP^ix+(q^i)^Tx+e^i$ for $i=1,\dots,N$. The elements of $P^i$, $q^i$ and $e^i$ are drawn from $U(0,1)$, $U(0,1)$ and $U(0,10)$, respectively.

The optimization problem is solved using Algorithm \ref{alg:algNPDexact}, Algorithm \ref{alg:algNPD} and ADMM (Algorithm \ref{alg:alg2}). For comparison, all three algorithms are terminated using the stop criteria of  Algorithm \ref{alg:algNPDexact}. The settings of each specific method are displayed in tables \ref{alg1_param}, \ref{alg2_param} and \ref{alg3_param}. 

We choose a relative value of the upper bounds on the stop criteria, $\epsilon$ and $\epsilon_{\textrm{feas}}$. That is, we choose them as

\begin{align*}
\epsilon=\epsilon_{\textrm{feas}}=10^{-6}\times\textrm{max}\{&1,\|\blkdiag(P^1, \dots,P^N)\|,\|\blkdiag(A_{in}^1, \dots, A_{in}^N)\|,\\
&\|\blkdiag(A^1, \dots,A^N)\|,\|(b_{in}^1,\dots,b_{in}^N)\|,\\
&\|(b^1,\dots,b^N)\|,\|(q^1,\dots,q^N)\| \}.
\end{align*}
For our specific problem generation we get $\epsilon=\epsilon_{\textrm{feas}}=0.0051$, that is we have a scaling factor equal to 5089.

We initialize the methods at the same point. The initial values of the global primal variable are drawn from $U(-10,10)$, the initial values of each local primal variable are given by the consistency constraints. The dual variables are all set to 10, except the dual variables that correspond to the consistency constraint, they are all set to zero. The search directions are initialized to zero if they are not warm started. For the approach using ADMM, we use indicator functions to represent all constraints except consensus. Consequently, we only have dual variables for the consistency constraints. When using ADMM on Problem \eqref{eq:DDEOP}, we get an inequality constrained optimization problem in the first primal variable update which we use \texttt{cvx}, a package for specifying and solving convex programs \citep{cvx}, to solve.

\begin{table}[h!]
\caption{Settings of Algorithm \ref{alg:algNPDexact}.}
\begin{center}
{\begin{tabular}{llllllll}
\hline
Parameter&$\mu$   &$\alpha$ &$\beta$ &$\epsilon_{\text{pri}}$&$\epsilon_{\text{dual}}$&$\rho$ &$\alpha_{\text{OR}}$ \\
 Value & 15 & 0.01 &0.5 &$50/2\times10^{-20}$&$50/2\times10^{-20}$ & 0.5 &1\\
\hline
\end{tabular}}
\end{center}
\label{alg1_param}
\end{table}

\begin{table}[h!]
\caption{Settings of Algorithm \ref{alg:algNPD}.}
\begin{center}
{\begin{tabular}{llllllll}
\hline
Parameter&$\eta_{\textrm{max}}^i$   &$\gamma^{i,(0)}$ &$\beta$ & $\theta$&$\rho$ &$\alpha_{\text{OR}}$ & $\epsilon_{\sigma}$ \\
 Value & 0.9 & 0.9 &0.1 &0.95& 0.5 &1 & 0.1\\
\hline
\end{tabular}}
\end{center}
\label{alg2_param}
\end{table}

\begin{table}[h!]
\caption{Settings of ADMM.}
\begin{center}
{\begin{tabular}{lll}
\hline
Parameter&$\rho$ &$\alpha_{\text{OR}}$ \\
 Value& 0.5 &1\\
\hline
\end{tabular}}
\end{center}
\label{alg3_param}
\end{table}

\subsubsection{Algorithm \ref{alg:algNPD} with and without Warm Starting}
We first compare Algorithm \ref{alg:algNPD} with and without warm starting. The total number of ADMM iterations is 17453 with warm starting and 21181 without warm starting. In Figure \ref{fig:iterADMMEx1}, the number of ADMM iterations in each instance of the primal-dual method is displayed. In Figure \ref{fig:pertADMMEx1}, the value of the perturbation parameter $\mu$ in each instance of the primal-dual method is shown. We benefit from warm staring ADMM with respect to the total number of ADMM iterations necessary for the primal-dual method to converge. However, for a specific iteration of the primal-dual method the warm started approach can require more iterations than the approach without warm starting, and the number of saved ADMM iterations fluctuate over the primal-dual method iterations. A smaller saving, or even no saving at all, occur when the new system of equations to solve is changed considerably from the previous one. Note, also that the number of iterations of the primal-dual method is affected. The number of iterations is 80 with warm starting and 81 without warm starting.

\begin{figure}[h!]
\begin{center}
%\begin{subfigure}
 \centering
\figuretikz{numberOfIterationsADMMEx1.tikz}{0.5}{0.5}
\caption{ADMM iterations. The number of iterations of ADMM with warm starting \eqref{solid} and without warm starting \eqref{dashed} in each instance of the primal-dual method are shown.} \label{fig:iterADMMEx1}
% \end{subfigure}
\end{center}
\end{figure}
\begin{figure}[h!]
\begin{center}
%\begin{subfigure}
 \centering
\figuretikz{perturbationADMMDistEx1.tikz}{0.5}{0.5}
\caption{Perturbation of KKT conditions.  The value of the perturbation parameter $\mu$ with warm starting ADMM \eqref{solid} and without warm starting \eqref{dashed} ADMM in each instance of the primal-dual method are shown.} \label{fig:pertADMMEx1}
%\end{subfigure}
\end{center}
\end{figure}

\subsubsection{Savings in ADMM Iterations Using Algorithm \ref{alg:algNPD} Compared to Algorithm \ref{alg:algNPDexact}}
We compare Algorithm \ref{alg:algNPDexact} and Algorithm \ref{alg:algNPD}, using warm starting in both methods. The total number of ADMM iterations is 38796 with Algorithm \ref{alg:algNPDexact} and 17453 with Algorithm \ref{alg:algNPD}. The number of iterations of the primal-dual method is 21 using Algorithm \ref{alg:algNPDexact} and 80 using Algorithm \ref{alg:algNPD}. We save 55\% in ADMM iterations, using Algorithm \ref{alg:algNPD}. This is due to the adaptive stop criteria of the search direction calculations. In particular, it dramatically decreases the amount of ADMM iterations necessary in the first couple of iterations of the primal-dual method. The convergence rate of Algorithm \ref{alg:algNPD} is highly dependent on the setting of the parameters in the method. With a different setting, we could get a much worse convergence rate.

\subsubsection{Stop Criteria}
We compare how the residuals and surrogate duality gap evolves between the three methods. We see in figures \ref{fig:primalResADMMEx1}, \ref{fig:dualResADMMEx1} and \ref{fig:surADMMEx1}, that the averaged, over all subproblems, value of the norm of local primal residuals, dual residuals and surrogate duality gaps are constantly decreasing as each method iterates. The same behavior is obtained for the global residuals and surrogate duality gap, see figures \ref{fig:globprimalResADMMEx1}, \ref{fig:dualResADMMEx1} and \ref{fig:surADMMEx1}. In addition, we can see that the stop criteria are fulfilled for all subproblems at approximately the same iteration. When comparing the local stop criteria with the global, we see that the local criteria can be conservative. In fact, if the global stop criteria were to be used instead of the local one, ADMM would terminate at an earlier iteration than with the local criteria, see Table \ref{alg_residuals} and figures \ref{fig:primalResADMMEx1}--\ref{fig:globsurADMMEx1}.

\begin{table}[h!]
\caption{Fulfillment of stop criteria.}
\begin{center}
{\begin{tabular}{llll}
\hline
&iteration $l$ at which&iteration $l$ at which&iteration $l$ at which\\
& local (global) primal& local (global) dual&local (global) gap\\
&constraint is ful-&constraint is ful-&constraint is ful-\\
&filled&filled&filled\\
\hline
 Alg. \ref{alg:algNPDexact} & 19 (18) &19 (20)& 21 (21)\\
 Alg. \ref{alg:algNPD} & 72 (68) & 75 (75)& 80 (80)\\
 ADMM & 304 (275)& 187 (169)& 0 (0)\\
\hline
\end{tabular}}
\end{center}
\label{alg_residuals}
\end{table}

\begin{figure}[h!]
\begin{center}
\begin{subfigure}
 \centering
\figuretikz{primalresidualADMMOrigEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{primalresidualADMMDistEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{primalresidualADMMFirstEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\caption{Norm of primal local residual.  The averaged, over all subproblems, value of the norm of the primal local residual along with its standard deviation \eqref{solid}  and upper bound, $\textrm{log}\left(\epsilon_{\textrm{feas}}^2/N\right)$, \eqref{dashed} are shown. In figures (a), (b) and (c), algorithms \ref{alg:algNPDexact}, \ref{alg:algNPD} and ADMM have been used, respectively.} \label{fig:primalResADMMEx1}\end{center}
\end{figure}
\begin{figure}[h!]
\begin{center}
\begin{subfigure}
 \centering
\figuretikz{globprimalresidualADMMOrigEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{globprimalresidualADMMDistEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{globprimalresidualADMMFirstEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\caption{Norm of primal global residual.  The value of the norm of the primal global residual \eqref{solid} and its upper bound $\textrm{log}\left(\epsilon_{\textrm{feas}}\right)$, \eqref{dashed} are shown. In figures (a), (b) and (c), algorithms \ref{alg:algNPDexact}, \ref{alg:algNPD} and ADMM have been used, respectively.} \label{fig:globprimalResADMMEx1}\end{center}
\end{figure}

\begin{figure}[h!]
\begin{center}
\begin{subfigure}
 \centering
\figuretikz{dualresidualADMMOrigEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{dualresidualADMMDistEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{dualresidualADMMFirstEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\caption{Norm of dual local residual.  The averaged, over all subproblems, value of the norm of the dual local residual along with its standard deviation \eqref{solid}  and upper bound, $\textrm{log}\left(\epsilon_{\textrm{feas}}^2/N\right)$, \eqref{dashed} are shown. In figures (a), (b) and (c), algorithms \ref{alg:algNPDexact}, \ref{alg:algNPD} and ADMM have been used, respectively.} \label{fig:dualResADMMEx1}\end{center}
\end{figure}
\begin{figure}[h!]
\begin{center}
\begin{subfigure}
 \centering
\figuretikz{globdualresidualADMMOrigEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{globdualresidualADMMDistEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{globdualresidualADMMFirstEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\caption{Norm of dual global residual.  The value of the norm of the dual global residual \eqref{solid} and its upper bound $\textrm{log}\left(\epsilon_{\textrm{feas}}\right)$, \eqref{dashed} are shown. In figures (a), (b) and (c), algorithms \ref{alg:algNPDexact}, \ref{alg:algNPD} and ADMM have been used, respectively.} \label{fig:globdualResADMMEx1}\end{center}
\end{figure}

\begin{figure}[h!]
\begin{center}
\begin{subfigure}
 \centering
\figuretikz{surADMMOrigEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{surADMMDistEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{surADMMFirstEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\caption{Local surrogate duality gap.  The averaged, over all subproblems, value of the surrogate duality gap along with its standard deviation \eqref{solid}  and upper bound, $\textrm{log}\left(\epsilon/N\right)$, \eqref{dashed} are shown. In figures (a), (b) and (c), algorithms \ref{alg:algNPDexact}, \ref{alg:algNPD} and ADMM have been used, respectively.} \label{fig:surADMMEx1}\end{center}
\end{figure}
\begin{figure}[h!]
\begin{center}
\begin{subfigure}
 \centering
\figuretikz{globsurADMMOrigEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{globsurADMMDistEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{globsurADMMOnTopEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\caption{Global surrogate duality gap.  The value of the surrogate duality gap \eqref{solid} and its upper bound, $\textrm{log}\left(\epsilon/N\right)$, \eqref{dashed} are shown. In figures (a), (b) and (c), algorithms \ref{alg:algNPDexact}, \ref{alg:algNPD} and ADMM have been used, respectively.} \label{fig:globsurADMMEx1}\end{center}
\end{figure}

\subsubsection{Optimal Value}
We compare how the value of the objective function evolves for each iteration of the three methods. In Figure \ref{fig:optvalADMMEx1}, the relative error of the objective function is shown. The so-called true optimal value is obtained using \texttt{cvx}. We see that the relative error evolves to approximately the same level for all three methods. However, ADMM gives a slightly higher value (the relative error is $1.2\times 10^{-5}$) than Algorithm \ref{alg:algNPDexact} (the relative error is $8.4\times 10^{-8}$) and Algorithm \ref{alg:algNPD} (the relative error is $6.3\times 10^{-8}$). Note that since the iterates are not necessarily feasible in each iteration, it is possible to obtain a value of the objective function that is very close to or even smaller than the true optimal value, while the algorithm has not yet terminated. This explains the dip of the relative error which can be seen in Figure \ref{fig:optvalADMMEx1}. 

\begin{figure}[h!]
\begin{center}
\begin{subfigure}
 \centering
\figuretikz{optValADMMOrigEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{optValADMMDistEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\begin{subfigure}
 \centering
\figuretikz{optValADMMEx1.tikz}{0.5}{0.5}
 \end{subfigure}
\caption{Relative error of objective function.  The relative error of the objective function \eqref{solid} is shown. In figures (a), (b) and (c), algorithms \ref{alg:algNPDexact}, \ref{alg:algNPD} and ADMM have been used, respectively. The true optimal value is obtained using \texttt{cvx} and is denoted $p^*$. The value of the objective function is denoted $p$.} \label{fig:optvalADMMEx1}\end{center}
\end{figure}

\subsubsection{Total Number of Iterations}
The total number of iterations, that is ADMM iterations for Algorithm \ref{alg:algNPDexact} and Algorithm \ref{alg:algNPD}, and interior-point iterations in \texttt{cvx} for ADMM are compared. We get the total number of iterations equal to 38796, 17453 and 8194, for each method respectively. Thus, for this specific simulation set-up, ADMM beats the proposed algorithms in terms of the number of iterations. Note, however, that Algorithm \ref{alg:algNPDexact} and Algorithm \ref{alg:algNPD} are computationally very cheap. The major computational effort takes place when calculating the search directions of the local primal variables, that is Step 9 in Algorithm \eqref{alg:alg4PD}, see Remark 3. If $\rho$ is kept constant, we only have to factorize $H^{i,(l)}_{\text{pd}} + \rho\left( I + (\bar A^i)^T\bar A^i\right)$ for $i = 1, \dots, N,$ once in each primal-dual iteration. However, if we use ADMM  (Algorithm \ref{alg:alg2}) we have to factorize the corresponding matrices, which are of the same sizes as in the other algorithms, for each interior-point iteration in each ADMM iteration. That is, instead of factorizing the matrices 21 or 80 times with Algorithm \ref{alg:algNPDexact} or Algorithm \ref{alg:algNPD}, respectively, we have to factorize the matrices 8194 times with ADMM. We pay a price for the savings in the number of factorizations necessary. In Algorithm \ref{alg:algNPDexact} and Algorithm \ref{alg:algNPD} the nodes have to communicate in each inner iteration (ADMM iteration) whilst for ADMM applied to the original problem the nodes only have to communicate in each outer iteration (ADMM iteration). 

The convergence rate of Algorithm \ref{alg:algNPD} is sensitive to the setting of its parameters. How to choose them optimally, or even wisely, is an open question. Note though, that we can make use of an ad-hoc adaptive stop criteria for the search direction calculation in Algorithm \ref{alg:algNPDexact} which does not guarantee global convergence. Such ad-hoc criteria, can yield much faster convergence rate without much tuning of the settings. For example, with $\epsilon_{pri}=2\epsilon_{dual}=50/4\times10^{-3}$ for $l\leq 5$, $\epsilon_{pri}=2\epsilon_{dual}=50/4\times10^{-5}$ for $6\leq l\leq10$ and $\epsilon_{pri}=2\epsilon_{dual}=50/4\times10^{-8}$ for $l\geq 11$, we get 5469 number of iterations, a 33\% saving compared to ADMM. This result provides an incitement to investigate further the convergence rate properties of Algorithm \ref{alg:algNPD}.

\subsection{Simulation Set-Up 2}
We consider ten subproblems ($N=10$). The rest of the problem set-up coincides with that of Section \ref{setup1}.

The optimization problem is solved using Algorithm \ref{alg:algNPD} and ADMM (Algorithm \ref{alg:alg2}). For comparison, both algorithms are terminated using the stop criteria of  Algorithm \ref{alg:algNPDexact}. We initialize the methods and use the same settings as in Section \ref{setup1}. We perform 50 Monte-Carlo runs of the simulation set-up.

\subsubsection{Total Number of Iterations}
The total number of iterations, that is ADMM iterations for Algorithm \ref{alg:algNPD} and inner-point iterations in \texttt{cvx} for ADMM, are compared. We get the averaged, over all Monte-Carlo runs, total number of iterations equal to 569 (the standard deviation is 62) and 1257 (the standard deviation is 29), for each method respectively. Thus, for this specific simulation set-up, Algorithm \ref{alg:algNPD} beats ADMM in terms of the number of iterations. In fact, Algorithm \ref{alg:algNPD} terminates after 55\% less iterations than ADMM (with respect to the averaged value).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}\label{sec:C}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We have proposed two efficient distributed primal-dual interior-point method for loosely coupled problems using ADMM (Algorithm \ref{alg:algNPDexact} and Algorithm \ref{alg:algNPD}). Due to the nature of the interior-point method, the loosely coupled structure of the problem is preserved in the linear system of equations that provides the primal-dual directions. ADMM takes advantage of this structure and makes the direction calculations highly parallellizable. Consequently, the proposed method has superior computational properties with respect to other distributed algorithms. Of course, we can use Algorithm \ref{alg:algNPDexact} and Algorithm \ref{alg:algNPD} on problems with completely coupled structure as well, but we can not expect the same superior properties as for the loosely coupled structure.

The latter of the methods (Algorithm \ref{alg:algNPD}) adaptively chooses the required accuracy in the termination condition of the inner iterations (ADMM iterations) with respect to the accuracy obtained in the outer iterations (interior-point method iterations). This is to avoid unnecessary inner iterations when the accuracy of the current outer iteration is low, as elaborated in \cite{bon:05}. We have stated under which assumptions the method converges to the optimal solution. In addition, we have provided comparisons between Algorithm~ \ref{alg:algNPDexact}, Algorithm \ref{alg:algNPD} and ADMM in simulation.
\vspace{20pt}
\section*{Funding}
This work was partially supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007–2013) / ERC Grant Agreement No. 267381, the Swedish Research Council and the Linnaeus Center ACCESS at KTH.

\vspace{94pt}
\begin{thebibliography}{51}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{DOI~\discretionary{}{}{}#1}\else
  \providecommand{\doi}{DOI~\discretionary{}{}{}\begingroup
  \urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{Al-Jeiroudi and Gondzio(2009)}]{al:09}
Al-Jeiroudi G, Gondzio J (2009) Convergence analysis of the inexact infeasible
  interior-point method for linear optimization. Journal of Optimization Theory
  and Applications 141(2):231--247

\bibitem[{Arrow(1964)}]{arr:64}
Arrow KJ (1964) Studies in linear and non-linear programming. Stanford
  University Press Stanford, Calif.

\bibitem[{Bellavia(1998)}]{bel:98}
Bellavia S (1998) Inexact interior-point method. Journal of Optimization Theory
  and Applications 96(1):109--121

\bibitem[{Bellavia and Pieraccini(2004)}]{bel:04}
Bellavia S, Pieraccini S (2004) Convergence analysis of an inexact infeasible
  interior point method for semidefinite programming. Computational
  Optimization and Applications 29(3):289--313

\bibitem[{Benzi et~al(2005)Benzi, Golub, and Liesen}]{ben:05}
Benzi M, Golub GH, Liesen J (2005) {Numerical solution of saddle point
  problems}. Acta Numerica 14(-1):1--137

\bibitem[{Bertsekas and Tsitsiklis(1997)}]{ber:97}
Bertsekas DP, Tsitsiklis JN (1997) Parallel and Distributed Computation:
  Numerical Methods. Athena Scientific

\bibitem[{Bonettini and Ruggiero(2007)}]{bon:07}
Bonettini S, Ruggiero V (2007) Some iterative methods for the solution of a
  symmetric indefinite {KKT} system. Computational Optimization and
  Applications 38(1):3--25

\bibitem[{Bonettini et~al(2005)Bonettini, Galligani, and Ruggiero}]{bon:05}
Bonettini S, Galligani E, Ruggiero V (2005) An inexact {N}ewton method combined
  with {H}estenes multipliers scheme for the solution of {Karush-Kuhn-Tucker}
  systems. Applied Mathematics and Computation 168(1):651--676

\bibitem[{Boyd and Vandenberghe(2004)}]{boyd:04}
Boyd S, Vandenberghe L (2004) Convex Optimization. Cambridge University Press

\bibitem[{Boyd et~al(2011)Boyd, Parikh, Chu, Peleato, and Eckstein}]{boyd:11}
Boyd S, Parikh N, Chu E, Peleato B, Eckstein J (2011) Distributed optimization
  and statistical learning via the alternating direction method of multipliers.
  Foundations and Trends in Machine Learning 3(1):1--122

\bibitem[{{Chen} et~al(2013){Chen}, {He}, {Ye}, and {Yuan}}]{che:13}
{Chen} C, {He} B, {Ye} Y, {Yuan} X (2013) {The Direct Extension of ADMM for
  Multi-block Convex Minimization Problems is Not Necessarily Convergent}.
  Optimization Online e-prints

\bibitem[{Chu et~al(2011)Chu, Gorinevsky, and Boyd}]{Chu:11}
Chu E, Gorinevsky D, Boyd S (2011) Scalable statistical monitoring of fleet
  data. In: Proceedings of the 18th IFAC World Congress, Milan, Italy, pp
  13,227--13,232

\bibitem[{Combettes and Pesquet(2011)}]{com:11}
Combettes PL, Pesquet JC (2011) Proximal splitting methods in signal
  processing. In: Fixed-Point Algorithms for Inverse Problems in Science and
  Engineering, Springer Optimization and Its Applications, vol~49, Springer New
  York, pp 185--212

\bibitem[{Curtis et~al(2010)Curtis, Schenk, and W\"achter}]{cur:10}
Curtis F, Schenk O, W\"achter A (2010) An interior-point algorithm for
  large-scale nonlinear optimization with inexact step computations. SIAM
  Journal on Scientific Computing 32(6):3447--3475

\bibitem[{Curtis et~al(2012)Curtis, Huber, Schenk, and W\"achter}]{cur:12}
Curtis FE, Huber J, Schenk O, W\"achter A (2012) A note on the implementation
  of an interior-point algorithm for nonlinear optimization with inexact step
  computations. Mathematical Programming 136(1):209--227

\bibitem[{Douglas and Rachford(1956)}]{dou:56}
Douglas J, Rachford HH (1956) {On the numerical solution of heat conduction
  problems in two and three space variables}. Transaction of the American
  Mathematical Society 82:421--489

\bibitem[{Durazzi and Ruggiero(2003)}]{dur:03}
Durazzi C, Ruggiero V (2003) A {N}ewton inexact interior-point method for large
  scale nonlinear optimization problems. {Annali dell'Universit\'a di Ferrara}
  49(1):333--357

\bibitem[{Eckstein(1989)}]{eck:89}
Eckstein J (1989) Splitting methods for monotone operators with application to
  parallel optimization. Ph{D} dissertation, {Massachussets Intitute of
  Technology}

\bibitem[{Eckstein and Bertsekas(1992)}]{eck:92}
Eckstein J, Bertsekas DP (1992) {On the Douglas-Rachford splitting method and
  the proximal point algorithm for maximal monotone operators}. Mathematical
  Programming 55:293--318

\bibitem[{Eisenstat and Walker(1994)}]{eis:94}
Eisenstat SC, Walker HF (1994) {Globally Convergent Inexact Newton Methods}.
  SIAM Journal on Optimization 4(2):393--422

\bibitem[{El-Bakry et~al(1996)El-Bakry, Tapia, Tsuchiya, and Zhang}]{bakry:96}
El-Bakry A, Tapia R, Tsuchiya T, Zhang Y (1996) On the formulation and theory
  of the newton interior-point method for nonlinear programming. Journal of
  Optimization Theory and Applications 89(3):507--541,
  \doi{10.1007/BF02275347},
  \urlprefix\url{http://dx.doi.org/10.1007/BF02275347}

\bibitem[{Fortin and Glowinski(1983)}]{for:83}
Fortin M, Glowinski R (1983) Chapter 1 augmented lagrangian methods in
  quadratic programming. In: Fortin M, Glowinski R (eds) Augmented Lagrangian
  Methods: Applications to the Numerical Solution of Boundary-Value Problems,
  Studies in Mathematics and Its Applications, vol~15, Elsevier, pp 1 -- 46

\bibitem[{Freund et~al(1999)Freund, Jarre, and Mizuno}]{fre:99}
Freund RW, Jarre F, Mizuno S (1999) Convergence of a class of inexact
  interior-point algorithms for linear programs. Mathematics of Operations
  Research 24(1):50--71

\bibitem[{Gabay(1983)}]{gab:83}
Gabay D (1983) Applications of the method of multipliers to variational
  inequalities. In: Fortin M, Glowinski R (eds) Augmented {L}agrangian Methods:
  {A}pplications to the Solution of Boundary-Value Problems, North-Holland

\bibitem[{Gabay and Mercier(1976)}]{gab:76}
Gabay D, Mercier B (1976) A dual algorithm for the solution of nonlinear
  variational problems via finite element approximation. Computers \&
  Mathematics with Applications 2(1):17--40

\bibitem[{{Glowinski} and {Marroco}(1975)}]{mar:75}
{Glowinski} R, {Marroco} A (1975) {Sur l'approximation, par \'el\'ements finis
  d'ordre un, et la r\'esolution, par penalisation-dualit\'e, d'une classe de
  probl\`emes de Dirichlet non lin\'eaires.} {Rev Franc Automat Inform Rech
  Operat, R} 9(2):41--76

\bibitem[{Goldfarb et~al(2012)Goldfarb, Ma, and Scheinberg}]{gol:12}
Goldfarb D, Ma S, Scheinberg K (2012) Fast alternating linearization methods
  for minimizing the sum of two convex functions. Mathematical Programming pp
  1--34

\bibitem[{Goldstein et~al(2012)Goldstein, ODonoghue, and Setzer}]{gols:12}
Goldstein T, O'Donoghue B, Setzer S (2012) Fast alternating direction
  optimization methods. Tech. Rep. CAM report 12-35, UCLA

\bibitem[{Grant and Boyd(2014)}]{cvx}
Grant M, Boyd S (2014) {CVX}: Matlab software for disciplined convex
  programming, version 2.1. \url{http://cvxr.com/cvx}

\bibitem[{Han and Yuan(2012)}]{han:12}
Han D, Yuan X (2012) A note on the alternating direction method of multipliers.
  Journal of Optimization Theory and Applications 155(1):227--238

\bibitem[{Hansson(2000)}]{han:00}
Hansson A (2000) A primal-dual interior-point method for robust optimal control
  of linear discrete-time systems. IEEE Transactions on Automatic Control
  45(9):1639--1655

\bibitem[{Hestenes(1969)}]{hes:69}
Hestenes MR (1969) Multiplier and gradient methods. Journal of Optimization
  Theory and Applications 4(5):303--320, \doi{10.1007/BF00927673}

\bibitem[{{Hong} and {Luo}(2012)}]{hon:12}
{Hong} M, {Luo} ZQ (2012) {On the Linear Convergence of the Alternating
  Direction Method of Multipliers}. ArXiv e-prints

\bibitem[{Iutzeler et~al(2012)Iutzeler, Ciblat, and Jakubowicz}]{Iut:12}
Iutzeler F, Ciblat P, Jakubowicz J (2012) Analysis of {Max-Consensus}
  algorithms in wireless channels. IEEE Transactions on Signal Processing
  60(11):6103--6107

\bibitem[{Korzak(2000)}]{kor:00}
Korzak J (2000) Convergence analysis of inexact infeasible-interior-point
  algorithms for solving linear programming problems. SIAM Journal on
  Optimization 11(1):133--148

\bibitem[{Liu et~al(2013)Liu, Hansson, and Vandenberghe}]{liu:13}
Liu Z, Hansson A, Vandenberghe L (2013) Nuclear norm system identification with
  missing inputs and outputs. Systems \& Control Letters 62(8):605--612

\bibitem[{Lu et~al(2009)Lu, Monteiro, and O'Neal}]{lu:09}
Lu Z, Monteiro RDC, O'Neal JW (2009) An iterative solver-based long-step
  infeasible primal-dual path-following algorithm for convex qp based on a
  class of preconditioners. Optimization Methods Software 24(1):123--143

\bibitem[{Mizuno and Jarre(1999)}]{miz:99}
Mizuno S, Jarre F (1999) Global and polynomial-time convergence of an
  infeasible-interior-point algorithm using inexact computation. Mathematical
  Programming 84(1):105--122

\bibitem[{Necoara and Suykens(2009)}]{nec:09}
Necoara I, Suykens JAK (2009) Interior-point lagrangian decomposition method
  for separable convex optimization. Journal of Optimization Theory and
  Applications 143(3):567--588

\bibitem[{Nedic and Ozdaglar(2009)}]{ned:09}
Nedic A, Ozdaglar A (2009) Distributed subgradient methods for multi-agent
  optimization. IEEE Transactions on Automatic Control 54(1):48--61

\bibitem[{Nedic et~al(2010)Nedic, Ozdaglar, and Parrilo}]{ned:10}
Nedic A, Ozdaglar A, Parrilo P (2010) Constrained consensus and optimization in
  multi-agent networks. IEEE Transactions on Automatic Control 55(4):922--938

\bibitem[{{Ohlsson} et~al(2013){Ohlsson}, {Chen}, {Khoshfetrat Pakazad},
  {Ljung}, and {Shankar Sastry}}]{ohl:13}
{Ohlsson} H, {Chen} T, {Khoshfetrat Pakazad} S, {Ljung} L, {Shankar Sastry} S
  (2013) {Scalable Anomaly Detection in Large Homogenous Populations}. ArXiv
  e-prints

\bibitem[{Parikh and Boyd(2014)}]{boy:14}
Parikh N, Boyd S (2014) Proximal algorithms. Foundations and Trends in
  Optimization 1(3):127--239

\bibitem[{Peaceman and Rachford(1955)}]{pea:55}
Peaceman DW, Rachford J H~H (1955) The numerical solution of parabolic and
  elliptic differential equations. Journal of the Society for Industrial and
  Applied Mathematics 3(1):28--41

\bibitem[{Powell(1969)}]{pow:69}
Powell MJD (1969) {A method for nonlinear constraints in minimization
  problems}. In: Fletcher R (ed) Optimization, Academic Press, New York, pp
  283--298

\bibitem[{Saad(2003)}]{saa:03}
Saad Y (2003) Iterative Methods for Sparse Linear Systems, 2nd edn. Society for
  Industrial and Applied Mathematics, Philadelphia, PA, USA

\bibitem[{Summers and Lygeros(2012)}]{sum:12}
Summers T, Lygeros J (2012) Distributed model predictive consensus via the
  alternating direction method of multipliers. In: 50th Annual Allerton
  Conference on Communication, Control, and Computing, pp 79--84

\bibitem[{{Teixeira} et~al(2013){Teixeira}, {Ghadimi}, {Shames}, {Sandberg},
  and {Johansson}}]{2013arXiv1303.6680T}
{Teixeira} A, {Ghadimi} E, {Shames} I, {Sandberg} H, {Johansson} M (2013)
  {Optimal scaling of the ADMM algorithm for distributed quadratic
  programming}. ArXiv e-prints \eprint{1303.6680}

\bibitem[{Toh(2008)}]{toh:08}
Toh K (2008) An inexact primal-dual path following algorithm for convex
  quadratic {SDP}. Mathematical Programming 112(1):221--254

\bibitem[{Wei et~al(2013)Wei, Ozdaglar, and Jadbabaie}]{wei:13}
Wei E, Ozdaglar A, Jadbabaie A (2013) A distributed {N}ewton method for network
  utility maximization--{I}: {A}lgorithm. IEEE Transactions on Automatic
  Control 58(9):2162--2175

\bibitem[{Zhou and Toh(2004)}]{zho:04}
Zhou G, Toh K (2004) Polynomiality of an inexact infeasible interior point
  algorithm for semidefinite programming. Mathematical Programming
  99(2):261--282

\end{thebibliography}

\appendix

\section{Global Convergence of Distributed Inexact Primal-Dual Interior-Point Method}\label{app:global}
Here we state the proof for global convergence of the proposed distributed primal-dual inexact interior-point method. To this end, we use the definitions of $\Omega^i(\epsilon)$ for $i=1,\dots,N$ and assumptions B1--B4. The lemmas, theorems and proofs are adapted with very minor changes from \cite{eis:94} and \cite{bel:98}, more detailed references are given before each lemma and theorem. They are all included here for the sake of completeness.

The collection of lemmas and theorems in this appendix enables us to show that, through the run of the algorithm, the iterates are persistently updated and guaranteed to converge to the optimal solution. We first address the concept of \emph{break down} of the algorithm. That is, when we are not able to find a suitable step direction and step size \cite{eis:94}. We show that such a break down of the algorithm will not occur under our assumptions as defined in Section \ref{distConRes}. Second, we show that the algorithm is convergent towards an optimal solution. 

Particularly, Theorem \ref{notZero} and  Lemma \ref{breakdown} assure a persistent update of the iterates. The theorem illustrates that the upper bound of the step size is bounded away from zero. Lemma \ref{breakdown} states that, given a suitable search direction, it is always possible to find a step size which yields a satisfactory decrease in the merit function $\|H^i(z^i)\|$ for $i=1,\dots,N$.

The subsequent theorems are related to the convergence of the algorithm, that is, that the generated $z^{(l)}$ converges and $\|H(z^{(l)})\|$ converges to zero. Particularly, Theorem \eqref{conLimitPointTheorem} states that the sequence of iterates  $\{z^{(l)}\}$ generated by our method converges to a point $z_*$, and Theorem \eqref{conLimitPointTheoremPre} is a result necessary for the proof of Theorem \eqref{conLimitPointTheorem}.  Finally, in Theorem \eqref{globalCon}, we state under which assumptions the proposed algorithm is global convergent, that is, generates a sequence of iterates such that $\|H(z^{(l)})\|\rightarrow 0$. Theorem \eqref{divergent} is used in the proof of Theorem \eqref{globalCon}.
\vspace{20pt}
\subsection{Break Down}
There are four steps in the algorithm where a possible break down can occur. They are when calculating the step direction (Step 9 in Algorithm \ref{alg:algNPD}), choosing the intermediate step sizes $\alpha_1^i$ and $\alpha_2^i$ (steps 10-13 in Algorithm \ref{alg:algNPD}) and calculating the actual step size $\alpha$ (steps 10-13 in Algorithm \ref{alg:algNPD}). The hypothesis under which our proposed algorithm\emph{ does not} break down is addressed in the rest of this section.

A step direction can always be calculated provided $H'(z)$ is invertible. Hence, if assumption B3 is fulfilled, the algorithm will not break down at Step 9.  Next,  Theorem \ref{notZero} shows that the intermediate step sizes for any such step direction is always bounded away from zero. This theorem is based on Theorem 3.2 in \cite{bel:98}.
\vspace{10pt}
\begin{theorem}\label{notZero}
Assume $\{z^{(l)}\}$ is generated by Algorithm \ref{alg:algNPD} and assumptions B1-- B4 are fulfilled, then the sequence $\{\alpha^{(l)}\}$ with $\alpha^{(l)}=\underset{i}{\text{min}}(\alpha^{i,(l)})$ and $\alpha^{i,(l)}=\text{min}(\alpha^{i}_1,\alpha^{i}_2)$ is bounded away from zero.
\end{theorem}
\vspace{10pt}
\begin{proof}
The proof follows closely that of Theorem 3.2 in \cite{bel:98}. We first show that the sum of the complementary KKT condition is bounded away from zero in $\Omega^i(\epsilon)$. First, note that $(s^{i})^T\lambda^{i}\geq\|S^{i}\Lambda^{i}\hat{e}\|$ and $(s^{i})^T\lambda^{i}\geq\bar{\tau}_2^i\gamma^{i}\|R^i(z^{i})\|$. Using the quadratic mean we get
\begin{align*}
(s^{i})^T\lambda^{i}&\geq \sqrt{\frac{\|S^{i}\Lambda^{i}\hat{e}\|^2+(\bar{\tau}_2^i\gamma^{i}\|R^i(z^{i})\|)^2}{2}}\\
&\geq \textrm{min}\left(1,\frac{\bar{\tau}_2^i}{2}\right)\|H^i(z^{i,(l)})\|\frac{1}{\sqrt{2}}\\
&\geq\textrm{min}\left(1,\frac{\bar{\tau}_2^i}{2}\right)\frac{\epsilon}{N}\frac{1}{\sqrt{2}}.
\end{align*}
The right hand side of the last inequality is independent of the iteration $l$ and strictly larger than zero, thus $(s^{i})^T\lambda^{i}$ is bounded away from zero in $\Omega^i(\epsilon)$ for $i=1,\dots,N$. Due to assumptions B1, B3 and the fact that $\|\hat{r}^{(l)}\|$ is bounded, we have that $\|\Delta z^{i}\|$ is bounded and there exists $M_1^{i}>0$ and $M_2^i>0$ such that

$$|\Delta s_j^i\Delta \lambda_j^i-\frac{\tau_1^i\gamma^i}{m_i}(\Delta s^i)^T\Delta \lambda^i|\leq M_1^i, \textrm{ for all } i=1,\dots,N$$

\noindent and

$$|(\Delta s^i)^T\Delta \lambda^i-\tau_2^i\gamma^iL^i\|\Delta z^{i,(l)}\|^2|\leq M_2^i, \textrm{ for all } i=1,\dots,N,$$

\noindent respectively. Furthermore, we have $\hat{r}^{(l)}=(\hat{r}_1^{(l)},0)$, where the zero corresponds to the complementary KKT condition which is solved exactly in Algorithm \ref{alg:algNPD}.

To determine $\alpha^{i}_1$ at iteration $l$, we consider the expression (where we omit the superscript $l$)
\begin{align*}
s_j^{i}(\alpha)\lambda_j^{i}(\alpha)-\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(s^{i}(\alpha))^T\lambda^{i}(\alpha)&=\\
(s_j^{i}+\alpha \Delta s_j^{i})(\lambda_j^{i}+\alpha \Delta\lambda_j^{i})-\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(s^{i}+\alpha \Delta s^{i})^T(\lambda^{i}+\alpha \Delta\lambda^{i})&=\\
\left(\!s_j^{i}\lambda_j^{i}\!-\!\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(s^{i})^T\lambda^{i}\!\right)\!+\!\alpha\left(\!\Delta\lambda_j^{i} s_j^{i}+ \Delta s_j^{i}\lambda_j^{i}\! -\!\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(s^{i})^T\Delta\lambda^{i}\! -\!\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(\lambda^{i})^T\Delta s^{i}\!\right)&+\\
\alpha^2\left(\Delta s_j^{i}\Delta\lambda_j^{i}-\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(\Delta s^{i})^T\Delta\lambda^{i}\right)&=\\
\left(\!s_j^{i}\lambda_j^{i}\!-\!\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(s^{i})^T\lambda^{i}\!\right)+\alpha\left(\!-s_j^{i}\lambda_j^{i}+\mu-\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}\hat{e}^T(-S^{i}\Lambda^{i}\hat{e}+\mu\hat{e})\right)&+\\
\alpha^2\left(\Delta s_j^{i}\Delta\lambda_j^{i}-\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(\Delta s^{i})^T\Delta\lambda^{i}\right)&=\\
(1-\alpha)\!\left(\!s_j^{i}\!\lambda_j^{i}\!-\!\frac{\bar{\tau}_1^i\gamma^{i}}{m_i}(s^{i})^T\lambda^{i}\!\right)\!+\alpha(1-\bar{\tau}_1^i\gamma^{i})\mu\!+\!\alpha^2\!\left(\!\Delta s_j^i\Delta \lambda_j^i-\frac{\tau_1^i\gamma^i}{m_i}(\Delta s^i)^T\Delta \lambda^i\!\right)\!&\geq\\
\alpha(1-\bar{\tau}_1^i\gamma^{i})\mu+\alpha^2\left(\Delta s_j^i\Delta \lambda_j^i-\frac{\tau_1^i\gamma^i}{m_i}(\Delta s^i)^T\Delta \lambda^i\right)&\geq\\
\alpha(1-\bar{\tau}_1^i\gamma^{i})\mu-\alpha^2\left |\Delta s_j^i\Delta \lambda_j^i-\frac{\tau_1^i\gamma^i}{m_i}(\Delta s^i)^T\Delta \lambda^i\right |&\geq\\
\alpha(1-\bar{\tau}_1^i\gamma^{i})\mu-\alpha^2M_1^i.
\end{align*}
Since $\alpha_1^i$ is defined as
\begin{align*}
\alpha_1^i=\underset{\alpha \in [0,1]}{\textrm{max}}\left\{\alpha|f_1(\alpha')\geq 0, \textrm{ for all }\alpha'\leq\alpha\right\},
\end{align*}
with
\begin{align*}f_1(\alpha)=\textrm{min}\left(S^{i,(l)}(\alpha)\Lambda^{i,(l)}(\alpha)\hat{e}\right)-\bar{\tau}_1^i\gamma^{i,(l)}\left(s^{i,(l)}(\alpha)\right)^T\lambda^{i,(l)}(\alpha)/m_i,\end{align*}
we see that
$$\alpha^{i}_1\geq (1-\bar{\tau}_1^i\gamma^{i})\mu\frac{1}{M^i_1}>0.$$
We have that $\mu$ is bounded away from zero (since both $\sigma$ and $(s^i)^T\lambda^i$ are bounded away from zero) and $\bar{\tau}_1^i\gamma^{i}$ is bounded away from one (since  $\bar{\tau}_1^i$ is at most one and $\gamma^{i}$ is bounded away from one). Hence, $\alpha^{i}_1$ is bounded away from zero in in $\Omega^i(\epsilon)$ for $i=1,\dots,N$.

To determine $\alpha^{i}_2$ at iteration $l$, we consider the expression (where we once again omit the superscript $l$)
$$(s^{i}(\alpha))^T\lambda^{i}(\alpha)-\bar{\tau}_2^i\gamma^{i}\|R^i(\alpha)\|.$$
In analogy to page 117 in \cite{bel:98}, we use the mean value theorem and the Lipschitz continuity of $(R^{i})'$. The mean value theorem gives
\begin{align*}
R^i(\alpha)&=R^i(z^{i})+\alpha (R^{i})'(z^{i})\Delta z^{i}+\alpha\left(\int_{0}^1 \left(  (R^{i})'(z^{i}+t\alpha\Delta z^{i})-(R^{i})'(z^{i}) \right)\Delta z^{i}\right)\\
&=R^i(z^{i})(1-\alpha)+\alpha\hat{r}^i_1+\alpha\left(\int_{0}^1 \left(  (R^{i})'(z^{i}+t\alpha\Delta z^{i})-(R^{i})'(z^{i}) \right)\Delta z^{i}\right),
\end{align*}
where $\hat{r}^i$ denotes the contributing elements of the $i$th term in $\|\hat{r}\|$, see \eqref{residualhat}.The Lipschitz continuity of $(R^{i})'$ gives
\begin{align}
\label{LipCon}
\begin{split}
\|R^i(\alpha)\|&\leq \|R^i(z^{i})\|(1-\alpha)+\alpha\|\hat{r}^i_1\|+\alpha^2L^i\|\Delta z^{i}\|^2.
\end{split}
\end{align}
From inequality \eqref{LipCon} and the stop criteria of the step direction calculation Algorithm~\ref{alg:algNPD}, we get
\begin{align*}
\|R^i(\alpha)\|\leq(1-\alpha)\|R^i(z^i)\|+\alpha\hat{\eta}\frac{(s^i)^T\lambda^i}{m}+\alpha^2L^i\|\Delta z^{i}\|^2.
\end{align*}
Thus,
\begin{align*}
(s^i(\alpha))^T\lambda^i(\alpha)-\bar{\tau}_2^i\gamma^i\|R^i(\alpha)\|&\geq\\
(\!s^{i}\!+\alpha \Delta s^{i})^T(\lambda^{i}\!+\alpha \Delta\lambda^{i})\!-\!\bar{\tau}_2^i\gamma^i\left(\!(1-\alpha)\|R^i(z^i)\|\!+\!\alpha\hat{\eta}\frac{(s^i)^T\lambda^i}{m}\!+\!\alpha^2L^i\|\Delta z^{i}\|^2\right)\!&=\\
(s^{i})^T\lambda^{i}+\alpha(s^{i})^T\Delta\lambda^{i}+\alpha(\Delta s^{i})^T\lambda^{i}+\alpha^2(\Delta s^{i})^T\Delta\lambda^{i}&-\\
\bar{\tau}_2^i\gamma^i\left((1-\alpha)\|R^i(z^i)\|+\alpha\hat{\eta}\frac{(s^i)^T\lambda^i}{m}+\alpha^2L^i\|\Delta z^{i}\|^2\right)&=\\
(1-\alpha)\left((s^i)^T\lambda^i-\bar{\tau}_2^i\gamma^i\|R^i(z^i)\|\right)+\alpha \left(m_i\mu-\frac{\bar{\tau}_2^i\gamma^i\hat{\eta}(s^i)^T\lambda^i}{m}\right)&+\\
\alpha^2\left((\Delta s^i)^T\Delta \lambda^i-\tau_2^i\gamma^iL^i\|\Delta z^{i,(l)}\|^2\right)&=\\
(1-\alpha)\left((s^i)^T\lambda^i-\bar{\tau}_2^i\gamma^i\|R^i(z^i)\|\right)+\alpha \left(m_i\mu-\frac{\bar{\tau}_2^i\gamma^i\hat{\eta}(s^i)^T\lambda^i}{m}\right)&+\\
\alpha^2\left((\Delta s^i)^T\Delta \lambda^i-\tau_2^i\gamma^iL^i\|\Delta z^{i,(l)}\|^2\right)&\geq\\
(1-\alpha)\left((s^i)^T\lambda^i-\bar{\tau}_2^i\gamma^i\|R^i(z^i)\|\right)+\alpha \left(\mu-\frac{\bar{\tau}_2^i\gamma^i\hat{\eta}(s^i)^T\lambda^i}{m}\right)&-\\
\alpha^2\left|(\Delta s^i)^T\Delta \lambda^i-\tau_2^i\gamma^iL^i\|\Delta z^{i,(l)}\|^2\right|&\geq\\
(1-\alpha)\!\left((s^i)^T\!\lambda^i\!\!-\!\bar{\tau}_2^i\gamma^i\|R^i(z^i)\|\right)\!+\!\alpha\! \left(\!\frac{\sigma}{m}\underset{i}{\textrm{min}}\left((s^i)^T\!\lambda^i\right)\!-\!\frac{\bar{\tau}_2^i\gamma^i\hat{\eta}(s^i)^T\!\lambda^i}{m}\right)\!-\!\alpha^2M_2^i&\geq\\
\alpha \left(\frac{\sigma}{m}\underset{i}{\textrm{min}}\left((s^i)^T\lambda^i\right)-\frac{\bar{\tau}_2^i\gamma^i\hat{\eta}(s^i)^T\lambda^i}{m}\right)-\alpha^2M_2^i.
\end{align*}
Since $\alpha_2^i$ is defined as
\begin{align*}
\alpha_2^i=\underset{\alpha \in [0,1]}{\textrm{max}}\left\{\alpha|f_2(\alpha')\geq 0, \textrm{ for all }\alpha'\leq\alpha\right\},
\end{align*}
with
$$f_2(\alpha)=\left(s^{i,(l)}(\alpha)\right)^T\lambda^{i,(l)}(\alpha)-\bar{\tau}_2^i\gamma^{i,(l)}\|R^i(\alpha)\|$$
and $\sigma^i>\bar{\tau}_2^i\gamma^{i}\hat{\eta}^i(s^i)^T\lambda^i/\underset{i}{\textrm{min}}\left((s^i)^T\lambda^i\right)+\epsilon_{\sigma}$,
we see that
$$\alpha^{i}_2\geq  \left(\sigma\underset{i}{\textrm{min}}\left((s^i)^T\lambda^i\right)-\bar{\tau}_2^i\gamma^i\hat{\eta}(s^i)^T\lambda^i\right)\frac{1}{mM_2^i}>0.$$
Hence, $\alpha^{i}_2$ is bounded away from zero in $\Omega^i(\epsilon)$ for $i=1,\dots,N$.
\end{proof}
Thus, we can always find a step size bounded away from zero. We continue by showing that, given a step direction $\tilde{p}^{i}$ calculated in accordance to Algorithm \ref{alg:algNPD}, there is a step $p^{i}$ such that 
\begin{align*}
\| H^i(z^{i}+p^{i}) \|\leq (1-\beta(1-\eta^{i}))\| H^i(z^{i}) \|.
\end{align*}
given $z^{i}$ and $\beta \in (0,1)$. That is, there exists an actual step size $\alpha$ such that the inequality above is fulfilled. To do this we state Lemma \eqref{breakdown} which is based on Lemma 3.1 in \cite{eis:94}.

\begin{lemma}\label{breakdown}
Given $z^{i}$ and $\beta \in (0,1)$, assume that $\tilde{p}^{i}$ is calculated in accordance to Algorithm \ref{alg:algNPD}. Then
\begin{align*}
\| H^i(z^{i})+(H^{i})'(z^{i})\tilde{p}^{i} \|<  \| H^i(z^{i}) \|,
\end{align*}
and there exists an $\eta_{\textrm{min}}^i\in [0,1)$ such that, for any $\eta^{i} \in [\eta_{\textrm{min}}^i,1)$, we can find a $p^{i}$ satisfying
\begin{align*}
\| H^i(z^{i}+p^{i}) \|\leq (1-\beta(1-\eta^{i}))\| H^i(z^{i}) \|.
\end{align*}
\end{lemma}
\begin{proof}
Constraint $\| H^i(z^{i})+(H^{i})'(z^{i})\tilde{p}^{i} \|< \| H^i(z^{i}) \|$ implies that $\|H^i(z^{i})\|\neq0$ and $\tilde{p}^{i}\neq0$. However, note that if $\|H^i(z^{i})\|=0$ we would have terminated the algorithm in the previous iteration. We start by showing the first inequality. The step direction calculations are terminated when the $i$th term of $\|\hat{r}\|$ (denoted $\|\hat{r}^i\|$) is less than $\hat{\eta}(s^i)^T\lambda^i/m$, see \eqref{residualhat}, \eqref{epsADMMInexact1} and \eqref{epsADMMInexact2}. Consequently, we get
\begin{align*}
\| H^i(z^{i})+(H^{i})'(z^{i})\tilde{p}^{i} \| &\leq \left\|\mu\begin{bmatrix}0 \\ 0\\ 0  \\\hat{e} \end{bmatrix}+\hat{r}^i\right\|\\
&\leq \sqrt{m_i}\mu+\left\|\hat{r}^i\right\|\\
&\leq\frac{\sqrt{m_i}\sigma}{m}\underset{i}{\textrm{min}}\left( (s^i)^T\lambda^i\right)+ \frac{\hat{\eta}}{m}\left( (s^i)^T\lambda^i\right)\\
&\leq \left(\sigma+ \hat{\eta}\right)\left( (s^i)^T\lambda^i\right)\\
&\leq\left(\sigma+ \hat{\eta}\right)\|H^i(z^i)\|\\
&=\bar{\eta}\|H^i(z^i)\|\\
&<\|H^i(z^i)\|,
\end{align*}
 since $\bar{\eta}<1$. The rest of the proof follows that of Lemma 3.1 in \cite{eis:94}. A slight modification has been made to incorporate $i=1,\dots,N$. Define
\begin{align*}
\tilde{\eta} &= \frac{\| H^i(z^{i})+(H^{i})'(z^{i})\tilde{p}^{i} \|}{\| H^i(z^{i})\|},\\
\epsilon ^i &=\frac{(1-\beta)(1-\tilde{\eta}^{i})\|H^i(z^{i})\|}{\|\tilde{p}^{i}\|},\\
\eta_{\textrm{min}}^i&=\textrm{max}\left\{ \tilde{\eta}, 1-\frac{(1-\tilde{\eta})\delta^i}{\|\tilde{p}^{i}\|}\right\},
\end{align*}
where $\delta ^i>0$ is chosen such that
\begin{align*}
\| H^i(z^{i}+p^{i})-H^i(z^{i})-(H^{i})'(z^{i})p^{i} \|\leq \epsilon^i \| p^{i} \|
\end{align*}
 whenever $\| p^{i} \|\leq \delta^i$. We set
\begin{align*}
 p^{i}= \frac{1-\eta^{i}}{1-\tilde{\eta}}  \tilde{p}^{i},
\end{align*}
for any $\eta^{i} \in [\eta_{\textrm{min}}^i,1)$.  Then
\begin{align*}
\| H^i(z^{i})+(H^{i})'(z^{i})p^{i} \|&\leq \frac{\eta^{i}-\tilde{\eta}}{1-\tilde{\eta}}  \| H^i(z^{i})\|+\frac{1-\eta^{i}}{1-\tilde{\eta}} \| H^i(z^{i})+(H^{i})'(z^{i})\tilde{p}^{i} \|\\
&=\frac{\eta^{i}-\tilde{\eta}}{1-\tilde{\eta}}  \| H^i(z^{i})\|+\frac{1-\eta^{i}}{1-\tilde{\eta}}\tilde{\eta} \| H^i(z^{i})\|\\
&=\eta^{i}\| H^i(z^{i})\|,
\end{align*}
and, due to
\begin{align*}
\| p^{i}\|= \frac{1-\eta^{i}}{1-\tilde{\eta}} \| \tilde{p}^{i}\|\leq \frac{1-\eta_{\textrm{min}}^i}{1-\tilde{\eta}} \| \tilde{p}^{i}\|\leq \delta^i,
\end{align*}
we get
\begin{align*}
\| H^i(z^{i}+p^{i}) \|&=\| H^i(z^{i}+p^{i})-H^i(z^{i})-(H^{i})'(z^{i})p^{i} +H^i(z^{i})+(H^{i})'(z^{i})p^{i} \|\\
& \leq\| H^i(z^{i}+p^{i})-H^i(z^{i})-(H^{i})'(z^{i})p^{i}\|+\|H^i(z^{i})+(H^{i})'(z^{i})p^{i} \|\\
&\leq \| H^i(z^{i}+p^{i})-H^i(z^{i})-(H^{i})'(z^{i})p^{i} \|\\
&\leq\epsilon^i\frac{1-\eta^{i}}{1-\tilde{\eta}}\| \tilde{p}^{i}\|+\eta^{i}\| H^i(z^{i})\|\\
&=(1-\beta )(1-\eta^{i})\| H^i(z^{i})\|+\eta^{i}\| H^i(z^{i})\|\\
&=(1-\beta(1-\eta^{i}))\| H^i(z^{i})\|.
\end{align*}
\end{proof}
Recall that the line search in Algorithm \ref{alg:algNPD} is initialized with $\tilde{p}^{i,(l)}=\alpha^{i,(l)}\Delta z^{i,(l)}$ and $\eta^{i,(l)}=1-\alpha^{i,(l)}(1-\bar{\eta}^{(l)})$, and evolves as $\tilde{p}^{i,(l)}=\theta\tilde{p}^{i,(l)}$ and $\eta^{i,(l)}=1-\theta(1-\eta^{i,(l)})$. Assume that $H^i(z^{i,(l)})\neq 0$ and $\|\Delta z^{i,(l)}\|\neq0$. Consequently, from the proof of Lemma~\ref{breakdown}, we conclude that the while loop terminates with
\begin{align}\label{terminate}1-\eta^{i,(l)}\geq\textrm{min}\left\{ \alpha^{i,(l)}(1-\bar{\eta}^{(l)}), \frac{\theta(1-\bar{\eta}^{(l)})\delta^{i}}{\|\Delta z^{i,(l)}\|}\right\},\end{align}
which is equivalent to the result on page 114 in \cite{bel:98} for $N=1$ and Lemma 5.1 in \cite{eis:94} for $N=1$ and an additional $\alpha$-update. To see this, we follow the reasoning made in the proof of Lemma 5.1 in \cite{eis:94}. If $\eta^{i,(l)} \in [\eta_{\textrm{min}}^{i,(l)},1)$ such that $1-\eta^{i,(l)}<\delta^i(1-\bar{\eta})/\|\Delta z^{i,(l)}\|$ we get, from the proof above, that
\begin{align*}
\| H^i(z^{i}+p^{i}) \|\leq (1-\beta(1-\eta^{i}))\| H^i(z^{i}) \|.
\end{align*}
Also, note that $1-\eta^{i,(l)}$ decreases with a factor $\theta\in(0,1)$ for each run of the while loop. If no iteration of the while loop is necessary we have $1-\eta^{i,(l)}=\alpha^{i,(l)}(1-\bar{\eta}^{(l)})$. Suppose instead that $1-\eta^{i,(l)}=\theta\delta^i(1-\bar{\eta})/\|\Delta z^{i,(l)}\|$ which is less than $\delta^i(1-\bar{\eta})/\|\Delta z^{i,(l)}\|$, then the loop terminates. Thus, \eqref{terminate} holds.

The value of $1-\eta^{i,(l)}$ is bounded away from zero (since $\alpha^{i,(l)}$ is bounded away from zero, $\bar{\eta}^{(l)}$ is bounded away from one, $\theta>0$ and independent of $l$, $\delta^i>0$ and independent of $l$, and $\|\Delta z^{i,(l)}\|$ is bounded). The fact that $1-\eta^{i,(l)}$ is bounded away from zero will be used in the proof of global convergence. So at this point we observe that the combination of Lemma \ref{breakdown} and Theorem \ref{notZero}, and the complementary discussion assure that the iterates are persistently updated  through the run of the algorithm. As a corollary, we can then state that the algorithm, provided assumptions B1-B4 are fulfilled, can only break down at some iteration point $z^{i,(l)}$ if and only if $\|H^i(z^{i,(l)})\|=0$. Notice that for convex problems the break down can only happen when we have arrived at an optimal solution, hence there will always exist a suitable search direction for updating the iterates. We will now focus on convergence properties of the algorithm.

\subsection{Convergence Properties}
We first discuss some results needed for the proof of global convergence. The following theorem is used in the proof of Theorem \ref{conLimitPointTheorem}, and is therefor included here. It is based on Theorem 3.5 in \cite{eis:94}.

\begin{theorem}\label{conLimitPointTheoremPre}
Assume that Algorithm \ref{alg:algNPD} does not break down. If $z_*$ is a limit point of $\{z^{(l)}\}$ such that there exists a $\Gamma$ independent of $l$ for which
\begin{align}\label{conLimitPoint}
\|p^{(l)} \|\leq\Gamma(1-\eta^{(l)}_*)\|H(z^{(l)})\|
\end{align}
when $z^{(l)}$ is sufficiently close to $z_*$ and $l$ is sufficiently large, then $z^{(l)}\rightarrow z_*$. Here, $(1-\eta^{(l)}_*)$ corresponds to the actual step size used by all subproblems $i=1,\dots, N$.
\end{theorem}
\begin{proof}
See the proof of Theorem 3.5 in \cite{eis:94}.
\end{proof}

We are now ready to state the following theorem, which is based on Theorem 3.1 in \cite{bel:98}.
\begin{theorem}\label{conLimitPointTheorem}
If $z_*$ is a limit point of $\{z^{(l)}\}$ such that $\|H'(z_*)\|$ is nonsingular, then the sequence $\{z^{(l)}\}$ generated by Algorithm \ref{alg:algNPD} converges to $z_*$ .
\end{theorem}
\begin{proof}
The proof follows closely that of Theorem 3.1 in \cite{bel:98}. We define $K=\|(H'(z_*))^{-1}\|$ and choose $\delta>0$ such that $(H'(z))^{-1}$ exists and $\|(H'(z))^{-1}\|\leq2K$ whenever $z\in N_{\delta}(z_*)$. The actual step used in subproblem $i$ is
\begin{align*}
 p^{i,(l)}= \underset{i}{min}\left\{1-\eta^{i,(l)}\right\}\frac{1}{1-\bar{\eta}^{(l)}}  \Delta z^{i}=\frac{1-\eta^{(l)}_*}{1-\bar{\eta}^{(l)}} \Delta z^{i},
\end{align*}
that is, we choose the minimum step size over all subproblems. Assume $z^{(l)}\in N_{\delta}(z_*)$, then
\begin{align*}
\| p^{(l)}\|&=\frac{1-\eta^{(l)}_*}{1-\bar{\eta}^{(l)}} \|\Delta z\|\\
&\leq \frac{1-\eta^{(l)}_*}{1-\bar{\eta}^{(l)}}\|(H'(z^{(l)}))^{-1}\|\|-H(z^{(l)})+\tilde{r}^{(l)}\|\\
&\leq \frac{1-\eta^{(l)}_*}{1-\bar{\eta}^{(l)}}\|(H'(z^{(l)}))^{-1}\|(\|H(z^{(l)})\|+\bar{\eta}^{(l)}\|H(z^{(l)})\|)\\
&\leq \frac{1-\eta^{(l)}_*}{1-\bar{\eta}^{(l)}}2K(1+\bar{\eta}^{(l)})\|H(z^{(l)})\|\\
&\leq \Gamma(1-\eta^{(l)}_*)\|H(z^{(l)})\|\\
\end{align*}
where $$\Gamma=2K\frac{1+\bar{\eta}_{\textrm{max}}}{1-\bar{\eta}_{\textrm{max}}}$$
and $\bar{\eta}_{\textrm{max}}$ is the maximum value that $\bar{\eta}^{(l)}$ can be set to. Thus, there exists a $\Gamma$ independent of $l$ for which
inequality \eqref{conLimitPoint} holds when $z^{(l)}$ is sufficiently close to $z_*$ and $l$ is sufficiently large. Consequently, by Theorem \ref{conLimitPointTheoremPre}, $z^{(l)}\rightarrow z_*$.
\end{proof}

The following theorem is needed in the proof of global convergence. It is based on Theorem 3.4 in \cite{eis:94}.

\begin{theorem}\label{divergent}
Assume that Algorithm \ref{alg:algNPD} does not break down. If $\sum_{l\geq 0}(1-\eta^{(l)}_*)$ is divergent then $H(z^{(l)}) \rightarrow 0$.
\end{theorem}
\begin{proof}
The proof follows closely that of Theorem 3.4 in \cite{eis:94}. From the line search in Algorithm \ref{alg:algNPD}, we get
\begin{align*}
\|H^i(z^{i,(l)})\|&\leq(1-\beta(1-\eta^{(l-1)}_*))\|H^i(z^{i,(l-1)})\|\\
&\leq \Pi_{0\leq j <l}(1-\beta(1-\eta^{(j)}_*))\|H^i(z^{i,(0)})\|\\
&\leq \textrm{exp}\left(-\beta\sum_{0\leq j <l}(1-\eta^{(j)}_*)\right)\|H^i(z^{i,(0)})\|.
\end{align*}
Thus
\begin{align*}
\|H(z^{(l)})\|\leq \textrm{exp}\left(-\beta\sum_{0\leq j <l}(1-\eta^{(j)}_*)\right)\|H(z^{(0)})\|.
\end{align*}
The divergence of $\sum_{l\geq 0}(1-\eta^{{(l)}_*})$ implies $H(z^{(l)}) \rightarrow 0$ since $\beta>0$ and $1-\eta^{(l)}_*\geq 0$.
\end{proof}

We are now ready to state the theorem of global convergence of the proposed method. It is based on Theorem 3.3 in \cite{bel:98}.

\begin{theorem}\label{globalCon}
Assume $\{z^{(l)}\}$ is generated by Algorithm \ref{alg:algNPD} and assumptions B1-B4 are fulfilled, then $\{\|H(z^{(l)})\|\}$ converges to zero.
\end{theorem}
\begin{proof}
The proof follows closely that of Theorem 3.3 in \cite{bel:98} and Theorem 5.2 in \cite{eis:94}. Algorithm \ref{alg:algNPD} does not break down (unless at the optimal solution) and has step sizes bounded away from zero, see Theorem \ref{notZero}, Lemma \ref{breakdown} and the complementary discussion. Furthermore, it follows from Theorem \ref{conLimitPointTheorem} that $z^{(l)}\rightarrow z_*$. 

The sequence $\{\|H(z^{(l)})\|\}$ is, by construction of the algorithm, decreasing and bounded, and consequently convergent. Assume that it converges to $\kappa >0$. So, for sufficiently large $l$ we have $z^{i,(l)}\in N_{\delta^i}(z_*)$ and the while loop of Algorithm \ref{alg:algNPD} terminates with inequality \eqref{terminate} fulfilled. Since $z^{(l)}\rightarrow z_*$, all but finitely many $l$ satisfy that $z^{i,(l)}\in N_{\delta^i}(z_*)$,  consequently $\sum_{l\geq 0}(1-\eta^{(l)}_*)$ is divergent. Hence, from Theorem \ref{divergent} we get that $H(z^{(l)}) \rightarrow 0$, which contradicts our assumption that $H(z^{(l)}) \rightarrow \kappa >0$. Thus, $\{\|H(z^{(l)})\|\}$ converges to zero.
\end{proof}




\section{ADMM, Fixed Point Iterations and Uzawa's Method}\label{sec:Appendix2}
As stated in Section \ref{sec:UMAFPI}, ADMM can be viewed as fixed point iterations and as a modified version of Uzawa's method. In this appendix, we explore these relations in detail for our specific problem formulation.
%=================================================================
\subsection{ADMM and Fixed Point Iterations}
%=================================================================
Consider functions $F_1(\Delta W)$ and $F_2(\Delta x)$ in \eqref{eq:ADMM}, and let us rewrite them in the form \[F_1(\Delta W)= \frac{1}{2}\Delta W^T\tilde{F}_1\Delta W+\tilde{f}_1^T\Delta W \text{ and }  F_2(\Delta x)=\tilde{f}_2^T\Delta x.\]
We can then rewrite the optimality conditions for \eqref{eq:ADMM}, in \eqref{sysOfEq}, as
%
\begin{equation}\label{kkt}
\underbrace{\begin{bmatrix} \tilde{F}_1 & 0 & \rho A^T  \\ 0& 0& \rho B^T\\\rho A & \rho B & 0 \end{bmatrix}}_{A_{\text{KKT}}}
\begin{bmatrix} \Delta W  \\ \Delta x \\ \Delta \bar{u}\end{bmatrix} =
\underbrace{\begin{bmatrix}-\tilde{f}_1  \\ -\tilde{f}_2\\\rho c\end{bmatrix}}_{b_{\text{KKT}}},
\end{equation}
%
where $\Delta \bar{u} = (\Delta \bar{v},\Delta \bar{v}_c)$. With the newly defined notation, the ADMM iterations for the problem in \eqref{eq:ADMM} can be written in closed form as
%
\begin{equation}\label{eq:ADMMIter}
\begin{split}
\Delta W^{(k+1)} &= M_1(-\tilde{f}_1+\rho A^Tc -\rho A^TB\Delta x^{(k)}-\rho A^T \Delta \bar{u}^{(k)}),\\
\Delta x^{(k+1)} &=M_2(-\tilde{f}_2+\!\rho B^Tc\! -\!\rho B^T\!A\Delta W^{(k+1)}\!-\!\rho B^T\! \Delta \bar{u}^{(k)}),\\
\Delta \bar{u}^{(k+1)} &=  \Delta \bar{u}^{(k)} + A \Delta W^{(k+1)} + B \Delta x^{(k+1)} - c,
\end{split}
\end{equation}
%
with
%
\begin{equation*}
M_1=(\tilde{F}_1+\rho A^TA)^{-1} \text{ and } M_2=(\rho B^TB)^{-1}.
\end{equation*}
%
We can rewrite the equations in \eqref{eq:ADMMIter} in a more compact manner as
%
\begin{equation}\label{eq:ADMMIterSS}
\begin{split}
 \begin{bmatrix} \Delta W^{(k+1)}   \\ \Delta x^{(k+1)}  \\  \Delta \bar{u}^{(k+1)} \end{bmatrix}=G \begin{bmatrix} \Delta W^{(k)}   \\ \Delta x^{(k)} \\  \Delta \bar{u}^{(k)}\end{bmatrix}+f,
\end{split}
\end{equation}
%
with
%
\begin{equation*}%\label{eq:Gmatrix}
G\!=\!
\begin{bmatrix}\begin{array}{c c c}
0 & -\rho M_1 A^TB &-\rho M_1 A^T\\
0 & \rho^2M_2 B^T\!AM_1A^TB& \rho M_2 B^T\!(\rho AM_1A^T\!- I) \\
0&\rho  (\rho B M_2 B^T\!\!-I)AM_1A^T\!B  & \rho (\rho B M_2 B^T\!\!-I)AM_1 A^T\!-\rho B M_2 B^T\!\!+\!I\!\end{array}\end{bmatrix}
\end{equation*}
%
and
%
\begin{equation*}%\label{eq:fvector}
\begin{split}
f=
\begin{bmatrix}M_1m_1\\
M_2m_2-\rho M_2B^TAM_1m_1\\
-c+AM_1m_1+B(M_2m_2-\rho M_2 B^TAM_1m_1)\\ \end{bmatrix}\!,
\end{split}
\end{equation*}
where
\begin{equation*}
m_1=-\tilde{f}_1+\rho A^Tc \text{ and } m_2=-\tilde{f}_2+\rho B^Tc.
\end{equation*}
%
The iterations in \eqref{eq:ADMMIterSS} clearly show that $\Delta x$ and $\Delta \bar{u}$ make up the state of the algorithm, whereas $\Delta W$ only is an intermediate result, see \cite{boyd:11}. Notice that we can view \eqref{eq:ADMMIterSS} as an iterative solver for a pre-conditioned version of the system of equations in \eqref{kkt}. That is, the iteration matrix $G$ and the vector $f$ can be expressed as
%
\begin{equation*}
G=I-M_{\text{PRE,1}}^{-1}A_{\text{KKT}}\textrm{ and }f = M_{\text{PRE,1}}^{-1}b_{\text{KKT}},
\end{equation*}
%
respectively, where $M_{\text{PRE,1}}$ is a pre-conditioner defined as
%
\begin{equation*}%\label{eq:precond}
\begin{split}
M_{\text{PRE,1}}=
\begin{bmatrix} \tilde{F}_1&-\rho A^TB&\rho A^T \\ 0&0&\rho B^T \\\rho A&\rho B &-\rho I\\ \end{bmatrix}.
\end{split}
\end{equation*}
%
For details regarding iterative solvers and pre-conditioners, see \cite{saa:03}. When ADMM converges, we get the fixed point iterations
%
\begin{equation*}%\label{eq:ADMMIterSSfixed}
\begin{split}
 \begin{bmatrix} \Delta W \\ \Delta x \\\Delta \bar{u} \end{bmatrix}=G \begin{bmatrix} \Delta W  \\\Delta  x \\ \Delta \bar{u}\end{bmatrix}+f \Leftrightarrow
A_{\text{KKT}}\begin{bmatrix} \Delta W  \\ \Delta x \\ \Delta \bar{u}\end{bmatrix}=b_{\text{KKT}},
\end{split}
\end{equation*}
%
which is equal to the system of equations that we would like to solve, namely \eqref{kkt}.

%=================================================================
\subsection{ADMM and Uzawa's Method}\label{admom}
%=================================================================

Notice that solving \eqref{kkt} is equivalent to finding a saddle point of the Lagrangian function
%
\begin{align*}%\label{lagrangian}
\begin{split}
\mathcal{L}(\Delta W,\Delta x,\Delta \bar{u}) =  F_1(\Delta W)+F_2(\Delta x)+\rho \Delta \bar u^T(AW+Bx-c),
\end{split}
\end{align*}
%
which also is a saddle point of the augmented Lagrangian
%
\begin{align*}%\label{auglagrangian}
\begin{split}
\mathcal{L}_{\rho}(\Delta W,\Delta x,\Delta \bar{u}) =   F_1(\Delta W)+F_2(\Delta x)+\frac{\rho}{2}\|A\Delta W+B\Delta x-c+\Delta \bar{u}\|_2^2.
\end{split}
\end{align*}
%
This augmented Lagrangian function is in turn the Lagrangian function of the optimization problem
%
\begin{equation}\label{eq:ADMMmod}
\begin{split}
\minimize_{\Delta S, \Delta x}  & \quad  F_1 (\Delta W) + F_2(\Delta x)+\frac{\rho}{2}\|A\Delta W+B\Delta x-c\|_2^2, \\
\subject & \quad  A \Delta W + B\Delta x = c,
\end{split}
\end{equation}
%
which is equivalent to \eqref{eq:ADMM}. ADMM applied to \eqref{eq:ADMM} is equivalent to Uzawa's method applied to \eqref{eq:ADMMmod} with one Gauss-Seidel iteration \citep{saa:03} in the update of the primal variables \citep{mar:75,gab:76} and the relaxation parameter equal to the penalty parameter $\rho$. This can be seen by first noting that the optimality conditions of \eqref{eq:ADMMmod}, i.e.,
%
\begin{equation}\label{augkkt}
\underbrace{\begin{bmatrix} M_1^{-1} & \rho A^TB & \rho A^T  \\ \rho B^TA& M_2^{-1}& \rho B^T\\\rho A & \rho B & 0 \end{bmatrix}}_{A_{\text{KKT}}^{\text{PRE,2}}}
\begin{bmatrix} \Delta W  \\ \Delta x \\ \Delta \bar{u}\end{bmatrix} =
\underbrace{\begin{bmatrix}m_1  \\ m_2\\\rho c\end{bmatrix}}_{b_{\text{KKT}}^{\text{PRE,2}}},
\end{equation}
%
are equivalent to those of Problem \eqref{eq:ADMM}. In fact, the system of equations \eqref{augkkt} is a preconditioned version of \eqref{kkt}, that is, $A_{\text{KKT}}^{\text{PRE,2}}=M_{\text{PRE,2}}^{-1}A_{\text{KKT}}$ and $b_{\text{KKT}}^{\text{PRE,2}}=M_{\text{PRE,2}}^{-1}b_{\text{KKT}}$,
with
%
\begin{align*}%\label{precond1}
M_{\text{PRE,2}} = \begin{bmatrix} I & 0 & -A^T  \\ 0& I& -B^T\\0 & 0 & I \end{bmatrix}.
\end{align*}
%
Uzawa's method minimizes the Lagrangian function by iteratively first minimizing with respect to primal variables and then with respect to dual variables. In ADMM, however, the minimization with respect to primal variables are performed sequentially over $\Delta W$ and $\Delta x$, that is, by performing one Gauss-Seidel iteration with respect to $\Delta W$ and $\Delta x$  \citep{saa:03}. The first step of Uzawa's method applied to our problem then requires solving
%
\begin{equation}\label{1uzawa}
\begin{bmatrix} M_1^{-1} & \rho A^TB\\ \rho B^TA& M_2^{-1}\end{bmatrix}
\begin{bmatrix} \Delta W  \\ \Delta x\end{bmatrix} =
\begin{bmatrix}m_1-\rho A^T\Delta \bar{u}^k  \\ m_2-\rho B^T\Delta \bar{u}^k\end{bmatrix},
\end{equation}
%
which is the closed form solution of
%
\begin{align*}
(\Delta W^{(k+1)},\Delta x^{(k+1)}) =& \argmin_{\Delta W,\Delta x} \left\{L_{\rho}(\Delta W,\Delta x,\Delta \bar{u}^k)  \right\}\\
=&\argmin_{\Delta W,\Delta x} \left\{\!\! F_1(\Delta W)\!\!+\!\!F_2(\Delta x)\!\!+\!\!\frac{\rho}{2}\|A\Delta W\!\!+\!\!B\Delta x\!\!-c\!\!+\!\!\Delta \bar{u}^{k}\|_2^2\right\}\!\!.
\end{align*}
%
The system of equations in \eqref{1uzawa} can be solved approximately using Gauss-Seidel, that is,
%
\begin{align}\label{gaussseidel}
\begin{split}
\begin{bmatrix} \Delta W^{k+1}  \\ \Delta x^{k+1}\end{bmatrix} =&
\begin{bmatrix} M_1^{-1} & 0\\ \rho B^T\!\!A& M_2^{-1}\end{bmatrix}^{\!-1}
\begin{bmatrix} 0 & -\rho A^TB\\ 0& 0\end{bmatrix}\!\!
\begin{bmatrix} \Delta W^{k}  \\ \Delta x^{k}\end{bmatrix}\!\!
+\\&
\begin{bmatrix} M_1^{-1} & 0\\ \rho B^T\!\!A& M_2^{-1}\end{bmatrix}^{\!-1}
\begin{bmatrix}m_1-\rho A^T\Delta \bar{u}^k  \\ m_2-\rho B^T\Delta \bar{u}^k\end{bmatrix}\\=&
\begin{bmatrix}0 & -\rho M_1 A^TB  \\
0 & \rho^2M_2 B^T\!\!AM_1 A^TB  \end{bmatrix}\!\!
\begin{bmatrix} \Delta W^{k}  \\ \Delta x^{k}\end{bmatrix}
+\\&
\begin{bmatrix}M_1(m_1-\rho A^T\Delta \bar{u}^k)\\
M_2(m_2\!-\!\rho B^T\!\Delta \bar{u}^k)-\rho M_2B^T\!\!AM_1(m_1\!-\!\rho A^T\!\Delta \bar{u}^k)\end{bmatrix}.
\end{split}
\end{align}
%
Notice that the iteration in \eqref{gaussseidel} is equivalent to the primal updates \eqref{eq:ADMMIterSS} in ADMM. The fact that the dual update in Uzawa's method is equivalent to that of ADMM then shows the equivalence. \\
%



\end{document}

1406.2192.extracted_bib
[{Al-Jeiroudi and Gondzio(2009)]{al:09 Al-Jeiroudi G, Gondzio J (2009) Convergence analysis of the inexact infeasible interior-point method for linear optimization. Journal of Optimization Theory and Applications 141(2):231--247
1406.2192.extracted_bib
[{Arrow(1964)]{arr:64 Arrow KJ (1964) Studies in linear and non-linear programming. Stanford University Press Stanford, Calif.
1406.2192.extracted_bib
[{Bellavia(1998)]{bel:98 Bellavia S (1998) Inexact interior-point method. Journal of Optimization Theory and Applications 96(1):109--121
1406.2192.extracted_bib
[{Bellavia and Pieraccini(2004)]{bel:04 Bellavia S, Pieraccini S (2004) Convergence analysis of an inexact infeasible interior point method for semidefinite programming. Computational Optimization and Applications 29(3):289--313
1406.2192.extracted_bib
[{Benzi et~al(2005)Benzi, Golub, and Liesen]{ben:05 Benzi M, Golub GH, Liesen J (2005) {Numerical solution of saddle point problems. Acta Numerica 14(-1):1--137
1406.2192.extracted_bib
[{Bertsekas and Tsitsiklis(1997)]{ber:97 Bertsekas DP, Tsitsiklis JN (1997) Parallel and Distributed Computation: Numerical Methods. Athena Scientific
1406.2192.extracted_bib
[{Bonettini and Ruggiero(2007)]{bon:07 Bonettini S, Ruggiero V (2007) Some iterative methods for the solution of a symmetric indefinite {KKT system. Computational Optimization and Applications 38(1):3--25
1406.2192.extracted_bib
[{Bonettini et~al(2005)Bonettini, Galligani, and Ruggiero]{bon:05 Bonettini S, Galligani E, Ruggiero V (2005) An inexact {Newton method combined with {Hestenes multipliers scheme for the solution of {Karush-Kuhn-Tucker systems. Applied Mathematics and Computation 168(1):651--676
1406.2192.extracted_bib
[{Boyd and Vandenberghe(2004)]{boyd:04 Boyd S, Vandenberghe L (2004) Convex Optimization. Cambridge University Press
1406.2192.extracted_bib
[{Boyd et~al(2011)Boyd, Parikh, Chu, Peleato, and Eckstein]{boyd:11 Boyd S, Parikh N, Chu E, Peleato B, Eckstein J (2011) Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning 3(1):1--122
1406.2192.extracted_bib
[{{Chen et~al(2013){Chen, {He, {Ye, and {Yuan]{che:13 {Chen C, {He B, {Ye Y, {Yuan X (2013) {The Direct Extension of ADMM for Multi-block Convex Minimization Problems is Not Necessarily Convergent. Optimization Online e-prints
1406.2192.extracted_bib
[{Chu et~al(2011)Chu, Gorinevsky, and Boyd]{Chu:11 Chu E, Gorinevsky D, Boyd S (2011) Scalable statistical monitoring of fleet data. In: Proceedings of the 18th IFAC World Congress, Milan, Italy, pp 13,227--13,232
1406.2192.extracted_bib
[{Combettes and Pesquet(2011)]{com:11 Combettes PL, Pesquet JC (2011) Proximal splitting methods in signal processing. In: Fixed-Point Algorithms for Inverse Problems in Science and Engineering, Springer Optimization and Its Applications, vol~49, Springer New York, pp 185--212
1406.2192.extracted_bib
[{Curtis et~al(2010)Curtis, Schenk, and W\"achter]{cur:10 Curtis F, Schenk O, W\"achter A (2010) An interior-point algorithm for large-scale nonlinear optimization with inexact step computations. SIAM Journal on Scientific Computing 32(6):3447--3475
1406.2192.extracted_bib
[{Curtis et~al(2012)Curtis, Huber, Schenk, and W\"achter]{cur:12 Curtis FE, Huber J, Schenk O, W\"achter A (2012) A note on the implementation of an interior-point algorithm for nonlinear optimization with inexact step computations. Mathematical Programming 136(1):209--227
1406.2192.extracted_bib
[{Douglas and Rachford(1956)]{dou:56 Douglas J, Rachford HH (1956) {On the numerical solution of heat conduction problems in two and three space variables. Transaction of the American Mathematical Society 82:421--489
1406.2192.extracted_bib
[{Durazzi and Ruggiero(2003)]{dur:03 Durazzi C, Ruggiero V (2003) A {Newton inexact interior-point method for large scale nonlinear optimization problems. {Annali dell'Universit\'a di Ferrara 49(1):333--357
1406.2192.extracted_bib
[{Eckstein(1989)]{eck:89 Eckstein J (1989) Splitting methods for monotone operators with application to parallel optimization. Ph{D dissertation, {Massachussets Intitute of Technology
1406.2192.extracted_bib
[{Eckstein and Bertsekas(1992)]{eck:92 Eckstein J, Bertsekas DP (1992) {On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming 55:293--318
1406.2192.extracted_bib
[{Eisenstat and Walker(1994)]{eis:94 Eisenstat SC, Walker HF (1994) {Globally Convergent Inexact Newton Methods. SIAM Journal on Optimization 4(2):393--422
1406.2192.extracted_bib
[{El-Bakry et~al(1996)El-Bakry, Tapia, Tsuchiya, and Zhang]{bakry:96 El-Bakry A, Tapia R, Tsuchiya T, Zhang Y (1996) On the formulation and theory of the newton interior-point method for nonlinear programming. Journal of Optimization Theory and Applications 89(3):507--541, \doi{10.1007/BF02275347, \urlprefix\url{http://dx.doi.org/10.1007/BF02275347
1406.2192.extracted_bib
[{Fortin and Glowinski(1983)]{for:83 Fortin M, Glowinski R (1983) Chapter 1 augmented lagrangian methods in quadratic programming. In: Fortin M, Glowinski R (eds) Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems, Studies in Mathematics and Its Applications, vol~15, Elsevier, pp 1 -- 46
1406.2192.extracted_bib
[{Freund et~al(1999)Freund, Jarre, and Mizuno]{fre:99 Freund RW, Jarre F, Mizuno S (1999) Convergence of a class of inexact interior-point algorithms for linear programs. Mathematics of Operations Research 24(1):50--71
1406.2192.extracted_bib
[{Gabay(1983)]{gab:83 Gabay D (1983) Applications of the method of multipliers to variational inequalities. In: Fortin M, Glowinski R (eds) Augmented {Lagrangian Methods: {Applications to the Solution of Boundary-Value Problems, North-Holland
1406.2192.extracted_bib
[{Gabay and Mercier(1976)]{gab:76 Gabay D, Mercier B (1976) A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Computers \& Mathematics with Applications 2(1):17--40
1406.2192.extracted_bib
[{{Glowinski and {Marroco(1975)]{mar:75 {Glowinski R, {Marroco A (1975) {Sur l'approximation, par \'el\'ements finis d'ordre un, et la r\'esolution, par penalisation-dualit\'e, d'une classe de probl\`emes de Dirichlet non lin\'eaires. {Rev Franc Automat Inform Rech Operat, R 9(2):41--76
1406.2192.extracted_bib
[{Goldfarb et~al(2012)Goldfarb, Ma, and Scheinberg]{gol:12 Goldfarb D, Ma S, Scheinberg K (2012) Fast alternating linearization methods for minimizing the sum of two convex functions. Mathematical Programming pp 1--34
1406.2192.extracted_bib
[{Goldstein et~al(2012)Goldstein, ODonoghue, and Setzer]{gols:12 Goldstein T, O'Donoghue B, Setzer S (2012) Fast alternating direction optimization methods. Tech. Rep. CAM report 12-35, UCLA
1406.2192.extracted_bib
[{Grant and Boyd(2014)]{cvx Grant M, Boyd S (2014) {CVX: Matlab software for disciplined convex programming, version 2.1. \url{http://cvxr.com/cvx
1406.2192.extracted_bib
[{Han and Yuan(2012)]{han:12 Han D, Yuan X (2012) A note on the alternating direction method of multipliers. Journal of Optimization Theory and Applications 155(1):227--238
1406.2192.extracted_bib
[{Hansson(2000)]{han:00 Hansson A (2000) A primal-dual interior-point method for robust optimal control of linear discrete-time systems. IEEE Transactions on Automatic Control 45(9):1639--1655
1406.2192.extracted_bib
[{Hestenes(1969)]{hes:69 Hestenes MR (1969) Multiplier and gradient methods. Journal of Optimization Theory and Applications 4(5):303--320, \doi{10.1007/BF00927673
1406.2192.extracted_bib
[{{Hong and {Luo(2012)]{hon:12 {Hong M, {Luo ZQ (2012) {On the Linear Convergence of the Alternating Direction Method of Multipliers. ArXiv e-prints
1406.2192.extracted_bib
[{Iutzeler et~al(2012)Iutzeler, Ciblat, and Jakubowicz]{Iut:12 Iutzeler F, Ciblat P, Jakubowicz J (2012) Analysis of {Max-Consensus algorithms in wireless channels. IEEE Transactions on Signal Processing 60(11):6103--6107
1406.2192.extracted_bib
[{Korzak(2000)]{kor:00 Korzak J (2000) Convergence analysis of inexact infeasible-interior-point algorithms for solving linear programming problems. SIAM Journal on Optimization 11(1):133--148
1406.2192.extracted_bib
[{Liu et~al(2013)Liu, Hansson, and Vandenberghe]{liu:13 Liu Z, Hansson A, Vandenberghe L (2013) Nuclear norm system identification with missing inputs and outputs. Systems \& Control Letters 62(8):605--612
1406.2192.extracted_bib
[{Lu et~al(2009)Lu, Monteiro, and O'Neal]{lu:09 Lu Z, Monteiro RDC, O'Neal JW (2009) An iterative solver-based long-step infeasible primal-dual path-following algorithm for convex qp based on a class of preconditioners. Optimization Methods Software 24(1):123--143
1406.2192.extracted_bib
[{Mizuno and Jarre(1999)]{miz:99 Mizuno S, Jarre F (1999) Global and polynomial-time convergence of an infeasible-interior-point algorithm using inexact computation. Mathematical Programming 84(1):105--122
1406.2192.extracted_bib
[{Necoara and Suykens(2009)]{nec:09 Necoara I, Suykens JAK (2009) Interior-point lagrangian decomposition method for separable convex optimization. Journal of Optimization Theory and Applications 143(3):567--588
1406.2192.extracted_bib
[{Nedic and Ozdaglar(2009)]{ned:09 Nedic A, Ozdaglar A (2009) Distributed subgradient methods for multi-agent optimization. IEEE Transactions on Automatic Control 54(1):48--61
1406.2192.extracted_bib
[{Nedic et~al(2010)Nedic, Ozdaglar, and Parrilo]{ned:10 Nedic A, Ozdaglar A, Parrilo P (2010) Constrained consensus and optimization in multi-agent networks. IEEE Transactions on Automatic Control 55(4):922--938
1406.2192.extracted_bib
[{{Ohlsson et~al(2013){Ohlsson, {Chen, {Khoshfetrat Pakazad, {Ljung, and {Shankar Sastry]{ohl:13 {Ohlsson H, {Chen T, {Khoshfetrat Pakazad S, {Ljung L, {Shankar Sastry S (2013) {Scalable Anomaly Detection in Large Homogenous Populations. ArXiv e-prints
1406.2192.extracted_bib
[{Parikh and Boyd(2014)]{boy:14 Parikh N, Boyd S (2014) Proximal algorithms. Foundations and Trends in Optimization 1(3):127--239
1406.2192.extracted_bib
[{Peaceman and Rachford(1955)]{pea:55 Peaceman DW, Rachford J H~H (1955) The numerical solution of parabolic and elliptic differential equations. Journal of the Society for Industrial and Applied Mathematics 3(1):28--41
1406.2192.extracted_bib
[{Powell(1969)]{pow:69 Powell MJD (1969) {A method for nonlinear constraints in minimization problems. In: Fletcher R (ed) Optimization, Academic Press, New York, pp 283--298
1406.2192.extracted_bib
[{Saad(2003)]{saa:03 Saad Y (2003) Iterative Methods for Sparse Linear Systems, 2nd edn. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA
1406.2192.extracted_bib
[{Summers and Lygeros(2012)]{sum:12 Summers T, Lygeros J (2012) Distributed model predictive consensus via the alternating direction method of multipliers. In: 50th Annual Allerton Conference on Communication, Control, and Computing, pp 79--84
1406.2192.extracted_bib
[{{Teixeira et~al(2013){Teixeira, {Ghadimi, {Shames, {Sandberg, and {Johansson]{2013arXiv1303.6680T {Teixeira A, {Ghadimi E, {Shames I, {Sandberg H, {Johansson M (2013) {Optimal scaling of the ADMM algorithm for distributed quadratic programming. ArXiv e-prints \eprint{1303.6680
1406.2192.extracted_bib
[{Toh(2008)]{toh:08 Toh K (2008) An inexact primal-dual path following algorithm for convex quadratic {SDP. Mathematical Programming 112(1):221--254
1406.2192.extracted_bib
[{Wei et~al(2013)Wei, Ozdaglar, and Jadbabaie]{wei:13 Wei E, Ozdaglar A, Jadbabaie A (2013) A distributed {Newton method for network utility maximization--{I: {Algorithm. IEEE Transactions on Automatic Control 58(9):2162--2175
1406.2192.extracted_bib
[{Zhou and Toh(2004)]{zho:04 Zhou G, Toh K (2004) Polynomiality of an inexact infeasible interior point algorithm for semidefinite programming. Mathematical Programming 99(2):261--282