BFGS-ADMM for Large-Scale Distributed Optimization

Distributed optimization acts as the computational engine for a wide range of applications in modern technology, including distributed control [1], power grid management [2], distributed machine learning [3], and resource allocation in sensor networks [4]. We consider a classical problem of minimizing a sum of cost functions and a convex regularizer, aiming at promoting sparsity in the solution structure. In specific:

where $f^{i}:\mathbb{R}^{d}\to\mathbb{R}$ captures the $i$-th agent’s objective of interest and $g:\mathbb{R}^{d}\to\mathbb{R}$ is a possibly nonsmooth regularizer, such as the weighted $\ell_{1}-$norm. In a network of agents, each agent aims to minimize its local cost function $f^{i}(\cdot)$ while communicating with its neighbors to cooperatively find the solution of the global problem. Distributed solutions are often pursued to fully utilize the computational power of the agents and reduce the amount of message passing in the network. Efficient communication is deemed especially desirable in scenarios where the number of participants $m$, and the dimension of the decision variable $d$ is large, such as in emerging applications of Cyber-Physical systems (CPS) [5]. An archetypal framework of distributed optimization introduces local decision variables at each agent, which are updated locally using private data and exchange variables with neighboring agents. A network-wide solution is achieved by means of asymptotic elimination of the consensus error, i.e., the difference between the values of neighboring agents converges to zero. Such a framework is termed distributed consensus optimization and offers significant flexibility for each agent to select appropriate updating schemes that suit its local hardware environment.

First order methods [6]-[9] constitute popular choices for solving (1) in a distributed fashion. The presence of the nonsmooth regularizer prohibits gradient descent from being implemented directly as $g(\cdot)$ is not differentiable. Subgradient methods [6] invoke the notion of subdifferential set to compute descent directions while letting agents exchange information over a time-varying topology. Proximal gradient [8] accommodates the nonsmoothness by splitting the objective function and evaluating the proximal mapping associated with the nondifferentiable part. Various acceleration schemes [7], [9] exist for first order methods, that involve storing past gradient and iterate information so as to form a momentum correction term. However, as the size of data grow the condition number of problem (1) tends to get large which causes first order methods to exhibit extremely slow convergence rate. A natural remedy for the aforementioned issues is to consider second order methods [10]-[12]. Through use of function curvature information, second order methods compute descent directions in the objective function level sets that are effective at accelerating the convergence over first order methods. Proximal Newton [12] can be considered as the second order counterpart of the proximal gradient but the associated proximal mapping increases computational burden drastically and renders distributed implementation infeasible as global information is required to compute the Newton step. Moreover, second order methods require solving a linear system to compute the Newton step which induces a computation cost of $\mathcal{O}(d^{3})$. Quasi-Newton methods [13] circumvent this procedure by using a finite difference of gradients to approximate the Hessian. The reduced computation costs along with competitive performance have rendered quasi-Newton methods a desirable alternative to second order methods, especially for Large-scale optimization problems.

Primal-dual methods [14]-[15] provide a different perspective for problem (1) within the framework of constrained optimization, where the consensus constraint is explicitly enforced. The Alternating Direction Method of Multipliers (ADMM) [15] falls into this catergory where the smooth and nonsmooth parts of the objective are considered separately and iterates are computed sequentially. However, at each iteration of the ADMM, a sub-optimization problem must be solved, which typically induces an expensive computational burden.

Contributions: (i) We propose BFGS-ADMM for convex composite optimization where the primal update (typically computational expensive) is replaced with a quasi-Newton step that does not require solving a linear system of equations. Moreover, by introducing an intermediate consensus variable, we eliminate the need for inner loops within the network when computing the update direction. (ii) We establish global linear convergence rate for the proposed algorithm without backtracking, under the assumption that private cost functions are strongly convex with Lipschitz continuous gradients. (iii) The advantage of BFGS-ADMM over first order methods is analytically established and experimentally demonstrated with numerical simulations on real datasets.

Notation: We denote column vectors $x\in\mathbb{R}^{d}$ with lower case letters and matrices $A\in\mathbb{R}^{n\times m}$ with capital letters. Superscripts are used to denote partitioned vector components and subscripts are used to denote iteration steps, e.g., $x^{i}_{t}$ denotes the value of subvector $x^{i}$ at step $t$. Matrix entries are denoted as $[A]^{ij}$. Unless specified otherwise, $\norm{x}$ and $\norm{A}$ denote the Euclidean norm of a vector and the corresponding induced norm of a matrix, respectively. We define the norm of a vector associated with a positive definite matrix $P\succ 0$ as $\norm{x}_{P}:=\sqrt{x^{\top}Px}$, and the set $\{1,\dots,m\}$ is abbreviated as $[m]$. The proximal mapping associated with the function $g(\cdot):\mathbb{R}^{d}\to\mathbb{R}$ is defined as $\textbf{prox}_{\mu g}(v)=\underset{\theta}{\mathrm{argmin}}\left\{g(\theta)+\tfrac{1}{2\mu}\norm{\theta-v}^{2}\right\}$.

We build a local model of the Augmented Lagrangian at step $t$ with respect to the primal variable $x$ as:

where $\mathcal{L}_{t}:=\mathcal{L}(x_{t},z_{t},\theta_{t};y_{t},\lambda_{t})$. Various algorithms can be designed by choosing different $H_{t}$. In this paper, we opt to use regularized BFGS approximation of the Augmented Lagrangian Hessian as follows:

