Hadrien Hendrikx

Lin Xiao

Sébastien Bubeck

Francis Bach

Laurent Massoulié

Statistically Preconditioned Accelerated Gradient Method for Distributed Optimization

We consider the setting of distributed empirical risk minimization where multiple machines compute the gradients in parallel and a centralized server updates the model parameters. In order to reduce the number of communications required to reach a given accuracy, we propose a preconditioned accelerated gradient method where the preconditioning is done by solving a local optimization problem over a subsampled dataset at the server. The convergence rate of the method depends on the square root of the relative condition number between the global and local loss functions. We estimate the relative condition number for linear prediction models by studying uniform concentration of the Hessians over a bounded domain, which allows us to derive improved convergence rates for existing preconditioned gradient methods and our accelerated method. Experiments on real-world datasets illustrate the benefits of acceleration in the ill-conditioned regime.

Introduction

We consider empirical risk minimization problems of the form minimize

x∈R d Φ(x) F (x) + ψ(x), (1)
where F is the empirical risk over a dataset {z 1 , . . . , z N }:

F (x) = 1 N N i=1 (x, z i ), (2)
and ψ is a convex regularization function. We incorporate smooth regularizations such as squared Euclidean norms (λ/2) x 2 into the individual loss functions (x, z i ), and 1 INRIA, DIENS, PSL Research University, Paris, France. 2 Microsoft Research, Redmond, WA, USA.. Correspondence to: Hadrien Hendrikx <hadrien.hendrikx@inria.fr>.

Proceedings of the 37 th International Conference on Machine Learning, Online, PMLR 119, 2020. Copyright 2020 by the author(s). leave ψ mainly for non-smooth regularizations such as the 1 -norm or the indicator function of a constraint set. In modern machine learning applications, the dataset is often very large and has to be stored at multiple machines. For simplicity of presentation, we assume N = mn, where m is the number of machines and n is the number of samples stored at each machine. Let D j = {z (j) 1 , . . . , z (j) n } denote the dataset at machine j and define the local empirical risk

f j (x) = 1 n n i=1
x, z (j) i

, j = 1, . . . , m.

(3)

The overall empirical risk of Equation ( 2) can then be written as

F (x) = 1 m m j=1 f j (x) = 1 nm m j=1 n i=1
x, z

.

We assume that F is L F -smooth and σ F -strongly convex over dom ψ, in other words,

σ F I d ∇ 2 F (x) L F I d , ∀x ∈ dom ψ, (4)
where I d is the d × d identity matrix. The condition number of F is defined as κ F = L F /σ F .

We focus on a basic setting of distributed optimization where the m machines (workers) compute the gradients in parallel and a centralized server updates the variable x. Specifically, during each iteration t = 0, 1, 2, . . ., (i) the server broadcasts x t to all m machines;

(ii) each machine j computes the gradient ∇f j (x t ) and sends it back to the server;

(iii) the server forms ∇F (x t ) = 1 m m j=1 ∇f j (x t ) and uses it to compute the next iterate x t+1 .

A standard way for solving problem (1) in this setting is to implement the proximal gradient method at the server:

x t+1 = argmin x∈R d ∇F (x t ) x+ψ(x)+ 1 2η t x-x t 2 , ( 5
)
where • denotes the Euclidean norm and η t > 0 is the step size. Setting η t = 1/L F leads to linear convergence:

Φ(x t ) -Φ(x * ) ≤ 1 -κ -1 F t L F 2 x * -x 0 2 , ( 6
)
where x * = arg min Φ(x) (e.g., Beck, 2017, Section 10.6).

In other words, in order to reach Φ(x t ) -Φ(x * ) ≤ , we need O(κ F log(1/ )) iterations, which is also the number of communication rounds between the workers and the server.

