A Distributed Newton Method for Large Scale Consensus Optimization

Introduction

Data analysis through machine and statistical learning has become an important tool in a variety of fields including artificial intelligence, biology, medicine, finance, and marketing. Though arising in diverse applications, these problems share key characteristics, such as an extremely large number (in the order of tens of millions) of training examples typically residing in high-dimensional spaces. With this unprecedented growth in data, the need for distributed computation across multiple processing units is ever-pressing. This direction holds the promise for algorithms that are both rich-enough to capture the complexity of modern data, and scalable-enough to handle “Big Data” efficiently.

In the distributed setting, central problems are split across multiple processors each having access to local objectives. The goal is then to minimize a sum of local costs while ensuring consensus (agreement) across all processors. To clarify, consider the example of linear regression in which the goal is to find a latent model for a given dataset. Rather than searching for a centralized solution, one can distribute the optimization across multiple processors each having access to local costs defined over random subsets of the full dataset. In such a case, each processor learns a separate “chunk” of the latent model, which is then unified by incorporating consensus constraints.

Generally, there are two popular classes of algorithms for distributed optimization. The first is sub-gradient based, while the second relies on a decomposition-coordination procedure. Sub-gradient algorithms proceed by taking a gradient related step then followed by an averaging with neighbors at each iteration. The computation of each step is relatively cheap and can be implemented in a distributed fashion . Though cheap to compute, the best known convergence rate of sub-gradient methods is relatively slow given by $\mathcal{O}\left(\frac{1}{\sqrt{t}}\right)$ with $t$ being the total number of iterations . The second class of algorithms solve constrained problems by relying on dual methods. One of the well-know methods (state-of-the-art) from this class is the Alternating Direction Method of Multipliers (ADMM) . ADMM decomposes the original problem to two subproblems which are then solved sequentially leading to updates of dual variables. In , the authors show that ADMM can be fully distributed over a network leading to improved convergence rates in the order of $\mathcal{O}\left(\frac{1}{t}\right)$.

Apart from accuracy problems inherent to ADMM-based methods , much rate improvements can be gained from adopting second-order (Newton) methods. Though a variety of techniques have been proposed , less progress has been made at leveraging ADMM’s accuracy and convergence rate issues. In a recent attempt , the authors propose a distributed second-order method for general consensus by using the approach in to compute the Newton direction. As detailed in Section 6, this method suffers from two problems. First, it fails to outperform ADMM and second, faces storage and computational defficiencies for large data sets, thus ADMM retains state-of-the-art statusThe approach involves power of matrices of sizes $np\times np$ with $n$ being the total number of nodes and $p$ the number of features..

Contributions: In this paper, we contribute to the above problems and propose a distributed Newton method for general consensus with the following characteristics: i) approximating the exact newton direction up-to any arbitrary $\epsilon>0$, ii) exhibiting super-linear convergence within a neighborhood of the optimal solution similar to exact Newton, and iii) outperforming ADMM and others in terms of iteration count, running times, and total message complexity on a set of benchmark datasets, including one on a real-world application of fMRI imaging. One can argue that our improvements arrive at increased communication costs compared to other techniques. In a set of experiments, we show that such an increase is relatively small for low accuracy requirements and demonstrate a growth proportional to the condition number of the processors’ graph as accuracy demands improve. Of course, as shown in our results (see Section 6), this increase is slower compared to other methods which can be exponential.

SDD Linear Systems

Symmetric Diagonally Dominant Matrices (SDD) play a vital role in the computation of the Newton direction in a distributed fashion. In this section, we briefly review SDD systems and present a summary of efficient methods for solving them. SDD systems are linear equations of the form:

The story of computing an approximation to the exact inverse, $\bm{M}^{-1}$, starts from standard splittings of symmetric matrices. Here, $\bm{M}$ is decomposed as:

where $\bm{D}_{0}$ is a diagonal matrix consisting of diagonal elements in $\bm{M}$, while $\bm{A}_{0}$ is a non-negative symmetric matrix collecting all off-diagonal components, i.e., $[\bm{A}_{0}]_{ij}=-[\bm{M}]_{ij}$ if $i\neq j$ and otherwise. Based on this splitting, the authors in recognize that the inverse of $\bm{M}_{0}$ can be written as:

Since $\bm{D}_{0}-\bm{A}_{0}\bm{D}_{0}^{-1}\bm{A}_{0}$ is also symmetric diagonally dominant, Spielman and Peng recurse the above for a length of $d=\mathcal{O}(\log n)$ to arrive at the inverse approximated chain, $\mathcal{C}=\{\bm{D}_{i},\bm{A}_{i}\}_{i=1}^{d}$ with:

Rewriting Equation 2 in terms of the approximated chain, an $\epsilon$-close solution to $\bm{x}^{\star}$ can be determined using a two step procedure. In the first a “crude” solution to $\bm{x}^{\star}$ is returned (Algorithm 1). The procedure operates in two loops, both running to an order $d=\mathcal{O}(\log n)$. In the forward loop, intermediate vectors are constructed which are then used in the backward loop to determine the solution $\bm{x}_{0}=\bm{Z}_{0}\bm{b}_{0}$, where $\bm{Z}_{0}\approx_{\epsilon_{d}}\bm{M}^{-1}$.

Since $\bm{Z}_{0}$ incurs an $\epsilon_{d}$ error (a constant error) to the real inverse $\bm{M}^{-1}$, Spielman and Peng introduce a Richardson pre-conditioning scheme, detailed in Algorithm 2, to arrive at any arbitrary precision.

