Network Newton-Part I: Algorithm and Convergence

I Introduction

Problems of this form arise often in, e.g., decentralized control systems , wireless systems , sensor networks , and large scale machine learning . In the latter case, distributed formulations are efficient in dealing with very large datasets where it is desirable to split training sets into smaller subsamples that are assigned to different servers . In this paper we assume that the local costs $f_{i}$ are twice differentiable and strongly convex. Therefore, the aggregate cost function $f$ is also twice differentiable and strongly convex.

There are different algorithms to solve (1) in a distributed manner. The most popular choices are decentralized gradient descent (DGD) , distributed implementations of the alternating direction method of multipliers , and decentralized dual averaging (DDA) . Although there are substantial differences between them, these methods can be generically abstracted as combinations of local descent steps followed by variable exchanges and averaging of information among neighbors. A feature common to all of these algorithms is the slow convergence rate in ill-conditioned problems since they operate on first order information only. This is not surprising because gradient descent methods in centralized settings where the aggregate function gradient is available at a single server have the same difficulties in problems with skewed curvature [see Chapter 9 of .]

This issue is addressed in centralized optimization by Newton’s method that uses second order information to determine a descent direction adapted to the objective’s curvature [see Chapter 9 of ]. In general, second order methods are not available in distributed settings because distributed approximations of Newton steps are difficult to devise. In the particular case of flow optimization problems, these approximations are possible when operating in the dual domain . As would be expected, these methods result in large reductions of convergence times.

Our goal here is to develop approximate Newton’s methods to solve (1) in distributed settings where agents have access to their local functions only and exchange variables with neighboring agents. We do so by introducing Network Newton (NN), a method that relies on distributed approximations of Newton steps for the global cost function $f$ to accelerate convergence of DGD. We begin the paper with an alternative formulation of (1) and a brief discussion of DGD (Section II). We then introduce a reinterpretation of DGD as an algorithm that utilizes gradient descent to solve a penalized version of (1) in lieu of the original optimization problem (Section II-A). This reinterpretation explains convergence of DGD to a neighborhood of the optimal solution. The volume of this neighborhood is given by the relative weight of the penalty function and the original objective which is controlled by a penalty coefficient.

If gradient descent on the penalized function finds an approximate solution to the original problem, the same solution can be found with a much smaller number of iterations by using Newton’s method. Alas, distributed computation of Newton steps requires global communication between all nodes in the network and is therefore impractical (Section III). To resolve this issue we approximate the Newton step of the penalized objective function by truncating the Taylor series expansion of the exact Newton step (Section III-A). This results in a family of methods indexed by the number of terms of the Taylor expansion that are kept in the approximation. The method that results from keeping $K$ of these terms is termed NN-$K$. A fundamental observation here is that the Hessian of the penalized function has a sparsity structure that is the same sparsity pattern of the graph. Thus, when computing terms in the Hessian inverse expansion, the first order term is as sparse as the graph, the second term is as sparse as the two hop neighborhood, and, in general, the $k$-th term is as sparse as the $k$-hop neighborhood of the graph. Thus, implementation of the NN-$K$ method requires aggregating information from $K$ hops away. Increasing $K$ makes NN-$K$ arbitrarily close to Newton’s method at the cost of increasing the communication overhead of each iteration.

Convergence of NN-$K$ to the optimal argument of the penalized objective is established (Section IV). We do so by establishing several auxiliary bounds on the eigenvalues of the matrices involved in the definition of the method (Propositions 1-3 and Lemma 2). Of particular note, we show that a measure of the error between the Hessian inverse approximation utilized by NN-$K$ and the actual inverse Hessian decays exponentially with the method index $K$. This exponential decrease hints that using a small value of $K$ should suffice in practice. Convergence is formally claimed in Theorem 1 that shows the convergence rate is at least linear. It follows from this convergence analysis that larger penalty coefficients result in faster convergence that comes at the cost of increasing the distance between the optimal solutions of the original and penalized objectives. The convergence guarantees established in this paper are not better than the corresponding guarantees for DGD. These advantages are established in a companion paper where we further show that the sequence of penalized objective function values generated by NN-$K$ has a convergence rate that is quadratic in a specific interval. This quadratic phase holds for all $K$ and can be made arbitrarily large by increasing $K$ . Numerical results in establish the advantages of NN-$K$ in terms of number of iterations and communications steps relative to DGD and establish that using $K=1$ or $K=2$ tends to work best in practice.

II Distributed Gradient Descent

Since the network is connected, the constraints ${\mathbf{x}}_{i}={\mathbf{x}}_{j}$ for all $i$ and $j\in{\mathcal{N}}_{i}$ imply that (1) and (II) are equivalent in the sense that we have ${\mathbf{x}}_{i}^{*}={\mathbf{x}}^{*}$ for all $i$. This must be the case because for a connected network the constraints ${\mathbf{x}}_{i}={\mathbf{x}}_{j}$ for all $i$ and $j\in{\mathcal{N}}_{i}$ collapse the feasible space of (II) to a hyperplane in which all local variables are equal. When all local variables are equal, the objectives in (1) and (II) coincide and, in particular, so do their optima.

Since $w_{ij}=0$ when $j\neq i$ and $j\notin\mathcal{N}_{i}$, it follows from (3) that each agent $i$ updates its estimate ${\mathbf{x}}_{i}$ of the optimal vector ${\mathbf{x}}^{*}$ by performing an average over the estimates ${\mathbf{x}}_{j,t}$ of its neighbors $j\in\mathcal{N}_{i}$ and its own estimate ${\mathbf{x}}_{i,t}$, and descending through the negative local gradient $-\nabla f_{i}({\mathbf{x}}_{i,t})$. DGD is a distributed method because to implement (3), node $i$ exchanges variables with neighboring nodes only.

If the considions in (4) are true, it is possible to show that (3) approaches the solution of (1) in the sense that ${\mathbf{x}}_{i,t}\approx{\mathbf{x}}^{*}$ for all $i$ and large $t$, . The accepted interpretation of why (3) converges is that nodes are gradient descending towards their local minima because of the term $-\alpha\nabla f_{i}({{\mathbf{x}}_{i,t}})$ but also perform an average of neighboring variables $\sum_{j=1}^{n}w_{ij}{\mathbf{x}}_{j,t}$. This latter consensus operation drives the agents to agreement. In the following section we show that (3) can be alternatively interpreted as a penalty method.

where in the second equality we added and subtracted ${\mathbf{y}}_{t}$ and regrouped terms. Inspection of (5) reveals that the DGD update formula at step $t$ is equivalent to a (regular) gradient descent algorithm being used to solve the program

Indeed, given the definition of the function $F({\mathbf{y}}):=(1/2){\mathbf{y}}^{T}({\mathbf{I}}-{\mathbf{Z}})\ {\mathbf{y}}+\alpha\sum_{i=1}^{n}f_{i}({\mathbf{x}}_{i})$ it follows that the gradient of $F({\mathbf{y}})$ at ${\mathbf{y}}={\mathbf{y}}_{t}$ is given by

Using (7) we rewrite (5) as ${\mathbf{y}}_{t+1}={\mathbf{y}}_{t}-{\mathbf{g}}_{t}$ and conclude that DGD descends along the negative gradient of $F({\mathbf{y}})$ with unit stepsize. The expression in (3) is just a distributed implementation of gradient descent that uses the gradient in (7). To confirm that this is true, observe that the $i$th element of the gradient ${\mathbf{g}}_{t}=[{\mathbf{g}}_{i,t};\ldots;{\mathbf{g}}_{i,t}]$ is given by