where $\epsilon>0$ and $\Delta=A^{\top}A$ is the Kronecker product of the graph degree matrix $\hat{\Delta}$ and $I_{d}$. Note that by setting $B_{t}=\nabla^{2}F(x_{t})$, we recover the Hessian of the Augmented Lagrangian with respect to $x$, i.e., $H_{t}=\nabla_{xx}^{2}\mathcal{L}$ in such case. The quadratic form of (12) admits a closed form solution when minimizing over $x$, i.e.,

Step (11b) reduces to solving the following linear system of equations:

Moreover, by completion of squares, step (11c) is equivalent to evaluating the following proximal mapping associated with $g(\cdot)$ with parameter $\mu_{2}$:

Dual variables are updated in verbatim as in (11d) and (11e). We proceed to present Lemma 1 which states that with an appropriate initialization, the intermediate variable $\{z^{ij}\}$ evolves on the manifold defined by the column space of $E_{u}$. Matrices without the hat denotes the Kronecker product of the corresponding matrix with the identity matrix, e.g., $E_{u}=\hat{E}_{u}\otimes I_{d}$ where $\hat{E}_{u}$ is defined in (4).

where $\textbf{diag}[s_{i}]\in\mathbb{R}^{md}$ is the diagonal matrix formed by $s_{i}I_{md}$, and the dual variables $(\theta_{t+1},\alpha_{t+1},\lambda_{t+1})$ are updated in the exact same way as in the BFGS-ADMM using (17b)-(17c). Note that our first order variant encapsulates existing first order methods EXTRA [21] as a special case by setting $s_{i}=\tfrac{1}{\epsilon}$. Indeed, when $g(\cdot)=0$, the proximal mapping and the associated dual variables in (17) are not needed and we obtain: $x_{t+1}=x_{t}-\tfrac{1}{\epsilon}\left[\nabla F(x_{t})+E_{s}^{\top}\alpha_{t}+\tfrac{1}{2\mu_{1}}L_{s}x_{t}\right].$ Taking the difference between the $x_{t+1}$ update and the $x_{t}$ update defined this way, and substituting the dual update for $\alpha_{t}$ using (17c), we obtain:

With appropriate choices of $\mu_{1},\mu_{2},\epsilon$, the update defined in (22) is equivalent to EXTRA. We present the comparison between BFGS-ADMM and its first order variant in terms of approximating the exact ADMM in the next section.
We proceed to present the distributed implementation of BFGS-ADMM in Algorithm 1. Note that in the primal update (17a), dual variables are invoked in the form of $E_{s}^{\top}\alpha_{t}$. By defining $\phi_{t}=E_{s}^{\top}\alpha_{t}\in\mathbb{R}^{md}$ and pre-multiplying both sides of (17c) with $E_{s}^{\top}$, we obtain an equivalent algorithm that allows for efficient distributed implementation. We let each agent hold the corresponding pair $(x^{i},\phi^{i})$ while the $l$-th agent additionally holds $(\theta,\lambda)$. At each iteration, each agent $i$ begins with calculating $h_{t}^{i}$ using $x_{t}^{j}$ from its neighbors as in step 3, where the $l$-th agent performs an additional update in step 5. Note that step 7 does not require solving a linear system since $(H_{t}^{-1})^{ii}$ is directly approximated using (18). Primal and dual updates are performed at each agent as in step 8 and step 9, respectively. The $l$-th agent evaluates an additional proximal mapping in step 11 to update $\theta_{t+1}$, and updates $\lambda_{t+1}$ in step 12. Finally, each agent $i$ updates its $(B_{t+1}^{-1})^{ii}$ and $(H_{t+1}^{-1})^{ii}$ using (10) and (18), respectively.

Linear convergence rate for ADMM is well established [15]. Since our goal here is to reduce the computational cost of ADMM, the best one can hope for is to reduce the gap between the approximation and the standard ADMM, while maintaining the linear convergence rate. In Lemma 4, we demonstrate the advantage of using BFGS approximation over first order methods by showing that the aforementioned gap decreases faster for BFGS-ADMM. We first state an additional assumption on the nonsmooth regularizer $g(\cdot)$ and the KKT conditions associated with (3) expressed in the primal optimal pair $(x_{\star},\theta_{\star})$ and the dual optimal pair $(\alpha_{\star},\lambda_{\star})$. Note that we have eliminated $z_{\star}$, as it is shown in Lemma 1 that it is a function of $x_{\star}$.

Assumption 2: The function $g(\cdot):\mathbb{R}^{d}\to\mathbb{R}$ is proper, closed, and convex. Equivalently, $\forall\,x,y\in\mathbb{R}^{d}$, the following inequality holds:

where $\partial g(\cdot)$ denotes the subdifferential set. We proceed to state a non-restrictive assumption that can be easily achieved by regularization.

We proceed to state a technical lemma that describes an inclusion relationship between the subdifferential set of $\partial g(\theta_{t+1})$ and the dual variable $\lambda_{t+1}$.

Recall the first order variant of the BFGS-ADMM where the primal variable is updated as in (21) while BFGS-ADMM updates $x_{t+1}$ as in (17a), and both update dual variables as in (17b)-(17c). We proceed to present results which capture their similarities and differences in the following.