To be used in determining the Newton direction, a distributed version of the above SDD solver has to be developed. In recent work, the authors in distribute the parallel SDD solver across multiple processors. To do so, the system in Equation 1 is interpreted as represented by an undirected weighted graph, $\mathcal{G}$, with $\bm{M}$ being its Laplacian. The authors then introduce an inverse approximate chain which can be computed in a distributed fashion. Consequently, both the “crude” and exact solutions can be determined using only local communication between nodes on the undirected graph $\mathcal{G}$. In the consensus problem, we adapt this solver for distributing the computation of the Newton direction (see Section 4).

Distributed Global Consensus

Clearly, the collection of $\bm{y}_{1},\dots,\bm{y}_{p}$ is locally distributed among the nodes of graph $\mathcal{G}$, since each node $i\in\mathcal{V}$ need only to have access to the $i^{th}$ components of such vectors. Consequently, we can rewrite the problem of Equation 3 in an equivalent distributed form:

where $\bm{M}=\bm{I}_{p\times p}\otimes\mathcal{L}$ is a block-diagonal matrix with Laplacian diagonal entries, and $\bm{y}$ is a vector concatenating $\bm{y}_{1}:\bm{y}_{p}$. At this stage, our aim is to solve the problem in Equation 5 using dual techniques. Before presenting properties of the dual problem, we next introduce a standard assumption on the associated functions $f_{i}$’s:

The cost functions $\left\{f_{i}\right\}_{i=1}^{n}$ are convex and:

twice continuously differentiable with: $\gamma\preceq\nabla^{2}f_{i}\preceq\Gamma$

2 Dual Problem

Having determined the dual variables, we still require a procedure which allows us to infer about the primal. Using the above, the primal variables are determined as the solution to the following system of differential equations:

Clearly, Equation 6 is locally defined for each node $i\in\mathcal{V}$, where for each $r=1,\dots,p$:

Hence, each node $i$ can construct its own system of equations by collecting $\{\lambda_{1}(j),\dots,\lambda_{p}(j)\}$ from its neighbors $j\in\mathcal{N}(i)$ without the need for full communication. Denoting the solution of the PDE as: $y_{1}(i)=\phi_{1}^{(i)}\left((\mathcal{L}\bm{\lambda}_{1})_{i},\dots,(\mathcal{L}\bm{\lambda}_{p})_{i}\right),\ \ \ \dots,\ \ \ y_{p}(i)=\phi_{p}^{(i)}\left((\mathcal{L}\bm{\lambda}_{1})_{i},\dots,(\mathcal{L}\bm{\lambda}_{p})_{i}\right)$, we can show the following essential theoretical guarantee on the partial derivatives:

Our method for computing the Newton direction relies on the fact that the Hessian of the dual problem is an SDD matrix. We prove this property in the following lemma:

The dual function $q(\bm{\lambda})=q(\bm{\lambda}_{1},\dots,\bm{\lambda}_{p})$ shares the following characteristics:

The dual Hessian $\bm{H}(\bm{\lambda})$ and gradient $\nabla q(\bm{\lambda})$ are given by:

Distributed Newton for General Consensus

we notice that Equation 7 can be split to two SDD linear systems of the form:

for $r=1,\dots,n$. Interestingly, the above computations can be performed completely locally by noting that each node $r\in\mathcal{V}$ can compute the $r^{th}$ component of each vector $\bm{b}_{1}^{[k]},\dots,\bm{b}_{p}^{[k]}$. This is true as such a node stores $f_{r}$ as well as the variables $z^{[k]}_{1}(r),\dots,z^{[k]}_{p}(r)$. Before commencing to the convergence analysis, the final step needed is to establish the connection between the approximate solutions:

Convergence Guarantees

We next analyze the convergence properties of our approximate Newton method showing similar three convergence phases to approximate Newton methods. We start with the following lemma:

Let $\bm{g}^{[k]}=\nabla q\left(\bm{\lambda}^{[k]}\right)$ be the dual gradient at the $k^{th}$ iteration. Then:

Now, we are ready to provide the theorem summarizing the three convergence phases:

Strict Decrease Phase: while $\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}\geq\eta_{1}$: $q\left(\bm{\lambda}^{[k+1]}\right)-q\left(\bm{\lambda}^{[k]}\right)\leq-\frac{\gamma^{3}}{\Gamma^{2}}\left(\frac{1-\epsilon}{1+\epsilon}\right)^{2}\frac{\mu_{2}^{4}(\mathcal{L})}{\mu_{n}^{7}(\mathcal{L})}\eta_{1}^{2}$,

Quadratic Decrease Phase: while $\eta_{0}\leq\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}\leq\eta_{1}$: $\left|\left|\bm{g}^{[k+1]}\right|\right|_{\bm{M}}\leq\frac{1}{\eta_{1}}\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}^{2}$,

Terminal Phase: while $\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}\leq\eta_{0}$: $\left|\left|\bm{g}^{[k+1]}\right|\right|_{\bm{M}}\leq\zeta\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}$, where $\eta_{0}=\frac{\zeta(1-\zeta)}{\xi}$, $\eta_{1}=\frac{1-\zeta}{\xi}$, and $\zeta=\sqrt{\left[1-\alpha_{k}+\epsilon\alpha_{k}\sqrt{\frac{\Gamma}{\gamma}\frac{\mu_{n}^{3}(\mathcal{L})}{\mu_{2}^{3}(\mathcal{L})}}\right]},\ \ \ \xi=\frac{B(\alpha_{k}\Gamma(1+\epsilon))^{2}}{2\mu_{2}^{4}(\mathcal{L})}$.

Evaluation