The gradient descent iteration ${\mathbf{y}}_{t+1}={\mathbf{y}}_{t}-{\mathbf{g}}_{t}$ is then equivalent to (3) if we entrust node $i$ with the implementation of the descent ${\mathbf{x}}_{i,t+1}={\mathbf{x}}_{i,t}-{\mathbf{g}}_{i,t}$, where, we recall, ${\mathbf{x}}_{i,t}$ and ${\mathbf{x}}_{i,t+1}$ are the $i$th components of the vectors ${\mathbf{y}}_{t}$ and ${\mathbf{y}}_{t+1}$. Observe that the local gradient component ${\mathbf{g}}_{i,t}$ can be computed using local information and the ${\mathbf{x}}_{j,t}$ iterates of its neighbors $j\in{\mathcal{N}}_{i}$. This is as it should be, because the descent ${\mathbf{x}}_{i,t+1}={\mathbf{x}}_{i,t}-{\mathbf{g}}_{i,t}$ is equivalent to (3).

Is it a good idea to descend on $F({\mathbf{y}})$ to solve (1)? To some extent. Since we know that the null space of ${\mathbf{I}}-{\mathbf{W}}$ is $\text{null}({\mathbf{I}}-{\mathbf{W}})=\text{span}({\mathbf{1}})$ and that ${\mathbf{Z}}={\mathbf{W}}\otimes{\mathbf{I}}$ we know that the span of ${\mathbf{I}}-{\mathbf{Z}}$ is $\text{null}({\mathbf{I}}-{\mathbf{Z}})=\text{span}({\mathbf{1}}\otimes{\mathbf{I}})$. Thus, we have that $({\mathbf{I}}-{\mathbf{Z}}){\mathbf{y}}={\mathbf{0}}$ holds if and only if ${\mathbf{x}}_{1}=\dots={\mathbf{x}}_{n}$. Since the matrix ${\mathbf{I}}-{\mathbf{Z}}$ is positive semidefinite – because it is stochastic and symmetric –, the same is true of the square root matrix $({{\mathbf{I}}-{\mathbf{Z}}})^{1/2}$. Therefore, we have that the optimization problem in (II) is equivalent to the optimization problem

III Network Newton

Instead of solving (6) with a gradient descent algorithm as in DGD, we can solve (6) using Newton’s method. To implement Newton’s method we need to compute the Hessian ${\mathbf{H}}_{t}:=\nabla^{2}F({\mathbf{y}}_{t})$ of $F$ evaluated at ${\mathbf{y}}_{t}$ so as to determine the Newton step ${\mathbf{d}}_{t}:=-{\mathbf{H}}_{t}^{-1}{\mathbf{g}}_{t}$. Start by differentiating twice in (6) in order to write ${\mathbf{H}}_{t}$ as

While the Hessian ${\mathbf{H}}_{t}$ is sparse, the inverse ${\mathbf{H}}_{t}$ is not. It is the latter that we need to compute the Newton step ${\mathbf{d}}_{t}:={\mathbf{H}}_{t}^{-1}{\mathbf{g}}_{t}$. To overcome this problem we split the diagonal and off diagonal blocks of ${\mathbf{H}}_{t}$ and rely on a Taylor’s expansion of the inverse. To be precise, write ${\mathbf{H}}_{t}={\mathbf{D}}_{t}-{\mathbf{B}}$ where the matrix ${\mathbf{D}}_{t}$ is defined as

Proceed now to factor ${\mathbf{D}}_{t}^{1/2}$ from both sides of the splitting relationship to write ${\mathbf{H}}_{t}={\mathbf{D}}_{t}^{{1}/{2}}({\mathbf{I}}-{\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}){\mathbf{D}}_{t}^{{1}/{2}}$. When we consider the Hessian inverse ${\mathbf{H}}^{-1}$, we can use the Taylor series $({\mathbf{I}}-{\mathbf{X}})^{-1}=\sum_{j=0}^{\infty}{\mathbf{X}}^{j}$ with ${\mathbf{X}}={\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ to write

Observe that the sum in (14) converges if the absolute value of all the eigenvalues of the matrix ${\mathbf{D}}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}^{-{1}/{2}}$ are strictly less than 1. For the time being we assume this to be the case but we will prove that this is true in Section IV. When the series converge, we can use truncations of this series to define approximations to the Newton step as we explain in the following section.

Network Newton (NN) is defined as a family of algorithms that rely on truncations of the series in (14). The $K$th member of this family, NN-$K$, considers the first $K+1$ terms of the series to define the approximate Hessian inverse

NN-$K$ uses the approximate Hessian ${\hat{\mathbf{H}}}_{t}^{(K)^{-1}}$ as a curvature correction matrix that is used in lieu of the exact Hessian inverse ${\mathbf{H}}^{-1}$ to estimate the Newton step. I.e., instead of descending along the Newton step ${\mathbf{d}}_{t}:=-{\mathbf{H}}_{t}^{-1}{\mathbf{g}}_{t}$ we descend along the NN-$K$ step ${\mathbf{d}}_{t}^{(K)}:=-{\hat{\mathbf{H}}}_{t}^{(K)^{-1}}{\mathbf{g}}_{t}$, which we intend as an approximation of ${\mathbf{d}}_{t}$. Using the explicit expression for ${\hat{\mathbf{H}}}_{t}^{(K)^{-1}}$ in (15) we write the NN-$K$ step as

where, we recall, the vector ${\mathbf{g}}_{t}$ is the gradient of objective function $F({\mathbf{y}})$ defined in (7). The NN-$K$ update formula can then be written as

Then observe that since the matrix ${\mathbf{B}}$ has the sparsity pattern of the graph, this recursion can be decomposed into local components

The matrix ${\mathbf{D}}_{ii,t}=\alpha\nabla^{2}f_{i}({\mathbf{x}}_{i,t})+2(1-w_{ii}){\mathbf{I}}$ is stored and computed at node $i$. The gradient component ${\mathbf{g}}_{i,t}=(1-w_{ii}){\mathbf{x}}_{i,t}-\sum_{j\in\mathcal{N}_{i}}w_{ij}{\mathbf{x}}_{j,t}+\alpha\nabla f_{i}({\mathbf{x}}_{i,t})$ is also stored and computed at $i$. Node $i$ can also evaluate the values of the matrix blocks ${\mathbf{B}}_{ii}=(1-w_{ii}){\mathbf{I}}$ and ${\mathbf{B}}_{ij}=w_{ij}{\mathbf{I}}$. Thus, if the NN-$k$ step components ${\mathbf{d}}_{j,t}^{(k)}$ are available at neighboring nodes $j$, node $i$ can then determine the NN-$(k+1)$ step component ${\mathbf{d}}_{i,t}^{(k+1)}$ upon being communicated that information.

The expression in (19) represents an iterative computation embedded inside the NN-$K$ recursion in (17). For each time index $t$, we compute the local component of the NN- step ${\mathbf{d}}_{i,t}^{(0)}=-{\mathbf{D}}_{ii,t}^{-1}{\mathbf{g}}_{i,t}$. Upon exchanging this information with neighbors we use (19) to determine the NN-$1$ step components ${\mathbf{d}}_{i,t}^{(1)}$. These can be exchanged and plugged in (19) to compute ${\mathbf{d}}_{i,t}^{(2)}$. Repeating this procedure $K$ times, nodes ends up having determined their NN-$K$ step component ${\mathbf{d}}_{i,t}^{(K)}$ .