If we use accelerated proximal gradient methods (e.g., [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF][START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF] at the server, then the iteration/communication complexity can be improved to O( √ κ F log(1/ )).

Statistical Preconditioning

In general, for minimizing F (x) = (1/m) m j=1 f j (x) with first-order methods, the communication complexity of O( √ κ F log(1/ )) cannot be improved [START_REF] Arjevani | Communication complexity of distributed convex learning and optimization[END_REF][START_REF] Scaman | Optimal algorithms for smooth and strongly convex distributed optimization in networks[END_REF]. However, for distributed empirical risk minimization (ERM), the additional finite-sum structure of each f j in (3) allows further improvement. A key insight here is that if the datasets D j at different workers are i.i.d. samples fr (3om the same source distribution, then the local empirical losses f j are statistically very similar to each other and to their average F , especially when n is large. Statistical preconditioning is a technique to further reduce communication complexity based on this insight.

An essential tool for preconditioning in first-order methods is the Bregman divergence. The Bregman divergence of a strictly convex and differentiable function φ is defined as D φ (x, y) φ(x) -φ(y) -∇φ(y) (x -y).

(7)

We also need the following concepts of relative smoothness and strong convexity [START_REF] Heinz | A descent lemma beyond Lipschitz gradient continuity: firstorder methods revisited and applications[END_REF][START_REF] Lu | Relatively smooth convex optimization by first-order methods, and applications[END_REF].

Definition 1. Suppose φ : R d → R is convex and twice differentiable. The function F is said to be L F/φ -smooth and σ F/φ -strongly convex with respect to φ if for all x ∈ R d ,

σ F/φ ∇ 2 φ(x) ∇ 2 F (x) L F/φ ∇ 2 φ(x). (8)
The classical definition in (4) can be viewed as relative smoothness and strong convexity where φ(x) = (1/2) x 2 . Moreover, it can be shown that (8) holds if and only if for all x, y ∈ R d σ F/φ D φ (x, y) ≤ D F (x, y) ≤ L F/φ D φ (x, y). (9)

Consequently, we define the relative condition number of F with respect to φ as κ F/φ = L F/φ /σ F/φ .

Following the Distributed Approximate Newton (DANE) method by [START_REF] Shamir | Communication-efficient distributed optimization using an approximate Newton-type method[END_REF], we construct the reference function φ by adding some extra regularization to one of the local loss functions (say f 1 , without loss of generality):

φ(x) = f 1 (x) + µ 2 x 2 . ( 10
)
Then we replace (1/2) x -x t 2 in the proximal gradient method (5) with the Bregman divergence of φ, i.e.,

x t+1 = argmin x∈R d ∇F (x t ) x+ψ(x)+ 1 η t D φ (x, x t ) . (11)
In this case, worker 1 acts as the server to compute x t+1 , which requires solving a nontrivial optimization problem involving the local loss function f 1 .

According to [START_REF] Shamir | Communication-efficient distributed optimization using an approximate Newton-type method[END_REF] and [START_REF] Lu | Relatively smooth convex optimization by first-order methods, and applications[END_REF], with η t = 1/L F/φ , the sequence {x t } generated by (11) satisfies

Φ(x t ) -Φ(x * ) ≤ 1 -κ -1 F/φ t L F/φ D φ (x * , x 0 ), (12)
which is a direct extension of (6). Therefore, the effectiveness of preconditioning hinges on how much smaller κ F/φ is compared to κ F . Roughly speaking, the better f 1 or φ approximates F , the smaller κ F/φ (≥ 1) is. In the extreme case of f 1 ≡ F (with only one machine m = 1), we can choose µ = 0 and thus φ ≡ F , which leads to κ F/φ = 1, and we obtain the solution within one step.

In general, we choose µ to be an upper bound on the spectral norm of the matrix difference ∇ 2 f 1 -∇ 2 F . Specifically, we assume that with high probability, for the operator norm between matrices (i.e., the largest singular value),

∇ 2 f 1 (x) -∇ 2 F (x) ≤ µ, ∀ x ∈ dom ψ, (13)
which implies (Zhang and Xiao, 2018, Lemma 3),

σ F σ F + 2µ ∇ 2 φ(x) ∇ 2 F (x) ∇ 2 φ(x). (14)
Now we invoke a statistical argument based on the empirical average structure in (3). Without loss of generality, we assume that D 1 contains the first n samples of {z 1 , . . . , z N } and thus

∇ 2 f 1 (x) = 1 n n i=1 ∇ 2 (x, z i ).
For any fixed x, we can use Hoeffding's inequality for matrices [START_REF] Tropp | An introduction to matrix concentration inequalities[END_REF] to obtain, with probability 1 -δ,

1 n n i=1 ∇ 2 (x, z i )-∇ 2 F (x) ≤ 32L 2 log(d/δ) n , ( 15
)
where L is the uniform upper bound on ∇ 2 (x, z i ) .

If the losses (x, z i ) are quadratic in x, then the Hessians are constant and (13) holds with µ = O(L / √ n), hiding the factor log(d/δ). In this case, we derive from (14) that

κ F/φ = 1 + 2µ σ F = 1 + O κ √ n , (16)
where we assume σ F ≈ σ , where ∇ 2 (x, z i ) σ I d for all x. Therefore, for large n, whenever we have κ F/φ < κ F , the communication complexity O(κ F/φ log(1/ )) is better than without preconditioning.

For non-quadratic loss functions, we need to ensure that (13) holds uniformly over a compact domain with high probability. Standard ball-packing arguments encounter an additional factor of √ d (e.g., Zhang and Xiao, 2018, Lemma 6). In this case, we have µ = O(L d/n) and

κ F/φ = 1 + 2µ σ F = 1 + O κ √ d √ n , (17)
which suggests that the benefit of preconditioning may degrade or disappear in high dimension.

Contributions and Outline

In this paper, we make the following two contributions.

First, we propose a Statistically Preconditioned Accelerated Gradient (SPAG) method that can further reduce the communication complexity. Accelerated methods with O( √ κ F/φ log(1/ )) complexity have been developed for quadratic loss functions (see related works in Section 2). However, [START_REF] Dragomir | Optimal complexity and certification of Bregman first-order methods[END_REF] have shown that acceleration is not possible in general in the relatively smooth and strongly convex setting, and that more assumptions are needed. Here, by leveraging the fact the reference function φ itself is smooth and strongly convex, we obtain

Φ(x t ) -Φ(x * ) ≤ t τ =1 1 - 1 κ F/φ G τ L F/φ D φ (x * , x 0 ),
where 1 ≤ G t ≤ κ φ and G t → 1 geometrically. Moreover, G t can be calculated at each iteration and serve as numerical certificate of the actual convergence rate. In all of our experiments, we observe G t ≈ 1 even in early iterations, which results in O( √ κ F/φ log(1/ )) iterations empirically.

Second, we derive refined bounds on the relative condition number for linear prediction models. Linear models such as logistic regression have the form (x, z i ) = i (a i x) + (λ/2) x 2 . Assume that i (a i x) ≤ 1 and a i ≤ R for all i, which implies L = R 2 and κ = R 2 /λ. Then the Hoeffding bounds in ( 16) for quadratics becomes κ

F/φ = 1 + O R 2 √
nλ , and for nonquadratics, the bound in (17

) (from previous work) becomes κ F/φ = 1 + O R 2 √ d √ nλ .
We show that:

• For quadratic losses, the bound on relative condition number can be improved by a factor of √ n, i.e.,

κ F/φ = 3 2 + O R 2 nλ log d δ .
• For non-quadratic losses, we derive a uniform concentration bound to remove the dependence of κ F/φ on d,

κ F/φ = 1 + O R 2 √ nλ RD + log(1/δ) ,
where D is the diameter of dom φ (bounded domain).

We also give a refined bound when the inputs a i are sub-Gaussian.

These new bounds on κ F/φ improve the convergence rates for all existing accelerated and non-accelerated preconditioned gradient methods (see related work in Section 2).

We start by discussing related work in Section 2. In Section 3, we introduce SPAG and give its convergence analysis. In Section 4, we derive sharp bounds on the relative condition number, and discuss their implications on the convergence rates of SPAG and other preconditioned gradient methods. We present experimental results in Section 5.

Related Work

Shamir et al. ( 2014) considered the case ψ ≡ 0 and introduced the statistical preconditioner (10) in DANE. Yet, they define a separate φ j (x) = f j (x) + (µ/2) x 2 for each worker j, compute m separate local updates using (11), and then use their average as x t+1 . For quadratic losses, they obtain the communication complexity O((κ 2 /n) log(1/ )), which is roughly O(κ 2 F/φ log(1/ )) in our notation, which is much worse than their result without averaging of O(κ F/φ log(1/ )) given in Section 1.1. We further improve this to O( √ κ F/φ log(1/ )) using acceleration. [START_REF] Zhang | DiSCO: Distributed optimization for self-concordant empirical loss[END_REF] proposed DiSCO, an inexact damped Newton method, where the Newton steps are computed by a distributed conjugate gradient method with a similar preconditioner as (10). They obtain a communication complexity of O(( √ κ /n (Lin et al., 2015) to accelerate DANE; their method, called AIDE, achieves the same improved complexity for quadratic functions. We obtain similar results for smooth convex functions using direct acceleration. [START_REF] Wang | Utilizing second order information in minibatch stochastic variance reduced proximal iterations[END_REF] apply statistical preconditioning to speed up a mini-batch variant of SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF], but they rely on generic Catalyst acceleration and their convergence results only hold for a very small ball around the optimum.

Distributed optimization methods that use dual variables to coordinate solutions to local subproblems include ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] and CoCoA [START_REF] Jaggi | Communication-efficient distributed dual coordinate ascent[END_REF][START_REF] Ma | Adding vs. averaging in distributed primal-dual optimization[END_REF][START_REF] Beck | First-Order Methods in Optimization[END_REF]. Numerical experiments demonstrate that they benefit from statistical similarities of local functions in the early iterations [START_REF] Yuan | On convergence of distributed approximate Newton methods: Globalization, sharper bounds and beyond[END_REF], but their established communication complexity is no better than O(κ F log(1/ )).

The SPAG Algorithm

Although our main motivation in this paper is distributed optimization, the SPAG algorithm works in the general setting of minimizing relatively smooth and strongly convex functions. In this section, we first present SPAG in the more general setting (Algorithm 1), then explain how to run it for distributed empirical risk minimization.

In the general setting, we consider convex optimization problems of the form (1), where ψ is a closed convex function and F satisfies the following assumption. Assumption 1. F is L F -smooth and σ F -strongly convex.

In addition, it is L F/φ -smooth and σ F/φ -strongly convex with respect to a differentiable convex function φ, and φ itself is L φ -smooth and σ φ -strongly convex.

Algorithm 1 requires an initial point x 0 ∈ dom ψ and two parameters L F/φ and σ F/φ . During each iteration, Line 6 finds a t+1 > 0 by solving a quadratic equation, then Line 7 calculates three scalars α t , β t and η t , which are used in the later updates for the three vectors y t , v t+1 and x t+1 . The function V t (•) being minimized in Line 10 is defined as

V t (x) = η t ∇F (y t ) x + ψ(x) + (1 -β t )D φ (x, v t ) + β t D φ (x, y t ). ( 18
)
The inequality that needs to be satisfied in Line 12 is

D φ (x t+1 , y t ) (19) ≤ α 2 t G t (1 -β t )D φ (v t+1 , v t ) + β t D φ (v t+1 , y t ) ,
where G t is a scaling parameter depending on the properties of D φ . It is a more flexible version of the triangle scaling gain introduced by [START_REF] Hanzely | Accelerated Bregman proximal gradient methods for relatively smooth convex optimization[END_REF].

Algorithm 1 SPAG(L F/φ , σ F/φ , x 0 )

1: v 0 = x 0 , A 0 = 0, B 0 = 1, G -1 = 1 2: for t = 0, 1, 2, . . . do 3: G t = max{1, G t-1 /2}/2 4: repeat 5: G t ← 2G t 6: Find a t+1 such that a 2 t+1 L F/φ G t = A t+1 B t+1 where A t+1 = A t + a t+1 , B t+1 = B t + a t+1 σ F/φ 7: α t = at+1 At+1 , β t = at+1 Bt+1 σ F/φ , η t = at+1
Bt+1 8:

y t = 1 1-αtβt (1 -α t )x t + α t (1 -β t )v t 9:
Compute ∇F (y t ) (requires communication)

10:

v t+1 = arg min x V t (x)
11:

x t+1 = (1 -α t )x t + α t v
G t = κ φ = L φ /σ φ .
Therefore, if φ = (1/2) • 2 , then we can set G t = 1 and there is no need to check (19). In general, Algorithm 1 always produces G t < 2κ φ for all t ≥ 0. Following the argument from Nesterov (2013, Lemma 4), the total number of gain-searches performed up to iteration t is bounded by

2(t + 1) + log 2 (G t ),
which also bounds the total number of gradient evaluations. Thus the overhead is roughly twice as if there were no gainsearch. Next we present a convergence theorem for SPAG.

Theorem 1. Suppose Assumption 1 holds. Then the sequences generated by SPAG satisfy for all t ≥ 0,

Φ(x t ) -Φ(x * ) + σ F/φ D φ (x * , v t ) ≤ 1 A t D φ (x * , v 0 ),
where

A t = 1 4σ F /φ t-1 τ =0 (1 + γ τ ) - t-1 τ =0 (1 -γ τ ) 2 ,
and

γ t = 1 2 √ κ F /φ Gt .
The proof of Theorem 1 relies on the techniques of [START_REF] Nesterov | Efficiency of the accelerated coordinate descent method on structured optimization problems[END_REF], and the details are given in Appendix A. We can estimate the convergence rate as follows:

1 A t = O t τ =0 1- 1 √ κ F /φ Gτ = O 1- 1 κ F /φ Gt t , where G t is such that G -1/2 t = (1/t) t τ =0 G -1/2 t , that is, G 1/2 t is the harmonic mean of G 1/2 0 , . . . , G 1/2
t-1 . In addition, it can be shown that A t ≥ t 2 /(4L F/φ G t ). Therefore, as σ F/φ → 0, Theorem 1 gives an accelerated sublinear rate:

Φ(x t ) -Φ(x * ) ≤ 4L F/φ G t t 2 D φ (x * , x 0 ).
To estimate the worst case when σ F/φ > 0, we replace G t by κ φ to obtain the iteration complexity

O √ κ F/φ κ φ log(1/ ) . Since κ F/φ κ φ ≈ κ F , this is roughly O √ κ F log(1/ )
, the same as without preconditioning. However, the next lemma shows that under a mild condition, we always have G t → 1 geometrically.

Lemma 2. Suppose Assumption 1 holds and in addition, ∇ 2 φ is M -Lipschitz-continuous, i.e., for all x, y ∈ dom ψ,

∇ 2 φ(x) -∇ 2 φ(y) ≤ M x -y .
Then the inequality (19) holds with

G t = min κ φ , 1 + (M/σ φ )d t , (20)
where

d t = v t+1 -v t + v t+1 -y t + x t+1 -y t .
In particular, if φ is quadratic, then we have M = 0 and G t = 1 always satisfies (19). In this case, the convergence rate in Theorem 1 satisfies

1/A t = O 1 -1/ √ κ F/φ t .
In general, M = 0, but it can be shown that the sequences generated by Algorithm

O 1 -1/ √ κ F/φ t .
The asymptotic nature of the preconditioned convergence rate in the nonquadratic case (G t converges to 1 instead of being a small constant) seems to be unavoidable, given the recent work on lower bounds for mirror descent methods by [START_REF] Dragomir | Optimal complexity and certification of Bregman first-order methods[END_REF].

Implementation for Distributed Optimization

In distributed optimization, Algorithm 1 is implemented at the server. During each iteration, communication between the server and the workers only happens when computing ∇F (y t ). Checking if the inequality ( 19) holds locally requires that the server has access to the preconditioner φ.

If the datasets on different workers are i.i.d. samples from the same source distribution, then we can use any f j in the definition of φ in (10) and assign worker j as the server. However, this is often not the case in practice and obtaining i.i.d. datasets on different workers may involve expensive shuffling and exchanging large amount of data among the workers. In this case, a better alternative is to randomly sample small portions of the data on each worker and send them to a dedicated server. We call this sub-sampled dataset D 0 and the local loss at the server f 0 , which is defined the same way as in (3). Then the server implements Algorithm 1 with φ(x) = f 0 (x) + (µ/2) x 2 . Here we only need D 0 be a uniform sub-sample of ∪ m j=1 D j , which is critical for effective preconditioning. On the other hand, it is not a problem at all if the datasets at the workers, D 1 , . . . , D m , are not shuffled to to be i.i.d., because it does not change the average gradients ∇F (y t ). In the rest of the paper, we omit the subscript to simply use f to represent the local empirical loss function. As discussed in Section 1.1, if

∇ 2 f (x) -∇ 2 F (x) ≤ µ, ∀ x ∈ dom ψ (21)
with high probability, then according to ( 14), we can choose

L F/φ = 1, σ F/φ = σ F σ F + 2µ
as the input to Algorithm 1. In the next section, we leverage matrix concentration bounds to estimate how µ varies with the number of subsamples n. With sufficiently large n, we can make µ small so that the relative condition number κ F/φ = 1 + 2µ/σ F is much smaller than κ F .

Bounding the Relative Condition Number

In this section, we derive refined matrix concentration bounds for linear prediction models. Suppose the overall dataset consists of N samples {z 1 , . . . , z N }, where each z i = (a i , b i ) with a i ∈ R d being a feature vector and b i the corresponding label or regression target. Linear models (including logistic and ridge regression) have the form (x,

z i ) = i (a i x) + λ 2 x 2
, where i is twice differentiable and may depend on b i , and λ > 0. We further assume that i = j for all i and j, which is valid for logistic and ridge regression as well. Since

f (x) = (1/n) n i=1 (x, z i ), we have ∇ 2 f (x) = 1 n n i=1 i (a i x)a i a i + λI d . ( 22
)
Here we omit the subscript j in f j since we only need one subsampled dataset at the server, as explained in Section 3.1.

For the overall loss function defined in (2), the Hessian ∇ 2 F (x) is defined similarly by replacing n with N .

We assume for simplicity that the strong convexity of F mainly comes from regularization, that is, σ F = σ = λ, but the results can be easily extended to account for the strong convexity from data. We start by showing tight results for quadratics, and then provide uniform concentration bounds of Hessians for more general loss functions. Finally, we give a refined bound when the a i 's are sub-Gaussian.

Quadratic Case

We assume in this section that i (a i x) = (a i x -b i ) 2 /2, and that there exists a constant R such that a i ≤ R for all i = 1, . . . , N . In this case we have L = R 2 and κ = R 2 /λ. Since the Hessians do not depend on x, we use the notation

H F = ∇ 2 F (x), H f = ∇ 2 f (x).
Previous works [START_REF] Shamir | Communication-efficient distributed optimization using an approximate Newton-type method[END_REF][START_REF] Sashank | AIDE: Fast and communication efficient distributed optimization[END_REF][START_REF] Yuan | On convergence of distributed approximate Newton methods: Globalization, sharper bounds and beyond[END_REF] use the Hoeffding bound ( 15) to obtain

1 + 2µ λ -1 (H f + µI d ) H F H f + µI d , (23) with µ = R 2 √ n 32 log(d/δ). ( 24
)
Our result is given in the following theorem.

Theorem 2. Suppose i is quadratic and a i ≤ R for all i.

For a fixed δ > 0, if n > 28 3 log 2d δ , then the following inequality holds with probability at least 1 -δ:

3 2 + 2µ λ -1 (H f + µI d ) H F 2 (H f + µI d ) , ( 25
)
with µ = 1 2 28R 2 3n log 2d δ -λ + . ( 26
)
Thus, for this choice of µ, σ

F/φ = 3 2 + 2µ λ -1 , L F/φ = 2 and so κ F/φ = O 1 + κ n log d δ with probability 1 -δ.
Theorem 2 improves on the result in ( 24) by a factor of √ n. The reason is that matrix inequality ( 23) is derived from the additive bound H f -H F ≤ µ (e.g., [START_REF] Shamir | Communication-efficient distributed optimization using an approximate Newton-type method[END_REF][START_REF] Yuan | On convergence of distributed approximate Newton methods: Globalization, sharper bounds and beyond[END_REF]. We derive the matrix inequality (25) directly from a multiplicative bound using the matrix Bernstein inequality (see proof in Appendix B.1). Note that by using matrix Bernstein instead of matrix Hoeffding inequality [START_REF] Tropp | An introduction to matrix concentration inequalities[END_REF], one can refine the bound for µ in (23) from L / √ n to L L F /n, which can be as small as L /n in the extreme case when all the a i 's are orthogonal. Our bound in (26) states that µ = O(L /n) in general for quadratic problems, leading to κ

F/φ = O(1 + κ /n).
Remark 1. Theorem 2 is proved by assuming random sampling with replacement. In practice, we mostly use random sampling without replacement, which usually concentrates even more than with replacement [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF].

Remark 2. In terms of reducing κ F/φ , there is not much benefit to having µ < λ. Indeed, higher values of µ regularize the inner problem of minimizing V t (x) in ( 18), because the condition number of

D φ (x, y) = D f (x, y)+(µ/2) x- y 2 is (L f + µ)/(λ + µ).
Increasing µ can thus lead to substantially easier subproblems when µ > λ, which reduces the computation cost at the server, although this may sometimes affect the rate of convergence.

Non-quadratic Case

For non-quadratic loss functions, we need ∇ 2 f (x) to be a good approximation of ∇ 2 F (x) for all iterations of the SPAG algorithm. It is tempting to argue that concentration only needs to hold for the iterates of SPAG, and a union bound would then give an extra log T factors for T iterations. Yet this only works for one step since x t depends on the points chosen to build f for t > 0, so the (a i x t )a i a i are not independent for different i (because of x t ). Therefore, the concentration bounds need to be written at points that do not depend on f . In order to achieve this, we restrict the optimization variable within a bounded convex set and prove uniform concentration of Hessians over the set. Without loss of generality, we consider optimization problems constrained in B(0, D), the ball of radius D centered at 0. Correspondingly, we set the nonsmooth regularization function as ψ(x) = 0 if x ∈ B(0, D) and infinity otherwise.

If the radius D is small, it is then possible to leverage the quadratic bound by using the inequality

H f (x) -H F (x) ≤ H f (x) -H f (y) + H f (y) -H F (y) + H F (x) -H F (y) .
Thus, under a Lipschitz-continuous Hessian assumption (which we have), only concentration at point y matters. Yet, such bounding is only meaningful when x is close to y, thus leading to the very small convergence radius of Wang and Zhang (2019, Theorem 13), in which they use concentration at the optimal point x * . Using this argument for several y's that pave B(0, D) leads to an extra √ d multiplicative factor since concentration needs to hold at exponentially (in d) many points, as discussed in Section 1.1. We take a different approach in this work, and proceed by directly bounding the supremum for all x ∈ B(0, D), thus looking for the smallest µ that satisfies:

sup x∈B(0,D) H f (x) -H F (x) op ≤ µ. ( 27
)
Equation ( 23) can then be used with this specific µ. We now introduce Assumption 2, which is for example verified for logistic regression with B = 1/4 and M = 1.

Assumption 2. There exist B and M such that i is M -Lipschitz continuous and 0 ≤ i (a x) ≤ B almost surely for all x ∈ B(0, D).

Theorem 3. If i satisfies Assumption 2, then Equation ( 27) is satisfied with probability at least 1 -δ for

µ = √ 4π R 2 √ n B 2 + 1 2π log(δ -1 ) + RM D .
Sketch of proof. The high probability bound on the supremum is obtained using Mc Diarmid inequality [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF]. This requires a bound on its expectation, which is obtained using symmetrization and the Sudakov-Fernique Lemma [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF]. The complete proof can be found in Appendix B.2.

The bound of Theorem 3 is relatively tight as long as RM D < B log(δ -1 ). Indeed, using the matrix Bernstein inequality for a fixed x ∈ B(0, D) would yield

µ = O R √ L F B log(d/δ)/ √ n .
Therefore, Theorem 3 is tight up to a factor R/ √ L F in this case. Remark 3. We consider D to be fixed in this work, although obtaining a meaningful solution to the ERM problem may require D to depend on the dimension d. Yet, D would actually depend on the intrinsic dimension of the data, which can be much smaller than d, especially when features are sparse.

Sub-Gaussian Bound

We show in this section that the bound of Theorem 3 can be improved under a stronger sub-Gaussian assumption on a. Definition 2. The random variable a ∈ R d is sub-Gaussian with parameter ρ > 0 if one has for all > 0, x ∈ B(0, D):

P(|a i x)| ≥ ρ ) ≤ 2e - 2 2 x 2 . ( 28
)
Theorem 4. If i satisfies Assumption 2 and the a i are sub-Gaussian with constant ρ, then denoting B = B /(M D), there exists C > 0 such that Equation ( 27) is satisfied with probability 1 -δ for

µ = C ρ 2 M D √ n (d + log(δ -1 )) ρ + B √ d + ρ + (R 2 B) 1 3 √ n .
Recall that this value of µ can be plugged into Equation ( 23) to bound the relative condition number.

Sketch of proof. This bound is a specific instantiation of a more general result based on chaining, which is a standard argument for proving results on suprema of empirical processes [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF]. The complete proof can be found in Appendix B.3.

The sub-Gaussian assumption (28) always holds with ρ = R, the almost sure bound on a i . However Theorem 4 improves over Theorem 3 only with a stronger sub-Gaussian assumption, i.e., when ρ < R. In particular for a i uniform over B(0, R),

one has ρ = R/ √ d. Assuming fur- ther that the (R 2 B) 1/3 / √ n term dominates yields µ = O(R 2 (R 2 B) 1/3 /n), a
√ n improvement over Theorem 3. We expect tighter versions of Theorem 4, involving the effective dimension d eff of vectors a i instead of the full dimension d, to hold.

Experiments

We have seen in the previous section that preconditioned gradient methods can outperform gradient descent by a large margin in terms of communication rounds, which was already observed empirically [START_REF] Shamir | Communication-efficient distributed optimization using an approximate Newton-type method[END_REF][START_REF] Sashank | AIDE: Fast and communication efficient distributed optimization[END_REF][START_REF] Yuan | On convergence of distributed approximate Newton methods: Globalization, sharper bounds and beyond[END_REF]. We compare in this section the performances of SPAG with those of DANE and its heavyball acceleration, HB-DANE [START_REF] Yuan | On convergence of distributed approximate Newton methods: Globalization, sharper bounds and beyond[END_REF], as well as accelerated gradient descent (AGD). Due to its better convergence guarantees [START_REF] Shamir | Communication-efficient distributed optimization using an approximate Newton-type method[END_REF][START_REF] Yuan | On convergence of distributed approximate Newton methods: Globalization, sharper bounds and beyond[END_REF], DANE refers in this section to the proximal gradient method with the Bregman divergence associated to φ = f 1 + (µ/2) • 2 (without averaging over m workers).

We apply these algorithms to train linear prediction models over a dataset

{(a i , b i )} N i=1
, where each a i ∈ R d is a feature vector and b i is the corresponding label or regression target. Specifically, we solve the problem

minimize x∈R d F (x) = 1 N n i=1 i (a i x) + λ 2 x 2 , where i (a i x) = log(1 + exp(-b i (a i x))) for logistic re- gression with b i ∈ {-1, +1} and i (a i x) = (a i x -b i ) 2
for ridge regression with b i ∈ R. We use two datasets from LibSVM 1 , RCV1 [START_REF] David D Lewis | RCV1: A new benchmark collection for text categorization research[END_REF] and the preprocessed version of KDD2010 (algebra) [START_REF] Yu | Feature engineering and classifier ensemble for KDD cup 2010[END_REF].

1 Accessible at https://www.csie.ntu.edu.tw/ ˜cjlin/libsvmtools/datasets/binary.html Note that, as mentioned in Section 3.1, the number of nodes used by SPAG does not affect its iteration complexity (but change the parallelism of computing ∇F (x t )). Only the size n of the dataset used for preconditioning matters. We initialize all algorithms at the same point, which is the minimizer of the server's entire local loss (with 10 5 samples regardless of how many samples are used for preconditioning). Note that this is not mandatory but provides a good communication-free initialization, although it requires partially solving the server's local loss. Not performing this would simply yield a worse initialization.

Tuning µ. Although µ can be estimated using concentrations results, as done in Section 4, these bounds are too loose to be used in practice. Yet, they show that µ depends very weakly on λ. This is verified experimentally, and we therefore use the same value for µ regardless of λ. To test the impact of µ on the iteration complexity, we fix a step-size of 1 and plot the convergence speed of SPAG for several values of µ. We see on Figure 1 that the value of µ drastically affects convergence, actually playing a role similar to the inverse of a step-size. Indeed, the smaller the µ the faster the convergence, up to a point at which the algorithm is not stable anymore. Convergence could be obtained for smaller values of µ by taking a smaller step-size. Yet, the step-size needs to be tuned for each value of µ, and we observed that this does not lead to significant improvements in practice. Thus, unless explicitly stated, we stick to the guidelines for DANE by [START_REF] Shamir | Communication-efficient distributed optimization using an approximate Newton-type method[END_REF], i.e., we choose L F/φ = 1 and tune µ. The tuning strategy for µ is the following: we tune the base value of µ by starting from 0.1/n for the smallest n and then decreasing it as long as it is stable, or increasing it as long as it is unstable. Then, we keep this base value and obtain µ for other values of n by relying on the fact that µ should be proportional to 1/n (slightly adjusting when necessary if the algorithm becomes unstable). Therefore, even though some tuning is required to set µ, this tuning is guided by the insight provided by Section 4.

Setting acceleration parameters. SPAG and HB-DANE require additional parameters compared with DANE. Yet, we use in our experiments the values given by the theory, i.e., we use SPAG with σ -1 F/φ = 1 + 2µ/λ and HB-DANE with β = (1 -(1 + 2µ/λ) -1/2 ) 2 . Fine tuning these parameters only leads to small improvements for both algorithms, as described in Appendix C. Therefore, SPAG and HB-DANE do not require more parameter tuning than DANE. We tune both the learning rate and the momentum of AGD.

Line search for G t . As explained in Section 3, the optimal G t is obtained through a line search. Yet, we observed in all our experiments that G t = 1 most of the time. This is due to the fact that we start at the minimizer of the local cost function, which can be close to the global solution. In addition, Equation (20) can actually be verified for G t < 1, even in the quadratic. Therefore, the line search generally has no added cost (apart from checking that G t = 1 works) and the effective rate in our experiments is κ -1/2 F/φ . Experiments for Figures 2 and4 use G t = 1 for simplicity.

Local subproblems. Local problems are solved using a sparse implementation of dual-free SDCA [START_REF] Shalev-Shwartz | Sdca without duality, regularization, and individual convexity[END_REF]. In practice, the ill-conditionned regime is very hard, especially when µ is small. Indeed, the local subproblems are very hard to solve, and it should be beneficial to use accelerated algorithms to solve the inner problems. In our experiments, we warm-start the local problems (initializing F/φ = 0.9 as it proved to be more stable.

on the solution of the previous one), and keep doing passes over the preconditioning dataset until ∇V t (x t ) ≤ 10 -9 (checked at each epoch).

RCV1. Figure 2 presents results for logistic regression on the RCV1 dataset with different regularization weights. All algorithms are run with N = 677399 (split over 4 nodes) and d = 47236. We see that in Figure 2(a), the curves can be clustered by values of n, meaning that when regularization is relatively high (λ = 10 -5 ), increasing the preconditioning sample size has a greater effect than acceleration since the problem is already well-conditioned. In particular, acceleration does not improve the convergence rate when n = 10 5 and λ = 10 -5 . When regularization is smaller (λ = 10 -7 ), SPAG and HB-DANE outperform DANE even when ten times less samples are used for preconditioning, as shown in Figure 2(b). As discussed in Appendix C, finer tuning (without using the theoretical parameters) of the momentum marginally improves the performances of SPAG and HB-DANE, at the cost of a grid search. SPAG generally outperforms HB-DANE in our experiments, but both methods have comparable asymptotic rates.

Figure 3 presents results for ridge regression on the same RCV1 dataset, using the class labels as the regression targets. In this case, φ is quadratic so we do not need line search on G t and always set G t = 1. The results are very similar to the logistic regression case. samples at the server for a given convergence speed. AGD competes with DANE when λ and n are small, but it is outperformed by SPAG in all our experiments. More experiments investigating the impact of line search, tuning and inaccurate local solutions are presented in Appendix C.

KDD2010.

Conclusion

We have introduced SPAG, an accelerated algorithm that performs statistical preconditioning for large-scale distributed optimization. Although our motivation in this paper is distributed empirical risk minimization, SPAG applies to much more general settings that can benefit from statistical preconditioning. We have given tight bounds on the relative condition number, a crucial quantity to understand the convergence rate of preconditioned algorithms. We have also shown, both in theory and in experiments, that acceleration allows SPAG to efficiently leverage rough preconditioning when only limited number of local samples are available. Preliminary experiments suggest that SPAG is more robust to inaccurate solution of the inner problems than HB-DANE. Characterizing the effects of inaccurate inner solutions in the preconditioning setting would be an interesting extension of this work.

Note that SPAG (Algorithm 1) can be considered as an accelerated variant of the general mirror descent method considered by Bauschke et al. ( 2017) and [START_REF] Lu | Relatively smooth convex optimization by first-order methods, and applications[END_REF]. Specifically, we can replace D φ by the Bregman divergence of any convex function of Legendre type (Rockafellar, 1970, Section 26). Recently, [START_REF] Dragomir | Optimal complexity and certification of Bregman first-order methods[END_REF] show that fully accelerated convergence rates, as those for Euclidean mirror-maps achieved by [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF], may not be attainable in the general setting. However, this negative result does not prevent us from obtaining better accelerated rates in the preconditioned setting. Indeed, we choose a smooth and strongly convex mirror map and further assume Lipschitz continuity of its Hessian. For smooth and strongly convex cost functions, the convergence rates of SPAG are almost always better than those obtained by standard accelerated algorithms (without preconditioning) as long as n is not too small, and can be much better with a good preconditioner.

A.1. Proof of Lemma 1

Using the second-order Taylor expansion (mean-value theorem), we have

D φ (x, y) = φ(x) -φ(y) -∇φ(y), x -y = 1 2 (x -y) ∇ 2 φ y + t(x -y) (x -y),
for some scalar t ∈ [0, 1]. We define

H(x, y) = ∇ 2 φ y + t(x -y) ,
where the dependence on t ∈ [0, 1] is made implicit with the ordered pair (x, y). Then we can write

D φ (x, y) = 1 2 x -y 2 H(x,y) .
By Assumption 1, φ is L φ -smooth and σ φ -strongly convex, which implies that for all x, y ∈ R d ,

σ φ x -y 2 ≤ x -y 2 H(x,y) ≤ L φ x -y 2 .
Let w t = (1 -β t )v t + β t y t . Then we have x t+1 -y t = α t v t+1 -w t and

D φ (x t+1 , y t ) = 1 2 x t+1 -y t 2 H(xt+1,yt) ≤ L φ 2 x t+1 -y t 2 = α 2 t L φ 2 v t+1 -w t 2 . Next we use v t+1 -w t = (1 -β t )(v t+1 -v t ) + β t (v t+1 -y t ) and convexity of • 2 to obtain D φ (x t+1 , y t ) ≤ α 2 t L φ 2 (1 -β t ) v t+1 -v t 2 + β t v t+1 -y t 2 ≤ α 2 t L φ 2σ φ (1 -β t ) v t+1 -v t 2 H(vt+1,vt) + β t v t+1 -y t 2 H(vt+1,yt) = α 2 t κ φ (1 -β t )D φ (v t+1 , v t ) + β t D φ (v t+1 , y t ) .
This finishes the proof of Lemma 1.

A.2. Proof of Theorem 1

Theorem 2 is a direct consequence of the following result, which is adapted from [START_REF] Nesterov | Efficiency of the accelerated coordinate descent method on structured optimization problems[END_REF].

Theorem 5 (Smooth and strongly convex mirror map φ). Suppose Assumption 1 holds. Then the sequences generated by Algorithm 1 satisfy for all t ≥ 0,

A t Φ(x t ) -Φ(x * ) + B t D(x * , v t ) ≤ A 0 F (x 0 ) -F (x * ) + B 0 D(x * , v 0 ).
Moreover, if we set A 0 = 0 and B 0 = 1 then for t ≥ 0,

A t ≥ 1 4σ F/φ π + t -π - t 2 , B t = 1 + σ F/φ A t ≥ 1 4 π + t + π - t 2 ,
where

π + t = t-1 i=0 1 + σ F/φ L F/φ G t , π - t = t-1 i=0 1 - σ F/φ L F/φ G t .
We first state an equivalent definition of relative smoothness and relative strong convexity [START_REF] Lu | Relatively smooth convex optimization by first-order methods, and applications[END_REF]. The function F is said to be L F/φ -smooth and σ F/φ -strongly convex with respect to φ if for all x, y ∈ R d ,

F (y) + ∇F (y) (x -y) + σ L/φ D φ (x, y) ≤ F (x) ≤ F (y) + ∇F (y) (x -y) + L L/φ D φ (x, y). ( 29
)
Obviously this is the same as ( 9). We also need the following lemma, which is an extension of a result from Chen and Teboulle (1993, Lemma 3.2), whose proof we omit.

Lemma 3 (Descent property of Bregman proximal point). Suppose g is a convex function defined over dom φ and

v t+1 = argmin x g(x) + (1 -β t )D φ (x, v t ) + β t D φ (x, y t ) , then for any x ∈ dom h, g(v t+1 ) + (1 -β t )D φ (v t+1 , v t ) + β t D φ (v t+1 , y t ) ≤ g(x) + (1 -β t )D φ (x, v t ) + β t D φ (x, y t ) -D φ (x, v t+1 ).
Proof of Theorem 5. The proof follows the same lines as [START_REF] Nesterov | Efficiency of the accelerated coordinate descent method on structured optimization problems[END_REF], with adaptations to use general Bregman divergences. Applying Lemma 3 with g(x) = η t ∇f (y t ) x + ψ(x) , we have for any x ∈ dom φ,

D(x, v t+1 ) + (1 -β t )D(v t+1 , v t ) + β t D(v t+1 , y t ) -(1 -β t )D(x, v t ) -β t D(x, y t ) ≤ η t ∇f (y t ) (x -v t+1 ) + η t ψ(x) -ψ(v t+1 ) .
Since by definition η t = at+1 Bt+1 , multiplying both sides of the above inequality by B t+1 yields

B t+1 D(x, v t+1 ) + B t+1 (1 -β t )D(v t+1 , v t ) + β t D(v t+1 , y t ) -B t+1 (1 -β t )D(x, v t ) -B t+1 β t D(x, y t ) ≤ a t+1 ∇f (y t ) (x -v t+1 ) + a t+1 ψ(x) -ψ(v t+1 ) .
Using the scaling property ( 19) and the relationships α t = at+1 At+1 and a 2 t+1 L f /φ G t = A t+1 B t+1 , we obtain

B t+1 (1 -β t )D(v t+1 , v t ) + β t D(v t+1 , y t ) ≥ B t+1 α 2 t G t D(x t+1 , y t ) = A 2 t+1 B t+1 a 2 t+1 G t D(x t+1 , y t ) = A t+1 L f /φ D(x t+1 , y t ).
Combining the last two inequalities and using the facts B t+1 (1 -β t ) = B t and B t+1 β t = a t+1 σ f /φ , we arrive at

B t+1 D(x, v t+1 ) + A t+1 L f /φ D(x t+1 , y t ) -B t D(x, v t ) -a t+1 σ f /φ D(x, y t ) ≤ a t+1 ∇f (y t ) (x -v t+1 ) + a t+1 ψ(x) -ψ(v t+1 ) . ( 30
)
We then expand the gradient term on the right-hand side of (30) into two parts:

a t+1 ∇f (y t ) (x -v t+1 ) = a t+1 ∇f (y t ) (x -w t ) + a t+1 ∇f (y t ) (w t -v t+1 ), (31)
where w t = (1 -β t )v t + β t y t . For the first part,

a t+1 ∇f (y t ) (x -w t ) = a t+1 ∇f (y t ) (x -y t ) + a t+1 (1 -α t ) α t ∇f (y t ) (x t -y t ) ≤ a t+1 f (x) -f (y t ) -σ f /φ D(x, y t ) + a t+1 (1 -α t ) α t (f (x t ) -f (y t )) . (32)
Notice that

a t+1 1 -α t α t = a t+1 1 α t -1 = a t+1 A t+1 a t+1 -1 = A t+1 -a t+1 = A t .
Therefore, Equation (32) becomes

a t+1 ∇f (y t ) (x -w t ) ≤ a t+1 f (x) -A t+1 f (y t ) + A t f (x t ) -a t+1 σ f /φ D(x, y t ). ( 33
)
For the second part on the right-hand side of (31),

a t+1 ∇f (y t ) (w t -v t+1 ) = - a t+1 α t ∇f (y t ) (x t+1 -y t ) = -A t+1 ∇f (y t ) (x t+1 -y t ) ≤ -A t+1 f (x t+1 ) -f (y t ) -L f /φ D(x t+1 , y t ) , (34)
where in the last inequality we used the relative smoothness assumption in (29).

Summing the inequalities ( 30), ( 32) and ( 34), we have

B t+1 D(x, v t+1 ) -B t D(x, v t ) ≤ a t+1 f (x) -A t+1 f (x t+1 ) + A t f (x t ) + a t+1 (ψ(v t+1 ) -ψ(x)) ≤ -A t+1 f (x t+1 ) -f (x) + A t f (x t ) -f (x) + a t+1 ψ(x) -ψ(v t+1 ) ,
which is the same as

A t+1 f (x t+1 ) -f (x) + B t+1 D(x, v t+1 ) ≤ A t f (x t ) -f (x) + B t D(x, v t ) + a t+1 ψ(x) -ψ(v t+1 ) . ( 35
)
Finally we consider the term a t+1 ψ(x) -ψ(v t+1 ) . Using x t+1 = (1 -α t )x t + α t v t+1 and convexity of ψ, we have

ψ(x t+1 ) ≤ (1 -α t )ψ(x t ) + α t ψ(v t+1 ).
Since by definition α t = at+1 At+1 and (1 -α t ) = At At+1 , the above inequality is equivalent to

A t+1 ψ(x t+1 ) ≤ A t ψ(x t ) + a t+1 ψ(v t+1 ),
which implies (using

A t+1 = A t + a t+1 ) that for any x ∈ dom φ, A t+1 ψ(x t+1 ) -ψ(x) ≤ A t ψ(x t ) -ψ(x) + a t+1 ψ(v t+1 ) -ψ(x) . ( 36
)
Summing the inequalities ( 35) and ( 36) and using Φ = f + ψ, we have

A t+1 Φ(x t+1 ) -Φ(x) + B t+1 D(x, v t+1 ) ≤ A t Φ(x t ) -Φ(x) + B t D(x, v t ).
This can then be unrolled, and we obtain the desired result by setting x = x * .

Finally, the estimates of A t and B t follow from a direct adaptation of the techniques in [START_REF] Nesterov | Efficiency of the accelerated coordinate descent method on structured optimization problems[END_REF]. The only difference is the use of time-varying γ t = σ F/φ /(L F/φ G t ) instead of a constant γ = σ F/φ /L F/φ , which does not impact the derivations.

A.3. Proof of Lemma 2

The analysis in Lemma 1 is very pessimistic, since we use uniform lower and upper bounds for the Hessian of φ, whereas what we actually want is to bound is the differences between Hessians. If the Hessian is well-behaved (typically Lipschitz, or if φ is self-concordant), we can prove Lemma 2, which leads to a finer asymptotic convergence rate.

We start with the local quadratic representation of Bregman divergence:

D φ (x t+1 , y t ) = 1 2 x t+1 -y t 2 H(xt+1,yt) = α 2 t 2 v t+1 -w t 2 H(xt+1,yt) ≤ α 2 t 2 (1 -β t ) v t+1 -v t 2 H(xt+1,yt) + β t v t+1 -y t 2 H(xt+1,yt) ≤ α 2 t 2 (1 -β t ) v t+1 -v t 2 H(vt+1,vt) + β t v t+1 -y t 2 H(vt+1,yt) + α 2 t 2 (1 -β t ) H(x t+1 , y t ) -H(v t+1 , v t ) • v t+1 -v t 2 + α 2 t 2 β t H(x t+1 , y t ) -H(v t+1 , y t ) • v t+1 -y t 2 .
Now we use the Lipschitz property of ∇ 2 φ to bound the spectral norms of differences of Hessians:

H(x t+1 , y t ) -H(v t+1 , v t ) ≤ M z xy -z vv , H(x t+1 , y t ) -H(v t+1 , y t ) ≤ M z xy -z vy ,
where

z vv ∈ [v t+1 , v t ], z xy ∈ [y t , x t+1 ] and z vy ∈ [y t , v t+1 ].
Using the triangle inequality of norms, we have

z xy -z vy = z xy -y t + y t -z vy ≤ z xy -y t + y t -z vy ≤ x t+1 -y t + y t -v t+1 , and
z vv -z xy ≤ z vv -v t+1 + v t+1 -y t + y t -z xy ≤ v t -v t+1 + v t+1 -y t + y t -x t+1 .
Therefore, we have

d t max z xy -z vv , z xy -z vy ≤ v t -v t+1 + v t+1 -y t + y t -x t+1 ,
and consequently,

D φ (x t+1 , y t ) ≤ α 2 t 2 (1 -β t ) v t+1 -v t 2 H(vt+1,vt) + β t v t+1 -y t 2 H(vt+1,yt) + M d t α 2 t 2 (1 -β t ) v t+1 -v t 2 + β t v t+1 -y t 2 ≤ α 2 t 2 (1 -β t ) v t+1 -v t 2 H(vt+1,vt) + β t v t+1 -y t 2 H(vt+1,yt) + M d t α 2 t 2σ φ (1 -β t ) v t+1 -v t 2 H(vt+1,vt) + β t v t+1 -y t 2 H(vt+1,yt) = α 2 t 2 1 + M d t σ φ (1 -β t ) v t+1 -v t 2 H(vt+1,vt) + β t v t+1 -y t 2 H(vt+1,yt) = α 2 t 1 + M d t σ φ (1 -β t )D(v t+1 , v t ) + β t D(v t+1 , y t ) .
Combining with Lemma 1, we see that G t = min{κ σ , 1 + (M/σ φ )d t } satisfies the inequality ( 19). This finishes the proof of Lemma 2.

Note that this condition is not directly useful. Indeed, x t+1 and v t+1 depend on G t . Yet, under the uniform choice of G t ≤ κ φ , it can be shown that d t → 0 at rate (1 -1/ √ κ φ κ F/φ ) t because the sequences v t , x t and y t all converge to x * at this rate in the strongly convex case (Lin and Xiao, 2015, Theorem 1). As a consequence, Algorithm 1 will eventually use G t ≤ 2, leading to an asymptotic rate of (1 -1/ √ κ F/φ ) t .

B. Concentration of Hessians

In practice, preconditioned gradient methods such as DANE are often used with a step-size of 1. This implies the assumption of L F/φ = 1, which holds if n is sufficiently large with a given µ or if µ is sufficiently large for a given n (but µ ≤ L F always). Otherwise convergence is not guaranteed (which is why it is sometimes considered as "rather unstable"). If µ is such that H f (x) -H F (x) ≤ µ for all x ∈ B(0, D) then L F/φ = 1 can safely be chosen since

H F (x) -H f (x) µI d .
Note that this choice of µ is completely independent of λ. In this case, we use that H F (x) -H f (x) -µI to write that

H f (x) + µ H F (x) + 2µ (1 + 2µH -1 F (x))H F (x) 1 + 2µ λ H F (x).
These derivations are similar to the ones of Zhang and Xiao (2018, Lemma 3), and so we obtain σ F/φ = 1 + 2µ λ -1 and the corresponding relative condition number κ F/φ = 1 + 2µ λ , as explained in Section 1.1. We see that µ is independent of λ, but the problem is still very ill-conditioned for small values of λ, meaning that acceleration makes a lot of sense. In the quadratic case, tighter relative bounds can be derived.

B.1. The quadratic case

This section is focusing on proving Theorem 2.

Proof of Theorem 2. We consider the random variable a, and (a i ) i∈{1,...,n} are n i.i.d. variables with the same law as a. We introduce matrices Ĥ and H such that H f = Ĥ + λI d and H F = H + λI d . In particular, H = E aa = E Ĥ. We define for α ≥ 0, β > 0, H α,β = αH + βI d , and

S i = 1 n H -1 2 α,β (a i a i -H)H -1 2 α,β ,
which is such that E [S i ] = 0. This allows to have bounds of the form i S i ≤ t with probability 1 -δ and a spectral bound µ that depends on α, β, δ (and other quantities related to H and a i a i ). We note that n i=1

S i = H -1 2 α,β ( Ĥ -H)H -1 2 α,β ,
and write the concentration bounds on the S i as -tH α,β Ĥ -H tH α,β for some t > 0, which can be rearranged as:

Ĥ + tβI d (1 -tα)H Ĥ -tβI d (1 + tα)H.
Using H f = Ĥ + λI d and H F = H + λI d , the first equation can be rearranged as:

H F 1 1 -tα (H f + t(β -αλ)I d ) . ( 37
)
The second equation can be written

H f (1 + tα)I d + t(β -αλ)H -1 F H F ,
which, by adding t(β -αλ)I d on both sides, leads to

H f + t(β -αλ)I d (1 + tα)I d + 2t(β -αλ)H -1 F H F .
We let µ = t(β -αλ) and use H -1 F λ -1 I d to write that:

1 + αt + 2µ λ -1 (H f + µI d ) H F 1 1 -αt (H f + µI d ) . (38)
We then use the fact that a i a i and H are positive semidefinite and upper bounded by R 2 I to write that:

S i ≤ 1 n H -1 α,β max aa , H ≤ R 2 βn . (39)
Using the fact that H = E aa , we bound the variance as:

i E S i S i = 1 n E H -1 2 α,β (aa -H)H -1 α,β (aa -H)H -1 2 α,β = 1 n H -1 2 α,β (E aa H -1 α,β aa -HH -1 α,β H)H -1 2 α,β ≤ 1 n max R2 H -1 2 α,β E aa H -1 2 α,β , H -1 2 α,β HH -1 α,β HH -1 2 α,β ≤ 1 n H -1 2 α,β HH -1 2 α,β max R2 , H -1 2 α,β HH -1 2 α,β , with R2 ≥ a H -1 α,β a almost surely. We first notice that a i H -1 α,β a i ≤ R 2 β .
Then, we use the positive definiteness of H α,β and H and the fact that βH -1 α,β I d to show that for α > 0:

H -1 2 α,β HH -1 2 α,β = α -1 H -1 2 α,β (αH + β -β)H -1 2 α,β = α -1 I d -βH -1 α,β ≤ α -1 1 - β αL + β = L αL + β ,
where L is the spectral norm of H, i.e., L = H . A quick calculation shows that this formula is also true for α = 0. In the case α = 0 and β = 1 (absolute bounds), H α,β = I d and we recover that we can bound the variance by LR 2 n , leading to the usual additive bounds.

For α > 0, we use the simpler bound H

-1 2 α,β HH -1 2 α,β ≤ α -1 and R2 ≤ β -1 R 2 , leading to i E S i S i ≤ max(β -1 R 2 , α -1 ) nα .
For any 1 > δ > 0, we note c δ = 28 3 log 2d δ . We now set α = βn c δ R 2 , and assume that n > c δ (otherwise concentration bounds will be very loose anyway). In this case, β -1 R 2 ≥ α -1 , meaning that the bound on the variance becomes:

i E S i S i ≤ 1 α 2 c δ .
Similarly, according to (39), every S i is almost surely bounded as: S i ≤ 1 αc δ . We can now use Matrix Bernstein Inequality (Tropp, 2015, Theorem (6.1.1)) to get that with probability 1 -p δ and for t ≥ 0,

n i=1 S i ≤ t, with p δ = 2d • exp(- t 2 /2 (α 2 c δ ) -1 + (αc δ ) -1 t/3
).

We choose t = (2α) -1 , which leads to p δ = δ. By substituting the expressions of αt = 1 2 and βt = R 2 c δ n αt into Equation (38), we obtain:

3 2 + 2µ λ -1 ( Ĥλ + µI d ) H λ 2 Ĥλ + µI d , with µ = t(β -αλ) = 1 2 28R 2 3n log 2d δ -λ .
In case β is very small so that µ < 0 then it is always possible to choose δ < δ so that µ > 0. This means that the same bound on µ holds with probability 1 -δ > 1 -δ.

B.2. Almost surely bounded a

We first introduce Theorem 6, which proves a general concentration result that implies Theorem 3 as a special case. Theorem 6. We consider functions ϕ 1 , ϕ 2 , which are respectively L 1 and L 2 Lipschitz-continuous. We consider two sets X and Y which are contained in balls of center 0 and radius D 1 and D 2 . We assume that |ϕ 1 (a i x)| ≤ B 1 and |ϕ 2 (a i y)| ≤ B 2 almost surely for all x ∈ X and y ∈ Y. We consider

Y = sup x∈X , y∈Y 1 n n i=1 ϕ 1 (a i x)ϕ 2 (a i y) -Eϕ 1 (a x)ϕ 2 (a y) .
Then, for all 1 ≥ δ > 0, with probability greater than 1 -δ:

Y ≤ √ 4π (E a 2 ) 1 2 √ n (B 2 L 1 D 1 + B 1 L 2 D 2 ) + 2B 1 B 2 √ 2n log 1 δ .
Theorem 3 is then a direct corollary of Theorem 6, as shown below:

Proof of Theorem 3. The result is obtained by applying Theorem 6 with ϕ 1 = and ϕ 2 = 1 2 (•) 2 . This implies that with probability at least 1 -δ, sup x∈B(0,D), y∈B(0,1)

y 1 n n i=1 (a i x)a i a i -E (a x)aa y ≤ µ,
where the value of µ can be obtained by letting

B 1 = B , L 1 = M , D 1 = D, D 2 = 1, B 2 = sup y: y ≤1 y a i a i y ≤ R 2 and L 2 = sup y: y ≤1 2 y a i = 2R.
Proof of Theorem 6. If changing any a i to some a i , then the deviation in Y is at most (almost surely):

1 n sup x∈X , y∈Y ϕ 1 (a i x)ϕ 2 (a i y) + sup x∈X , y∈Y ϕ 1 (a i x)ϕ 2 (a i y) ≤ 2 n B 1 B 2 .
Mac-Diarmid's inequality (see, e.g., Vershynin, 2019, Theorem 2.9.1) thus implies that with probability greater than 1 -δ,

Y ≤ EY + 2B 1 B 2 √ 2n log 1 δ . (40)
In order to bound EY , we first use classical symmetrization property (see, e.g., Vershynin, 2019, Section 6.4)

EY ≤ √ 2π • E sup x∈X , y∈Y 1 n n i=1 ε i ϕ 1 (a i x)ϕ 2 (a i y),
where each ε i is an independent standard normal variable.

Denoting Z x,y = 1 n n i=1 ε i ϕ 1 (a i x)ϕ 2 (a i y),
we have, for any x, y, x , y , assuming the a i are fixed,

E(Z x,y -Z x ,y ) 2 = 1 n 2 n i=1 ϕ 1 (a i x)ϕ 2 (a i y) -ϕ 1 (a i x )ϕ 2 (a i y ) 2 = 1 n 2 n i=1 ϕ 1 (a i x) ϕ 2 (a i y) -ϕ 2 (a i y ) + ϕ 1 (a i x) -ϕ 1 (a i x ) ϕ 2 (a i y ) 2 ≤ 1 n 2 n i=1 2 ϕ 1 (a i x) 2 ϕ 2 (a i y) -ϕ 2 (a i y ) 2 + 2 ϕ 2 (a i y ) 2 ϕ 1 (a i x) -ϕ 1 (a i x ) 2 ≤ 1 n 2 n i=1 2B 2 1 ϕ 2 (a i y) -ϕ 2 (a i y ) 2 + 2B 2 2 ϕ 1 (a i x) -ϕ 1 (a i x ) 2 .
We then have, using Lipschitz-continuity:

E(Z x,y -Z x ,y ) 2 ≤ 1 n 2 n i=1 2B 2 1 L 2 2 a i y -a i y 2 + 2B 2 2 L 2 1 a i x -a i x 2 = E( Zx,y -Zx ,y ) 2 , for Zx,y = 1 n n i=1 √ 2B 2 L 1 ε1i a i x + √ 2B 1 L 2 ε2i a i y ,
with all ε1i and ε2i independent standard random variables.

Using Sudakov-Fernique inequality (Vershynin, 2019, Theorem 7.2.11), we get

EY = √ 2πE sup x∈X , y∈Y Z x,y ≤ √ 2πE sup x∈X , y∈Y Zx,y = √ 4πB 2 L 1 E sup x∈X 1 n n i=1 ε1i a i x + √ 4πB 1 L 2 E sup y∈Y 1 n n i=1 ε2i a i y ≤ √ 4πB 2 L 1 D 1 E 1 n n i=1 ε1i a i + √ 4πB 1 L 2 D 2 E 1 n n i=1 ε2i a i ≤ √ 4πB 2 L 1 D 1 E 1 n n i=1 ε1i a i 2 + √ 4πB 1 L 2 D 2 E 1 n n i=1 ε2i a i 2 ≤ √ 4πB 2 L 1 D 1 (E a 2 ) 1 2 √ n + √ 4πB 1 L 2 D 2 (E a 2 ) 1 2 √ n .
Plugging this into Equation ( 40), we obtain that with probability greater than 1 -δ,

Y ≤ √ 4π (E a 2 ) 1 2 √ n (B 2 L 1 D 1 + B 1 L 2 D 2 ) + 2B 1 B 2 √ 2n log 1 δ .
Remark 4 (Relative bounds). In the quadratic case, considering relative bounds allowed to choose smaller values of µ and to tighten the bounds on the relative condition number by a √ n factor. Theorem 3 consists in bounding (using the definition of the operator norm) sup x∈B(0,D),y∈B(0,1)

1 n n i=1 (a i x)(a i y) 2 -y H(x)y ,
and heavily relies on the fact that (a i y) 2 is independent of x. The proof needs to be adapted in the case of the relative bounds since this term becomes (a i H -1 2 α,β (x)y) 2 , which now depends on x as well, and thus requires a different control.

B.3. Subgaussian a

We considered in the previous section a splitting of the summands of the Hessians as a product of 2 functions. We now present a different bound that is designed for a product of an arbitrary number of functions ϕ 1 , . . . , ϕ r : R → R. This section is devoted to proving Theorem 7, which is based on the chaining argument (Boucheron et al., 2013, Chapter 13), and from which Theorem 4 can be derived directly. Theorem 7. Assume that for all i, ϕ i (0) = 0 and ϕ i is 1-Lipschitz. Assume that a is ρ-subgaussian, and that for all k,

sup x∈B(0,1) |ϕ k (a i x)| ≤ B k . Denote B =
1 n i∈[n] r k=1 ϕ k (a i x k ) -E r k=1 ϕ k (a x k ) ≥ ρ r C r (d + γ) 1 √ dn + (ρ -r B) 1-2/r n   ≤ r π 2 6 e -γ .
We are primarily interested in the case r = 3, ϕ 1 = ϕ 2 = id (the identity mapping) to control distances between Hessians.

Proof. We look for bounds on

Y := sup x1,...,xr∈S1 1 n i∈[n] r k=1 ϕ k (a i x k ) -E a r k=1 ϕ k (a x k ) . (41)
For all j ≥ 0, let N j be an -net of S 1 that approximates S 1 to distance 2 -j . Then, N j can be chosen as |N j | ≤ (1 + 2 j+1 ) d (see, e.g., Vershynin, 2019, Section 4.2). For all x ∈ S 1 , let Π j (x) be some point in N j such that x -Π j (x) ≤ 2 -j . By convention we take Π 0 (x) = 0.

Then for all (x 1 , . . . , x r ) ∈ S r , using the chaining approach [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF], we write

1 n i∈[n] k∈[r] ϕ k (a i x k ) = j≥0 1 n i∈[n]    k∈[r] ϕ k (a i Π j+1 (x k )) - k∈[r] ϕ k (a i Π j (x k ))    = j≥0 k∈[r] 1 n i∈[n] k-1 =1 ϕ (a i Π j+1 (x )) ϕ k (a i Π j+1 (x k )) -ϕ k (a i Π j (x k )) r =k+1 ϕ (a i Π j (x )).
Let j ≥ 0 and k ∈ [r] be fixed. Consider a term of the form Z = 1 n i∈[n] Z i , with

Z i = k-1 =1 ϕ (a i u ) ϕ k (a i u k ) -ϕ k (a i v k ) r =k+1 ϕ (a i v ), (42)
where u ∈ N j , v ∈ N j+1 , and u k -v k ≤ j := 2 -j+1 . By the triangle inequality, for all x ∈ S 1 , letting u = Π j (x ) and v = Π j+1 (x ), these assumptions are satisfied. For each Z i and t > 0, we have:

P(Z i ≥ j ρ r t) ≤ P |ϕ (a i u )| ≥ ρt 1/r for some < r, or |ϕ k (a i u k ) -ϕ k (a i v k )| ≥ ρ j t 1/r ,
or |ϕ (a i v )| ≥ ρt 1/r for some > k .

Therefore, we have

P(Z i ≥ j ρ r t) ≤ 2re -t 2/r /2 if t ≤ ρ -r P j,k and
P(Z i ≥ j ρ r t) = 0 if t > ρ -r P j,k ,
where we noted P j,k := min 2B/ j , 2(B/B k )R . We will also make use of notation j * (k) := log 2 (R/B k ) , so that

j ≤ j * (k) ⇒ P j,k = 2B/ j , j > j * (k) ⇒ P j,k = 2(B/B k )R. Fixing j ≥ 0, k ∈ [r],
we write for any θ > 0 (a specific θ will be chosen later):

E e (θ/n)ρ -r [Zi-EZi]/ j = 1 + θ n 2 E ( -1 j ρ -r (Z i -EZ i )) 2 F ((θ/n)( -1 j ρ -r (Z i -EZ i ))) ,
where

F (x) := x -2 [e x -x -1] ≤ e |x| .
Thus using this bound and the inequality xy ≤ x 2 + y 2 :

E e (θ/n)ρ -r [Zi-EZi]/ j ≤ 1 + θ n 2 E(( -1 j ρ -r (Z i -EZ i )) 4 + Ee 2(θ/n)ρ -r |Zi-EZi|/ j . (43)
By the sub-gaussian tail assumption, E( -1 j ρ -r (Z i -EZ i )) 4 is bounded by a constant κ r dependent on r. We now assume that θ is such that θ n ≤ min (ρ -r P j,k ) 2/r-1 8 , 1 , which is equivalent to having (θ/n)y ≤ y 2/r /8 for y ∈ [0, ρ -r P j,k ] and r ≥ 2. Then, Ee 2(θ/n)ρ -r |Zi-EZi|/ j is also bounded by another constant κ r dependent on r. Indeed, by the sub-gaussian tail assumption, |EZ i | ≤ ρ r j s r for some r-dependent constant, and we can then use the fact that:

E e 2(θ/n)ρ -r |Zi-EZi|/ j ≤ E e 2(θ/n)ρ -r (|Zi|+|EZi|)/ j ≤ e 2(θ/n)sr Ee 2θ/nρ -r |Zi|/ j ≤ e 2(θ/n)sr [1 + 2θ n ∞ 0 e 2(θ/n)y [P(Z i ≥ yρ r j ) + P(-Z i ≥ yρ r j )]dy] ≤ e 2(θ/n)sr [1 + 2θ n 2r ρ -r P j,k 0 e 2(θ/n)y-y 2/r /2 dy], ≤ e 2(θ/n)sr [1 + 2θ n 2r ∞ 0 e -y 2/r /4 dy] = e 2(θ/n)sr [1 + θ n c r ].
We finally use the fact that θ/n ≤ 1 to write Ee 2(θ/n)ρ -r |Zi-EZi|/ j ≤ κ r , with κ r = e 2sr [1 + c r ]. We write κ r = κ r + κ r and use Equation ( 43) together with the independence of the Z i to obtain:

E e θρ -r [ 1 n n i=1 Zi-EZi]/ j ≤ 1 + θ n 2 κ r n ≤ e θ 2
n κ r .

Thus, using that P(X ≥ y) = P(e X ≥ e y ) ≤ e -y Ee X (Markov Inequality), we have that for fixed u u , v , ∈ [r] in the suitable -nets is upper bounded for all θ ∈ [0, min(n, n(ρ -r P j,k ) 2/r-1 /8)] as:

P   1 n i∈[n] Z i -EZ i ≥ ρ r j t j,k   ≤ exp (r + 1)d ln(1 + 2 j+2 ) -θt j,k + κ r θ 2 /n . (44)
We see in Equation ( 42) that the variables Z i are built by fixing a specific either u for < k, v for > k, and u k and v k , meaning that there are actually r + 1 variables to be fixed in nets of resolution either 2 -j or 2 -j-1 . Note that all Z i for i ∈ {1, • • • , n} are constructed with the same choice of u and v . Therefore, the number of possible choices for u ∈ N j and v ∈ N j+1 involved in the definition of Z i is upper-bounded by

|N j+1 | r+1 ≤ e d(r+1) ln(1+2 j+2 ) .
Combining this with Equation (44), we obtain using a union bound that:

P sup u ,v Z -EZ ≥ ρ r j t j,k = P ∪ u ,v Z -EZ ≥ ρ r j t j,k ≤ u ,v P (Z -EZ ≥ ρ r j t j,k ) ≤ exp (r + 1)d ln(1 + 2 j+2 ) -θt j,k + κ r θ 2 /n .
Let now θ j,k = min(n, n(ρ -r P j,k ) 2/r-1 /8, √ nd), and

t j,k = κ r θ j,k n + 1 θ j,k [d(r + 1) ln(1 + 2 j+2 ) + γ + 2 ln(j + 1)],
where γ > 0 is a free parameter. We then use the chaining decomposition of Y = sup x j≥0,k∈[r] Z -EZ , and another union bound on j and k to write that:

P   Y ≥ ρ r j≥0,k∈[r] j t j,k   = P   sup x j≥0,k∈[r] Z -EZ ≥ ρ r j≥0,k∈[r] j t j,k   ≤ P   j≥0,k∈[r] sup x {Z -EZ} ≥ ρ r j≥0,k∈[r] j t j,k   ≤ j≥0,k∈[r] P sup u ,v {Z -EZ} ≥ ρ r j t j,k ≤ j≥0,k∈[r]
e -γ-2 ln(1+j) .

In the end, using that j≥1 j -2 = π 2 /6, we obtain:

P   Y ≥ ρ r j≥0,k∈[r] j t j,k   ≤ r π 2 6 e -γ .
Moreover, one has

j t j,k ≤ j A r (1 + j)(d + γ) (ρ -r P j,k ) 1-2/r n + 1 √ nd + 1 n ,
for some suitable constant A r dependent only on r. Fix some k ∈ [r]. Write:

1 A r j≥0 j t j,k ≤ 4 d + γ √ nd + d + γ n j * (k) j=0 j (2ρ -r B/ j ) 1-2/r (1 + j) + d + γ n j>j * (k) j (2ρ -r BR/B k ) 1-2/r (1 + j) ≤ 4 d + γ √ nd + d + γ n A r (ρ -r B) 1-2/r 1 + (B k /R) 2/r ln(R/B k ) ,
where A r is another constant depending only on r.

Since B k ≤ R, then (B k /R) 2/r ln(R/B k ) is bounded by a (r-dependent) constant.
We know present Corollary 1, which is a consequence of Theorem 7. We consider again i.i.d. a i , bounded by R, satisfying the subgaussian tail assumption with parameter ρ, and some function ϕ that is 1-Lipschitz, and uniformly bounded by B ϕ . Writing

H(x) = 1 n n i=1 a i a i ϕ(a i x), (45)
We have the following corollary.

Corollary 1. Thus for 1 > δ > 0, with probability at least 1 -δ, it holds for some C > 0 that

sup x∈S1 H(x) op ≤ C ρ 3 (d + ln(1/δ) + ln(5π 2 /6))[(1 + ρ -1 B ϕ )/ √ dn + (1 + {ρ -3 R 2 B ϕ } 1-2/3 )/n]. (46)
Proof. Let us write ϕ 3 (u) = ϕ(u) -ϕ(0). Then ϕ 3 satisfies our assumptions (1-Lipschitz, ϕ 3 (0) = 0). Moreover, we can decompose matrix H(x) -EH(x) into M (x) + N , where

M (x) = 1 n n i=1 a i a i ϕ 3 (a i x) -Ea 1 a 1 ϕ 3 (a 1 x) , N = 1 n n i=1 a i a i ϕ(0) -Ea 1 a 1 ϕ(0) .
Taking r = 2, ϕ 1 = ϕ 2 = Id, the Theorem 7 gives us that

P( N op ≥ C 2 |ϕ(0)|ρ 2 (γ + d)(1/ √ dn + 1/n)) ≤ 2 π 2 6 e -γ . (47)
Taking next r = 3, and B = R 2 B ϕ , we obtain

P( sup x∈S1 M (x) op ≥ C 3 ρ 3 (d + γ)(1/ √ dn + [ρ -3 R 2 B ϕ ] 1-2/3 /n)) ≤ 3 π 2 6 e -γ . (48)
these two bounds give us that for all γ > 0, we have, setting C = C 2 + C 3 :

P( sup x∈S1 H(x) op ≥ C ρ 3 (d + γ)[(1 + ρ -1 B ϕ )/ √ dn + (1 + {ρ -3 R 2 B ϕ } 1-2/3 )/n]) ≤ 5 π 2 6 e -γ . (49)
We finally take γ = -ln 6δ 5π 2 .

The last step required to prove Theorem 4 is to consider the supremum over B(0, D) with an arbitrary M -Lipschitz function, which can be done by direct reduction:

Proof of Theorem 4. To apply this to , defined on B(0, D) and M Lipschitz, we apply Corollary 1 to ϕ(x) = 1 M D (Dx) (which is 1-Lipschitz on B(0, 1)). Then, B ϕ = B /M D and the right hand side must be multiplied by M D.

Remark 5. Note that there is a difference in the way Theorem 6 and Theorem 7 are applied to our linear models problem. In particular, Theorem 6 considers ϕ 1 = • 2 and ϕ 2 = , whereas Theorem 7 uses ϕ 1 = ϕ 2 = Id and ϕ 3 = . Theorem 6 can be adapted to work with r = 3, but the bound does not improve when splitting • 2 into Id × Id. Similarly, Theorem 7 could be used with r = 2 and ϕ 1 = • 2 /(2R) (to respect the 1-Lipschitz assumption), but in this case the bound can only be worse since the main difference is that the ρ 3 factor becomes Rρ 2 , and ρ is generally smaller than R.

B.4. Tightness of Theorem 7

Consider that the a i uniformly distributed on the sphere with radius R = √ d, and take for f k the identity. Such vectors can be constructed by taking vectors A i with coordinates i.i.d. standard gaussian, and setting a i = √ d A i -1 A i . The subgaussianity parameter ρ can then be taken equal to 1.

Then, using known results about maximal correlation between variables with fixed marginals (e.g., Vershynin, 2019, Section 3), the expectation E r k=1 a 1 x k is maximized, over choices x k ∈ S 1 , by taking x 1 = • • • = x r . We may choose x 1 = e 1 , the first unit vector, by rotational invariance, and thus the expectation is upper-bounded as:

E r k=1 a 1 x k ≤ Ed r/2 E |A i (1)| r A i r .
This is of order 1, as can be shown using concentration inequalities on the deviations of A i from

√ d. Consider then the empirical sum 1 n i∈[n] k∈[r] a i x k . Choose x k = d -1/2 a 1 for all k ∈ [r].
Then this empirical sum evaluates to

1 n i∈[n] k∈[r] a i x k = 1 n d r/2 + 1 n n i=2 k∈[r]
a i a 1 .

The second sum can be shown to be of order 1 (conditioning on a 1 , and then using, e.g., Bienaymé-Tchebitchev inequality). Thus, one cannot hope to establish concentration without extra assumptions on the data distribution unless d r/2 = O(n).

Contrast this with the result of Theorem 7: for R

= √ d, B = R r and ρ = 1, it gives that Y ≤ O(d (r/2)(1-2/r) d/n) = O(d r/2 /n).
Thus the result is sharp for the particular example we just considered.

C. Experimental setting

Some implementation details are omitted in the main text due to lack of space. To ease the reader's understanding, we provide these details here, along with some additional experimental results. We also provide code for SPAG in supplementary material.

Optimization problem. We used the logistic loss with quadratic regularization, meaning that the function at node i is:

f i : x → 1 m m j=1 log 1 + exp(-y i,j x a (i) j ) + λ 2 x 2 ,
where y i,j ∈ {-1, 1} is the label associated with a (i) j , the j-th sample of node i. The local datasets are constructed by shuffling the LibSVM datasets, and then assigning a fixed portion to each worker. Then, the server subsamples n points from its local dataset to construct the preconditioning dataset. To assess the suboptimality, we let the best algorithm run for more time in order to get a good approximation of the minimum error. Then, we subtract it to the running error of an algorithm to get the suboptimality at each step.

Adjusting α t and β t . We found that choosing A 0 = 0 and B 0 = 1 for SPAG is usually not the best choice. Indeed, rates are asymptotic and sequences α t and β t converge very slowly when σ F/φ is small, whereas we typically rarely use more than about 100 iterations of SPAG. Therefore, we start the algorithm with A t0 and B t0 with t 0 > 0 instead. We used t 0 = 50, but SPAG is not very sensitive to this choice.

Tuning the momentum. Figure 5(a) evaluates the relevance of tuning the parameters controlling the momentum of SPAG and HB-DANE. To do so, we compare the default values of β = (1 -(1 + 2µ/λ) -1/2 ) 2 (for HB-DANE) and σ F/φ = 1/(1 + 2µ/λ) (for SPAG) to values obtained through a grid search on the KDD2010 dataset with λ = 10 -7 . We tune HB-DANE by using a grid-search of resolution 0.05 to test the values between 0.5 and 1. For n = 10 3 , theory predicts a momentum of β = 0.86 and the grid search gives β = 0.85. For n = 10 4 , theory predicts β = 0.81 and the grid search gives β = 0.8. For SPAG, we test σ F/φ = 10 -2 , 3 × 10 -3 , 10 -3 and so on until σ F/φ = 10 -5 (so roughly divided by 3 at each step). For n = 10 3 , theory predicts σ F/φ = 0.005 and the tuning yields σ F/φ = 0.006. For n = 10 4 , theory predicts σ F/φ = 0.0099 and the grid-search leads to σ F/φ = 0.01. We do not display the curves in this case (n = 10 4 ) since they are nearly identical. Therefore, the grid-search always obtains the value on the grid that is closest to the theoretical value of the parameter, and the difference in practice is rather small, as can be seen in Figure 5(a). This is why we use default values in the main text.

Local subproblems. Local problems are solved using a sparse implementation of SDCA [START_REF] Shalev-Shwartz | Sdca without duality, regularization, and individual convexity[END_REF]. In practice, the ill-conditionned regime is very hard, especially when µ is small. Indeed, the local subproblems are very hard to solve, and it should be beneficial to use accelerated algorithms to solve the inner problems. In our experiments, we warm-start the local problems (initializing on the solution of the previous one), and keep doing passes over the preconditioning dataset until ∇V t (x t ) ≤ 10 -9 (checked at each epoch). This threshold is important because it greatly affects the performances of preconditioned gradient methods. Figure 5(b) compares the performances of SPAG, DANE and HB-DANE for different number of passes on the inner problems for the RCV1 dataset for n = 10 4 and λ = 10 -5 . We use µ = 2 × 10 -5 and a step-size of 1 for all algorithms. We first see that increasing the number of passes significantly improves the convergence speed of all algorithms. Besides, heavy-ball acceleration does not seem very efficient when local problems are not solved accurately enough. On the contrary, SPAG seems to enjoy faster rates than DANE nevertheless. It would be interesting to understand these different behaviours more in details.

Gain far from the optimum. So far, we have presented experiments with good initializations (solution for the local dataset), and argued why G t was very small in this case. Because of Lemma 2, one would expect that G t could be large when x t is very far from x * . Yet, We see in the proof of Lemma 2 that the Lipschitz constant of the Hessian only needs to be considered for any convex set that contains x t+1 , v t+1 , y t and v t . In the case of logistic regression, the third derivative decreases very fast when far from 0, meaning that the local Lipschitz constant of the Hessian is small when the iterates are far from 0. In other words, the Hessian changes slowly when far from the optimum (at least for logistic regression).

We believe that this is the reason why G t can always be chosen of order 1 (smaller than 2) in our experiments, and that this holds regardless of the initialization. To support this claim, we plot in Figure 5(c) the values of the gain for the RCV1 dataset with λ = 10 -7 and 5 different x 0 sampled from N (0, 10 3 ), the normal law centered at 0 with variance 10 3 . We use a step-size of 0.9 and µ = 2 × 10 -5 . We first see that for G min = 1, the gain is always very low, and actually increases at some point instead of becoming lower and lower, so the fact that we were able to choose G t of order 1 in the other experiments is not linked to the good initialization. We had to choose a slightly higher µ than in the other experiments in order to satisfy the relative smoothness condition, which was not satisfied at each iteration otherwise. Since G t is small in practice and the smaller the G t the better the rate, we test SPAG with no minimum value for the gain G t . The curve for the gain in this case is shown by G min = 0, and we see that the true gain stabilizes to a higher value, since updates are more agressive. We discuss the efficiency of this version in the next paragraph. Note that the oscillations are not due to numerical instability or inaccurate solving of the inner problems, but rather to the fact that the step-size is slightly too big so sometimes the smoothness inequality is not verified. Yet, this does not affect the convergence of SPAG, as shown in Figure 5(d).

Line Search with no minimum value. Since the gain is almost always smaller than 1, the line-search in SPAG generally only consists in checking that G t = 1 works, which can be done locally. Therefore, there is no added communication cost. As discussed earlier, it is possible to allow G t < 1 when performing line search, which makes SPAG slightly more adaptative at the cost of a few more line-search loops. Figure 5(d) presents the difference between SPAG using a line search with G min = 0 and G min = 1. The curves show the suboptimality for the runs used to generate Figure 5(c). Note that we omit the cost of line search in the iteration cost (we still count in terms of number of iterations, even though more communication rounds are actually needed when G min = 0). We see that setting G min = 0 is initially slightly faster but that the rate is very similar, so that using G min = 0 may slightly improve iteration complexity but is not worth doing in this case. Note that suboptimality curves for different initializations are almost indistinguishable, which can be explained by the fact that the quadratic penalty term dominates and that all initializations have roughly the same norm (since d is high).

Figure 1 .

 1 Figure 1. Effect of µ on the convergence speed of SPAG on RCV1 with λ = 10 -7 and n = 10 4 .

Figure 2 .

 2 Figure2. Logistic regression on RCV1 dataset. The legend is the same for all figures. We use µ = 0.1/n, except for n = 10 5 for which we use µ = 2 • 10 -6 .

Figure 3 .

 3 Figure3. Ridge regression on RCV1 dataset with λ = 10 -5 . The legend is the same as in Figure4. We use µ = 4/n, and L F/φ = 1, except for SPAG with n = 10 5 for which we use L -1 F/φ = 0.9 as it proved to be more stable.

Figure 4 .

 4 Figure4. Logistic regression on KDD2010 with λ = 10 -7 . We use µ = 0.1/(2n), except for n = 10 3 where µ = 10 -5 , and L F/φ = 2 for SPAG.

  r k=1 B k . For suitable constant C r , for all γ > 0, one has that P   sup x1,...,xr∈B(0,1)

e

  αy p(X ≥ y)dy, with α = 2θ/n and X = ρ -r |Z i | j to get:

  Impact of Gmin on the suboptimality.

Figure 5 .

 5 Figure 5. Impact of several implementation details.

  approximates the overall Hessian by the harmonic mean of the local Hessians. It is equivalent to DANE in the quadratic case. They obtain a communication complexity that has logarithmic dependence on the condition number but requires local sample size n > d.[START_REF] Mahajan | An efficient distributed learning algorithm based on effective local functional approximations[END_REF] proposed a distributed algorithm based on local function approximation, which is related to the preconditioning idea of DANE.

	Yuan and Li (2019) revisited the analysis of DANE and
	found that the worse complexity of O((κ 2 /n) log(1/ )) is
	due to the lost statistical efficiency when averaging m differ-
	ent updates computed by (11). They propose to use a single
	local preconditioner at the server and obtain a communica-tion complexity of O((1 + κ / √ n) log(1/ )) for quadratic
	functions. In addition, they propose a variant of DANE
	with heavy-ball momentum (DANE-HB), and show that it has communication complexity O(( √ κ /n 1/4 ) log(1/ ))
	for quadratic loss functions, matching that of DiSCO and

AIDE. For non-quadratic functions, they show DANE-HB has accelerated local convergence rate near the solution.

[START_REF] Wang | GIANT: Globally improved approximate Newton method for distributed optimization[END_REF]

proposed GIANT, an approximate New-ton method that

Acknowledgements

This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). We also acknowledge support from the European Research Council (grant SEQUOIA 724063) and from the MSR-INRIA joint centre.

Appendix

A. Convergence Analysis of SPAG

This section provides proofs for Lemma 1, Theorem 1 and Lemma 2 presented in Section 3. Before getting to the proofs, we first comment on the nature of the accelerated convergence rate obtained in Theorem 1.