We evaluate our method against five other approaches: 1) distributed Newton ADD, an adaptation of ADD that we introduce to compute the Newton direction of general consensus, 2) distributed ADMM , 3) distributed averaging , an algorithm solving general consensus using local averaging, 4) network newton 1 and 2 , and 5) distributed gradients . We are chiefly interested in the convergence speeds of both the objective value and the consensus error. The comparison against ADMM allows us to understand whether we surpass state-of-the-art, while comparisons against ADD and network newton sheds-the-light on the accuracy of our Newton’s direction approximation. Real-World Distributed Implementation: To simulate a real-world distributed environment, we used the Matlab parallel pool running on an 8 core server. After generating the processors’ graph structure with random edge assignment (see below for specifics on the node-edge configuration), we split nodes equally across the 8 cores. Hence, each processor was assigned a collection of nodes for performing computations. Communication between these nodes was handled using the Matlab Message Passing Interface (MatlabMPI) of , which allows for efficient scripting. As for bandwidth, it has been shown in that MatlabMPI can match C-MPI for large messages, and can maintain high-bandwidth even when multiple processors are communicating.

We performed three sets of experiments on standard machine learning problems: 1) linear regression, 2) logistic regression, and 3) reinforcement learning. We transformed centralized problems to fit within the distributed consensus framework. This can be easily achieved by factoring the summation running over all the available training examples to partial summations across multiple processors while introducing consensus (see Appendix H). Due to space constraints, we report two of these in this section and leave the others to the supplementary material (see Appendix G). We considered both synthetic as well as real-world data sets: Synthetic Regression Task: We created a dataset for regression with $10^{8}$ data points each being an $80$ dimensional vector. The task parameter vector $\bm{\theta}$ was generated as a linear combination of these features. The training data set $\bm{X}$ was generated from a standard normal distribution in $80$ dimensions. The training labels were given as $\bm{y}=\bm{X}\bm{\theta}+\bm{\zeta}$, where each element in $\bm{\zeta}$ was an independent univariate Gaussian noise. MNIST Data The MNIST data set is a large database of handwritten digits which has been used as a benchmark for classification algorithms . The goal is to classify among 10 different digits amounting from 1 to 10. After reading each image, we perform dimensionality reduction to reduce the number of features of each instance image to $150$ features using principle component analysis and follow a one-versus-all classification scheme.

2 Linear Regression Results

Synthetic Data: We randomly distributed the regression objective over a network of 100 nodes and 250 edges. The edges were chosen uniformly at random. An $\epsilon$ of $1/10$ was provided to the SDD solver for determining the approximate Newton direction. Step-sizes were determined separately for each algorithm using a grid-search-like-technique over {0.01, 0.1, 0.2, 0.3, 0.5, 0.6, 0.9, 1} to ensure best operating conditions. We used the local objective and the consensus error as performance metrics. Results shown in Figures 1(a) and 1(b) demonstrate that our method (titled Distributed SDD Newton) significantly outperforms all other techniques in both objective value and consensus error. Namely, distributed SDD Newton converges to the optimal value in about 40 iterations compared to about 200 for the second-best performing algorithm. It is also interesting to recognize that the worst performing algorithms were distributed gradients and network newton 1 and 2 from .

3 Logistic Regression Results

We chose the most successful algorithms from previous experiments to perform image classification. We considered both smooth (L2 norms) and non-smooth (L1 norms) regularization forms on latent parameters. The processor graph was set to 10 nodes and 20 edges generated uniformly at random. Results depicted in Figures 1(e)-1(f) demonstrate that our algorithm is again capable of outperforming state-of-the-art methods.

4 fMRI Experiment

Having shown that our approach outperforms others on relatively dense benchmark datasets, we are now interested in the performance on sparse datasets where the number of features is much larger than the number of inputs. To do so, we used the functional Magnetic Resonance Imaging (fMRI) dataset from . The goal in these experiments is to classify the cognitive state (i.e., wether looking at a picture or a sentence) of a subject based on fMRI data. Six subjects were considered in total. Each had 40 trials that lasted for 27 seconds attaining in total 54 images per-subject. After preprocessing as described in , we acquired a sparse data-set with 240 input data points, each having 43,720 features. We then performed logistic regression with an L1 regularization and reported objective values and consensus errors. Figures 2(a) and 2(b) demonstrate the objective value and consensus errors on the fMRI dataset. First, it is clear that our approach outperforms others on both criteria. It is worth noting that the second-best performing algorithm to ours is Distributed ADD-Newton; an approach we proposed in this paper for computing the Newton direction. Distributed ADMM and Distributed Averaging perform the worst on such a sparse problem. Second, Figure 2(b) clearly manifests the drawback of ADMM which requires substantial amounts of iterations for converging to the optimal feasible point. Due to the size of the feature set (i.e., 43,720) even small deviations from the optimal model can lead to significant errors in the value of the objective function. This motivates the need for the accurate solutions as acquired by our method.

5 Communication Overhead & Running Times

It can be argued that our results arrive at a high communication cost between processors. This can be true as our method relies on an SDD-solver while others allow only for few messages per iteration. We conducted a final experiment measuring local communication exchange with respect to accuracy requirements. For that, we chose the London Schools data set as all algorithms performed relatively well. Results reported in Figure 2(c) demonstrate that this increase is negligible compared to other methods. Clearly, as accuracy improves so does the communication overhead of all other algorithms. Distributed SDD-Newton has a growth rate proportional to the condition number of the graph being much slower compared to the exponential growth observed by other techniques. Finally, Figure 2(d) reports running times till convergence on the same dataset. Clearly, our method is the fastest when compared with others. The worst performing algorithms were Network Newton, distributed averaging and sub-gradients.

Conclusions & Future Work