The resulting NN-$K$ method is summarized in Algorithm 1. The descent iteration in (17) is implemented in Step 11. Implementation of this descent requires access to the NN-$K$ descent direction ${\mathbf{d}}_{i,t}^{(K)}$ which is computed by the loop in steps 6-10. Step $6$ initializes the loop by computing the NN-0 step ${\mathbf{d}}_{i,t}^{(0)}=-{\mathbf{D}}_{ii,t}^{-1}{\mathbf{g}}_{i,t}$. The core of the loop is in Step 9 which corresponds to the recursion in (19). Step 8 stands for the variable exchange that is necessary to implement Step 9. After $K$ iterations through this loop, the NN-$K$ descent direction ${\mathbf{d}}_{i,t}^{(K)}$ is computed and can be used in Step 11. Both, steps 6 and 9, require access to the local gradient component ${\mathbf{g}}_{i,t}$. This is evaluated in Step 5 after receiving the prerequisite information from neighbors in Step 4. Steps 1 and 3 compute the blocks ${\mathbf{B}}_{ii,t}$, ${\mathbf{B}}_{ij,t}$, and ${\mathbf{D}}_{ii,t}$ that are also necessary in steps 6 and 9.

IV Convergence Analysis

In this section we show that as time progresses the sequence of objective function values $F({\mathbf{y}}_{t})$ [cf. (6)] approaches the optimal objective function value $F({\mathbf{y}}^{*})$. In proving this claim we make the following assumptions.

There exists constants $0\leq\delta\leq\Delta<1$ that lower and upper bound the diagonal weights for all $i$,

The local objective functions $f_{i}({\mathbf{x}})$ are twice differentiable and the eigenvalues of the local objective function Hessians are bounded with positive constants $0<m\leq M<\infty$, i.e.

Notice that the lower bound in Assumption 1 is more a definition than a constraint since we may have $\delta=0$. This is not recommendable as it is implies that the weight $w_{ii}$ assigned to the local variable ${\mathbf{x}}_{i}$ in (3) is null, but nonetheless allowed. The upper bound $\Delta<1$ on the weights $w_{ii}$ is true for all connected networks as long as neighbors $j\in{\mathcal{N}}_{i}$ are assigned nonzero weights $w_{ij}>0$. This is because the matrix ${\mathbf{W}}$ is doubly stochastic [cf. (4)], which implies that $w_{ii}=1-\sum_{j\in{\mathcal{N}}_{i}}w_{ij}<1$ as long as $w_{ij}>0$.

The lower bound $m$ for the eigenvalues of local objective function Hessians $\nabla^{2}f_{i}({\mathbf{x}})$ is equivalent to the strong convexity of local objective functions $f_{i}({\mathbf{x}})$ with parameter $m$. The strong convexity assumption for the local objective functions $f_{i}({\mathbf{x}})$ stated in Assumption 2 is customary in convergence proofs of Newton-based methods, since the Hessian of objective function should be invertible to establish Newton’s method [Chapter 9 of ]. The upper bound $M$ for the eigenvalues of local objective function Hessians $\nabla^{2}f_{i}({\mathbf{x}})$ is similar to the condition that gradients $\nabla f_{i}({\mathbf{x}})$ are Lipschitz continuous with parameter $M$ for the case that functions are twice differentiable.

The restriction imposed by Assumption 3 is also typical of second order methods . Assumption 3 guarantees that the Hessian matrices of objective functions $F({\mathbf{y}})$ are also Lipschitz continuous as we show in the following lemma.

Consider the definition of objective function $F({\mathbf{y}})$ in (6). If Assumption 3 holds then the objective function Hessian ${\mathbf{H}}({\mathbf{y}})=:\nabla^{2}F({\mathbf{y}})$ is Lipschitz continuous with parameter $\alpha L$, i.e.

Lemma 1 states that the penalty objective function introduced in (6) has the property that the Hessians are Lipschitz continuous, while the Lipschitz constant is a function of the penalty coefficient $1/\alpha$. This observation implies that as we increase the penalty coefficient $1/\alpha$, or, equivalently, decrease $\alpha$, the objective function $F({\mathbf{y}})$ approaches a quadratic form because the curvature becomes constant.

To prove convergence properties of NN we need bounds for the eigenvalues of the block diagonal matrix ${\mathbf{D}}_{t}$, the block sparse matrix ${\mathbf{B}}$, and the Hessian ${\mathbf{H}}_{t}$. These eigenvalue bounds are established in the following proposition using the conditions imposed by Assumptions 1 and 2.

Consider the definitions of matrices ${\mathbf{H}}_{t}$, ${\mathbf{D}}_{t}$, and ${\mathbf{B}}$ in (10), (12), and (13), respectively. If Assumptions 1 and 2 hold true, then the eigenvalues of matrices ${\mathbf{H}}_{t}$, ${\mathbf{D}}_{t}$, and ${\mathbf{B}}$ are uniformly bounded as

Proposition 1 states that Hessian matrix ${\mathbf{H}}_{t}$ and block diagonal matrix ${\mathbf{D}}_{t}$ are positive definite, while matrix ${\mathbf{B}}$ is positive semidefinite.

As we noted in Section III, for the expansion in (14) to be valid the eigenvalues of the matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ must be nonnegative and strictly smaller than $1$. The following proposition states that this is true for all times $t$.

Consider the definitions of the matrices ${\mathbf{D}}_{t}$ in (12) and ${\mathbf{B}}$ in (13). If Assumptions 1 and 2 hold true, the matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ is positive semidefinite and its eigenvalues are bounded above by a constant $\rho<1$

where $\rho:=2(1-\delta)/(2(1-\delta)+{\alpha m})$.

Error matrix ${\mathbf{E}}_{t}$ measures closeness of the Hessian inverse approximation matrix ${\hat{\mathbf{H}}}_{t}^{(K)^{-1}}$ and the exact Hessian inverse ${\mathbf{H}}^{-1}_{t}$ at time $t$. Based on the definition of error matrix ${\mathbf{E}}_{t}$, if the Hessian inverse approximation ${\hat{\mathbf{H}}}_{t}^{(K)^{-1}}$ approaches the exact Hessian inverse ${\mathbf{H}}_{t}^{-1}$ the error matrix ${\mathbf{E}}_{t}$ approaches the zero matrix ${\mathbf{0}}$. We therefore bound the error of the Hessian inverse approximation by developing a bound for the eigenvalues of the error matrix ${\mathbf{E}}_{t}$. This bound is provided in the following proposition where we further show that the error of the Hessian inverse approximation for NN-$K$ decreases exponentially as we increases $K$.

Consider the NN-$K$ method as introduced in (12)-(17) and the definition of error matrix ${\mathbf{E}}_{t}$ in (28). Further, recall the definition of the constant $\rho:=2(1-\delta)/(\alpha+2(1-\delta))<1$ in Proposition 2. The error matrix ${\mathbf{E}}_{t}$ is positive semidefinite and all its eigenvalues are upper bounded by $\rho^{K+1}$,

Proposition 3 asserts that the error in the approximation of the Hessian inverse, thereby on the approximation of the Newton step, is bounded by $\rho^{K+1}$. This result corroborates the intuition that the larger $K$ is, the closer that ${\mathbf{d}}_{i,t}^{(K)}$ approximates the Newton step. This closer approximation comes at the cost of increasing the communication cost of each descent iteration. The decrease of this error being proportional to $\rho^{K+1}$ hints that using a small value of $K$ should suffice in practice. This has been corroborated in numerical experiments where $K=1$ and $K=2$ tend to work best – see . Further note that to decrease $\rho$ we can increase $\delta$ or increase $\alpha$. Increasing $\delta$ calls for assigning substantial weight to $w_{ii}$. Increasing $\alpha$ comes at the cost of moving the solution of (6) away from the solution of (II-A) and its equivalent (1).

Bounds on the eigenvalues of the objective function Hessian ${\mathbf{H}}_{t}$ are central to the convergence analysis of Newton’s method [Chapter 9 of]. Lower bounds for the Hessian eigenvalues guarantee that the matrix is nonsingular. Upper bounds imply that the minimum eigenvalue of the Hessian inverse ${\mathbf{H}}^{-1}$ is strictly larger than zero, which, in turn, implies a strict decrement in each Newton step. Analogous bounds for the eigenvalues of the NN approximate Hessian inverses ${{\hat{\mathbf{H}}}_{t}^{(K)^{-1}}}$ are required. These bounds are studied in the following lemma.

Consider the NN-$K$ method as defined in (12)-(17). If Assumptions 1 and 2 hold true, the eigenvalues of the approximate Hessian inverse ${\hat{\mathbf{H}}}_{t}^{(K)^{-1}}$ are bounded as

where constants $\lambda$ and $\Lambda$ are defined as

According to the result of Lemma 2, the NN-$K$ approximate Hessian inverses ${\hat{\mathbf{H}}}_{t}^{(K)^{-1}}$ are strictly positive definite and have all of their eigenvalues bounded between the positive and finite constants $\lambda$ and $\Lambda$. This is true for all $K$ and uniform across all iteration indexes $t$. Considering these eigenvalue bounds and the fact that $-{\mathbf{g}}_{t}$ is a descent direction, the approximate Newton step $-{\hat{\mathbf{H}}}_{t}^{(K)^{-1}}{\mathbf{g}}_{t}$ enforces convergence of the iterate ${\mathbf{y}}_{t}$ to the optimal argument ${\mathbf{y}}^{*}$ of the penalized objective function $F({\mathbf{y}})$ in (6). In the following theorem we show that if the stepsize $\epsilon$ is properly chosen, the sequence of objective function values $F({\mathbf{y}}_{t})$ converges at least linearly to the optimal objective function value $F({\mathbf{y}}^{*})$.

Consider the NN-$K$ method as defined in (12)-(17) and the objective function $F({\mathbf{y}})$ as introduced in (6). Further, recall the definitions of the lower and upper bounds $\lambda$ and $\Lambda$, respectively, for the eigenvalues of the approximate Hessian inverse ${\hat{\mathbf{H}}}_{t}^{(K)^{-1}}$ in (31). If the stepsize $\epsilon$ is chosen as

and Assumptions 1, 2, and 3 hold true, the sequence $F({\mathbf{y}}_{t})$ converges to the optimal argument $F({\mathbf{y}}^{*})$ at least linearly with constant $0<1-\zeta<1$. I.e.,

where the constant $0<\zeta<1$ is explicitly given by

Theorem 1 shows that the objective function error sequence $F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*})$ asymptoticly converges to zero and that the rate of convergence is at least linear. Note that according to the definition of the convergence parameter $\zeta$ in Theorem 1 and the definitions of $\lambda$ and $\Lambda$ in (31), increasing $\alpha$ leads to faster convergence. This observation verifies existence of a tradeoff between rate and accuracy of convergence. For large values of $\alpha$ the sequence generated by Network Newton converges faster to the optimal solution of (6). These faster convergence comes at the cost of increasing the distance between the optimal solutions of (6) and (1). Conversely, smaller $\alpha$ implies smaller gap between the optimal solutions of (6) and (1), but the convergence rate of NN-$K$ is slower. This suggests value in the use of adaptive strategies for the selection of $\alpha$ that we develop in .

V Conclusions

This paper developed the network Newton method as an approximate Newton method for solving distributed optimization problems where the components of the objective function are available at different nodes of a network. The algorithm builds on a reinterpretation of distributed gradient descent as a penalty method and relies on an approximation of the Newton step of the corresponding penalized objective function. To approximate the Newton direction we truncate the Taylor series of the exact Newton step. This leads to a family of methods defined by the number $K$ of Taylor series terms kept in the approximation. When we keep $K$ terms of the Taylor series, the method is called NN-$K$ and can be implemented through the aggregation of information in $K$-hop neighborhoods. We showed that the proposed method converges at least linearly to the solution of the penalized objective, and, consequently, to a neighborhood of the optimal argument for the original optimization problem. It follows from this convergence analysis that larger penalty coefficients result in faster convergence that comes at the cost of increasing the distance between the optimal solutions of the original and penalized objectives.

This paper does not show any advantage of NN relative to distributed gradient descent, other than the expectation to see improved convergence times due to the attempt to approximate the Newton direction of the penalized objective. These advantages are shown in a companion paper where we: (i) Show that the convergence rate is quadratic in a specific interval that can be made arbitrarily large by increasing $K$ . (ii) Use numerical results to establish the advantages of NN-$K$ in terms of number of iterations and communications steps relative to DGD.

Appendix A Proof of Lemma 1

The result in (35) is implied by the fact that the matrix ${\mathbf{I}}-{\mathbf{Z}}$ does not depend on the argument ${\mathbf{y}}$ of the Hessian ${\mathbf{H}}({\mathbf{y}})$. The next step is to bound the norm of the difference for two ${\mathbf{G}}$ matrices $\|{\mathbf{G}}({\mathbf{y}})-{\mathbf{G}}({\hat{\mathbf{y}}})\|$ in terms of the difference between two vectors ${\mathbf{y}}$ and ${\hat{\mathbf{y}}}$, i.e. $\|{\mathbf{y}}-{\hat{\mathbf{y}}}\|$.

According to the definition of ${\mathbf{G}}({\mathbf{y}})$ in (11), the difference matrix ${\mathbf{G}}({\mathbf{y}})-{\mathbf{G}}({\hat{\mathbf{y}}})$ is block diagonal and the $i$th diagonal block is

Observe that each summand in (A) can be upper bounded by applying Cauchy-Schwarz inequality as

Substituting the upper bound in (38) into (A) implies that the squared norm $\left\|{\mathbf{G}}({\mathbf{y}})-{\mathbf{G}}({\hat{\mathbf{y}}})\right\|_{2}^{2}$ is bounded above as

Observe that Assumption 2 states that the local objective function Hessians $\nabla^{2}f_{i}({\mathbf{x}}_{i})$ are Lipschitz continuous with parameter $L$, i.e., $\|\nabla^{2}f_{i}({\mathbf{x}}_{i})-\nabla^{2}f_{i}({\hat{\mathbf{x}}}_{i})\|\leq L\|{\mathbf{x}}_{i}-{\hat{\mathbf{x}}}_{i}\|$. Considering this inequality the upper bound in (39) can be changed by replacing $\|\nabla^{2}f_{i}({\mathbf{x}}_{i})-\nabla^{2}f_{i}({\hat{\mathbf{x}}}_{i})\|$ by $L\|{\mathbf{x}}_{i}-{\hat{\mathbf{x}}}_{i}\|$ which yields

Note now that for any sequences of scalars $a_{i}$ and $b_{i}$, the inequality $\sum_{i=1}^{n}a_{i}^{2}b_{i}^{2}\leq(\sum_{i=1}^{n}a_{i}^{2})(\sum_{i=1}^{n}b_{i}^{2})$ holds. If we divide both sides of this relation by $\sum_{i=1}^{n}b_{i}^{2}$ and set $a_{i}=\|{\mathbf{x}}_{i}-{\hat{\mathbf{x}}}_{i}\|$ and $b_{i}=\|{\mathbf{v}}_{i}\|$, we obtain

Combining the two inequalities in (40) and (41) leads to

Since the right hand side of (42) does not depend on the vector ${\mathbf{v}}$ we can eliminate maximization with respect ${\mathbf{v}}$. Further, note that according to the structure of vectors ${\mathbf{y}}$ and ${\hat{\mathbf{y}}}$, we can write $\left\|{\mathbf{y}}-{\hat{\mathbf{y}}}\right\|_{2}^{2}=\sum_{i=1}^{n}\left\|{\mathbf{x}}_{i}-{\hat{\mathbf{x}}}_{i}\right\|_{2}^{2}$. These two observations in association with (42) imply that squared norm of the difference between matrices ${\mathbf{G}}({\mathbf{y}})$ and ${\mathbf{G}}({\hat{\mathbf{y}}})$ is bounded above by