In this paper, we proposed a distributed Newton method for solving general consensus optimization. Our method exploits the SDD property of the dual Hessian leading to an accurate computation of the Newton direction up-to-any arbitrary $\epsilon>0$. We showed that our method exhibits three phases of convergence with a quadratic phase in the neighborhood of optimal solution. In a set of experiments on standard machine learning benchmarks (including non-smooth cost functions) we demonstrated that our algorithm is capable of outperforming state-of-the-art methods, including ADMM. Finally, we empirically demonstrated that such an improvement arrives at a negligible increase in communication overhead between processors.

Our next step is to develop incremental versions of this algorithm, and use generalized Hessians to allow for non differentiable cost functions. We also plan on taking such a framework to the lifelong machine learning setting.

References

Appendix A Synopsis

We organized appendix as follows. The proofs of Lemmas 1, 2, 3, 4 are presented in sections B,C,D,E. Theorem 1 is proved is Section F. The experimental result for reinforcement learning and London Schools datasets are presented in section G. Finally, the reductions of standard machine learning problems (regression, classification, reinforcement learning) to global consensus are given in section H.

Appendix B Proof Primal-Dual Properties

Using the definition of $z_{1},\ldots z_{p}$, the primal-dual variable system can be written as:

Taking the derivative of the above system with respect to $z_{1}$ gives:

Denoting $\boldsymbol{u}_{1}=[\frac{\partial\phi^{(i)}_{1}}{\partial z_{1}},\frac{\partial\phi^{(i)}_{2}}{\partial z_{1}},\ldots,\frac{\partial\phi^{(i)}_{p}}{\partial z_{1}}]^{\mathsf{T}}$ we have:

with $\boldsymbol{u}_{r}=[\frac{\partial\phi^{(i)}_{1}}{\partial z_{r}},\frac{\partial\phi^{(i)}_{2}}{\partial z_{r}},\ldots,\frac{\partial\phi^{(i)}_{p}}{\partial z_{r}}]^{\mathsf{T}}$. For convenience, we rewrite the above systems as:

It can be clearly seen that Equation 12 implies:

Hence, using $||\boldsymbol{U}||_{2}\leq\frac{1}{\gamma}$ we have for each entry of $\boldsymbol{U}$ we have:

The above finalizes the proof of the Lemma.

Appendix C Proof Dual Function Properties

The dual function $q(\bm{\lambda})=q(\bm{\lambda}_{1},\dots,\bm{\lambda}_{p})$ shares the following characteristics:

The dual Hessian $\bm{H}(\bm{\lambda})$ and gradient $\nabla q(\bm{\lambda})$ are given by:

Recall that $\bm{y(\lambda)}$ minimizes the Lagrangian, given by

Using conjugate $f^{*}()$, the dual function can be written as:

Denote $\boldsymbol{u}=-\boldsymbol{M}\boldsymbol{\lambda}$, then the $k^{th}$ component of vector $\nabla f^{*}(-\boldsymbol{M}\boldsymbol{\lambda})$ is given by:

Hence, using $(\ref{eq_2})$ and the relation between gradients of a function and its conjugate in the expression for vector $\nabla f^{*}(-\boldsymbol{M}\boldsymbol{\lambda})$ gives:

In the next step we target matrix $\boldsymbol{F}(\boldsymbol{y}^{+})$. Using (15):

Taking the partial derivative $\frac{\partial}{\partial\lambda_{1}}$ from the both sides of the above equation gives the following for the left and right hand sides. Left hand side:

Repeating the same step for partial derivatives $\frac{\partial}{\partial\lambda_{2}},\ldots,\frac{\partial}{\partial\lambda_{np}}$ gives:

Finally, combining this result with (21) gives:

To commence, we consider each of the above statements separately. Noting that the proof of the first statement can be found in , we begin with proving the second statement.

with $B=\frac{\delta p}{\gamma}\mu_{n}^{2}(\mathcal{L})\sqrt{\mu_{n}(\mathcal{L})}$, where $\mu_{n}(\mathcal{L})$ is the largest eigenvalue of $\mathcal{L}$ and the constants $\gamma$ and $\delta$ are these given in 1.

We first remind the reader that weighted norm of a vector $\bm{v}$ and matrix $\bm{A}$ are given by:

Fixing a vector $\bm{v}$ such that $||\bm{v}||_{\bm{M}}\neq 0$, then we have:

The next step is to upper bound $\mu_{\max}(|[\nabla^{2}f(\boldsymbol{y}(\bar{\boldsymbol{\lambda}}))]^{-1}-[\nabla^{2}f(\boldsymbol{y}(\boldsymbol{\lambda}))]^{-1}|)$. For this purpose, we next study the properties of primal Hessian in details and recognize:

Claim: The primal Hessian $\nabla^{2}f(\boldsymbol{y}(\boldsymbol{\lambda}))$ shares the following properties:

Firstly, notice that for any $j\neq i$ and any $r=1\ldots,p$, we have:

Hence, the sparsity pattern of the primal hessian allows the symmetric reordering of rows and columns that transform $\nabla^{2}f(\boldsymbol{\boldsymbol{\lambda}})$ into the block diagonal matrix as:

Note that $\boldsymbol{W}(\boldsymbol{\lambda})$ preserves the important properties of $\nabla^{2}f(\boldsymbol{\boldsymbol{\lambda}})$. Particularly, the spectrum of these two matrices are the same. That can be easily seen by considering a matrix $\boldsymbol{A}$ and letting $\boldsymbol{T}_{ij}$ to be the operator that swaps $i^{th}$ and $j^{th}$ rows of $\boldsymbol{A}$. Now, consider the matrix $\bar{\boldsymbol{A}}$ resultant of the swapping of the $i^{th}$ and $j^{th}$ row and column of $\boldsymbol{A}$. Then, $\bar{\boldsymbol{A}}=\boldsymbol{T}_{ij}\boldsymbol{A}\boldsymbol{T}_{ij}$, and