Taking the square root of both sides of (43) implies

According to (44) we can conclude that the matrix ${\mathbf{G}}$ is Lipschitz continuos with parameter $L$. Considering the expression in (35) and the inequality in (44), the claim in (23) follows.

Appendix B Proof of Proposition 1

We first study the bounds for the eigenvalues of matrix ${\mathbf{I}}-{\mathbf{Z}}$. Notice that since ${\mathbf{I}}-{\mathbf{Z}}$ can be written as $({\mathbf{I}}_{n}-{\mathbf{W}})\otimes{\mathbf{I}}_{p}$, all the eigenvalues of matrix ${\mathbf{I}}-{\mathbf{Z}}$ are in the spectrum of matrix ${\mathbf{I}}-{\mathbf{W}}$. Therefore, we can study the bounds for the eigenvalues of matrix ${\mathbf{I}}-{\mathbf{W}}$ in lieu of matrix ${\mathbf{I}}-{\mathbf{Z}}$. The Gershgorin circle theorem states that each eigenvalue of a matrix ${\mathbf{A}}$ lies within at least one of the Gershgorin discs $D(a_{ii},R_{ii})$ where the center $a_{ii}$ is the $i$th diagonal element of $A$ and the radius $R_{ii}:=\sum_{j\neq i}|a_{ij}|$ is the sum of the absolute values of all the non-diagonal elements of the $i$th row. Note that matrix ${\mathbf{I}}-{\mathbf{W}}$ is symmetric and as a result all eigenvalues are real. Hence, Gershgorin discs can be considered as intervals of width $[a_{ii}-R_{ii},a_{ii}+R_{ii}]$ for matrix ${\mathbf{I}}-{\mathbf{W}}$, where $a_{ii}=1-w_{ii}$ and $R_{ii}=\sum_{j\neq i}|w_{ij}|$. Since all the elements of matrix ${\mathbf{W}}$ are non-negative, $|w_{ij}|$ can be substituted by $w_{ij}$. Therefore, all the eigenvalues of matrix ${\mathbf{I}}-{\mathbf{W}}$ in at least one of the intervals $[1-w_{i}-\sum_{j\neq i}w_{ij},1-w_{i}+\sum_{j\neq i}w_{ij}]$. Now observing that sum of the weights that a node assigns to itself and all the other nodes is one, i.e. $\sum_{j}w_{ij}=1$, it can be derived that $1-w_{ii}=\sum_{j\neq i}^{n}w_{ij}$. This observation implies that the Gershgorin intervals can be simplified as $[0,2(1-w_{ii})]$ for $i=1,\dots,n$. This observation in association with the fact that $2(1-w_{ii})\leq 2(1-\delta)$ implies that all the eigenvalues of matrix ${\mathbf{I}}-{\mathbf{W}}$ are in the interval $[0,2(1-\delta)]$ and consequently the eigenvalues of matrix ${\mathbf{I}}-{\mathbf{Z}}$ are bounded as

To prove bounds for the eigenvalues of Hessian ${\mathbf{H}}_{t}$, first we find lower and upper bounds for the eigenvalues of matrix ${\mathbf{G}}_{t}$. Since matrix ${\mathbf{G}}_{t}$ is block diagonal and the eigenvalues of each diagonal block ${\mathbf{G}}_{ii,t}=\nabla^{2}f_{i}({\mathbf{x}}_{i,t})$ are bounded by constants $0<m\leq M<\infty$ as mentioned in (21), we obtain that the eigenvalues of matrix ${\mathbf{G}}_{t}$ are bounded as

Considering the definition of Hessian ${\mathbf{H}}_{t}:={\mathbf{I}}-{\mathbf{Z}}+\alpha{\mathbf{G}}_{t}$ and the bounds in (45) and (46), the first claim follows.

We proceed now to prove bounds for the eigenvalues of block diagonal matrix ${\mathbf{D}}_{t}$. According to the definition of matrix ${\mathbf{D}}_{t}$ in (12) we can write

where ${\mathbf{W}}_{d}$ is defined as ${\mathbf{W}}_{d}:=\text{diag}({\mathbf{W}})$. Note that matrix ${\mathbf{I}}_{n}-{\mathbf{W}}_{d}$ is diagonal and the $i$-th diagonal component is $1-w_{ii}$. Since the local weights satisfy $\delta\leq w_{ii}\leq\Delta$, we obtain that eigenvalues of matrix ${\mathbf{I}}_{n}-{\mathbf{W}}_{d}$ are bounded below and above by $1-\Delta$ and $1-\delta$, respectively. Observe that the eigenvalue sets of matrices $({\mathbf{I}}_{n}-{\mathbf{W}}_{d})$ and $({\mathbf{I}}_{n}-{\mathbf{W}}_{d})\otimes{\mathbf{I}}_{p}$ are identical which implies

Considering the relation in (47) and bounds in (46) and (48), the second claim follows.

Based on the definition of matrix ${\mathbf{B}}$ in (13) and the relation that ${\mathbf{Z}}={\mathbf{W}}\otimes{\mathbf{I}}$ we can write

Hence, to bound eigenvalues of matrix ${\mathbf{B}}$ we study lower and upper bounds for the eigenvalues of matrix ${\mathbf{I}}-2{\mathbf{W}}_{d}+{\mathbf{W}}$. Observe that in the $i$-th row of matrix ${\mathbf{I}}-2{\mathbf{W}}_{d}+{\mathbf{W}}$, the diagonal component is $1-w_{ii}$ and the $j$th component is $w_{ij}$ for all $j\neq i$. Using Gershgorin theorem and the same argument that we established for the eigenvalues of ${\mathbf{I}}-{\mathbf{Z}}$, we can write

Considering (50) and the expression for matrix ${\mathbf{B}}$ in (49), the last claim follows.

Appendix C Proof of Proposition 2

According to the result of Proposition 1, the block diagonal matrix ${\mathbf{D}}_{t}$ is positive definite and matrix ${\mathbf{B}}$ is positive semidefinite which immediately implies that matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ is positive semidefinite and the lower bound in (27) follows.

Recall the definition of block diagonal matrix ${\mathbf{D}}_{t}$ in (12) and define matrix ${\hat{\mathbf{D}}}$ as a special case of matrix ${\mathbf{D}}_{t}$ for $\alpha=0$. I.e., ${\hat{\mathbf{D}}}:=2({\mathbf{I}}-{\mathbf{Z}}_{d})$. Notice that matrix ${\hat{\mathbf{D}}}$ is diagonal, only depends on the structure of the network, and that it is also time invariant. Since matrix ${\hat{\mathbf{D}}}$ is diagonal and each diagonal component $1-w_{ii}$ is strictly larger than 0, matrix ${\hat{\mathbf{D}}}$ is positive definite and invertible. Therefore, we can write ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ as

The next step is to find an upper bound for the eigenvalues of the symmetric term ${\hat{\mathbf{D}}}^{-{1/2}}{\mathbf{B}}{\hat{\mathbf{D}}}^{-{1/2}}$ in (51). Observing the fact that matrices ${\hat{\mathbf{D}}}^{-{1}/{2}}{\mathbf{B}}{\hat{\mathbf{D}}}^{-{1}/{2}}$ and ${\mathbf{B}}{\hat{\mathbf{D}}}^{-1}$ are similar, eigenvalues of these matrices are identical. Therefore, we proceed to characterize an upper bound for the eigenvalues of matrix ${\mathbf{B}}{\hat{\mathbf{D}}}^{-1}$. Based on the definitions of matrices ${\mathbf{B}}$ and ${\hat{\mathbf{D}}}$, the product ${\mathbf{B}}{\hat{\mathbf{D}}}^{-1}$ is given by

Therefore, the general form of matrix ${\mathbf{B}}{\hat{\mathbf{D}}}^{-1}$ is

Note that each diagonal component of matrix ${\mathbf{B}}{\hat{\mathbf{D}}}^{-1}$ is $1/2$ and that the sum of non-diagonal components of column $i$ is

Now, by considering the result in (54) and applying Gershgorin theorem we can conclude that eigenvalues of matrix ${\mathbf{B}}{\hat{\mathbf{D}}}^{-1}$ are bounded as

where $\mu_{i}({\mathbf{B}}{\hat{\mathbf{D}}}^{-1})$ indicates the $i$-th eigenvalue of matrix ${\mathbf{B}}{\hat{\mathbf{D}}}^{-1}$. The bounds in (55) and similarity of matrices ${\mathbf{B}}{\hat{\mathbf{D}}}^{-1}$ and ${\hat{\mathbf{D}}}^{-1/2}{\mathbf{B}}{\hat{\mathbf{D}}}^{-1/2}$ show that the eigenvalues of matrix ${\hat{\mathbf{D}}}^{-1/2}{\mathbf{B}}{\hat{\mathbf{D}}}^{-1/2}$ are uniformly bounded in the interval

Based on the decomposition in (51) to characterize the bounds for the eigenvalues of matrix ${\mathbf{D}}_{t}^{-1/2}{\mathbf{B}}{\mathbf{D}}_{t}^{-1/2}$, the bounds for the eigenvalues of matrix ${\hat{\mathbf{D}}}^{1/2}{\mathbf{D}}_{t}^{-1/2}$ should be studied as well. Notice that according to the definitions of matrices ${\hat{\mathbf{D}}}$ and ${\mathbf{D}}_{t}$, the product ${\hat{\mathbf{D}}}^{1/2}{\mathbf{D}}_{t}^{-1/2}$ is block diagonal and the $i$-th diagonal block is

Observe that according to Assumption 1, the eigenvalues of local Hessian matrices $\nabla^{2}f_{i}({\mathbf{x}}_{i})$ are bounded by $m$ and $M$. Further notice that the diagonal elements of weight matrix $w_{ii}$ are bounded by $\delta$ and $\Delta$, i.e. $\delta\leq w_{ii}\leq\Delta$. Considering these bounds we can show that the eigenvalues of matrices $(\alpha/2(1-w_{11}))\nabla^{2}f_{i}({\mathbf{x}}_{i,t})+{\mathbf{I}}$ are lower and upper bounded as

By considering the bounds in (58), the eigenvalues of each block of matrix ${\hat{\mathbf{D}}}^{1/2}{\mathbf{D}}_{t}^{-1/2}$ as introduced in (57) are bounded below and above as

Since (59) holds for all the diagonal blocks of matrix ${\hat{\mathbf{D}}}^{1/2}{\mathbf{D}}_{t}^{-1/2}$, the eigenvalues of this matrix also satisfy the bounds in (59) which implies that

for $i=1,\dots,n$. Observing the decomposition in (51), the norm of the matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ is upper bounded as

Considering the symmetry of matrices ${\hat{\mathbf{D}}}^{1/2}{\mathbf{D}}_{t}^{-1/2}$ and ${\hat{\mathbf{D}}}^{-{1}/{2}}{\mathbf{B}}{\hat{\mathbf{D}}}^{-{1}/{2}}$, and the upper bounds for their eigenvalues in (56) and (60), respectively, we can substitute the norm of these two matrices by the upper bounds of their eigenvalues and simplify the upper bound in (61) to

Based on the upper bound for the norm of the matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ in (62) and the fact that matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ is positive semidefinite, we can conclude that the eigenvalues of matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ are upper bounded by ${2(1-\delta)}/({2(1-\delta)+\alpha m})$ and the right hand side of (27) follows.

Appendix D Proof of Proposition 3

In this proof and the rest of the proofs we denote the Hessian approximation as ${\hat{\mathbf{H}}}_{t}^{-1}$ instead of ${\hat{\mathbf{H}}}_{t}^{{(K)^{-1}}}$ for simplification of equations. To prove lower and upper bounds for the eigenvalues of the error matrix ${\mathbf{E}}_{t}$ we first develop a simplification for the matrix ${\mathbf{I}}-{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1}$ in the following lemma.

Consider the NN-$K$ method as defined in (12)-(17). The matrix ${\mathbf{I}}-{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1}$ can be simplified as

Proof : Considering the definitions of the Hessian inverse approximation ${\hat{\mathbf{H}}}_{t}^{-1}$ in (15) and the matrix decomposition for the exact Hessian ${\mathbf{H}}_{t}={\mathbf{D}}_{t}-{\mathbf{B}}$, we obtain

By considering the result in (D), we simplify the expression ${\mathbf{I}}-{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1}$ as

In the right hand side of (D) the identity matrix cancels out the first term in the sum $\sum_{k=0}^{K}\left({\mathbf{B}}{\mathbf{D}}_{t}^{-1}\right)^{k}$. The remaining terms of this sum are cancelled out by the first $K$ terms of the sum $\sum_{k=0}^{K}\left({\mathbf{B}}{\mathbf{D}}_{t}^{-1}\right)^{k+1}$ so that the whole expression simplifies to $({\mathbf{B}}{\mathbf{D}}_{t}^{-1})^{K+1}$ as is claimed in (63). $\blacksquare$

Observing the fact that the error matrix ${\mathbf{E}}_{t}$ is a conjugate of the matrix ${\mathbf{I}}-{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1}$ and considering the simplification in Lemma 3 we show that the eigenvalues of error matrix ${\mathbf{E}}_{t}$ are bounded.

Proof of Proposition 3: Recall the result of Proposition 2 that all the eigenvalues of matrix ${\mathbf{D}}_{t}^{-1/2}{\mathbf{B}}{\mathbf{D}}_{t}^{-1/2}$ are uniformly bounded between and $\rho$. Since matrices ${\mathbf{D}}_{t}^{-1/2}{\mathbf{B}}{\mathbf{D}}_{t}^{-1/2}$ and ${\mathbf{B}}_{t}{\mathbf{D}}_{t}^{-1}$ are similar (conjugate) the sets of eigenvalues of these two matrices are identical. Therefore, eigenvalues of matrix ${\mathbf{B}}{\mathbf{D}}^{-1}$ are bounded as

for $i={1,2,\dots,np}$. The bounds for the eigenvalues of matrix ${\mathbf{B}}{\mathbf{D}}^{-1}$ in association with expression (63) leads to the following bounds for the eigenvalues of matrix ${\mathbf{I}}-{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1}$,

Observe that the error matrix ${\mathbf{E}}_{t}={\mathbf{I}}-{\hat{\mathbf{H}}}_{t}^{-1/2}{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1/2}$ is the conjugate of matrix ${\mathbf{I}}-{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1}$. Hence, the bounds for the eigenvalues of matrix ${\mathbf{I}}-{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1}$ also hold for the eigenvalues of error matrix ${\mathbf{E}}_{t}$ and the claim in (29) follows.

Appendix E Proof of Lemma 2

According to the Cauchy-Schwarz inequality, the product of the norms is larger than norm of the products. This observation and the definition of the approximate Hessian inverse ${\hat{\mathbf{H}}}_{t}^{-1}$ in (15) leads to