where in the last step, we used the fact that $\boldsymbol{T}^{2}_{ij}=\boldsymbol{I}$. Since $\boldsymbol{W}(\boldsymbol{\lambda})$ is constructed from $\nabla^{2}f(\boldsymbol{y}(\boldsymbol{\lambda}))$ by a symmetric reordering of rows and columns, we deduce that $\text{Spectrum}(\boldsymbol{W}(\boldsymbol{\lambda}))=\text{Spectrum}(\nabla^{2}f(\boldsymbol{y}(\boldsymbol{\lambda})))$. Therefore, we can write:

which implies the property in Equation 23. To prove the property in Equation 24, we notice that if $\bar{\boldsymbol{A}}=\boldsymbol{T}_{ij}\boldsymbol{A}\boldsymbol{T}_{ij}$ and $\boldsymbol{A}$ is invertible, then so is $\bar{\boldsymbol{A}}$ and:

where we used the fact that $\boldsymbol{T}^{-1}_{ij}=\boldsymbol{T}_{ij}$. Let us denote $\{\boldsymbol{T}_{1},\ldots,\boldsymbol{T}_{l}\}$ to be a collection of operators that swap the rows of $\nabla^{2}f(\boldsymbol{y}(\boldsymbol{\lambda}))$ to transforming it to $\boldsymbol{W}(\boldsymbol{\lambda})$, i.e.:

Then, $[\nabla^{2}f(\boldsymbol{y}(\boldsymbol{\lambda}))]^{-1}=\boldsymbol{T}_{l}\cdots\boldsymbol{T}_{1}\boldsymbol{W}^{-1}(\boldsymbol{\lambda})\boldsymbol{T}_{1}\cdots\boldsymbol{T}_{l}$, and:

The above proves the property in Equation 24. ∎

Now, consider the term $\left(y_{k}(i)(\bar{\boldsymbol{\lambda}})-y_{k}(i)(\boldsymbol{\lambda})\right)$. Using the result of Lipschitzness on the solution of the partial differential equations, we can write:

Combining this result with that from Equation 24 gives:

Applying the previous equation to that in Equation 22 gives:

The above finalizes the statement of the claim and consequently that of the lemma.

Appendix D Proof Approximation Accuracy

We next denote the exact solution of the system in Equation 8 as:

Second, we let $\boldsymbol{z}^{*[k]}=\left[\nabla^{2}f(\boldsymbol{y}^{[k]})\right]^{-1}\boldsymbol{M}\boldsymbol{d}^{*[k]}=[(\boldsymbol{z}^{*[k]}_{1})^{\mathsf{T}},\ldots,(\boldsymbol{z}^{*[k]}_{p})^{\mathsf{T}}]^{\mathsf{T}}$ be the corresponding vector $\boldsymbol{z}^{[k]}$. Hence:

The next step is to rewrite the right and left hand sides of Equation 28 in terms of $\boldsymbol{\hat{d}}^{[k]}$ and $\boldsymbol{d}^{*[k]}$:

where we use the fact that $\boldsymbol{M}[\nabla^{2}f(\boldsymbol{y}^{[k]})]^{-1}\boldsymbol{M}[\nabla^{2}f(\boldsymbol{y}^{[k]})]^{-1}\boldsymbol{M}\succeq\frac{1}{\mu_{n}(\boldsymbol{\mathcal{L}})}[\boldsymbol{H}^{[k]}]^{2}$. The last transition follows from the fact that $\ker(\boldsymbol{H}^{[k]})=\ker(\boldsymbol{M})$ and, therefore:

Combining the above results for (28) immediately gives:

with $\epsilon_{1}=\epsilon_{0}\frac{\mu_{n}(\boldsymbol{\mathcal{L}})}{\mu_{2}(\boldsymbol{\mathcal{L}})}\sqrt{\frac{\Gamma}{\gamma}}$.

where we used $\boldsymbol{H}^{[k]}\succeq\frac{\Gamma}{\mu_{2}(\boldsymbol{\mathcal{L}})}\boldsymbol{M}$. Hence:

Notice that from Equation 29, it follows that:

Therefore, combining Equations 29, 31, and 32, in Equation 30, yields:

The above finalizes the statement of the lemma. ∎

Appendix E Proof Dual Gradient Change

Let $\bm{g}^{[k]}=\nabla q\left(\bm{\lambda}^{[k]}\right)$ be the dual gradient at the $k^{th}$ iteration. Then:

We start with the following claim, which plays a crucial role in our analysis:

Claim: Let $\nabla q(\boldsymbol{\lambda})$ be the dual gradient and $\boldsymbol{H}(\boldsymbol{\lambda})$ its Hessian, then for any $\bar{\boldsymbol{\lambda}},\boldsymbol{\lambda}$:

We proceed by adding and subtracting $\boldsymbol{H}(\boldsymbol{\lambda})(\bar{\boldsymbol{\lambda}}-\boldsymbol{\lambda})$ to the integral in the right hand side of Equation 34:

Consequently, we can separate the integral in Equation 35 as:

The second integral on the right hand side of Equation 36 is independent of $t$. Therefore, we can simplify the integral as $\boldsymbol{H}(\boldsymbol{\lambda})(\bar{\boldsymbol{\lambda}}-\boldsymbol{\lambda})$, which implies:

By rearranging the terms in Equation 37 and taking the norm of both sides we obtain:

Considering the inequality in Equation 38 and the fact that norm of integral is less than the integral of the norms, we can write:

Now, let us consider the following three cases:

$\boldsymbol{v}\in\ker\{\boldsymbol{M}\}^{\perp}$: In this case Equation 40 follows immediately from the definition: $||\boldsymbol{A}||_{\boldsymbol{M}}=\sup_{\boldsymbol{v}:\boldsymbol{v}\notin\ker\{\boldsymbol{M}\}}\frac{||\boldsymbol{Av}||_{\boldsymbol{M}}}{||\boldsymbol{v}||_{\boldsymbol{M}}}$.

$\boldsymbol{v}\in\ker\{\boldsymbol{M}\}$: In this case $||\boldsymbol{v}||_{\boldsymbol{M}}=0$, and $[\boldsymbol{H}(\bar{\boldsymbol{\lambda}})-\boldsymbol{H}(\boldsymbol{\lambda})]\boldsymbol{v}=\boldsymbol{0}-\boldsymbol{0}=\boldsymbol{0}$.

$\boldsymbol{v}=\boldsymbol{u}_{1}+\boldsymbol{u}_{2}$, where $\boldsymbol{u}_{1}\in\ker\{\boldsymbol{M}\}^{\perp},\boldsymbol{u}_{2}\in\ker{\boldsymbol{M}}$. In this case $||\boldsymbol{v}||_{\boldsymbol{M}}=||\boldsymbol{u}_{1}||_{\boldsymbol{M}}$, and using the first case result for $\boldsymbol{u}_{1}\in\ker\{\boldsymbol{M}\}^{\perp}$, we have:

Applying the above result to Equation 39 gives:

Applying the result of claim E to $\boldsymbol{\lambda}^{[k+1]}$ and $\boldsymbol{\lambda}^{[k]}$ gives:

Applying the triangular inequality, we have:

Because $\boldsymbol{g}^{[k]}\in\ker\{\boldsymbol{H}^{[k]}\}=\ker\{\boldsymbol{M}\}$ and using the result in Equation 42, it follows that:

where $\mu_{\min}(\boldsymbol{M})$ is the second smallest eigenvalue of $\boldsymbol{M}$ which is equal to $\mu_{2}(\boldsymbol{\mathcal{L}})$.

Therefore, our goal now is to upper bound the term $||\boldsymbol{H}^{[k]}\boldsymbol{c}^{[k]}||_{\boldsymbol{M}}$. Using Equation 42 and the fact that $\boldsymbol{M}\preceq\mu_{n}(\boldsymbol{\mathcal{L}})\boldsymbol{I}_{np\times np}$, we have:

Applying this result in Equation 44 gives:

Applying the results of Equations 43 and 45, in Equation 41:

Appendix F Proof Convergence Phases

Strict Decrease Phase: while $\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}\geq\eta_{1}$: $q\left(\bm{\lambda}^{[k+1]}\right)-q\left(\bm{\lambda}^{[k]}\right)\leq-\frac{\gamma^{3}}{\Gamma^{2}}\left(\frac{1-\epsilon}{1+\epsilon}\right)^{2}\frac{\mu_{2}^{4}(\mathcal{L})}{\mu_{n}^{7}(\mathcal{L})}\eta_{1}^{2}$,

Quadratic Decrease Phase: while $\eta_{0}\leq\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}\leq\eta_{1}$: $\left|\left|\bm{g}^{[k+1]}\right|\right|_{\bm{M}}\leq\frac{1}{\eta_{1}}\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}^{2}$,

Terminal Phase: while $\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}\leq\eta_{0}$: $\left|\left|\bm{g}^{[k+1]}\right|\right|_{\bm{M}}\leq\zeta\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}$, where $\eta_{0}=\frac{\zeta(1-\zeta)}{\xi}$, $\eta_{1}=\frac{1-\zeta}{\xi}$, and $\zeta=\sqrt{\left[1-\alpha_{k}+\epsilon\alpha_{k}\sqrt{\frac{\Gamma}{\gamma}\frac{\mu_{n}^{3}(\mathcal{L})}{\mu_{2}^{3}(\mathcal{L})}}\right]},\ \ \ \xi=\frac{B(\alpha_{k}\Gamma(1+\epsilon))^{2}}{2\mu_{2}^{4}(\mathcal{L})}$.

We will proof each phase separately. We start with phase one when $\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}\geq\eta$. Taking the Taylor expansion of the dual function gives:

Hence, choosing ${\alpha}_{k}=\alpha^{\star}=\left(\frac{\gamma}{\Gamma}\right)^{2}\left(\frac{\mu_{2}(\mathcal{L})}{\mu_{n}(\mathcal{L})}\right)^{4}\frac{1-\epsilon}{(1+\epsilon)^{2}}$ and using $\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}^{2}\geq\eta_{1}$, we arrive at the strict decrease phase. To prove the quadratic decrease phase, we let $\eta_{0}\leq\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}^{2}<\eta$. It immediately follows that:

Finally for $\left|\left|\bm{g}^{[k]}\right|\right|_{\bm{M}}^{2}<\eta_{0}$, we have

The above finalizes the statement of the theorem. ∎

Appendix G London Schools & Reinforcement Learning Results

The London Schools data set consists of examination scores from 15,362 students in 139 schools. This is a benchmark regression task with a goal of predicting examination scores of each student. We use the same feature encoding used in , where four school-specific categorical variables along with three student-specific categorical variables are encoded as a collection of binary features. In addition, we use the examination year and a bias term as additional features, giving each data instance 27 features.

G.2 Reinforcement Learning