Observe that as a result of Proposition 1 the eigenvalues of matrix ${\mathbf{D}}_{t}$ are bounded below by ${2(1-\Delta)+\alpha m}$. Therefore, the maximum eigenvalue of its inverse ${\mathbf{D}}_{t}^{-1}$ is smaller than $1/(2(1-\Delta)+\alpha m)$. It then follows that the norm of the matrix ${\mathbf{D}}_{t}^{-1/2}$ is bounded above as

Based on the result in Proposition 2 the eigenvalues of the matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$ are smaller than $\rho$. Further using the symmetry and positive definiteness of the matrix ${\mathbf{D}}_{t}^{-{1}/{2}}{\mathbf{B}}{\mathbf{D}}_{t}^{-{1}/{2}}$, we obtain

Using the triangle inequality in (E) to claim that the norm of the sum is smaller than the sum of the norms and substituting the upper bounds in (69) and (70) in the resulting expression we obtain

By considering the fact that $\rho$ is smaller than $1$, the sum $\sum_{k=0}^{K}\rho^{k}$ can be simplified to $(1-\rho^{K+1})/(1-\rho)$. Considering this simplification for the sum in (71), the upper bound in (30) for the eigenvalues of the approximate Hessian inverse ${\hat{\mathbf{H}}}_{t}^{-1}$ follows.

The next step is to provide a lower bound for the eigenvalues of the Hessian inverse approximation matrix ${\hat{\mathbf{H}}}_{t}^{-1}$. In the Hessian inverse approximation formula (15), all the summands except the first one, ${\mathbf{D}}_{t}^{-1}$, are positive semidefinite. Hence, the approximate Hessian inverse ${\hat{\mathbf{H}}}_{t}^{-1}$ is the sum of matrix ${\mathbf{D}}_{t}^{-1}$ and $K$ positive semidefinite matrices and as a result we can conclude that

Proposition 1 shows that the eigenvalues of matrix ${\mathbf{D}}_{t}$ are bounded above by $2(1-\delta)+\alpha M$ which leads to the conclusion that there exits a lower bound for the eigenvalues of matrix ${\mathbf{D}}_{t}^{-1}$,

Observing the relation in (72) we realize that the lower bound for the eigenvalues of matrix ${\mathbf{D}}_{t}^{-1}$ in (73) holds for the eigenvalues of the Hessian inverse approximation ${\hat{\mathbf{H}}}_{t}^{-1}$. Therefore, all the eigenvalues of the Hessian inverse approximation ${\hat{\mathbf{H}}}_{t}^{-1}$ are greater than ${1}/({2(1-\delta)+\alpha M})$. This completes the proof of the claim in (30).

Appendix F Proof of Theorem 1

To prove global convergence of the Network Newton method we first introduce two technical lemmas. In the first lemma we use the result of Lemma 1, namely, that the objective function Hessian $\nabla^{2}F({\mathbf{y}})$ is Lipschitz continuous, to develop an upper bound for the objective function value $F({\mathbf{y}})$ using the first three terms of its Taylor expansion. In the second lemma we construct an upper bound for the objective function error at step $t+1$, namely $F({\mathbf{y}}_{t+1})-F({\mathbf{y}}^{*})$, in terms of the error at step $t$, namely $F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*})$.

Proof : Since objective function $F$ is twice differentiable, based on the Fundamental Theorem of Calculus we can write

where $\nabla F$ is the gradient of function $F$. We proceed by adding and subtracting the term $\nabla F({\mathbf{y}})^{T}({\hat{\mathbf{y}}}-{\mathbf{y}})$ to the right hand side of (75) which yields

where $\nabla^{2}F$ is the Hessian of function $F$. After setting ${\hat{\mathbf{z}}}={\mathbf{y}}+\omega({\hat{\mathbf{y}}}-{\mathbf{y}})$ and ${\mathbf{z}}={\mathbf{y}}$ in (77) and rearranging terms it follows that

Observing the fact that ${\mathbf{y}}+\omega({\hat{\mathbf{y}}}-{\mathbf{y}})-{\mathbf{y}}=\omega({\hat{\mathbf{y}}}-{\mathbf{y}})$, we can further simplify (78) to

Based on the relation for the difference of gradients $\nabla F({\mathbf{y}}+\omega({\hat{\mathbf{y}}}-{\mathbf{y}}))-\nabla F({\mathbf{y}})$ in (79), we can rewrite (76) by applying this substitution. This yields

We proceed by adding and subtracting the quadratic integral $\int_{0}^{1}\int_{0}^{1}\omega({\hat{\mathbf{y}}}-{\mathbf{y}})^{T}\nabla^{2}F({\mathbf{y}})({\hat{\mathbf{y}}}-{\mathbf{y}})dsd\omega$ to the right hand side of (80) to write

Observe that the term $({\hat{\mathbf{y}}}-{\mathbf{y}})^{T}\nabla^{2}F({\mathbf{y}})({\hat{\mathbf{y}}}-{\mathbf{y}})$ in the third summand of (F) is not a function of $\omega$ or $s$. Hence, we can move this term outside of the integral and simplify the integral to $\int_{0}^{1}\int_{0}^{1}\omega\ dsd\omega=1/2$. As a result of these observation the third summand of (F) can be replaced by $(1/2)({\hat{\mathbf{y}}}-{\mathbf{y}})^{T}\nabla^{2}F({\mathbf{y}})({\hat{\mathbf{y}}}-{\mathbf{y}})$ and we can rewrite (F) as

We proceed now to construct an upper bound for the integral in (F). Observe that according to the definition of the Euclidean norm of a matrix we have the inequality ({\hat{\mathbf{y}}}-{\mathbf{y}})^{T}\big{[}\nabla^{2}F({\mathbf{y}}+s\omega({\hat{\mathbf{y}}}-{\mathbf{y}}))-\nabla^{2}F({\mathbf{y}})\big{]}({\hat{\mathbf{y}}}-{\mathbf{y}})\leq\|\nabla^{2}F({\mathbf{y}}+s\omega({\hat{\mathbf{y}}}-{\mathbf{y}}))-\nabla^{2}F({\mathbf{y}})\|_{2}\|{\hat{\mathbf{y}}}-{\mathbf{y}}\|^{2}. By applying this inequality we have an upper bound for the integral in (F) that results in

The next step is to provide an upper bound for the term $\|\nabla^{2}F({\mathbf{y}}+s\omega({\hat{\mathbf{y}}}-{\mathbf{y}}))\!-\!\nabla^{2}F({\mathbf{y}})\|_{2}$ in the right hand side of (F). Lemma 1 shows that the penalized objective function Hessian $\nabla^{2}F$ is Lipschitz continuous with parameter $\alpha L$. Therefore, we can write

By considering (F) and substituting $\|\nabla^{2}F({\mathbf{y}}+s\omega({\hat{\mathbf{y}}}-{\mathbf{y}}))-\nabla^{2}F({\mathbf{y}})\|_{2}$ by the upper bound in (F), we obtain

Now observe that since $\|{\hat{\mathbf{y}}}-{\mathbf{y}}\|$ does not depend on $s$ or $\omega$, the integral in the last summand of (F) can be simplified as

The simplification in (F) for the last summand of (F) implies the claim in (74) is valid. $\blacksquare$

Lemma 4 shows an upper bound for the Taylor expansion of the objective function value $F({\hat{\mathbf{y}}})$. We use the result of Lemma 4 to establish an upper bound for the objective function error at step $t+1$ in terms of the error at step $t$. This result is proven in the following lemma.

Consider the NN-$K$ method as defined in (12)-(17) and the objective function $F({\mathbf{y}})$ as defined in (6). Further, recall the definition of ${\mathbf{y}}^{*}$ as the optimal argument of the objective function $F({\mathbf{y}})$. If assumptions 1, 2, and 3 hold true, the sequence of objective function value errors $\{F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*})\}$ satisfies

Proof : Recall the result of Lemma 4. By setting ${\hat{\mathbf{y}}}:={\mathbf{y}}_{t+1}$ and ${\mathbf{y}}:={\mathbf{y}}_{t}$ in (74) we obtain

where ${\mathbf{g}}_{t}:=\nabla F({\mathbf{y}}_{t})$ and ${\mathbf{H}}_{t}:=\nabla^{2}F({\mathbf{y}}_{t})$. From the definition of the NN-$K$ update formula in (16) we can write the difference of two consecutive variables as ${\mathbf{y}}_{t+1}-{\mathbf{y}}_{t}=-\epsilon{\hat{\mathbf{H}}}_{t}^{-1}{\mathbf{g}}_{t}$. Making this substitution in (88) implies

According to the definition of error matrix ${\mathbf{E}}_{t}$ in (28), we can substitute ${\hat{\mathbf{H}}}_{t}^{-1/2}{\mathbf{H}}_{t}{\hat{\mathbf{H}}}_{t}^{-1/2}$ by ${\mathbf{I}}-{\mathbf{E}}_{t}$. By making this substitution into the third summand of (F) we obtain

Proposition 3 shows that the error matrix ${\mathbf{E}}_{t}$ is always positive semidefinite. As a result, we conclude that the quadratic form ${\mathbf{g}}_{t}^{T}{\hat{\mathbf{H}}}_{t}^{-{1}/{2}}{\mathbf{E}}_{t}{\hat{\mathbf{H}}}_{t}^{-{1}/{2}}{\mathbf{g}}_{t}$ is always nonnegative. Considering this lower bound we can simplify (F) to

Note that since the stepsize is not larger than $1$, we obtain that $2\epsilon-\epsilon^{2}$ is positive. Moreover, recall the result of Lemma 2 that all the eigenvalues of the Hessian inverse approximation ${\hat{\mathbf{H}}}_{t}^{-1}$ are lower and upper bounded by $\lambda$ and $\Lambda$, respectively. These two observations imply that we can replace the term ${\mathbf{g}}_{t}^{T}{\hat{\mathbf{H}}}_{t}^{-1}{\mathbf{g}}_{t}$ by its lower bound $\lambda\|{\mathbf{g}}_{t}\|^{2}$. Moreover, existence of upper bound $\Lambda$ for the eigenvalues of Hessian inverse approximation ${\hat{\mathbf{H}}}_{t}^{-1}$ implies that the term $\|{\hat{\mathbf{H}}}_{t}^{-1}{\mathbf{g}}_{t}\|^{3}$ is upper bounded by $\Lambda^{3}\|{\mathbf{g}}_{t}\|^{3}$. Substituting these bounds for the second and third terms of (91) and subtracting the optimal objective function value $F({\mathbf{y}}^{*})$ from both sides of inequality (91) leads to

We now find lower and upper bounds for the norm of gradient $\|{\mathbf{g}}_{t}\|$ in terms of the objective function error $F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*})$. As it follows from Proposition 1, the eigenvalues of Hessian ${\mathbf{H}}_{t}$ are bounded by $\alpha m$ and $2+\alpha M$. Taking a Taylor expansion of the objective function $F({\mathbf{y}})$ around ${\mathbf{w}}$ and using the lower bound $\alpha m$ for the Hessian eigenvalues yields

For fixed ${\mathbf{w}}$, the right hand side of (93) is a quadratic function of ${\mathbf{y}}$ whose minimum argument we can find by setting its gradient to zero. Doing this yields the minimizing argument ${\hat{\mathbf{y}}}={\mathbf{w}}-(1/m)\nabla F({\mathbf{w}})$ implying that for all ${\mathbf{y}}$ we must have

The bound in (F) is true for all ${\mathbf{w}}$ and ${\mathbf{y}}$. In particular, for ${\mathbf{y}}={\mathbf{y}}^{*}$ and ${\mathbf{w}}={\mathbf{y}}_{t}$ (F) yields

Rearrange terms in (95) to obtain $2\alpha m(F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*}))$ as a lower bound for $\|\nabla F({\mathbf{y}}_{t})\|^{2}=\|{\mathbf{g}}_{t}\|^{2}$. Now substitute the lower bound $2\alpha m(F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*}))$ for squared norm of gradient $\|{\mathbf{g}}_{t}\|^{2}$ in the second summand of (F) to obtain

Notice that according to the definition of $\lambda$ in (31) we can substitute $2(1-\delta)+\alpha M$ by $1/\lambda$. Implementing this substitution and minimizing both sides of the equality with respect to ${\mathbf{y}}$ yields

Setting ${\hat{\mathbf{y}}}={\mathbf{y}}_{t}$, observing that by definition $\|\nabla F({\mathbf{y}}_{t})\|=\|{\mathbf{g}}_{t}\|$, rearranging terms, and taking the square root of both sides of the resulting inequality leads to

Replacing the upper bound in (99) for the norm of the gradient $\|{\mathbf{g}}_{t}\|$ in the last term of (F) yields the claim in (5). $\blacksquare$

We use the result of Lemma 5 to prove linear convergence of the sequence of objective function errors $F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*})$ to zero.

Proof of Theorem 1: To simplify upcoming derivations define the sequence $\beta_{t}$ as

Recall the result of Lemma 5. Factorizing $F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*})$ from the terms of the right hand side of (5) in association with the definition of $\beta_{t}$ in (100) implies that we can simplify (5) as

To prove global convergence of objective function error $F({\mathbf{y}}_{t})-F({\mathbf{y}}^{*})$ we need to show that for all time steps $t$, the constants $\beta_{t}$ are strictly smaller than $1$ and larger than , i.e., that $0<\beta_{t}<1$ for all times $t$.

We first show that $\beta_{t}$ is less than $1$ for all $t\geq 0$. To do so observe that the second term in the right and side of (100) is nonnegative. It is therefore true that

Considering the inequality $(\epsilon-1)^{2}\geq 0$ it is trivial to derive that $\epsilon(2-\epsilon)\leq 1$. Moreover, considering the facts that $m<M$ and $1-\delta>0$, we obtain $\alpha m<\alpha M+(1-\delta)$ which yields $\alpha m/(\alpha M+2(1-\delta))<1$. Considering the definition of $\lambda$ in (31) we can substitute $1/(2(1-\delta)+\alpha M)$ by $\lambda$ and write $\alpha m\lambda<1$. By multiplying these two ratios, both of which are smaller than $1$, we conclude that

That $\beta_{t}<1$ follows by combining (102) with (103).

To prove that $0<\beta_{t}$ for all $t\geq 0$ we prove that this is true for $t=0$ and then prove that the $\beta_{t}$ sequence is increasing. To show that $\beta_{0}$ is positive first note that since the stepsize $\epsilon$ satisfies the condition in (32) we can write

By computing the squares of both sides of (104), multiplying the right hand side of the resulting inequality by 2 to make the inequality strict, and factorizing ${\alpha m}\lambda$ from the term in the resulting right hand side we obtain

If we now divide both sides of the inequality in (105) by the first multiplicand in the right hand side of (105) we obtain

Observe that based on the hypothesis in (32) the step size $\epsilon$ is smaller than $1$ and it is then trivially true that $2-\epsilon\geq 1$. This observation shows that if we multiply the right hand side of (106) by $2(1-\epsilon/2)$ the inequality still holds,

Furhter multiplying both sides of inequality (107) by $\epsilon$ and rearranging terms leads to

According to the definition of $\beta_{t}$ in (100), the result in (108) implies that $\beta_{0}>0$.