We considered the policy search framework to control a double cart-pole system (DCP). As detailed in , the DCP adds a second inverted pendulum to the standard cart-pole system, with six parameters and six state features. The goal is to balance both poles upright. We generated 20,000 rollouts each with a length of 150 time steps.

Appendix H Reductions to Global Consensus

with $\mu_{i}$ being a regularization parameterIn our experiments $\mu_{i}$’s were fixed for all nodes. These were set for values between {0.01, 0.02, 0.05, 0.06, 0.1} depending on the size of input dataset., and $\sum_{i=1}^{n}m_{i}=m$. To simplify the computations, next we rewrite Equation 47 in an equivalent matrix-vector form:

Introducing $\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{p}$, we can write the linear regression problem as:

where $\boldsymbol{\theta}_{i}=[y_{1}(i),y_{2}(i),\ldots,y_{p}(i)]^{\mathsf{T}}$, and:

Hence, the Lagrangian of the problem can be written as:

Therefore, the primal variables can be recovered using:

where $(\boldsymbol{L\Lambda})(i,:)$ is $i^{th}$ row of $\mathcal{L}\bm{\Lambda}$. Note, that to apply the symmetric diagonally dominant the second time the Hessian of the local objective function $f_{r}(y_{1}(r),y_{2}(r),\ldots,y_{p}(r))$ must be computed. It is easy to see, that $\nabla^{2}f_{r}(y_{1}(r),y_{2}(r),\ldots,y_{p}(r))=2\boldsymbol{P}_{i}$.

In this section, we will describe distributed ADMM method for linear regression. Recall, that in distributed ADMM each node implements the following instructions:

Each agent $i$ updates its estimate of $\boldsymbol{\theta}^{[k]}_{i}$ in a sequential order with

Each agent updates $\boldsymbol{\lambda}_{ji}$ for $j\in P(i)$ as follows:

We can get the closed form solution for (49):

H.1.2 Linear Regression via Distributed Averaging

In distributed averaging, each node keeps three variables:

where $\beta$ is a step-size and $\boldsymbol{g}_{i}(t)$ is the sub-gradient of $f_{i}$ evaluated at $\boldsymbol{w}_{i}(t)$, i.e.:

H.2 Logistic Regression

Here, however, each local cost is given as a logistic loss defined as:

with $\mu_{i}$ being a regularization parameter and $\sum_{i=1}^{n}m_{i}=m$. $\Psi(\bm{\theta}_{i})$ is a regularization function. Next, we consider two such cases when $\Psi$ is both smooth and non-smooth.

When $\Psi(\bm{\theta}_{i})=||\bm{\theta}_{i}||_{2}^{2}$, the local objective in Equation 52 can be written as:

It is easy to see that Equation 53 can be simplified as:

where $\boldsymbol{\theta}_{i}=[y_{1}(i),y_{2}(i),\ldots,y_{p}(i)]^{\mathsf{T}}$. Using Equation 53:

Hence, the Lagrangian of the above problem can be written as:

Applying the standard Newton method for Equation 56:

where $t$ is the iteration count, and $\boldsymbol{y}^{[i]}_{Newton}|_{(t)}$ is Newton direction, evaluated at $[y_{1}(i),\ldots,y_{p}(i)]_{(t)}$, and given by the solution of the following system:

with $\boldsymbol{H}_{i}|_{(t)}$ and $\nabla\zeta_{i}|_{(t)}$ being the Hessian and noting that the gradient of $\zeta_{i}$ is evaluated at $[y_{1}(i),\ldots,y_{p}(i)]_{(t)}$. The components of the gradient $\nabla\zeta_{i}|_{(t)}$ can be computed as:

This can be written in the following matrix-vector form:

where $(\boldsymbol{L\Lambda})(i,:)$ is $i^{th}$ row of matrix $\boldsymbol{L\Lambda}$ and:

The Hessian of $\zeta_{i}$ can be immediately written as:

where $\boldsymbol{D}_{i}|_{(t)}$ is diagonal $m_{i}\times m_{i}$ matrix, such that

It is again easy to see that to apply the SDDM-solver the Hessian of the local objective is needed. This can be derived as:

where $\boldsymbol{D}_{r}$ is diagonal matrix, given by:

The derivations, so-far, are appropriate for the proposed distributed Newton method. To be able to compare against ADMM, corresponding mathematics needs to be developed. This section details such constructs. We start by recalling that each node in distributed ADMM implements the following instructions:

Each agent $i$ updates its estimate of $\boldsymbol{\theta}^{[k]}_{i}$ in a sequential order with:

Each agent updates $\boldsymbol{\lambda}_{ji}$ for $j\in P(i)$ as follows:

For solving the optimization problem in Equation 65, standard Newton is used for function $\xi_{i}$:

where $\alpha_{t}$ is a step-size and $\boldsymbol{\theta}^{[i]}_{Newton}|_{(t)}$ is the Newton Direction, evaluated at $\boldsymbol{\theta}^{[k+1]}_{i}(t)$, and given by the solution to the following system:

with $\boldsymbol{H}_{i}|_{(t)}$ and $\nabla\xi_{i}(\boldsymbol{\theta}^{[k+1]}_{i}(t))$ being the Hessian and the gradient of function $\xi_{i}$ evaluated at vector $\boldsymbol{\theta}^{[k+1]}_{i}(t)$. Clearly, the gradient can be written in a vector form as:

The Hessian of $\xi_{i}$ can be immediately written as:

where $\boldsymbol{D}_{i}|_{(t)}$ is diagonal $m_{i}\times m_{i}$ matrix, such that

In this section, we describe the necessary mathematics needed for distributed averaging. We recall that distributed averaging operates as:

where $\beta$ is a step-size and $\boldsymbol{g}_{i}(t)$ is the sub-gradient of $f_{i}$ evaluated at $\boldsymbol{w}_{i}(t)$, i.e.:

H.2.2 Non-Smooth Regularizers

Now, we consider the case when the local objective is defined by:

It is easy to see that (74) can be simplified as:

Similar to Equations 54 and 55, we introduce $\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{p}$ and $\boldsymbol{b}_{1},\ldots,\boldsymbol{b}_{m_{i}}$. Then, the problem can be written as:

where $\boldsymbol{\theta}_{i}=[y_{1}(i),y_{2}(i),\ldots,y_{p}(i)]^{\mathsf{T}}$, and using Equation 74 we have:

Clearly, $\zeta_{i}$ is not smooth. To proceed, we use the smooth approximations of the L1 norm, given by:

where $\alpha$ is parameter controlling the approximation’s quality. Therefore, the local objective can be written as:

Applying standard Newton for Equation 76 gives:

where $t$ is the iteration count, and $\boldsymbol{y}^{[i]}_{Newton}|_{(t)}$ is Newton direction, evaluated at vector $[y_{1}(i),\ldots,y_{p}(i)]_{(t)}$, and given by the solution to the following system:

with $\boldsymbol{H}_{i}|_{(t)}$ and $\nabla\zeta_{i}|_{(t)}$ being the Hessian and the gradient of $\zeta_{i}$ evaluated at vector $[y_{1}(i),\ldots,y_{p}(i)]_{(t)}$. The components of the gradient $\nabla\zeta_{i}|_{(t)}$ can be computed as:

The above can be written in a vector-matrix form as:

where $\boldsymbol{B}_{i},\boldsymbol{\Lambda},\boldsymbol{\delta}|_{(t)}$ are defined in (61),(62), (63) respectively, $(\boldsymbol{L\Lambda})(i,:)$ is the $i^{th}$ row of matrix $\boldsymbol{L\Lambda}$, and

The Hessian of $\zeta_{i}$ can be immediately written as:

where $\boldsymbol{D}_{i}|_{(t)}$ is a diagonal $m_{i}\times m_{i}$ matrix, such that

and $\boldsymbol{\Delta}_{i}|_{(t)}$ is diagonal $p\times p$ matrix, such that

To apply SDDM-solver, the Hessian of the local objective function $f_{r}(y_{1}(r),y_{2}(r),\ldots,y_{p}(r))$ must be computed. It is easy to see that:

where $\boldsymbol{D}_{r}$ and $\boldsymbol{\Delta}_{r}$ are diagonal matrices, given by:

In this section, we describe ADMM for logistic regression problems with $L1$ regularization. The approach is very similar to L2 case, with the following differences:

The gradient of function $\xi_{i}$ is given as:

where $[\boldsymbol{\theta}^{[k+1]}_{i}(t)]_{j}$ is $j^{th}$ component of vector $\boldsymbol{\theta}^{[k+1]}_{i}(t)$.

The Hessian of function $\xi_{i}$ can be written:

where $\boldsymbol{D}_{i}|_{(t)}$ is diagonal $m_{i}\times m_{i}$ matrix, given in (70), and $\boldsymbol{\Delta}_{i}|_{(t)}$ is diagonal $p\times p$ matrix such that

where $[\boldsymbol{\theta}^{[k+1]}_{i}(t)]_{j}$ is $j^{th}$ component of vector $\boldsymbol{\theta}^{[k+1]}_{i}(t)$.

In this section, we describe the mathematics needed for distributed averaging. Again, these are similar to the L2 regularization setting with the following differences:

The sub-gradient of $f_{i}$ can be written as :

where $\boldsymbol{\delta}|_{(t)}$ is given in (72) and

where $[\boldsymbol{w}_{i}(t)]_{j}$ is $j^{th}$ component of vector $\boldsymbol{w}_{i}(t)$.

H.3 Reinforcement Learning

We consider the policy search framework for reinforcement learning with a uni-variate Gaussian policy. Here, the data is represented as a collection of trajectories $\{\boldsymbol{\tau}_{i}\}^{m}_{i=1}$, where

with $\mu_{i}$ being a regularization parameter, and $\sum_{i=1}^{n}m_{i}=m$. Easy to see that (88) can be simplified as:

To further simplify expressions (89) let us introduce:

Similarly to (54) we introduce vectors $\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{p}$, then problem (87) can be written as:

where $\boldsymbol{\theta}_{i}=[y_{1}(i),y_{2}(i),\ldots,y_{p}(i)]^{\mathsf{T}}$, and:

and primal variables can be recovered from the dual by the following equation:

where $(\boldsymbol{L\Lambda})(i,:)$ is $i^{th}$ row of matrix $\boldsymbol{L\Lambda}$, and $\boldsymbol{\Lambda}$ is given as (62).

In this section we will describe distributed ADMM for the reinforcement learning problem in (87). Recall, that in distributed ADMM each nodes implements the following instructions:

Each agent $i$ updates its estimate of $\boldsymbol{\theta}^{[k]}_{i}$ in a sequential order with

Each agent updates $\boldsymbol{\lambda}_{ji}$ for $j\in P(i)$ as follows:

We can get the closed form solution for (65):

H.3.2 Reinforcement Learning via Distributed Averaging

In this section, we describe distributed averaging for reinforcement learning (87). Recall, that:

where $\beta$ is a step-size and $\boldsymbol{g}_{i}(t)$ is the sub-gradient of $f_{i}$ evaluated at $\boldsymbol{w}_{i}(t)$, i.e.: