Communication-Efficient Accurate Statistical Estimation

TL;DR

This paper introduces CEASE, a set of communication-efficient iterative algorithms for distributed statistical estimation that achieve rapid convergence and statistical efficiency with minimal communication, suitable for large-scale data scenarios.

Contribution

The paper proposes novel CEASE algorithms that adapt to loss function similarities, converge quickly without good initialization, and achieve statistical efficiency in finite steps.

Findings

Algorithms converge rapidly with linear rates.

Statistical efficiency achieved in finite steps.

Validated through extensive experiments on synthetic and real data.

Abstract

When the data are stored in a distributed manner, direct application of traditional statistical inference procedures is often prohibitive due to communication cost and privacy concerns. This paper develops and investigates two Communication-Efficient Accurate Statistical Estimators (CEASE), implemented through iterative algorithms for distributed optimization. In each iteration, node machines carry out computation in parallel and communicate with the central processor, which then broadcasts aggregated information to node machines for new updates. The algorithms adapt to the similarity among loss functions on node machines, and converge rapidly when each node machine has large enough sample size. Moreover, they do not require good initialization and enjoy linear converge guarantees under general conditions. The contraction rate of optimization errors is presented explicitly, with dependence on the local sample size unveiled. In addition, the improved statistical accuracy per iteration is derived. By regarding the proposed method as a multi-step statistical estimator, we show that statistical efficiency can be achieved in finite steps in typical statistical applications. In addition, we give the conditions under which the one-step CEASE estimator is statistically efficient. Extensive numerical experiments on both synthetic and real data validate the theoretical results and demonstrate the superior performance of our algorithms.

Equations

$$\underset{\mathbf{\theta }\in {\mathbb{R}}^p}{\min }\left\{f(\mathbf{\theta })+g(\mathbf{\theta })\right\},$$

$$\underset{\mathbf{\theta }\in {\mathbb{R}}^p}{\min }\left\{f(\mathbf{\theta })+g(\mathbf{\theta })\right\},f(\mathbf{\theta })=\frac{1}{m}\sum _{k=1}^m{f}_k(\mathbf{\theta }).$$

$$f(\mathbf{\theta })=f^{(1)}(\mathbf{\theta })+R(\mathbf{\theta }),\text{mbox}whereR(\mathbf{\theta })=f(\mathbf{\theta })-f^{(1)}(\mathbf{\theta }).$$

$$R_k(\mathbf{\theta })=f_k(\mathbf{\theta })-[f_k(\stackrel{ˉ}{}\mathbf{\theta })+⟨\nabla f_k(\stackrel{ˉ}{}\mathbf{\theta }),\mathbf{\theta }-\stackrel{ˉ}{}\mathbf{\theta }⟩],$$

$$f_k(\mathbf{\theta })-⟨\nabla f_k(\stackrel{ˉ}{}\mathbf{\theta })-\nabla f(\stackrel{ˉ}{}\mathbf{\theta }),\mathbf{\theta }⟩$$

$${\mathbf{\theta }}_{t+1}={\mathrm{argmin}}_{\mathbf{\theta }}\left\{f_1(\mathbf{\theta })+g(\mathbf{\theta })-⟨\nabla f_1({\mathbf{\theta }}_t)-\nabla f({\mathbf{\theta }}_t),\mathbf{\theta }⟩\right\}$$

$${\mathbf{\theta }}_{t,k}={\mathrm{argmin}}_{\mathbf{\theta }}\left\{f_k(\mathbf{\theta })+g(\mathbf{\theta })-⟨\nabla f_k({\mathbf{\theta }}_t)-\nabla f({\mathbf{\theta }}_t),\mathbf{\theta }⟩\right\}$$

$$f_k(\mathbf{\theta })=\frac{1}{2n}\sum _{i\in {\mathcal{I}}_k}(y_i-{\mathbf{x}}_i^{\top }\mathbf{\theta }{)}^2=\frac{1}{2}{\mathbf{\theta }}^{\top }{\hat{\mathbf{\Sigma }}}_k\mathbf{\theta }-{\hat{\mathbf{w}}}_k^{\top }\mathbf{\theta }+\frac{1}{2n}\sum _{i\in {\mathcal{I}}_k}{y}_i^2.$$

$${\mathbf{\theta }}_{t+1,k}$$

$${\mathbf{\theta }}_{t+1}$$

$${\mathbf{x}}_{t+1}={\mathrm{prox}}_{{\alpha }^{-1}h}({\mathbf{x}}_t)={\mathrm{argmin}}_{\mathbf{x}\in {\mathbb{R}}^p}\{h(\mathbf{x})+(\alpha /2)\Vert \mathbf{x}-{\mathbf{x}}_t{\Vert }_2^2\},\forall t\ge 0,$$

$${\mathbf{\theta }}_{t+1}={\mathrm{prox}}_{{\alpha }^{-1}(f+g)}({\mathbf{\theta }}_t)={\mathrm{argmin}}_{\mathbf{\theta }\in {\mathbb{R}}^p}\left\{f(\mathbf{\theta })+g(\mathbf{\theta })+\frac{\alpha }{2}\Vert \mathbf{\theta }-{\mathbf{\theta }}_t{\Vert }_2^2\right\}.$$

$$\stackrel{~}{}{\mathbf{\theta }}_{s,k}$$

$${\tilde{\mathbf{\theta }}}_{s+1}$$

$${\mathbf{\theta }}_{t+1}={\mathrm{argmin}}_{\mathbf{\theta }}\left\{f_1(\mathbf{\theta })+g(\mathbf{\theta })-⟨\nabla f_k({\mathbf{\theta }}_t)-\nabla f({\mathbf{\theta }}_t),\mathbf{\theta }⟩+\frac{\alpha }{2}\Vert \mathbf{\theta }-{\mathbf{\theta }}_t{\Vert }_2^2\right\},$$

$${\mathbf{\theta }}_{t,k}={\mathrm{argmin}}_{\mathbf{\theta }}\left\{f_k(\mathbf{\theta })+g(\mathbf{\theta })-⟨\nabla f_k({\mathbf{\theta }}_t)-\nabla f({\mathbf{\theta }}_t),\mathbf{\theta }⟩+\frac{\alpha }{2}\Vert \mathbf{\theta }-{\mathbf{\theta }}_t{\Vert }_2^2\right\}$$

$${\rho }_0=\sup \left\{c\in [0,\rho ]:\{f_k+g{\}}_{k=1}^m\text{ are }c\text{-strongly convex in }B(\hat{}\mathbf{\theta },R)\right\}.$$

$$\Vert {\mathbf{\theta }}_{t+1}-\hat{}\mathbf{\theta }{\Vert }_2\le \Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2\cdot \frac{\frac{\delta }{{\rho }_0+\alpha }\sqrt{{\rho }^2+2\alpha \rho }+\alpha }{\rho +\alpha },0\le t\le T-1;$$

$$\Vert {\mathbf{\theta }}_{t+1}-\hat{}\mathbf{\theta }{\Vert }_2\le \Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2\cdot \frac{{\gamma }_t\sqrt{{\rho }^2+2\alpha \rho }+\alpha }{\rho +\alpha },0\le t\le T-1,$$

$$\Vert {\mathbf{\theta }}_{t+1}-\hat{}\mathbf{\theta }{\Vert }_2\le C_1(\delta /\rho )\Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2,0\le t\le T-1;$$

$$\Vert {\mathbf{\theta }}_{t+1}-\hat{}\mathbf{\theta }{\Vert }_2\le C_2(\delta /\rho {)}^2\Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2.$$

$$\Vert {\mathbf{\theta }}_{t+1}-\hat{}\mathbf{\theta }{\Vert }_2\le (\delta /{\rho }_0)\Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2,\forall t\ge 0.$$

$$\Vert {\mathbf{\theta }}_{t+1}-\hat{}\mathbf{\theta }{\Vert }_2\le \frac{\delta }{{\rho }_0}\Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2\cdot \min \left\{1,\frac{\delta }{\rho }\left(1+\frac{M}{{\rho }_0}\Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2\right)\right\},\forall t\ge 0.$$

$$h(y_i;{\mathbf{x}}_i,{\mathbf{\theta }}^{\ast })=c({\mathbf{x}}_i,y_i)\exp \left({\mathbf{y}}_i({\mathbf{x}}_i^{\top }{\mathbf{\theta }}^{\ast })-b({\mathbf{x}}_i^{\top }{\mathbf{\theta }}^{\ast })\right).$$

$$f_k(\mathbf{\theta })=\frac{1}{n}\sum _{i\in {\mathcal{I}}_k}\left[b({\mathbf{x}}_i^{\top }\mathbf{\theta })-y_i({\mathbf{x}}_i^{\top }\mathbf{\theta })\right].$$

$$\nabla f_k(\mathbf{\theta })=\frac{1}{n}\sum _{i\in {\mathcal{I}}_k}[b^{\prime }({\mathbf{x}}_i^{\top }\mathbf{\theta })-y_i]{\mathbf{x}}_i\text{mbox}and{\nabla }^2{f}_k(\mathbf{\theta })=\frac{1}{n}\sum _{i\in {\mathcal{I}}_k}{b}^{\prime \prime }({\mathbf{x}}_i^{\top }\mathbf{\theta }){\mathbf{x}}_i{\mathbf{x}}_i^{\top }.$$

$$\underset{k\in [m]}{\max }\underset{\mathbf{\theta }\in B(\hat{}\mathbf{\theta },R)}{\max }\Vert {\nabla }^2{f}_k(\mathbf{\theta })-{\nabla }^{2}f(\mathbf{\theta }){\Vert }_2=O_{\mathbb{P}}\left(\Vert \mathbf{\Sigma }{\Vert }_2\sqrt{\frac{p(\log p+\log N)}{n}}\right),$$

$$\Vert {\mathbf{\theta }}_t-{\mathbf{\theta }}^{\ast }{\Vert }_2\le \Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2+\Vert \hat{}\mathbf{\theta }-{\mathbf{\theta }}^{\ast }{\Vert }_2.$$

$$\Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2\le {\left[1-\frac{\rho }{10(\alpha +\rho )}\right]}^t\Vert {\mathbf{\theta }}_0-\hat{}\mathbf{\theta }{\Vert }_2,\forall t\ge 0;$$

$$\Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2=O_{\mathbb{P}}({\eta }^{t/2}\Vert {\mathbf{\theta }}_0-\hat{}\mathbf{\theta }{\Vert }_2),\forall t\ge 0;$$

$$\Vert {\mathbf{\theta }}_t-\hat{}\mathbf{\theta }{\Vert }_2=O_{\mathbb{P}}({\eta }^{t-t_0}\Vert {\mathbf{\theta }}_{t_0}-\hat{}\mathbf{\theta }{\Vert }_2),\forall t\ge t_0,$$

Full text

1 Introduction

Statistical inference in modern era faces tremendous challenge on computation and storage. The exceedingly large size of data often makes it impossible to store all of them on a single machine. Moreover, many applications have individual agents (e.g. local governments, research labs, hospitals, smart phones) collecting data independently. Communication is prohibitively expensive due to the limited bandwidth, and direct data sharing raises concerns in privacy and loss of ownership. These constraints make it necessary to develop methodologies for distributed systems, solving statistical problems with divide-and-conquer procedures and communicating only certain summary statistics.

Distributed statistical inference has received considerable attention recently, covering a wide spectrum of topics including $M$-estimation (Zhang et al., 2013; Chen and Xie, 2014; Shamir et al., 2014; Rosenblatt and Nadler, 2016; Wang et al., 2017a; Lee et al., 2017b; Battey et al., 2018; Wang et al., 2018; Shi et al., 2018; Banerjee et al., 2019), principal component analysis (Fan et al., 2019; Garber et al., 2017), nonparametric regression (Shang and Cheng, 2017; Szabó and Van Zanten, 2019; Han et al., 2018), quantile regression (Volgushev et al., 2019; Chen et al., 2021), bootstrap (Kleiner et al., 2014), confidence intervals (Jordan et al., 2019; Chen et al., 2021), Bayesian methods (Wang and Dunson, 2013; Jordan et al., 2019), etc. In the commonly-used setting, the overall dataset is partitioned and stored on $m$ node machines connected to a central processor. Most of the approaches only require one round of communication: the node machines work in parallel and send their results to the central processor, which then aggregates the information to get a final result. As typical examples, Zhang et al. (2013) average the $M$-estimators on node machines; Battey et al. (2018) average debiased estimators; and Fan et al. (2019) average subspaces via eigen-decomposition. While these one-shot methods are communication-efficient, they only work with a small number of node machines (e.g. $m=o(\sqrt{N})$, where $N$ is the total sample size) and require large sample on each, as their theories heavily rely on asymptotic expansions of estimators. Such conditions are easily violated in practice.

Multi-round procedures come as a remedy, which alternate between local computations and global aggregations. It is possible to achieve optimal statistical precision after a few rounds of communications, under broader settings than those for one-shot procedures. Shamir et al. (2014) propose a Distributed Approximate NEwton (DANE) algorithm where, in each iteration, each node machine minimizes a modified loss function based on its own samples and the gradient information from all other machines obtained through communication. However, for non-quadratic losses, the analysis in Shamir et al. (2014) does not imply any advantage of DANE in terms of communication over distributed implementation of gradient descent. Other approximate Newton algorithms include Zhang and Xiao (2015), Wang et al. (2018), Chen et al. (2021) and Crane and Roosta (2019). Jordan et al. (2019) develop a Communication-efficient Surrogate Likelihood (CSL) framework for estimation and inference in regular parametric models, penalized regression, and Bayesian statistics. A similar method also appears in Wang et al. (2017a). These methods no longer have restrictions on the number of machines such as $m=o(\sqrt{N})$.

Due to the nature of Newton-type methods, existing theories for these algorithms heavily rely on good initialization or even self-concordance assumption on loss functions. They essentially focus on improving an initial estimator that is already consistent but not efficient, whose ideas coincide with the classical one-step estimator (Bickel, 1975). Such initialization itself needs additional efforts and assumptions. Moreover, current results still require each machine to have sufficiently many samples so that loss functions on different machines are similar to each other. These all make the proposed methods unreliable in practice.

Aside from distributed statistical inference, there has also been a vast literature in distributed optimization. The ADMM (Boyd et al., 2011) is a celebrated example among the numerous algorithms that handle deterministic optimization problems with minimum structural assumption. Yet, the convergence can be quite slow and it cannot fully utilize the similarity among loss functions on node machines.

In this paper, we develop and study two Communication-Efficient Accurate Statistical Estimators (CEASE) based on multi-round algorithms for distributed statistical estimation. Our new algorithms extend the DANE algorithm (Shamir et al., 2014) to regularized empirical risk minimization. Moreover, we provide sharp convergence guarantees for general scenarios, even if the local loss functions are dissimilar and regularization is nonsmooth.

We assume that all the $m$ node machines have the same sample size $n$. Each has a regularized empirical risk function $f_{k}+g$ defined by the samples stored there, and the goal is to compute the minimizer of the overall regularized risk function $\frac{1}{m}\sum_{k=1}^{m}f_{k}+g$ to statistical precision. When $n$ is sufficiently large, their rates of convergence are better than or comparable to existing methods designed for this large-sample regime. For moderate or small $n$, they are still guaranteed to converge linearly even without good initialization, while other statistical methods fail. In addition, our algorithms take advantage of the similarity among $\{f_{k}\}_{k=1}^{m}$ and thus improve over general-purpose algorithms like ADMM. They interpolate between distributed algorithms for statistical estimation and general deterministic problems. Theoretical findings are verified by extensive numerical experiments. From a technical point of view, our algorithms use the proximal point algorithm (Rockafellar, 1976) as the backbone and obtain inexact updates in a distributed manner. This turns out to be crucial for proving convergence under general conditions. Our techniques are potentially useful for studying other distributed algorithms.

The rest of this paper is organized as follows. Section 2 introduces the algorithms. Section 3 presents deterministic convergence results. Section 4 provides guarantees in statistical problems. Section 5 shows numerical results on both synthetic and real data. Section 6 concludes the paper and discusses possible future directions.

Here we list the notations used throughout the paper. We denote by $[n]$ the set $\{1,2,\cdots,n\}$. We write $a_{n}=O(b_{n})$ or $a_{n}\lesssim b_{n}$ if there exists a constant $C>0$ such that $a_{n}\leq Cb_{n}$ holds for sufficiently large $n$; and $a_{n}\asymp b_{n}$ if $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. Given $\mathbf{x},\mathbf{y}\in{\mathbb{R}}^{k}$ and $r>0$, we define $B(\mathbf{x},r)=\{\mathbf{z}\in{\mathbb{R}}^{k}:\|\mathbf{z}-\mathbf{x}\|_{2}\leq r\}$ and $\langle\mathbf{x},\mathbf{y}\rangle=\sum_{j=1}^{k}x_{j}y_{j}$. For a convex function $h$ on ${\mathbb{R}}^{k}$, we let $\partial h(\mathbf{x})$ be its sub-differential set at $\mathbf{x}\in{\mathbb{R}}^{k}$, and ${\rm argmin}_{\mathbf{x}\in{\mathbb{R}}^{k}}h(\mathbf{x})$ be the set of its minimizers if $\inf_{\mathbf{x}\in{\mathbb{R}}^{k}}h(\mathbf{x})>-\infty$. We use $\|\cdot\|_{2}$ to denote the $\ell_{2}$ norm of a vector or operator norm of a matrix. For two sequences of random variables $\{X_{n}\}_{n=1}^{\infty}$ and $\{Y_{n}\}_{n=1}^{\infty}$ where $Y_{n}\geq 0$, we write $X_{n}=O_{{\mathbb{P}}}(Y_{n})$ if for any $\varepsilon>0$ there exists $C>0$ such that ${\mathbb{P}}(|X_{n}|\geq CY_{n})\leq\varepsilon$ for sufficiently large $n$. We use $\|X\|_{\psi_{2}}=\sup_{p\geq 1}\frac{1}{\sqrt{p}}{\mathbb{E}}^{1/p}|X|^{p}$ to refer to the sub-Gaussian norm of random variable $X$, and $\|\mathbf{X}\|_{\psi_{2}}=\sup_{\|u\|_{2}=1}\|\langle u,\mathbf{X}\rangle\|_{\psi_{2}}$ to denote the sub-Gaussian norm of random vector $\mathbf{X}$.

2 The CEASE algorithm

2.1 Problem setup

Let $\mathcal{P}$ be an unknown probability distribution over some sample space $\mathcal{X}$. For any parameter $\boldsymbol{\theta}\in{\mathbb{R}}^{p}$, define its population risk $F(\boldsymbol{\theta})={\mathbb{E}}_{\mathbf{X}\sim\mathcal{P}}\ell(\boldsymbol{\theta};\mathbf{X})$ based on a loss function $\ell:{\mathbb{R}}^{p}\times\mathcal{X}\to{\mathbb{R}}$. In parametric inference problems, $\ell$ is often chosen as the negative log-likelihood function of some parametric family. Under mild conditions, $F$ is well-defined and has a unique minimizer $\boldsymbol{\theta}^{*}$. A ubiquitous problem in statistics and machine learning is to estimate $\boldsymbol{\theta}^{*}$ given i.i.d. samples $\{\mathbf{X}_{i}\}_{i=1}^{N}$ from $\mathcal{P}$, and the minimizer of the empirical risk $f(\boldsymbol{\theta})=\frac{1}{N}\sum_{i=1}^{N}\ell(\boldsymbol{\theta};\mathbf{X}_{i})$ becomes a natural candidate. To achieve desirable precision in high-dimensional problems, it is often necessary to incorporate prior knowledge of $\boldsymbol{\theta}^{*}$. A principled approach is the regularized empirical risk minimization

[TABLE]

where $g(\boldsymbol{\theta})$ is a deterministic penalty function. Common choices for $g(\boldsymbol{\theta})$ include the $\ell_{2}$ penalty $\lambda\|\boldsymbol{\theta}\|_{2}^{2}$ (Hoerl and Kennard, 1970), the $\ell_{1}$ penalty $\lambda\|\boldsymbol{\theta}\|_{1}$ (Tibshirani, 1996), and a family of folded concave penalty functions $\|p_{\lambda}(|\boldsymbol{\theta}|)\|_{1}$ such as SCAD (Fan and Li, 2001) and MCP (Zhang, 2010), where $\lambda>0$ is a regularization parameter. Throughout the paper, we assume that both $\ell$ and $g$ are convex in $\boldsymbol{\theta}$, and $\ell$ is twice continuously differentiable in $\boldsymbol{\theta}$. We allow $g$ to be non-smooth (e.g. the $\ell_{1}$ penalty).

Consider the distributed setting where the $N$ samples are stored on $m$ machines connected to a central processor. Denote by $\mathcal{I}_{k}$ the index set of samples on the $k$th machine and $f_{k}(\boldsymbol{\theta})=\frac{1}{|\mathcal{I}_{k}|}\sum_{i\in\mathcal{I}_{k}}\ell(\boldsymbol{\theta};\mathbf{X}_{i})$. For simplicity, we assume that $\{\mathcal{I}_{k}\}_{k=1}^{m}$ are disjoint, $N$ is a multiple of $m$, and $|\mathcal{I}_{k}|=n=N/m$ for all $k\in[m]$. Then (2.1) can be rewritten as

[TABLE]

Each machine $k$ only has access to its local data and hence local loss function $f_{k}$ and the penalty $g$. We aim to solve (2.2) in a distributed manner with both statistical efficiency and communication-efficiency.

2.2 Adaptive gradient enhancements and distributed algorithms in large-sample regimes

In the large-sample regime, we drop the regularization term for now and consider the empirical risk minimization problem $\min_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}f(\boldsymbol{\theta})$ for estimating $\boldsymbol{\theta}^{*}={\rm argmin}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}F(\boldsymbol{\theta})$. In some problems, direct minimization of $f$ is costly, while it is easy to obtain some rough estimate $\bar{}\boldsymbol{\theta}$ that is close to $\boldsymbol{\theta}^{*}$ but not as accurate as the global minimimizer $\widehat{}\boldsymbol{\theta}={\rm argmin}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}f(\boldsymbol{\theta})$. Bickel (1975) proposes the one-step estimator based on the local quadratic approximation and shows that it is as efficient as $\widehat{}\boldsymbol{\theta}$ if the initial estimator $\bar{}\boldsymbol{\theta}$ is accurate enough. Iterating this further results in multiple-step estimators that improve the optimization error and hence statistical errors when the initial estimator is not good enough (Robinson, 1988). This inspires us to refine an existing estimator using some proxy of $f$.

In the distributed environment, starting from an initial estimator $\bar{}\boldsymbol{\theta}$, the gradient vector $\nabla f(\bar{}\boldsymbol{\theta})$ can easily be communicated. Construct a linear function $f^{(1)}(\boldsymbol{\theta})=f(\bar{}\boldsymbol{\theta})+\langle\nabla f(\bar{}\boldsymbol{\theta}),\boldsymbol{\theta}-\bar{}\boldsymbol{\theta}\rangle$, the first-order Taylor expansion of $f$ around $\bar{}\boldsymbol{\theta}$. The object function to be minimized can be written as

[TABLE]

Since the linear function $f^{(1)}(\boldsymbol{\theta})$ can easily be communicated to each node machine whereas $R(\cdot)$ can not, the latter is naturally replaced by its subsampled version at node $k$:

[TABLE]

where $f_{k}(\boldsymbol{\theta})$ is the loss function based on the data at node $k$. With this replacement, the target of optimization at node $k$ becomes $f^{(1)}(\boldsymbol{\theta})+R_{k}(\boldsymbol{\theta})$, which equals to

[TABLE]

up to an additive constant. This function will be called gradient-enhanced loss function, in which the gradient at point $\bar{}\boldsymbol{\theta}$ based on the local data is replaced by the global one. This function has one very nice fixed point at the global minimum $\widehat{}\boldsymbol{\theta}$: the minimizer of the adaptive gradient-enhanced function at $\bar{}\boldsymbol{\theta}=\widehat{}\boldsymbol{\theta}$ is still $\widehat{\boldsymbol{\theta}}$. This can easily be seen by verifying that the gradient at the point $\widehat{\boldsymbol{\theta}}$ is zero.

The idea of using such an adaptive gradient-enhanced function has been proposed in Shamir et al. (2014) and Jordan et al. (2019), though the motivations are different. Jordan et al. (2019) develop a Commmunication-efficient Surrogate Likelihood (CSL) method using the gradient-enhanced loss function $f_{1}(\boldsymbol{\theta})-\langle\nabla f_{1}(\bar{}\boldsymbol{\theta})-\nabla f(\bar{}\boldsymbol{\theta}),\boldsymbol{\theta}\rangle$ on the first machine, uses the minimizer on that machine as a new estimate, and iterates these steps until convergence. In the presence of a regularizer $g$ in (2.1), one simply adds $g$ to the gradient-enhanced loss; see the Algorithm 1 below.

Note that in Algorithm 1, only the first machine solves optimization problems and others just evaluate gradients. These machines are idling while the first one is working hard. To fully utilize the computing power of machines and accelerate convergence, all the machines can optimize their corresponding gradient-enhanced loss functions in parallel and the central processor then aggregates the results. This is motivated by the Distributed Approximate NEwton (DANE) algorithm (Shamir et al., 2014). Algorithm 2 describes the procedure in detail. Intuitively, the averaging step requires little computation but helps reduce the variance of estimators on node machines and enhance the accuracy.

We now illustrate Algorithm 2 in the context of linear regression. Given samples $\{(\mathbf{x}_{i},y_{i})\}_{i\in[N]}$, the $k$th machine defines a quadratic loss function

[TABLE]

Here, $\widehat{}\boldsymbol{\Sigma}_{k}=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}$ and $\widehat{}\mathbf{w}_{k}=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}\mathbf{x}_{i}y_{i}$. The overall loss function is $f(\boldsymbol{\theta})=\frac{1}{m}\sum_{k=1}^{m}f_{k}(\boldsymbol{\theta})=\frac{1}{2}\boldsymbol{\theta}^{\top}\widehat{\boldsymbol{\Sigma}}\boldsymbol{\theta}-\widehat{\mathbf{w}}^{\top}\boldsymbol{\theta}$, where $\widehat{}\boldsymbol{\Sigma}=\frac{1}{m}\sum_{k=1}^{m}\widehat{}\boldsymbol{\Sigma}_{k}$, $\widehat{}\mathbf{w}=\frac{1}{m}\sum_{k=1}^{m}\widehat{}\mathbf{w}_{k}$. Then the update of Algorithm 2 in one iteration is

[TABLE]

Intuitively, this is a form of contraction towards the global minimizer $\widehat{}\boldsymbol{\theta}$. As for the logistic regression, we can also write out the corresponding enhanced losses and minimize them using Newton’s method. Due to space limitations, we refer to Appendix C for details.

2.3 The CEASE Algorithm in general regimes

Algorithms 1 and 2 are built upon large-sample regimes, with sufficiently strong convexity of $\{f_{k}+g\}_{k=1}^{m}$ and small discrepancy between them. This requires the local sample size $n$ to be large enough, which may not be the case in practice. Even worse, the required local sample size depends on structural parameters, making such a condition unverifiable. In fact, our numerical experiments confirm the instability of Algorithms 1 and 2 even for moderate $n$. A naive method of remedy is to add strict convex quadratic regularization $q(\boldsymbol{\theta})$. While this remedy can make the algorithm converge rapidly, the nonadaptive nature of $q(\boldsymbol{\theta})$ will lead to a wrong target. Instead of using a fixed $q$, we will adjust it according to current solutions. The idea stems from the proximal point algorithm (Rockafellar, 1976).

Definition 2.1**.**

For any convex function $h:~{}{\mathbb{R}}^{p}\to{\mathbb{R}}$, define the proximal mapping $\mathrm{prox}_{h}:~{}{\mathbb{R}}^{p}\to{\mathbb{R}}^{p}$, $\mathbf{x}\mapsto{\rm argmin}_{\mathbf{y}\in{\mathbb{R}}^{p}}\{h(\mathbf{y})+\|\mathbf{y}-\mathbf{x}\|_{2}^{2}/2\}$.

For a given $\alpha>0$, the proximal point algorithm for minimizing $h$ iteratively computes

[TABLE]

starting from some initial value $\mathbf{x}_{0}$. It is a strongly convex optimization, shrinking towards the current value $\mathbf{x}_{t}$. Under mild conditions, $\{\mathbf{x}_{t}\}_{t=0}^{\infty}$ converges linearly to some $\widehat{}\mathbf{x}\in{\rm argmin}_{{\mathbb{R}}^{p}}h(\mathbf{x})$ (Rockafellar, 1976).

Now we take $h=f+g$ and write the proximal point iteration for our problem (2.2):

[TABLE]

Each iteration (2.5) is a distributed optimization problem, whose object function is not available to node machines. But it can be solved by Algorithms 1 and 2. Specifically, suppose we have already obtained $\boldsymbol{\theta}_{t}$ and aim for $\boldsymbol{\theta}_{t+1}$ in (2.5). Letting $\tilde{g}(\boldsymbol{\theta})=g(\boldsymbol{\theta})+(\alpha/2)\|\boldsymbol{\theta}-\boldsymbol{\theta}_{t}\|_{2}^{2}$, Algorithm 2 starting from $\tilde{\boldsymbol{\theta}}_{0}=\boldsymbol{\theta}_{t}$ produces iterations over $s=0,1,\cdots$

[TABLE]

When $\alpha+\rho_{0}>\delta$, $\{\tilde{\boldsymbol{\theta}}_{s}\}_{s=0}^{\infty}$ converges $Q$-linearly ***According to Nocedal and Wright (2006), a sequence $\{\mathbf{x}_{n}\}_{n=1}^{\infty}$ in ${\mathbb{R}}^{p}$ is said to converge $Q$-linearly to $\mathbf{x}^{*}\in{\mathbb{R}}^{p}$ if there exists $r\in(0,1)$ such that $\|\mathbf{x}_{n+1}-\mathbf{x}^{*}\|_{2}\leq r\|\mathbf{x}_{n}-\mathbf{x}^{*}\|_{2}$ for $n$ sufficiently large. to $\boldsymbol{\theta}_{t+1}$. On the other hand, there is no need to solve (2.5) exactly, as $\mathrm{prox}_{\alpha^{-1}(f+g)}(\boldsymbol{\theta}_{t})$ is merely an intermediate quantity for computing $\widehat{}\boldsymbol{\theta}$. We therefore only run one iteration of Algorithm 2 and use the resulting approximate solution as $\boldsymbol{\theta}_{t+1}$. This considerably simplifies the algorithm, reducing double loops to a single loop, and enhances statistical interpretation of the method as a multi-step estimator. However, it makes technical analysis more challenging. Similarly, we may also use one step of Algorithm 1 to compute the inexact proximal update.

The above discussions lead us to propose two Communication-Efficient Accurate Statistical Estimators (CEASE) in Algorithms 3 and 4, which use the proximal point algorithm as the backbone and obtain inexact updates in a distributed manner. They are regularized versions of Algorithms 1 and 2, with an additional proximal term in the objective functions. That term reduces relative differences of the local loss functions on individual machines, and is crucial for convergence when $\{f_{k}\}_{k=1}^{m}$ are not similar enough. In Appendix A we introduce a variant of Algorithm 4 which also stablizes Algorithm 2.

In each iteration, Algorithm 3 has one round of communication and one optimization problem to solve. Although Algorithm 4 has two rounds of communication per iteration, only one round involves parallel optimization and the other is simply averaging. We will compare their theoretical guarantees as well as practical performances in the sequel.

Algorithm 4 is an extension of the DANE algorithm in Shamir et al. (2014) to regularized empirical risk minimization. While DANE is originally motivated by mirror descent, we view it as a distributed implementation of the proximal point algorithm. The new perspective helps us obtain stronger convergence guarantees. Ideas from the proximal point algorithm have appeared in the literature of distributed stochastic optimization for different purposes such as accelerating first-order algorithms (Lee et al., 2017a) and regularizing sizes of updates (Wang et al., 2017b).

3 Deterministic analysis

We first present in Section 3.1 the deterministic (almost sure) results for Algorithms 3 and 4 based on high-level structural assumptions. As special cases of these algorithms with $\alpha=0$, Algorithms 1 and 2 will be analyzed in Section 3.2.

3.1 Deterministic analysis of the CEASE algorithm

Definition 3.1**.**

Let $h:{\mathbb{R}}^{p}\to{\mathbb{R}}$ be a convex function, $\Omega\subseteq{\mathbb{R}}^{p}$ be a convex set, and $\rho\geq 0$. $h$ is $\rho$-strongly convex in $\Omega$ if $h(\mathbf{y})\geq h(\mathbf{x})+\langle\mathbf{g},\mathbf{y}-\mathbf{x}\rangle+\frac{\rho}{2}\|\mathbf{y}-\mathbf{x}\|_{2}^{2}$, $\forall\mathbf{x},\mathbf{y}\in\Omega$ and $\mathbf{g}\in\partial h(\mathbf{x})$.

Assumption 3.1** (Strong convexity).**

$f+g$* has a unique minimizer $\widehat{}\boldsymbol{\theta}\in{\mathbb{R}}^{p}$, and is $\rho$-strongly convex in $B(\widehat{}\boldsymbol{\theta},R)$ for some $R>0$ and $\rho>0$.*

Assumption 3.2** (Homogeneity).**

$\|\nabla^{2}f_{k}(\boldsymbol{\theta})-\nabla^{2}f(\boldsymbol{\theta})\|_{2}\leq\delta$, $\forall k\in[m],\boldsymbol{\theta}\in B(\widehat{}\boldsymbol{\theta},R)$.

We will refer to $\delta$ as a homogeneity parameter. Based on both assumptions, we define

[TABLE]

A simple but useful fact is $\max\{\rho-\delta,0\}\leq\rho_{0}\leq\rho$. In most interesting problems, the population risk $F$ is smooth and strongly convex on any compact set. When $\{\mathbf{X}_{i}\}_{i=1}^{N}$ are i.i.d. and the total sample size $N$ is large, the empirical risk $f$ concentrates around $F$ and inherits nice properties from the latter, making Assumption 3.1 hold easily.

On the other hand, since the empirical risk functions $\{f_{k}\}_{k=1}^{m}$ are i.i.d. stochastic approximations of the population risk $F$, they should not be too far away from their average $f$ provided that $n$ is not too small. Assumption 3.2 is a natural way of characterizing this similarity. It is a generalization of the concept “$\delta$-related functions” for quadratic losses in Arjevani and Shamir (2015). With high probability, it holds with reasonably small $\delta$ and large $R$ under general conditions. Large $n$ implies small homogeneity parameter $\delta$ and thus similar $\{f_{k}\}_{k=1}^{m}$. Assumption 3.2 always holds with $\delta=\max_{k\in[m]}\{\sup_{\boldsymbol{\theta}\in B(\widehat{}\boldsymbol{\theta},R)}\|\nabla^{2}f_{k}(\boldsymbol{\theta})-{\mathbb{E}}\nabla^{2}f_{k}(\boldsymbol{\theta})\|_{2}+\sup_{\boldsymbol{\theta}\in B(\widehat{}\boldsymbol{\theta},R)}\|\nabla^{2}f(\boldsymbol{\theta})-{\mathbb{E}}\nabla^{2}f_{k}(\boldsymbol{\theta})\|_{2}\}$.

The following additional assumption on smoothness of the Hessian matrix of $f+g$ is not necessary for contraction, but it helps us obtain a much stronger result on the contraction rate of Algorithm 2, justifying the power of the simple averaging step.

Assumption 3.3** (Smoothness of Hessian).**

$g\in C^{2}({\mathbb{R}}^{p})$, and there exists $M\geq 0$ such that $\|[\nabla^{2}f(\boldsymbol{\theta}^{\prime})+\nabla^{2}g(\boldsymbol{\theta}^{\prime})]-[\nabla^{2}f(\boldsymbol{\theta}^{\prime\prime})+\nabla^{2}g(\boldsymbol{\theta}^{\prime\prime})]\|_{2}\leq M\|\boldsymbol{\theta}^{\prime}-\boldsymbol{\theta}^{\prime\prime}\|_{2}$, $\forall\boldsymbol{\theta}^{\prime},\boldsymbol{\theta}^{\prime\prime}\in B(\widehat{}\boldsymbol{\theta},R)$.

Theorem 3.1 gives contraction guarantees for Algorithms 3 and 4. It is deterministic and non-asymptotic by nature.

Theorem 3.1**.**

Let Assumptions 3.1 and 3.2 hold. Consider the multi-step estimators $\{\boldsymbol{\theta}_{t}\}_{t=0}^{T}$ generated by Algorithm 3 or 4. Suppose that $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R/2)$ and $[\delta/(\rho_{0}+\alpha)]^{2}<\rho/(\rho+2\alpha)$.

1. (a)

For both Algorithms 3 and 4, we have

[TABLE] 2. (b)

If Assumption 3.3 also holds, then for Algorithm 4 we have

[TABLE]

where we define $\gamma_{t}=\frac{\delta}{\rho_{0}+\alpha}\cdot\min\{1,\frac{\delta}{\rho+\alpha}(1+\frac{M}{\rho_{0}+\alpha}\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2})\}$; 3. (c)

Both multiplicative factors in (3.2) and (3.3) are strictly less than 1.

In the contraction factor in (3.2), the two summands $\frac{\delta\sqrt{\rho^{2}+2\alpha\rho}}{(\rho_{0}+\alpha)^{2}}$ and $\frac{\alpha}{\rho+\alpha}$ come from bounding the inexact proximal update $\|\boldsymbol{\theta}_{t+1}-\mathrm{prox}_{\alpha^{-1}(f+g)}(\boldsymbol{\theta}_{t})\|_{2}$ and the residual $\|\mathrm{prox}_{\alpha^{-1}(f+g)}(\boldsymbol{\theta}_{t})-\widehat{}\boldsymbol{\theta}\|_{2}$, respectively. Similar results hold for (3.3). The condition $[\delta/(\rho_{0}+\alpha)]^{2}<\rho/(\rho+2\alpha)$ ensures that both contraction factors are less than 1. Note that (3.3) requires Assumption 3.3, which forces $g$ to be smooth.

Theorem 3.1 shows the $Q$-linear convergence of the sequence $\{\boldsymbol{\theta}_{t}\}_{t=0}^{\infty}$ generated by both Algorithms 3 and 4 under quite general settings. The contraction rate depends explicitly on the structural parameters and the choice of $\alpha$. The local loss functions $\{f_{k}\}_{k=1}^{m}$ just need to be convex and smooth, and the convex penalty $g$ is allowed to be non-smooth, e.g. the $\ell_{1}$ norm. On the contrary, most algorithms for distributed statistical estimation are only designed for smooth problems, and many of them are only rigorously studied when the loss functions are quadratic or self-concordant (Shamir et al., 2014; Zhang and Xiao, 2015; Wang et al., 2017b). This is another important aspect of our contributions.

We immediately see from Theorem 3.1 that Algorithms 3 and 4 converge linearly as long as $[\delta/(\rho_{0}+\alpha)]^{2}<\rho/(\rho+2\alpha)$, which is guaranteed to hold by choosing sufficiently large $\alpha$. By contrast, however, we’ll show in Section 3.2 that Algorithms 1 and 2 (corresponding to $\alpha=0$) hinges on the homogeneity assumption $\rho_{0}>\delta$ in Theorem 3.2, i.e. the functions $\{f_{k}\}_{k=1}^{m}$ must be similar enough. In the statistical setting, this requires the local sample size $n$ to be large. Therefore, proper regularization provides a safety net for the algorithms under general regimes with potentially insufficient local sample size. Corollary 3.1 below gives a guideline for choosing $\alpha$ to make Algorithms 3 and 4 converge in general.

Corollary 3.1**.**

Let Assumptions 3.1 and 3.2 hold, $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R/2)$, and $\{\boldsymbol{\theta}_{t}\}_{t=0}^{T}$ be the iterates of Algorithm 3 or 4. With any $\alpha\geq 4\delta^{2}/\rho$, both algorithms converge with contraction factors in (3.2) and (3.3) bounded by $(1-\frac{\rho}{10(\alpha+\rho)})$.

On the other hand, consider the case where the local loss functions have small relative difference $\delta/\rho$. In this case, Theorem 3.2 states that the contraction factors for unregularized versions ($\alpha=0$) of Algorithms 3 and 4 are in the same order of $\delta/\rho$ and $(\delta/\rho)^{2}$, respectively, which are smaller than the contracting factors with $\alpha>0$. The following corollary characterizes the upper bound for $\alpha$ so that the contraction factors remain at these small orders.

Corollary 3.2**.**

Let Assumptions 3.1 and 3.2 hold, $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R/2)$, and suppose $\alpha\leq C\delta^{2}/\rho$ for some constant $C$. There exist constants $C_{1}$ and $C_{2}$ such that the followings hold when $\delta/\rho$ is sufficiently small:

• (a)

Algorithms 3 and 4 have the contraction property

[TABLE]

• (b)

If Assumption 3.3 also holds and $\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2}\leq\rho/M$, then for Algorithm 4

[TABLE]

Note that the second result above only holds given Assumption 3.3, which requires a smooth regularization $g$. Corollary 3.2 reveals that by choosing $\alpha\asymp\delta^{2}/\rho$, the contraction factors are essentially the same as those of the unregularized ($\alpha=0$) algorithms when $\delta/\rho$ is small. By combining Corollaries 3.1 and 3.2, we use $\alpha\asymp\delta^{2}/\rho$ as a default choice for Algorithms 3 and 4 to become both fast and robust. They are reliable in general cases (Corollary 3.1) and efficient in nice cases (Corolary 3.2 and Theorem 3.1 with $\alpha=0$)

Algorithms 3 and 4 attain communication efficiency by utilizing similarity among local loss functions: The contraction factors in Corollary 3.2 go to zero if $\delta/\rho$ does. In fact, both algorithms achieve $\varepsilon$-accuracy within $O(\max\{1,(\delta/\rho)^{2}\}\log(\frac{\|\boldsymbol{\theta}_{0}-\widehat{}\boldsymbol{\theta}\|_{2}}{\varepsilon}))$ rounds of communication. In contrast, the distributed accelerated gradient descent requires $O(\sqrt{\kappa_{0}}\log(\frac{\|\boldsymbol{\theta}_{0}-\widehat{}\boldsymbol{\theta}\|_{2}}{\varepsilon}))$ rounds of communication to achieve $\varepsilon$-accuracy (Shamir et al., 2014), with $\kappa_{0}$ being the condition number of $(f+g)$, which does not take advantage of sample size $n$. As long as $\delta/\rho\ll\kappa_{0}^{1/4}$, Algorithms 3 and 4 communicate less than the distributed accelerated gradient descent. And again, our general results for Algorithms 3 and 4 also apply to the case with nonsmooth penalty functions while those for distributed accelerated gradient descent do not.

Moreover, if $(f+g)$ is smooth and $\boldsymbol{\theta}_{t}$ is reasonably close to $\widehat{}\boldsymbol{\theta}$, Corollary 3.2 shows that each iteration of Algorithm 4 is roughly equivalent to two iterations of Algorithm 3, although the former only has one round of optimization. The averaging step in Algorithm 4 reduces the error as much as the optimization step, while taking much less time. In this case, Algorithm 4 is preferable, and our numerical experiments also confirm this.

For unregularized empirical risk minimization, i.e. $g=0$ in (2.2), Algorithm 4 reduces to an extension or a useful case of the DANE algorithm (Shamir et al., 2014). In this case, Theorem 3.1 and its corollaries refine the analysis of DANE (Shamir et al., 2014) in several aspects. On the one hand, our analysis handle both smooth and nonsmooth problems, while in Shamir et al. (2014), the theoretical analysis beyond quadratic loss requires extremal choice of tuning parameters and does not show any advantage over distributed implementation of the gradient descent. On the other hand, as mentioned in Section 2, we derive Algorithm 4 from the proximal point algorithm with a new prospective, which leads to sharp convergence analysis along with suggestions on choosing the tuning parameter $\alpha$. As a by-product, we close a gap in the theory of DANE in non-quadratic settings. Our analysis techniques are potentially useful for other distributed optimization algorithms, especially when the loss is not quadratic.

3.2 Deterministic analysis in large-sample regimes

In this section, we restrict ourselves to large-sample regimes where the local sample size $n$ is sufficiently large such that $\rho_{0}>\delta\geq 0$, where $\rho_{0}$ is the strong convexity parameter in (3.1). The following theorem gives deterministic results for Algorithms 1 and 2.

Theorem 3.2**.**

Let Assumptions 3.1 and 3.2 hold, and $\rho_{0}>\delta\geq 0$. Consider the iterates $\{\boldsymbol{\theta}_{t}\}_{t=0}^{\infty}$ produced by Algorithm 1 or 2, with $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R)$. Then

[TABLE]

In addition, if Assumption 3.3 also holds, then for Algorithm 2 we have

[TABLE]

Note that the last inequality requires Assumption 3.3 and thus a smooth regularization $g$. The first part of Theorem 3.2 is a refinement of the analysis in Jordan et al. (2019), since we allow the initial estimator to be inaccurate and we have more explicit rates of contraction of optimization errors. This will be further demonstrated in Section 4.2.

The second part points out benefits of the averaging step, which is a novel result. Similar to the results on Algorithm 4, Theorem 3.2 shows that when $\boldsymbol{\theta}_{t}$ is close to $\widehat{}\boldsymbol{\theta}$, with an additional standard assumption on Hessian smoothness, the averaging step alone in Algorithm 2 is almost as powerful as an optimization step in terms of contraction: The contracting constant will eventually be $\frac{\delta}{\rho_{0}}\frac{\delta}{\rho}$. With negligible computational cost, averaging significantly improves upon individual solutions $\{\boldsymbol{\theta}_{t,k}\}_{k=1}^{m}$ by doubling the speed of convergence.

4 Statistical analysis

We further analyze the statistical properties of the above algorithms under a generalized linear model. Essentially, both the CSL methods and the CEASE algorithm are $T$-step estimators, starting from the initial estimator $\boldsymbol{\theta}_{0}$. The question here is the effect of iterations in the multiple step estimators and the role of the initial estimator. We start with statistical analysis of Algorithms 3 and 4 in Section 4.1, and then study Algorithms 1 and 2 in Section 4.2. In Section 4.3, we provide practical guidance when implementing the CEASE algorithm based on these analysis.

4.1 Multi-step estimators in general regimes

The deterministic analysis in Section 3.2 applies to a wide range of statistical models. Here we consider the generalized linear model with canonical link, where our samples are i.i.d. pairs $\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{N}$ of covariates and responses and the conditional density of $y_{i}$ given $\mathbf{x}_{i}$ is given by

[TABLE]

For simplicity, we let the dispersion parameter to be 1 as we do not consider the issue of over-dispersion; $b(\cdot)$ is some known convex function, and $c$ is a known function such that $h$ is a valid probability density function. The negative log-likelihood of the whole data is an affine transformation of $f(\boldsymbol{\theta})=\frac{1}{m}\sum_{k=1}^{m}f_{k}(\boldsymbol{\theta})$ with

[TABLE]

It is easy to verify that

[TABLE]

Assume that $\mathbf{x}_{i}=(1,\mathbf{u}_{i}^{\top})^{\top}\in{\mathbb{R}}^{p}$, where $\{\mathbf{u}_{i}\}_{i=1}^{N}\subseteq{\mathbb{R}}^{p-1}$ are i.i.d. random covariate vectors with zero mean and covariance matrix $\boldsymbol{\Sigma}$. Suppose there exist universal positive constants $A_{1}$, $A_{2}$ and $A_{3}$ such that $A_{1}\leq\lVert\boldsymbol{\Sigma}\rVert_{2}\leq A_{2}p^{A_{3}}$. Let $\boldsymbol{\Sigma}^{*}={\mathbb{E}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\top})=\begin{pmatrix}1&\mathbf{0}\\ \mathbf{0}&\boldsymbol{\Sigma}\end{pmatrix}$, $g$ be a deterministic penalty function, and $F(\boldsymbol{\theta})={\mathbb{E}}f(\boldsymbol{\theta})$ be the population risk function. Below we impose some standard regularity assumptions.

Assumption 4.1**.**

• •

$\{\boldsymbol{\Sigma}^{-1/2}\mathbf{u}_{i}\}_{i=1}^{N}$* are i.i.d. sub-Gaussian random vectors.*

• •

For all $x\in{\mathbb{R}}$, $|b^{\prime\prime}(x)|$ and $|b^{\prime\prime\prime}(x)|$ are bounded by some constant.

• •

$\|\boldsymbol{\theta}^{*}\|_{2}$* is bounded by some constant.*

As in Assumptions 3.1 and 3.3, the following general assumptions are also needed for our analysis. Here $R$ is some positive quantity that satisfies $R<A_{4}p^{A_{5}}$ for some universal constants $A_{4}$ and $A_{5}$.

Assumption 4.2**.**

There exists a universal constant $\rho>0$ such that $(F+g)$ is $\rho$-strongly convex in $B(\boldsymbol{\theta}^{*},2R)$.

The following smoothness assumption is only needed for a part of our theory; it is used to show that the averaging step in Algorithm 4 can significantly enhance the accuracy.

Assumption 4.3**.**

$g\in C^{2}({\mathbb{R}}^{p})$, and there exists a universal constant $M\geq 0$ such that $\|[\nabla^{2}F(\boldsymbol{\theta}^{\prime})+\nabla^{2}g(\boldsymbol{\theta}^{\prime})]-[\nabla^{2}F(\boldsymbol{\theta}^{\prime\prime})+\nabla^{2}g(\boldsymbol{\theta}^{\prime\prime})]\|_{2}\leq M\|\boldsymbol{\theta}^{\prime}-\boldsymbol{\theta}^{\prime\prime}\|_{2}$, $\forall\boldsymbol{\theta}^{\prime},\boldsymbol{\theta}^{\prime\prime}\in B(\boldsymbol{\theta}^{*},2R)$.

Under the model assumptions above, we can explicitly determine the rate of $\delta$ in Assumption 3.2. In particular, we will show in Lemma E.5 in Appendix E that

[TABLE]

provided that $n\geq cp$ for an arbitrary positive constant $c$. Therefore, with high probability, $\delta\asymp\lVert\boldsymbol{\Sigma}\rVert_{2}\sqrt{p(\log p+\log N)/{n}}$. Omitting the logarithmic terms, we see that the contraction factor is approximately $\kappa\sqrt{p/n}$, where $\kappa\triangleq\lVert\boldsymbol{\Sigma}\rVert_{2}/\rho$ can be viewed as a condition number. This rate is more explicit on $p$ and $\kappa$ than that in Jordan et al. (2019), where finite $p$ and $\kappa$ are assumed. In addition, with a smooth regularization, Algorithm 4 benefits from the averaging step in that it improves the contraction rate to approximately $\kappa^{2}p/n$.

Let $\boldsymbol{\theta}_{t}$ be the $t$-th iterate of one of the proposed algorithms. It is clear that the statistical error of the estimator $\boldsymbol{\theta}_{t}$ is bounded by its optimization error and the statistical error of $\widehat{}\boldsymbol{\theta}$:

[TABLE]

The second term is well-studied in statistics, which is of order $O_{{\mathbb{P}}}(\sqrt{p/N})$ under mild conditions. In the following theorem, we show that for the first term (i.e. the optimization error), with proper choice of $\alpha$, each iteration of Algorithm 3 or 4 makes $\boldsymbol{\theta}_{t}$ closer to the global minimum $\widehat{}\boldsymbol{\theta}$ by some order depending on local sample size. Thus, through finite steps, the optimization errors are eventually negligible in comparison with statistical errors (assuming $N$ is of order $(n/p)^{a}$ for a finite $a$ in typical applications), and the distributed multi-step estimator will work as well as the global minimum as if the data were aggregated in the central server.

Theorem 4.1**.**

Suppose that Assumptions 4.1 and 4.2 hold, and with high probability the initial value $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R/2)$. Let $\eta=\kappa^{2}(\log N)p/n$ and $\kappa=\lVert\boldsymbol{\Sigma}\rVert_{2}/\rho$. For any $c_{1},c_{2}>0$, there exists $C>0$ such that the followings hold with high probability:

• (a)

If $n\geq c_{1}p$ and $\alpha\geq C\rho\eta$, then both Algorithms 3 and 4 have linear convergence

[TABLE]

• (b)

If $\eta$ is sufficiently small and $\alpha\leq c_{2}\rho\eta$, then for both algorithms

[TABLE]

in addition, if Assumption 4.3 also holds, then for Algorithm 4 we have

[TABLE]

where $t_{0}=\lceil\frac{2\log(CMR/\rho)}{\log(1/\eta)}\rceil$.

In contrast to a fixed contraction derived by Shamir et al. (2014), Theorem 4.1 explains the significant benefits of large local sample size even in the presence of a non-smooth penalty: the optimization error shrinks by a factor that converges to zero explicitly in $n$. As a brief illustration, let us consider the case with smooth loss functions, sufficient local sample size and no regularization. Let $\boldsymbol{\theta}_{0}$ be the average of individual estimators on node machines. By Corollary 2 in Zhang et al. (2013), this simple divide-and-conquer estimator has accuracy $\lVert\boldsymbol{\theta}_{0}-\boldsymbol{\theta}^{*}\rVert_{2}=O_{{\mathbb{P}}}(\max\{\sqrt{\frac{p}{N}},\frac{p}{n},\frac{\kappa\sqrt{p\log p}}{n}\})$. Using the explicit expression of $\eta$, we can easily deduce from Theorem 4.1 (b) that the one-step estimator $\boldsymbol{\theta}_{1}$ obtained by Algorithm 3 behaves the same as the global minimizer $\widehat{}\boldsymbol{\theta}$ if the local sample size satisfies $n^{3}\gg N(\kappa^{2}p\log N)(p+\kappa^{2}\log p).$ In this case, the local optimization in Algorithm 3 can further be replaced by using the explicit one-step estimator as in Bickel (1975) and Jordan et al. (2019), since the initial estimator is in a consistent neighborhood. More generally, the $t$-step estimator $\boldsymbol{\theta}_{t}$ has negligible optimization error under even weaker local sample size requirement: $n^{t+2}\gg N(\kappa^{2}p\log N)^{t}(p+\kappa^{2}\log p).$ A similar remark applies to Algorithm 4.

Similar to the deterministic results, the averaging step is about as effective as the optimization step when $g$ is smooth and $\boldsymbol{\theta}_{t}$ is sufficiently close to $\widehat{}\boldsymbol{\theta}$ in that after a finite $t_{0}$ iterations. See a simplified example in Section 4.2 with $\alpha=0$.

As for the initialization, the condition $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R/2)$ is mild, since $\widehat{}\boldsymbol{\theta}$ is usually a consistent estimate and $\|\boldsymbol{\theta}^{*}\|_{2}$ is bounded (Assumption 4.1). In contrast with Jordan et al. (2019), we allow inaccurate initial value such as $\boldsymbol{\theta}_{0}=\mathbf{0}$ and give more explicit rates of contraction even when $p$ and $\kappa$ diverge. On the other hand, the accuracy of the initial estimator $\boldsymbol{\theta}_{0}$ does help reduce the number of iterations.

Combining the results to be presented in Theorem 4.2, we’ll see that by choosing $\alpha\asymp\rho\eta$, Algorithms 3 and 4 inherit all the merits of Algorithms 1 and 2 in the large-$n$ regime – fast linear contraction of rate $\sqrt{\eta}=\kappa\sqrt{p(\log N)/n}$, and for Algorithm 4, a even faster rate of $\eta=\kappa^{2}p(\log N)/n$ to $\widehat{}\boldsymbol{\theta}$ when the loss and the penalty functions are smooth. These facts also guarantee that Algorithms 3 and 4 reach the statistical efficiency in $O(\frac{\log\lVert\boldsymbol{\theta}_{0}-\widehat{}\boldsymbol{\theta}\rVert_{2}+\log(N/p)}{\log(1/\eta)})$ iterations. On the other hand, compared to Algorithms 1 and 2, Algorithms 3 and 4 overcome the difficulties with a small local sample size $n$ in that as long as $n/p$ is bounded away from some small constant (which is reasonable for many big-data problems of interest), shrinkage of optimization error is guaranteed. Moreover, while it is hard to check whether $n$ is sufficiently large in practice, proper choice of $\alpha$ always guarantees linear convergence, and the contraction rates adapt to the sample size $n$. In this way, Algorithms 3 and 4 perfectly resolve the main issue of their vanilla versions.

We can get stronger results in the specific case of distributed linear regression, where the contraction rate has nearly no dependence on the conditional number $\kappa$. Due to space constraints, we put all the details in Appendix B.

4.2 Multi-step estimators in large-sample regimes

We now present the contraction of optimization error of Algorithms 1 and 2.

Theorem 4.2**.**

Suppose that Assumptions 4.1 and 4.2 hold, and with probability tending to one, $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R)$ for some $R>\lVert\widehat{}\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\rVert_{2}$. For Algorithms 1 and 2, we have

[TABLE]

where $\eta=\kappa^{2}p(\log N)/n$. In addition, let Assumption 4.3 also hold. There exists some constant $C$ such that for Algorithm 2 we have

[TABLE]

where $t_{0}=\lceil\frac{2\log(CMR/\rho)}{\log(1/\eta)}\rceil$.

The strengthened result (4.1) requires Assumption 4.3 and thus smooth $g$. Theorem 4.2 shows that when $n$ is sufficiently large, Algorithms 1 and 2 behave similarly as Algorithms 3 and 4 – faster convergence with the larger $n$, mild restrictions on initialization, and averaging speeds up contraction given smooth loss. However, there is no convergence guarantee in general regimes. Section 5 further shows that with insufficient local sample size, the practical performance Algorithms 1 and 2 is less satisfactory.

Finally, to see how the averaging step reduces the statistical error (i.e. the distance between the estimator and $\boldsymbol{\theta}^{*}$), we continue to look at the linear regression example mentioned at the end of Section 2.2. For simplicity, assume $\{\mathbf{x}_{i}\}_{i\in[N]}$ are i.i.d. standard normal random vectors. (2.3) and (2.4) can be expressed as

[TABLE]

If further $n\gg p$, then the two contraction factors satisfy

[TABLE]

When the Algorithm 1 applies to this problem, the expression changes slightly to

[TABLE]

The smaller magnitude in (4.3) is due to the averaging.

Suppose that we initialize the algorithm using the one-shot average $\boldsymbol{\theta}_{0}=\frac{1}{m}\sum_{k=1}^{m}\widehat{}\boldsymbol{\theta}_{k}$, where $\widehat{}\boldsymbol{\theta}_{k}$ is the least squares solution on the $k$th machine. When $n=O(\sqrt{Np})$, we have $\sqrt{p/N}=O(p/n)$. Zhang et al. (2013) assert that $\|\widehat{}\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\|_{2}=O_{\mathbb{P}}(\sqrt{\frac{p}{N}})$, $\|\boldsymbol{\theta}_{0}-\widehat{}\boldsymbol{\theta}\|_{2}=O_{\mathbb{P}}(\frac{p}{n})$. Thus, the initial statistical error is $\|\boldsymbol{\theta}_{0}-\boldsymbol{\theta}^{*}\|_{2}=O_{\mathbb{P}}(\max\{\sqrt{\frac{p}{N}},\frac{p}{n}\})=O_{\mathbb{P}}(\frac{p}{n})$. By the contraction properties in (4.2) and (4.3), the optimization errors are

[TABLE]

We can see that when $p^{2}/n^{2}\ll\sqrt{p/N}\ll p^{3/2}/n^{3/2}$ (or equivalently, $p^{3/4}N^{1/4}\ll n\ll p^{3/2}N^{1/3}$), the optimization error is negligible for $\boldsymbol{\theta}_{1}$, but is not negligible for $\boldsymbol{\theta}_{1,k}$. When $n$ is smaller, even more iterations are needed. A refined analysis of distributed least squares is in Appendix B.

4.3 Guidance on practice

We now provide some general guidance on how to implement the CEASE algorithm in practice. First, we recommend choosing an initialization depending on the magnitude of the local sample size $n$. In particular, when $n$ is not very large compared to the dimension $p$, a zero initialization would be more robust. On the other hand, given a moderate or large $n$, the one-shot average estimator will lead to extremely fast convergence.

Second, according to Theorem 4.1, it suffices to take $\alpha$ of the order of $\rho\kappa^{2}p\log N/n$. In practice, setting $\alpha$ to be a small multiple of $p/n$ seems suitable in many occasions.

Finally, as is already shown, both Algorithm 3 and 4 reach statistical efficiency in $O(\frac{\log\lVert\boldsymbol{\theta}_{0}-\widehat{}\boldsymbol{\theta}\rVert_{2}+\log(N/p)}{\log(1/\eta)})$ iterations. In both our simulations and real data examples, the CEASE algorithms with a properly chosen $\alpha$ converge to the centralized estimator within 10 iterations. With a moderate $n$, a warm start further boosts the convergence speed.

5 Numerical experiments

5.1 Synthetic data

We first conduct distributed logistic regression to illustrate the effect of local sample size and initialization on convergence. We keep the total sample size $N=10000$ and the dimensionality $p=101$ fixed, and generate the i.i.d. data $\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{N}$ as follows: $\mathbf{x}_{i}=(1,\mathbf{u}_{i}^{\top})^{\top}$ with $\mathbf{u}_{i}\sim N(\mathbf{0}_{p-1},\boldsymbol{\Sigma})$ and $\boldsymbol{\Sigma}={\rm diag}(10,5,2,1\cdots 1)\in{\mathbb{R}}^{(p-1)\times(p-1)}$; ${\mathbb{P}}(y_{i}=1)=1-{\mathbb{P}}(y_{i}=0)=1/(1+e^{-\mathbf{x}_{i}^{\top}\boldsymbol{\theta}^{*}})$ where $\boldsymbol{\theta}^{*}\in{\mathbb{R}}^{p}$ is a random vector with norm 3 whose direction is chosen uniformly at random from the sphere. We use the natural logarithm of the estimation error $\|\boldsymbol{\theta}_{t}-\boldsymbol{\theta}^{*}\|_{2}$ to measure the performance of different algorithms, including multiple versions of the CEASE algorithms, GIANT (Wang et al., 2018), ADMM (Boyd et al., 2011) and accelerated gradient descent (Nesterov, 1983).

Figure 1 shows how the estimation errors evolve with iterations. The curves show the average values over 100 independent runs; the error bands correspond to one standard deviation. The regimes “large $n$”, “moderate $n$” and “small $n$” refer to $(n,m)=(2000,5)$, $(1000,10)$ and $(250,40)$; “zero initialization” and “good initialization” refer to $\boldsymbol{\theta}_{0}=\mathbf{0}$ (bottom panel) and $\bar{\boldsymbol{\theta}}$ (top panel), respectively. Here $\bar{\boldsymbol{\theta}}$ is the one-shot distributed estimator (Zhang et al., 2013) that averages the individual estimators on node machines. According to Figure 1, the standard deviation of each iterate is around 0.1. For the “large $n$, good initialization” regime, all of the iterates are unsurprisingly very close to the optimal solution. Their error bands will cover up the curves. So we omit the bands in that case for the sake of clarity.

With proper regularization, the two CEASE algorithms are the only ones that converge rapidly in all scenarios. The purely deterministic methods ADMM (Boyd et al., 2011) and accelerated gradient descent (Nesterov, 1983) are also reliable but slow. Other distributed algorithms like unregularized CEASE and GIANT (Wang et al., 2018) easily fail when the local sample size is small or the initialization is uninformative. In addition, the CEASE with averaging (Algorithm 4) is superior to the one without averaging (Algorithm 3). For example, when $(n,m)=(1000,10)$, the averaged CEASE with $\alpha=0$ converges while the one without averaging does not. Hence the averaging step leads to better performance.

We also test the efficacy of our algorithms in the distributed $\ell_{1}$-regularized logistic regression, where the penalty $g$ is nonsmooth (See Appendix D for details). To summarize, our simulations demonstrate several important properties of the CEASE Algorithms:

• •

In all scenarios, the CEASE Algorithms converge rapidly, usually within several steps, which is consistent with our theory;

• •

The CEASE Algorithms efficiently utilize statistical structures and similarities among local losses, and benefit from the averaging step with smooth loss functions;

• •

The CEASE Algorithms are also able to handle the most general situations (e.g. small local sample size, uninformative initialization) with convergence guarantees.

5.2 Real data

As a real data example, we choose the Fashion-MNIST dataset (Xiao et al., 2017) as a testbed for comparison of algorithms. The whole dataset consists of 70000 grayscale images of fashion products in 10 classes, each of which has 6000 training samples and 1000 testing samples. We choose the 7th and 9th classes (Sneakers and Ankle boots) and the goal is to train a classifier that distinguishes them. Each image has $28\times 28=784$ pixels, represented by a feature vector in $[0,1]^{784}$. The number of training (or testing) samples is $6000\times 2=12000$ (or $1000\times 2=2000$). We randomly partition the training set and conduct logistic regression in a distributed manner. The performance metric is the classification error on the testing set. Figure 2 shows the average performance of the CEASE algorithms, ADMM, GIANT and AGD based on 100 independent runs, together with error bars showing one standard deviation. Here “large $n$”, “moderate $n$” and “small $n$” refer to $(n,m)=(1200,10)$, $(480,25)$ and $(240,50)$, respectively. All of the iterations are initialized with the one-shot average (Zhang et al., 2013). The experiments on this real data example also support our theoretical findings.

6 Discussions

We have developed two CEASE distributed estimators (Algorithms 3 and 4) for statistical estimation, with theoretical guarantees and superior performance on real data. Several new directions are worth exploring. First, while we assumed exact computation for simplicity, finer analysis should allow for inexact updates in practice. Second, we hope to extend the algorithms to decentralized and asynchronous settings. Third, distributed versions of confidence regions and hypothesis tests are of great importance, and our point estimation strategies may serve as a starting point. Finally, it will be interesting to explore non-convex statistical optimization problems such as mixture models and deep learning. We believe that the idea of gradient-enhanced loss function still plays an important role.

Acknowledgement

We gratefully acknowledge NSF grants DMS-1662139 and DMS-1712591, NIH grant 2R01-GM072611-14, and ONR grant N00014-19-1-2120. We acknowledge computing resources from Columbia University’s Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement Grant 1G20-RR030893-01, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR) Contract C090171, both awarded April 15, 2010.

Appendices

Section A introduces a variant of Algorithm 4. Section B outlines the deferred results for distributed linear regression. Section C shows the CEASE iterates for distributed logistic regression. Section D presents the numerical results on the distributed $\ell_{1}$-regularized logistic regression. Section E presents the proofs of the main results. Section F lists a few technical lemmas that are used throughout the proofs.

Appendix A A variant of Algorithm 4

As mentioned in the main text, Algorithm 2 is unstable when the local sample size $n$ is not sufficiently large. The proximal gradient method Algorithm 5 was introduced to stabilize the solution path by shrinking towards the solution in the previous step. A variant of CEASE that stablizes Algorithm 2 is to take smaller step-sizes, which we now present. The idea is applicable to stabilize Algorithm 1 too, resulting a variant to Algorithm 3.

The only difference between Algorithms 2 and 5 lies in the aggregation step. From $\boldsymbol{\theta}_{t}$ the former directly jumps to the average of new individual estimators $\{\boldsymbol{\theta}_{t,k}\}_{k=1}^{m}$, while the latter proceeds more cautiously in that direction. Algorithm 2 is a special case of Algorithm 5 with $\alpha_{t}=1$ for all $t$. Choosing $\alpha_{t}\in(0,1)$ helps stablize the iterates especially when $n$ is small.

Algorithm 5 is conceptually simple and easy to implement. To see its performance, we conduct distributed logistic regression on the first set of synthetic data in Section 5.1. The numerical results there show that Algorithm 2 fails to converge when $(n,m)=(250,40)$. We run Algorithm 5 with constant step-size $\alpha_{0}=\alpha_{1}=\cdots$ under this setting, where $\alpha_{0}\in\{1,1/2,1/4,1/8\}$. We also run CEASE for comparison. Figure 3 summarizes all the results. Again, the curves show the average values over 100 independent runs; the error bands correspond to one standard deviation. In this experiment, the performance of Algorithm 5 with $\alpha_{t}=1/4$ is similar to that of CEASE; $\alpha_{t}=1/2$ leads to even faster convergence; and $\alpha_{t}=1/8$ slows it down. For simplicity, we take the step-size to be constant over time. It would be interesting to explore decaying schemes such as $\alpha_{t}\asymp t^{-\beta}$ for some $\beta>0$.

The new algorithm also has some drawbacks. For instance, it implicitly assumes that the local sample size is large enough so that each local machine can solve its optimization problems reliably. Consider the example of the distributed logistic regression with $N$ samples and $p$ variables in total. When $m$ is large so that $n=N/m<p$, the loss functions $\{f_{k}\}_{k=1}^{m}$ on local machines are no longer strongly convex and thus do not have unique minima. Then Algorithm 5 needs to be modified because

[TABLE]

is not uniquely defined. We need some pivoting rule to choose one optimal solution, such as the one with the minimum Euclidean norm. That complicates the algorithm and it is not clear how to establish theoretical guarantees then. In contrast, the quadratic proximity term in CEASE always make the objective function strongly convex. It ensures the uniqueness of $\{\boldsymbol{\theta}_{t,k}\}_{k=1}^{m}$ and facilitates computation.

Appendix B Distributed linear regression

In distributed linear regression, recall that the $k^{th}$ machine defines a quadratic loss function

[TABLE]

where $\widehat{}\boldsymbol{\Sigma}_{k}=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}$ and $\widehat{}\mathbf{w}_{k}=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}\mathbf{x}_{i}y_{i}$. Let $f(\boldsymbol{\theta})=\frac{1}{2}\boldsymbol{\theta}^{\top}\widehat{\boldsymbol{\Sigma}}\boldsymbol{\theta}-\widehat{\mathbf{w}}^{\top}\boldsymbol{\theta}$. Without loss of geneality we write $\mathbf{x}_{i}=(1,\mathbf{u}_{i}^{\top})^{\top}\in{\mathbb{R}}^{p}$.

Assumption B.1**.**

• •

${\mathbb{E}}\mathbf{u}_{i}=0$* and ${\mathbb{E}}(\mathbf{u}_{i}\mathbf{u}_{i}^{\top})=\boldsymbol{\Sigma}\succ 0$. $\{\boldsymbol{\Sigma}^{-1/2}\mathbf{u}_{i}\}_{i=1}^{N}$ are i.i.d. sub-Gaussian random vectors with bounded $\|\boldsymbol{\Sigma}^{-1/2}\mathbf{u}_{i}\|_{\psi_{2}}$.*

• •

The minimum eigenvalue $\lambda_{\min}(\boldsymbol{\Sigma})$ is bounded away from zero.

• •

$N/{\rm Tr}(\boldsymbol{\Sigma})\geq C>0$* and $n/\log m\geq c>0$ where $C$ and $c$ are constants.*

For the least-squares, Algorithm 3 admits a close-form:

[TABLE]

and so does Algorithm 4:

[TABLE]

Intuitively, the averaging step in Algorithm 4 reduces variance and accelerates convergence. Below we study Algorithm 4 with the help of these analytical expressions. In the large sample regime, we achieve a contraction factor of $O(p/n)$ without any condition number; in the general regime, linear convergence is still guaranteed.

Theorem B.1**.**

Suppose Assumption B.1 holds and $n/p$ is bounded away from zero. Then, there exist positive constants $C_{1},C_{2}$ and $C_{3}$ such that when (i) $n\geq C_{1}p$ and $\alpha\geq 0$ or (ii) $\alpha\geq C_{1}{\rm Tr}(\boldsymbol{\Sigma})/n$, with probability tending to 1,

[TABLE]

where $\kappa=\lambda_{\max}(\boldsymbol{\Sigma})/\lambda_{\min}(\boldsymbol{\Sigma})$ and $\eta=1-\frac{1-\min\{1/2,C_{2}p/n\}}{1+C_{3}\alpha}$.

Theorem B.1 reveals the following remarkable facts about Algorithm 4: No matter what relationship $n$ and $p$ have, proper regularization always guarantees linear convergence, and the rate exhibits a smooth transition as $p/n$ grows. Hence we can handle the distributed statistical estimation problem without assuming large enough $n$, overcoming the difficulty of other algorithms in literature (Zhang et al., 2013; Battey et al., 2018; Jordan et al., 2019).

If $n/p$ is large enough, the regularization is not necessary, but choosing $\alpha\asymp p/n$ does not hurt much. This is because we can control the contraction factor as:

[TABLE]

When $n/p$ is not that large, most distributed statistical estimation procedures fail. By choosing $\alpha=\tilde{C}{\rm Tr}(\boldsymbol{\Sigma})/n$ for $\tilde{C}>C_{1}$ (see Condition (ii) of Theorem B.1) we still have linear convergence with contraction factor at most

[TABLE]

In most situations of interest we have ${\rm Tr}(\boldsymbol{\Sigma})\asymp p$ (even for pervasive factor models). Therefore we see that $\alpha\asymp p/n$ is a universal and adaptive choice of regularization over all the possible relation between $n$ and $p$.

Another benefit of the Algorithms is that the condition number $\kappa$ has only logarithmic effect on the iteration complexity, and the contraction factor in Theorem B.1 does not depend on $\kappa$ at all. This is in stark contrast to the analysis under the same setting in Shamir et al. (2014), and helps relax the commonly used boundedness assumption on the condition number in Zhang et al. (2013), Battey et al. (2018), Jordan et al. (2019), among others. It is worth mentioning that Wang et al. (2018) derive similar results for distributed linear regression when the local sample size $n$ is sufficiently large.

Appendix C Distributed logistic regression

In this section, we demonstrate the iterates of CEASE in distributed logistic regression. For illustration purposes, we use Newton’s method on each local machine to solve the optimization problem. It is worth pointing out that Newton’s method is not the only choice. First-order methods such gradient descent can also efficiently do the job.

Given samples $\{(\mathbf{x}_{i},y_{i})\}_{i\in[N]}$, the loss on the $k$th machine is

[TABLE]

Thus for Algorithm 4, we have for each iteration

[TABLE]

Here $\widehat{}\mathbf{w}_{k}=\frac{1}{m}\sum_{i\in\mathcal{I}_{k}}y_{i}\mathbf{x}_{i}$, $\boldsymbol{\eta}_{t,k}=\nabla f_{k}(\boldsymbol{\theta}_{t})-\nabla f(\boldsymbol{\theta}_{t})$. Assume that we start with $\boldsymbol{\theta}_{t,k}^{(0)}=\boldsymbol{\theta}_{t}$. Then one Newton iteration of the above minimization problem is

[TABLE]

Here, $\mathbf{X}_{k}\in\mathbb{R}^{n\times p}$ is a matrix with rows consisting of $\{\mathbf{x}_{i}\}_{i\in\mathcal{I}_{k}}$, $\mathbf{Y}_{k}\in\mathbb{R}^{n}$ denotes the vector of $\{y_{i}\}_{i\in\mathcal{I}_{k}}$, $\mathbf{p}_{k}^{(j)}\in\mathbb{R}^{n}$ is the vector of fitted probabilities with entry $i$ equal to $e^{(\boldsymbol{\theta}_{t,k}^{(j)})^{\top}\mathbf{x}_{i}}/(1+e^{(\boldsymbol{\theta}_{t,k}^{(j)})^{\top}\mathbf{x}_{i}})$, and $\mathbf{W}_{k}^{(j)}\in\mathbb{R}^{n\times n}$ is a diagonal matrix with entries $\{e^{(\boldsymbol{\theta}_{t,k}^{(j)})^{\top}\mathbf{x}_{i}}/(1+e^{(\boldsymbol{\theta}_{t,k}^{(j)})^{\top}\mathbf{x}_{i}})^{2}\}_{i\in\mathcal{I}_{k}}$. Running one-step or multi-step Newton’s iteration will result in some $\boldsymbol{\theta}_{t,k}^{(j)}$ sufficiently close to $\boldsymbol{\theta}_{t,k}$. After that, we simply average them across all the machines and obtain $\boldsymbol{\theta}_{t+1}$.

For Algorithm 2, in each iteration, we simply set $\alpha=0$ in the above procedure. For Algorithm 3 and Algorithm 1, we simply omit the averaging step in Algorithm 4 and Algorithm 2 respectively.

Appendix D Numerical results on distributed $\ell_{1}$-regularized logistic regression

In this section, we present the performance of our algorithms in the distributed $\ell_{1}$-regularized logistic regression problem where nonsmooth penalty is present. We fix the total sample size $N=5000$ and the dimensionality $p=1001$, and generate the i.i.d. data $\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{N}$ as follows: $\mathbf{x}_{i}=(1,\mathbf{u}_{i}^{\top})^{\top}$ with $\mathbf{u}_{i}\sim N(\mathbf{0}_{p-1},\mathbf{I}_{p-1})$; ${\mathbb{P}}(y_{i}=1)=1-{\mathbb{P}}(y_{i}=0)=\frac{e^{\mathbf{x}_{i}^{\top}\boldsymbol{\theta}^{*}}}{1+e^{\mathbf{x}_{i}^{\top}\boldsymbol{\theta}^{*}}}$ where $\boldsymbol{\theta}^{*}=(\mathbf{1}_{10}^{\top},\mathbf{0}_{991}^{\top})^{\top}/\sqrt{2}\in{\mathbb{R}}^{p}$. We define the penalty function $g(\boldsymbol{\theta})=\lambda\|\boldsymbol{\theta}\|_{1}$ with $\lambda=0.5\sqrt{\frac{\log p}{N}}$, such that the regularized MLE over the whole dataset recovers the nonzeros of $\boldsymbol{\theta}^{*}$ accurately. Figure 4 shows the performance of CEASE algorithms and ADMM, where “large $n$”, “moderate $n$” and “small $n$” refer to $(n,m)=(1000,5)$, $(500,10)$ and $(250,20)$, and “zero initialization” and “good initialization” refer to $\boldsymbol{\theta}_{0}=\mathbf{0}$ and $\bar{\boldsymbol{\theta}}$, respectivley. Again, $\bar{\boldsymbol{\theta}}$ is the one-shot distributed estimator (Zhang et al., 2013). All the results are average values of 100 independent runs.

Similar to the distributed logistic regression case, the CEASE algorithms with proper regularization (Algorithms 3 and 4) work well in general; without regularization, the CEASE algorithm fails to converge when the local sample size $n$ is small and the initialization in uninformative. For this nonsmooth problem, the CEASE algorithm with averaging (Algorithm 4) does not seem to have advantage over the single version (Algorithm 3). The ADMM converges quickly to a region near the minimizer but then proceeds quite slowly, which appears to be a common phenomenon in many distributed optimization problems (Boyd et al., 2011).

Appendix E Proofs

E.1 Proof of Theorem 3.1

Lemma E.1**.**

Let Assumptions 3.1 and 3.2 hold. Consider the iterates $\{\boldsymbol{\theta}_{t}\}_{t=0}^{T}$ generated by Algorithm 4. Define

[TABLE]

If $0<\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2}<R/2$, $\delta<\rho_{0}+\alpha$ and $\gamma_{t}^{2}<\rho/(\rho+2\alpha)$, then

[TABLE]

Theorem 3.1 directly follows from Lemma E.1 and induction. Below we only prove Lemma E.1 with Assumption 3.3. The other part in Lemma E.1 without Assumption 3.3 can be derived by slightly modifying this proof.

Proof of Lemma E.1 with Assumption 3.3.

Let $\boldsymbol{\theta}_{t}^{+}=\mathrm{prox}_{\alpha^{-1}(f+g)}(\boldsymbol{\theta}_{t})$. By the triangle inequality,

[TABLE]

We first invoke Theorem 3.2 to bound the first term $\|\boldsymbol{\theta}_{t+1}-\boldsymbol{\theta}_{t}^{+}\|_{2}$ in (E.1). Define $\tilde{g}(\boldsymbol{\theta})=g(\boldsymbol{\theta})+(\alpha/2)\|\boldsymbol{\theta}-\boldsymbol{\theta}_{t}\|_{2}^{2}$ for $\boldsymbol{\theta}\in{\mathbb{R}}^{p}$. Then $\boldsymbol{\theta}_{t+1}$ is the first iterate of Algorithm 2 initialized at $\boldsymbol{\theta}_{t}$ for computing $\boldsymbol{\theta}^{+}_{t}={\rm argmin}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}\{\frac{1}{m}\sum_{k=1}^{m}f_{k}(\boldsymbol{\theta})+\tilde{g}(\boldsymbol{\theta})\}$.

From $\widehat{}\boldsymbol{\theta}=\mathrm{prox}_{\alpha^{-1}(f+g)}(\widehat{}\boldsymbol{\theta})$ and Lemma F.3 we obtain that $\|\boldsymbol{\theta}_{t}^{+}-\widehat{}\boldsymbol{\theta}\|_{2}\leq\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2}$. Then the condition $\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2}<R/2$ leads to $B(\boldsymbol{\theta}_{t}^{+},R/2)\subseteq B(\widehat{}\boldsymbol{\theta},R)$. By Assumptions 3.1 and 3.2,

• •

in $B(\boldsymbol{\theta}^{+}_{t},R/2)$, $\{f_{k}+\tilde{g}\}_{k=1}^{m}$ are $(\rho_{0}+\alpha)$-strongly convex and $(f+\tilde{g})$ is $(\rho+\alpha)$-strongly convex;

• •

$\|\nabla^{2}f_{k}(\boldsymbol{\theta})-\nabla^{2}f(\boldsymbol{\theta})\|_{2}\leq\delta$ holds for all $k\in[m]$ and $\boldsymbol{\theta}\in B(\boldsymbol{\theta}^{+}_{t},R/2)$.

Furthermore, Lemma F.3 also yields

[TABLE]

Then $\|\boldsymbol{\theta}_{t}-\boldsymbol{\theta}_{t}^{+}\|_{2}\leq\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2}<R/2$ and $\tilde{\boldsymbol{\theta}}_{0}\in B(\boldsymbol{\theta}^{+}_{t},R/2)$. Based on these conditions and $\alpha+\rho_{0}>\delta$, we use Theorem 3.2 to get

[TABLE]

From here, (E.1) and (E.2) we obtain that

[TABLE]

where $h(x)=\gamma_{t}\sqrt{1-x^{2}}+x$, $\forall x\in[0,1]$. From $h^{\prime}(x)=1-\gamma_{t}x/\sqrt{1-x^{2}}$ we see that $h^{\prime}\geq 0$ on $[0,1/\sqrt{1+\gamma_{t}^{2}}]$.

On the one hand, Lemma F.3 asserts that $\|\boldsymbol{\theta}^{+}-\widehat{}\boldsymbol{\theta}\|_{2}/\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2}\leq\alpha/(\rho+\alpha)$. On the other hand, the assumption $\gamma_{t}^{2}<\rho/(\rho+2\alpha)$ forces

[TABLE]

The proof is completed by computation:

[TABLE]

where we used the assumption $\gamma_{t}^{2}<\rho/(\rho+2\alpha)$ again. ∎

E.2 Proof of Corollary 3.1

We claim that $(\frac{\delta}{\rho_{0}+\alpha})^{2}\leq\frac{7}{9}\cdot\frac{\rho}{\rho+2\alpha}$. Given this, Corollary 3.1 follows from Theorem 3.1 and

[TABLE]

The claim trivially holds if $\delta=0$. When $\delta>0$, let us first assume $0<\delta\leq\rho$ and define $b=\alpha\rho/\delta^{2}$. Then $\alpha=b\delta^{2}/\rho$, $b\geq 4$ and $\rho_{0}\geq\max\{\rho-\delta,0\}$ force

[TABLE]

and $[\delta/(\rho_{0}+\alpha)]^{2}\leq[\delta/(\rho_{0}+\alpha)]^{2}\leq 1/h_{1}(\delta/\rho)$, where $h_{1}(x)=(bx+x^{-1}-1)^{2}$. On the other hand, $\rho/(\rho+2\alpha)=1/(1+2\alpha/\rho)=1/h_{2}(\delta/\rho)$, where $h_{2}(x)=1+2bx^{2}$.

We are going to show $h_{2}(x)\leq 7h_{1}(x)/9$, $\forall x\in(0,1]$, which leads to the desired result under $0<\delta\leq\rho$. If $0<x\leq\sqrt{3}/2$, then $h_{1}(x)\geq(2\sqrt{b}-1)^{2}\geq(2\sqrt{b}-\sqrt{b}/2)^{2}\geq 9b/4$ and

[TABLE]

If $\sqrt{3}/2<x\leq 1$, then $h_{1}(x)\geq b^{2}x^{2}\geq 3b^{2}/4$, $h_{2}(x)\leq 1+2b\leq(b/4)+2b=9b/4$, and $h_{2}(x)/h_{1}(x)=3/b\leq 3/4\leq 7/9$.

Suppose now that $\delta>\rho$, and define $b=\alpha\rho/\delta^{2}$. Then

[TABLE]

From $b\geq 4$ and $\delta/\rho>1$ we get $(\frac{\delta}{\rho_{0}+\alpha})^{2}\leq\frac{7}{9}\cdot\frac{\rho}{\rho+2\alpha}$ from

[TABLE]

E.3 Proof of Corollary 3.2

Throughout the proof we assume that $\delta/\rho$ is sufficiently small. The regularity conditions in Theorem 3.1 are easily verified as $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R/2)$ and $[\delta/(\rho_{0}+\alpha)]^{2}<\rho/(\rho+2\alpha)$. Here we used the fact $\rho_{0}\geq\rho-\delta$.

From $\rho_{0}\geq\rho-\delta$ we get $\rho_{0}+\alpha\geq\rho_{0}\geq\rho/2$ and $\delta/(\rho_{0}+\alpha)\leq 2\delta/\rho$. Also, $\sqrt{\rho^{2}+2\alpha\rho}\leq\rho\sqrt{1+2C(\delta/\rho)^{2}}\lesssim\rho$. We control the contraction factor in (3.2):

[TABLE]

Recall that $\gamma_{t}=\frac{\delta}{\rho_{0}+\alpha}\cdot\min\{1,\frac{\delta}{\rho+\alpha}(1+\frac{M}{\rho_{0}+\alpha}\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2})\}$ in Theorem 3.1. When $\delta/\rho$ is small and $\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2}\leq\rho/M$, we have $\rho_{0}+\alpha\geq\rho/2$, $\frac{M}{\rho_{0}+\alpha}\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2}\leq 2$, and $\gamma_{t}\leq(2\delta/\rho)^{2}$. This help bound the contraction factor in (3.2):

[TABLE]

E.4 Proof of Theorem 3.2

Theorem 3.2 is a direct summary of the following two lemmas.

Lemma E.2** (Contraction).**

Let Assumptions 3.1 and 3.2 hold, with $\rho_{0}>\delta\geq 0$. Then $\|\varphi_{k}(\boldsymbol{\theta})-\widehat{}\boldsymbol{\theta}\|_{2}\leq(\delta/\rho_{0})\|\boldsymbol{\theta}-\widehat{}\boldsymbol{\theta}\|_{2}$, $\forall\boldsymbol{\theta}\in B(\widehat{}\boldsymbol{\theta},R)$, $\forall k\in[m]$.

Proof.

Fix $\boldsymbol{\theta}\in B(\widehat{}\boldsymbol{\theta},R)$. By the first order condition of $\varphi_{k}(\boldsymbol{\theta})$, we have that

[TABLE]

Using the fixed point property $\varphi_{k}(\widehat{}\boldsymbol{\theta})=\widehat{}\boldsymbol{\theta}$, we have $\nabla f_{k}(\widehat{}\boldsymbol{\theta})-\nabla f(\widehat{}\boldsymbol{\theta})\in\partial[f_{k}(\widehat{}\boldsymbol{\theta})+g(\widehat{}\boldsymbol{\theta})]$. By the Taylor expansion and Assumption 3.2,

[TABLE]

From this, (E.3) and Lemma F.2, we obtain that $\|\varphi_{k}(\boldsymbol{\theta})-\widehat{}\boldsymbol{\theta}\|_{2}\leq(\delta/\rho_{0})\|\boldsymbol{\theta}-\widehat{}\boldsymbol{\theta}\|_{2}$.

∎

Lemma E.3** (Averaging).**

Let Assumptions 3.1, 3.2 and 3.3 hold, with $\rho_{0}>\delta\geq 0$. We have

[TABLE]

Proof.

Define $L_{k}(\boldsymbol{\theta})=f_{k}(\boldsymbol{\theta})+g(\boldsymbol{\theta})$ and $L(\boldsymbol{\theta})=f(\boldsymbol{\theta})+g(\boldsymbol{\theta})$ for $\boldsymbol{\theta}\in{\mathbb{R}}^{p}$. Then $\widehat{}\boldsymbol{\theta}={\rm argmin}_{\boldsymbol{\xi}\in{\mathbb{R}}^{p}}L(\boldsymbol{\xi})$ and $\varphi_{k}(\boldsymbol{\theta})={\rm argmin}_{\boldsymbol{\xi}\in{\mathbb{R}}^{p}}\{L_{k}(\boldsymbol{\xi})-\langle\nabla L_{k}(\boldsymbol{\theta})-\nabla L(\boldsymbol{\theta}),\boldsymbol{\xi}\rangle\}$. By the optimality conditions,

[TABLE]

After subtracting $\nabla L_{k}(\widehat{}\boldsymbol{\theta})$ from both sides and rearranging terms, we get

[TABLE]

Note that the average of the right hand side over $k\in[m]$ is $\mathbf{0}$.

Define $\mathbf{H}_{k}=\int_{0}^{1}\nabla^{2}L_{k}[(1-t)\widehat{\boldsymbol{\theta}}+t\varphi_{k}(\boldsymbol{\theta})]\mathrm{d}t$ for $k\in[m]$ and $\widehat{}\mathbf{H}=\nabla^{2}L(\widehat{\boldsymbol{\theta}})$. Then

[TABLE]

where we let $\bar{\varphi}(\boldsymbol{\theta})=\frac{1}{m}\sum_{k=1}^{m}\varphi_{k}(\boldsymbol{\theta})$. As a result,

[TABLE]

Lemma E.2 forces $\max_{k\in[m]}\|\varphi_{k}(\boldsymbol{\theta})-\widehat{}\boldsymbol{\theta}\|_{2}\leq(\delta/\rho_{0})\|\boldsymbol{\theta}-\widehat{}\boldsymbol{\theta}\|_{2}$, and Assumption 3.1 yields $\widehat{}\mathbf{H}\succeq\rho\mathbf{I}$ and $\|\widehat{}\mathbf{H}\|_{2}\leq 1/\rho$. Furthermore, we use Assumptions 3.2 and 3.3 to get

[TABLE]

The proof is finished by combining all the estimates above. ∎

E.5 Proof of Theorem 4.1

The proof is implied by combining proof of Corollary 3.2 with the results of the following two lemmas, the first of which is a direct counterpart of Theorem 3.1 in the stochastic setting, given an additional condition on similarity between local Hessians. The second lemma below specifies the order of Hessian difference in the generalized linear model, hence providing a contraction rate and guiding the choice of $\alpha$.

Lemma E.4**.**

Let Assumption 4.2 hold. Denote

[TABLE]

Consider the iterates $\{\boldsymbol{\theta}_{t}\}_{t=0}^{T}$ generated by Algorithm 3 or Algorithm 4. Suppose that $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R/2)$ and $[\widehat{\delta}/(\rho_{0}+\alpha)]^{2}<\rho/(\rho+2\alpha)$.

• •

For both Algorithms 3 and 4, we have

[TABLE]

• •

If in addition, Assumption 4.3 also holds, then for Algorithm 4 we have

[TABLE]

where we define $\gamma_{t}=\frac{\widehat{\delta}}{\rho_{0}+\alpha}\cdot\min\{1,\frac{\widehat{\delta}}{\rho+\alpha}(1+\frac{M}{\rho_{0}+\alpha}\|\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\theta}\|_{2})\}$;

• •

Both multiplicative factors in (E.4) and (E.5) are strictly less than 1.

Proof of Lemma E.4.

We first assume that $\rho_{0}>\widehat{\delta}$ and analyze the vanilla DANE algorithm under the new assumptions. Let Assumption 4.2 hold, and $\{\boldsymbol{\theta}_{t}\}_{t=0}^{\infty}$ be the iterates with $\boldsymbol{\theta}_{0}\in B(\widehat{}\boldsymbol{\theta},R)$. For any $\boldsymbol{\theta}\in B(\widehat{}\boldsymbol{\theta},R)$, $\|\nabla^{2}f(\boldsymbol{\theta})-\nabla^{2}F(\boldsymbol{\theta})\|_{2}\leq\widehat{\delta}/2$ and thus $\|\nabla^{2}f_{k}(\boldsymbol{\theta})-\nabla^{2}f(\boldsymbol{\theta})\|_{2}\leq\widehat{\delta}$ for $k\in[m]$. Hence it implies Assumption 3.2 with $\delta=\widehat{\delta}$, and Lemma E.2 continues to hold. We can also get the result in Lemma E.3 under Assumption 4.3, by replacing $\widehat{}\mathbf{H}=\nabla^{2}f(\widehat{}\boldsymbol{\theta})$ in the proof of Lemma E.3 by $\nabla^{2}F(\widehat{}\boldsymbol{\theta})+\nabla^{2}g(\boldsymbol{\theta})$. Then we drop the assumption $\rho_{0}>\widehat{\delta}$ can reproduce the results in Theorem 3.1 under the new setting, by following its original proof. ∎

Lemma E.5**.**

Under Assumption 4.1, for an arbitrarily small positive constant $c$, there exist universal positive constants $C_{1},C_{2}$ and $C_{3}$ depending only on $c$ such that as long as $n\geq cp$, with probability at least $1-2e^{-C_{2}n}-Ne^{-C_{3}p}$,

[TABLE]

Proof of Lemma E.5.

Let $\tilde{\mathbf{x}}_{i}=(\boldsymbol{\Sigma}^{*})^{-1/2}\mathbf{x}_{i}$, $\tilde{}\boldsymbol{\theta}=(\boldsymbol{\Sigma}^{*})^{1/2}\boldsymbol{\theta}$, and define a new loss function $\tilde{l}(\tilde{}\boldsymbol{\theta},\tilde{}\mathbf{x}_{i})=b(\tilde{}\mathbf{x}_{i}^{\top}\tilde{}\boldsymbol{\theta})-y_{i}(\tilde{}\mathbf{x}_{i}^{\top}\tilde{}\boldsymbol{\theta}).$ Let $\widehat{R}_{k}(\tilde{}\boldsymbol{\theta})=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}\tilde{l}(\tilde{}\boldsymbol{\theta},\tilde{}\mathbf{x}_{i})$ for $k\in[m]$. Then we have that $\nabla^{2}\widehat{R}_{k}(\tilde{}\boldsymbol{\theta})=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}b^{\prime\prime}(\tilde{}\boldsymbol{\theta}^{\top}\tilde{}\mathbf{x}_{i})\tilde{}\mathbf{x}_{i}\tilde{}\mathbf{x}_{i}^{\top}$ and $\nabla^{2}f_{k}(\boldsymbol{\theta})=(\boldsymbol{\Sigma}^{*})^{1/2}\nabla^{2}\widehat{R}_{k}(\tilde{}\boldsymbol{\theta})(\boldsymbol{\Sigma}^{*})^{1/2}$. Similarly we have $\nabla^{2}F(\boldsymbol{\theta})=(\boldsymbol{\Sigma}^{*})^{1/2}{\mathbb{E}}\nabla^{2}\widehat{R}_{k}(\tilde{}\boldsymbol{\theta})(\boldsymbol{\Sigma}^{*})^{1/2}$.

Therefore

[TABLE]

and we only need to control the quantity on the right hand side. Here $\tilde{R}=\lVert\boldsymbol{\Sigma}^{*}\rVert_{2}^{1/2}R$.

Define $\boldsymbol{\Delta}_{0}=\max_{k\in[m]}\max_{\tilde{}\boldsymbol{\theta}\in B((\boldsymbol{\Sigma}^{*})^{1/2}\widehat{}\boldsymbol{\theta},\tilde{R})}\|\nabla^{2}\widehat{R}_{k}(\tilde{}\boldsymbol{\theta})-\nabla^{2}{\mathbb{E}}\widehat{R}_{k}(\tilde{}\boldsymbol{\theta})\|_{2}$ and

[TABLE]

Here $\tilde{}\mathbf{X}$ shares the distribution with $\tilde{}\mathbf{x}_{i}$. Firstly we bound $\phi_{k}({\tilde{}\boldsymbol{\theta}})$ for any fixed ${\tilde{}\boldsymbol{\theta}}\in B({\tilde{}\boldsymbol{\theta}}^{*},2\tilde{R})$, where ${\tilde{}\boldsymbol{\theta}}^{*}:=(\boldsymbol{\Sigma}^{*})^{1/2}\widehat{}\boldsymbol{\theta}$. For any $k\in[m]$, under Assumption 4.1, there exist constants $c_{1},c_{2}$ such that for any $\epsilon\geq 0$

[TABLE]

To see this, notice that $\phi_{k}({\tilde{}\boldsymbol{\theta}})=\sup_{\mathbf{u}\in{\mathbb{S}}^{p-1}}g_{k}(\mathbf{u})$, where $g_{k}(\mathbf{u})=\mathbf{u}^{\top}\{\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}b^{\prime\prime}({\tilde{}\mathbf{x}}_{i}^{\top}{\tilde{}\boldsymbol{\theta}}){\tilde{}\mathbf{x}}_{i}{\tilde{}\mathbf{x}}_{i}^{\top}-{\mathbb{E}}b^{\prime\prime}({\tilde{}\mathbf{X}}^{\top}{\tilde{}\boldsymbol{\theta}}){\tilde{}\mathbf{X}}{\tilde{}\mathbf{X}}^{\top}\}\mathbf{u}.$ Let $\mathcal{N}$ be a $\frac{1}{4}$-covering of ${\mathbb{S}}^{p-1}$, and $|\mathcal{N}|\leq 9^{p}$. Denote $\widehat{}\mathbf{u}=\arg\max_{\mathbf{u}}g_{k}(\mathbf{u})$. Find $\tilde{}\mathbf{u}\in\mathcal{N}$ such that $\lVert\widehat{}\mathbf{u}-\tilde{}\mathbf{u}\rVert_{2}\leq\frac{1}{4}$. Then

[TABLE]

and thus

[TABLE]

On the other hand from Bernstein’s inequality we see that there exists a constant $c_{2}$ such that for any $\mathbf{u}\in\mathcal{N}$, $\epsilon\geq 0$, ${\mathbb{P}}(g_{k}(\mathbf{u})\geq\frac{\epsilon}{2})\leq 2e^{-c_{2}\min\{\epsilon,\epsilon^{2}\}n}$. Therefore

[TABLE]

where $c_{1}=\log 9$.

Now for $t\geq 1$, define the event $E_{t}\triangleq\left\{\max_{i=1}^{N}\lVert{\tilde{}\mathbf{x}}_{i}\rVert_{2}^{3}<(8t)^{3/2}{\mathbb{E}}\lVert{\tilde{}\mathbf{X}}\rVert_{2}^{3}\right\}.$ Then by Theorem 2.1 in Hsu et al. (2012), ${\mathbb{P}}(E_{t}^{c})={\mathbb{P}}(\max_{i=1}^{N}\lVert{\tilde{}\mathbf{x}}_{i}\rVert_{2}^{2}\geq 8t[{\mathbb{E}}\lVert{\tilde{}\mathbf{X}}\rVert_{2}^{3}]^{2/3})\leq{\mathbb{P}}(\max_{i=1}^{N}\lVert{\tilde{}\mathbf{x}}_{i}\rVert_{2}^{2}\geq 8t{\mathbb{E}}\lVert{\tilde{}\mathbf{X}}\rVert_{2}^{2})\leq N{\mathbb{P}}(\lVert{\tilde{}\mathbf{X}}\rVert_{2}^{2}\geq 8t{\mathbb{E}}\lVert{\tilde{}\mathbf{X}}\rVert_{2}^{2})\leq Ne^{-tp}$. Under the event $E_{t}$, for ${\tilde{}\boldsymbol{\theta}}_{1},{\tilde{}\boldsymbol{\theta}}_{2}\in B({\tilde{}\boldsymbol{\theta}}^{*},2\tilde{R})$, we have

[TABLE]

for some constant $c_{3}$ depending only on $B_{3}$.

Now let $\mathcal{N}_{\delta}$ be a $\delta$-covering of $B({\tilde{}\boldsymbol{\theta}}^{*},2\tilde{R})$, where $\delta=\frac{\epsilon}{c_{3}(tp)^{3/2}}$. We can also assume that $|{\mathcal{N}}_{\delta}|\leq(\frac{6\tilde{R}}{\delta})^{p}$. Therefore for any $k\in[m]$,

[TABLE]

Thus

[TABLE]

It is easily seen that the last expression is no more than $Ne^{-tp}+2e^{-C_{2}p}$ when

[TABLE]

which is satisfied if $t$ is chosen as a suitable constant depending on $c$, and that

[TABLE]

Here $C_{i}^{\prime}$ and $C$ is a constant depending only on $c$.

Finally, note that

[TABLE]

and that $\lVert\boldsymbol{\Sigma}\rVert_{2}\geq A_{1}$ for a universal $A_{1}>0$, we have $\lVert\boldsymbol{\Sigma}^{*}\rVert_{2}\leq\max\{1,1/A_{1}\}\lVert\boldsymbol{\Sigma}\rVert_{2}$. Thus combining (E.8) and (E.6) completes the proof. ∎

E.6 Proof of Theorem 4.2

Theorem 4.2 is a special of Theorem 4.1 by taking $\alpha=0$.

E.7 Proof of Theorem B.1

We first present three lemmas, based on which we build the proof of the main theorem.

Lemma E.6**.**

Suppose that $\widehat{}\boldsymbol{\Sigma}$ is positive-definite, i.e. $\lambda_{\min}(\widehat{}\boldsymbol{\Sigma})>0$. Define $\boldsymbol{\varepsilon}_{t}=\widehat{}\boldsymbol{\Sigma}^{1/2}(\boldsymbol{\theta}_{t}-\widehat{}\boldsymbol{\Sigma}^{-1}\widehat{}\mathbf{w})$ for $t\geq 0$ and

[TABLE]

If $\alpha\geq 0$ is appropriately chosen such that $\Delta=\max\limits_{k\in[m]}\|\boldsymbol{\Delta}_{k}\|_{2}\leq 1/2$. Then

[TABLE]

which guarantees linear convergence of $\left\{\boldsymbol{\varepsilon}_{t}\right\}_{t=0}^{\infty}$.

Proof of Lemma E.6.

Define $\boldsymbol{\varepsilon}_{t,k}=\widehat{}\boldsymbol{\Sigma}^{1/2}(\boldsymbol{\theta}_{t,k}-\widehat{}\boldsymbol{\Sigma}^{-1}\widehat{}\mathbf{w})$. Then

[TABLE]

Define $\widetilde{\boldsymbol{\Sigma}}_{k}^{(1)}=\widehat{}\boldsymbol{\Sigma}^{-1/2}\widehat{\boldsymbol{\Sigma}}_{k}\widehat{}\boldsymbol{\Sigma}^{-1/2}$ for $k\in[m]$. The fact

[TABLE]

gives $\widehat{}\boldsymbol{\Sigma}^{1/2}(\widehat{\boldsymbol{\Sigma}}_{k}+\alpha\mathbf{I})^{-1}\widehat{}\boldsymbol{\Sigma}^{1/2}=(\widetilde{\boldsymbol{\Sigma}}_{k}^{(1)}+\alpha\widehat{}\boldsymbol{\Sigma}^{-1})^{-1}$ and $\boldsymbol{\varepsilon}_{t+1,k}=[\mathbf{I}-(\widetilde{\boldsymbol{\Sigma}}_{k}^{(1)}+\alpha\widehat{}\boldsymbol{\Sigma}^{-1})^{-1}]\boldsymbol{\varepsilon}_{t}$.

Define $\widehat{}\mathbf{D}=(\mathbf{I}+\alpha\widehat{}\boldsymbol{\Sigma}^{-1})^{-1}$ and $\widetilde{{\boldsymbol{\Sigma}}}_{k}=\widehat{}\mathbf{D}^{1/2}\widetilde{\boldsymbol{\Sigma}}_{k}^{(1)}\widehat{}\mathbf{D}^{1/2}$. From

[TABLE]

and

[TABLE]

we get

[TABLE]

where $\mathbf{C}_{k}=(\mathbf{I}+\boldsymbol{\Delta}_{k})^{-1}$ and $\mathbf{C}=\frac{1}{m}\sum_{k=1}^{m}\mathbf{C}_{k}$. Let $\mathbf{R}_{k}=\mathbf{C}_{k}-(\mathbf{I}-\boldsymbol{\Delta}_{k})$ and $\mathbf{R}=\frac{1}{m}\sum_{k=1}^{m}\mathbf{R}_{k}=\mathbf{C}-\mathbf{I}$. We have $\|\mathbf{I}-\widehat{}\mathbf{D}^{1/2}\mathbf{C}\widehat{}\mathbf{D}^{1/2}\|_{2}=\|\mathbf{I}-\widehat{}\mathbf{D}^{1/2}(\mathbf{I}+\mathbf{R})\widehat{}\mathbf{D}^{1/2}\|_{2}$. Below we control the right-hand side.

By $\Delta=\max\limits_{k\in[m]}\|\boldsymbol{\Delta}_{k}\|_{2}\leq 1/2$ and Lemma F.4, we obtain that $\|\mathbf{C}_{k}\|_{2}\leq\frac{{1}}{1-\Delta}\leq 2$ and $\|\mathbf{R}_{k}\|_{2}\leq 2\Delta^{2}$. Consequently, $\|\mathbf{C}\|_{2}\leq 2$ and $\|\mathbf{R}\|_{2}\leq 2\Delta^{2}\leq 1/2$. Then we obtain that

[TABLE]

Consequently,

[TABLE]

Let $\{\widehat{\lambda}_{j}\}_{j=1}^{p}$ be the eigenvalues of $\widehat{\boldsymbol{\Sigma}}$ sorted in descending order. Since $\widehat{}\mathbf{D}$ has eigenvalues $\{(1+\alpha/\widehat{\lambda}_{j})^{-1}\}_{j=1}^{p}\subseteq(0,1]$, the eigenvalues of $\mathbf{I}-(1\pm 2\Delta^{2})\widehat{}\mathbf{D}$ are $\left\{1-\frac{1\pm 2\Delta^{2}}{1+\alpha/\widehat{\lambda}_{j}}\right\}_{j=1}^{p}$. Then

[TABLE]

By elementary calculation and the fact $2\Delta^{2}\leq 1/2<1$ we get

[TABLE]

Therefore,

[TABLE]

∎

Lemma E.7**.**

Suppose Assumption B.1 hold. Then there exists a constant $C$ determined by $\|\mathbf{u}_{i}\|_{\psi_{2}}$, such that $\mathbb{P}(\Delta\leq 1/2)\geq 1-2me^{-n/C}$ holds under either of the two conditions: (i) $n\geq Cp$ and $\alpha\geq 0$; (ii) $\alpha\geq C{\rm Tr}(\boldsymbol{\Sigma})/n$.

Proof of Lemma E.7.

Define $\widehat{}\mathbf{S}_{k}=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}$ and $\bar{\mathbf{u}}_{(k)}=\frac{1}{n}\sum_{\in\mathcal{I}_{k}}\mathbf{u}_{i}$ for $k\in[m]$. Then we have $\widehat{}\boldsymbol{\Sigma}_{k}=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}=\begin{pmatrix}1&\bar{\mathbf{u}}_{(k)}^{\top}\\ \bar{\mathbf{u}}_{(k)}&\widehat{}\mathbf{S}_{k}\end{pmatrix}$. Let $\boldsymbol{\Sigma}^{*}={\mathbb{E}}\widehat{}\boldsymbol{\Sigma}_{k}=\begin{pmatrix}1&\mathbf{0}\\ \mathbf{0}&\boldsymbol{\Sigma}\end{pmatrix}$ and observe that

[TABLE]

Let $\mathbf{B}_{k}=(\widehat{}\boldsymbol{\Sigma}+\alpha\mathbf{I})^{-1/2}(\widehat{\boldsymbol{\Sigma}}_{k}-\boldsymbol{\Sigma}^{*})(\widehat{}\boldsymbol{\Sigma}+\alpha\mathbf{I})^{-1/2}$. Since

[TABLE]

it boils down to bound $\|\mathbf{B}_{k}\|_{2}$. To do this, we let $\mathbf{A}_{0}=(\boldsymbol{\Sigma}^{*}+\alpha\mathbf{I})^{1/2}(\widehat{}\boldsymbol{\Sigma}+\alpha\mathbf{I})^{-1/2}$ and write

[TABLE]

Define $\mathbf{D}=(\mathbf{I}+\alpha(\boldsymbol{\Sigma}^{*})^{-1})^{-1}$ and $\widehat{\boldsymbol{\Sigma}}^{(1)}_{k}=(\boldsymbol{\Sigma}^{*}+\alpha\mathbf{I})^{-1/2}\widehat{\boldsymbol{\Sigma}}_{k}(\boldsymbol{\Sigma}^{*}+\alpha\mathbf{I})^{-1/2}$. On the one hand,

[TABLE]

On the other hand,

[TABLE]

where we used Lemma F.4. Based on these, we have

[TABLE]

and it suffices to prove under the given conditions that

[TABLE]

holds for some constant $C$.

By definition, we have

[TABLE]

and as a result,

[TABLE]

Here we used a simple fact that $\left\|\begin{pmatrix}\mathbf{0}&\mathbf{A}^{\top}\\ \mathbf{A}&\mathbf{0}\end{pmatrix}\right\|_{2}=\|\mathbf{A}\|_{2}$ for any matrix $\mathbf{A}$.

Observe that $\mathbf{v}_{i}=(\boldsymbol{\Sigma}+\alpha\mathbf{I})^{-1/2}\mathbf{u}_{i}$ is a sub-gaussian random variable with zero mean and covariance matrix $(1+\alpha\boldsymbol{\Sigma}^{-1})^{-1}$. On the other hand, $(\boldsymbol{\Sigma}+\alpha\mathbf{I})^{-1/2}(\widehat{}\mathbf{S}_{k}-\boldsymbol{\Sigma})(\boldsymbol{\Sigma}+\alpha\mathbf{I})^{-1/2}=\frac{1}{n}\sum_{i\in\mathcal{I}_{k}}\mathbf{v}_{i}\mathbf{v}_{i}^{\top}-(1+\alpha\boldsymbol{\Sigma}^{-1})^{-1}$. Lemma F.5 forces

[TABLE]

where $C$ is the constant therein. These estimates and (E.12) lead to (E.11). ∎

Following the similar idea in the proof above, we get the following results.

Lemma E.8**.**

Suppose Assumption B.1 holds with $C$ being the constant in Lemma E.7. Then

• •

$\mathbb{P}(\|(\boldsymbol{\Sigma}^{*})^{-1/2}(\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}^{*})(\boldsymbol{\Sigma}^{*})^{-1/2}\|_{2}\leq 1/2)\geq 1-2me^{-N/C}$;

• •

$\mathbb{P}(\Delta\leq C_{2}^{\prime}\sqrt{p/n})\geq 1-2e^{-C_{1}^{\prime}p}$* holds for some constants $C_{1}^{\prime}$ and $C_{2}^{\prime}$.*

Now we come back to the main proof. We will use the three lemmas above to show that with high probability, $\lambda_{\max}(\widehat{}\boldsymbol{\Sigma})/\lambda_{\min}(\widehat{}\boldsymbol{\Sigma})\leq 3\kappa$ and

[TABLE]

Then we conclude the proof by induction and some simple linear algebra.

First, Lemma E.8 asserts that with high probability, $\|(\boldsymbol{\Sigma}^{*})^{-1/2}(\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}^{*})(\boldsymbol{\Sigma}^{*})^{-1/2}\|_{2}\leq 1/2$. On this event, we have $-\frac{1}{2}\boldsymbol{\Sigma}^{*}\preceq\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}^{*}\preceq\frac{1}{2}\boldsymbol{\Sigma}^{*}$ and thus $\frac{1}{2}\boldsymbol{\Sigma}^{*}\preceq\widehat{\boldsymbol{\Sigma}}\preceq\frac{3}{2}\boldsymbol{\Sigma}^{*}$. Hence

[TABLE]

and $\lambda_{\max}(\widehat{}\boldsymbol{\Sigma})/\lambda_{\min}(\widehat{}\boldsymbol{\Sigma})\leq 3\lambda_{\max}(\boldsymbol{\Sigma}^{*})/\lambda_{\min}(\boldsymbol{\Sigma}^{*})=3\kappa$.

Second, Lemma E.6 forces

[TABLE]

Lemmas E.7 and E.8 imply that $\Delta\leq 1/2$, $\Delta\leq\tilde{C}\sqrt{p/n}$, and $\lambda_{\min}(\widehat{}\boldsymbol{\Sigma})\geq 1$ hold simultaneously with high probability, where $\tilde{C}$ is some constant. On this event, we have

[TABLE]

for some constants $C_{2}$ and $C_{3}$. Then we get (E.13) from the estimates above and complete the proof.

Appendix F Technical lemmas

The following lemma lists basic properties of strongly convex functions, which can be found in standard textbooks on convex optimization (Nesterov, 2013).

Lemma F.1**.**

Suppose $f$ is a convex function defined on some convex open set $\Omega\subseteq{\mathbb{R}}^{p}$, and $\partial f(\mathbf{x})$ denotes its subdifferential set at $\mathbf{x}\in\Omega$. The followings are equivalent:

• •

$f$* is $\rho$-strongly convex in $\Omega$;*

• •

$f[(1-t)\mathbf{x}+t\mathbf{y}]\leq(1-t)f(\mathbf{x})+tf(\mathbf{y})-(\rho/2)t(1-t)\|\mathbf{y}-\mathbf{x}\|_{2}^{2}$, $\forall\mathbf{x},\mathbf{y}\in\Omega$ and $t\in[0,1]$;

• •

$\langle\mathbf{h}-\mathbf{g},\mathbf{y}-\mathbf{x}\rangle\geq\rho\|\mathbf{y}-\mathbf{x}\|_{2}^{2}$, $\forall\mathbf{x},\mathbf{y}\in\Omega$, $\mathbf{g}\in\partial f(\mathbf{x})$ and $\mathbf{h}\in\partial f(\mathbf{y})$.

If any of the above holds, $f$ is said to be $\rho$-strongly convex.

Lemma F.2**.**

Let $f:~{}{\mathbb{R}}^{p}\to{\mathbb{R}}$ be a convex function. Suppose there exists $\mathbf{x}\in{\mathbb{R}}^{p}$ and $r>0$ such that $f$ is $\rho$-strongly convex in $B(\mathbf{x},r)$. If $\|\mathbf{h}-\mathbf{g}\|_{2}<\rho r$ holds for some $\mathbf{g}\in\partial f(\mathbf{x})$ and $\mathbf{h}\in\partial f(\mathbf{y})$, then $\|\mathbf{y}-\mathbf{x}\|_{2}\leq\|\mathbf{h}-\mathbf{g}\|_{2}/\rho\leq r$.

Proof.

If we know a priori that $\|\mathbf{y}-\mathbf{x}\|_{2}\leq r$, then we use the strong convexity of $f$ in $B(\mathbf{x},r)$ and Cauchy-Schwarz inequality to obtain

[TABLE]

and get the desired result. Suppose on the contrary that $\|\mathbf{y}-\mathbf{x}\|_{2}>r$, and define $\bar{\mathbf{y}}=\mathbf{x}+r(\mathbf{y}-\mathbf{x})/\|\mathbf{y}-\mathbf{x}\|_{2}$. Then $\|\bar{\mathbf{y}}-\mathbf{x}\|_{2}=r$. The strong convexity of $f$ in $B(\mathbf{x},r)$ and Lemma F.1 yield

[TABLE]

By the convexity of $f$, we always have

[TABLE]

Summing up the two inequalities above, we get

[TABLE]

where we also used the Cauchy-Schwarz inequality. Then $\|\mathbf{h}-\mathbf{g}\|_{2}\geq\rho\|\bar{\mathbf{y}}-\mathbf{x}\|_{2}=\rho r$ leads to contradiction. Hence, we must have only the case $\|\mathbf{y}-\mathbf{x}\|_{2}\leq r$. ∎

Lemma F.3**.**

Let $f:~{}{\mathbb{R}}^{p}\to{\mathbb{R}}$ be a convex function. For any $\mathbf{x},\mathbf{y}\in{\mathbb{R}}^{p}$, we have

[TABLE]

If $\inf_{\mathbf{x}\in{\mathbb{R}}^{p}}f(\mathbf{x})>-\infty$, then $\|\mathrm{prox}_{f}(\mathbf{x})-\mathbf{x}^{*}\|_{2}\leq\|\mathbf{x}-\mathbf{x}^{*}\|_{2}$ and

[TABLE]

hold for any $\mathbf{x}^{*}\in{\rm argmin}_{\mathbf{x}\in{\mathbb{R}}^{p}}f(\mathbf{x})$.

If $f$ is $\rho$-strongly convex in $B(\mathbf{x}^{*},r)$ for some $r>0$ and $\mathbf{x}^{*}={\rm argmin}_{\mathbf{x}\in{\mathbb{R}}^{p}}f(\mathbf{x})$, then $\|\mathrm{prox}_{\alpha^{-1}f}(\mathbf{x})-\mathbf{x}^{*}\|_{2}\leq\frac{\alpha}{\alpha+\rho}\|\mathbf{x}-\mathbf{x}^{*}\|_{2}$, $\forall\mathbf{x}\in B(\mathbf{x}^{*},r)$ and $\alpha>0$.

Proof of Lemma F.3.

The first claim is the well-known “firm non-expansiveness” property of the proximal mapping (Parikh and Boyd, 2014).

If $\inf_{\mathbf{x}\in{\mathbb{R}}^{p}}f(\mathbf{x})>-\infty$, then any $\mathbf{x}^{*}\in{\rm argmin}_{\mathbf{x}\in{\mathbb{R}}^{p}}f(\mathbf{x})$ is a fixed point of $\mathrm{prox}_{f}$. The firm non-expansiveness with $\mathbf{y}=\mathbf{x}^{*}$ yields

[TABLE]

and then $\|\mathrm{prox}_{f}(\mathbf{x})-\mathbf{x}^{*}\|_{2}\leq\|\mathbf{x}-\mathbf{x}^{*}\|_{2}$. The next claim is proved by

[TABLE]

where the inequality follows from (F.1).

For the last claim, we fix any $\alpha>0$ and $\mathbf{x}\in B(\mathbf{x}^{*},r)$ and define $\mathbf{x}^{+}=\mathrm{prox}_{\alpha^{-1}f}(\mathbf{x})$. Then $\|\mathbf{x}^{+}-\mathbf{x}^{*}\|_{2}\leq\|\mathbf{x}-\mathbf{x}^{*}\|_{2}<r$. The optimality conditions for $\mathbf{x}^{*}={\rm argmin}_{\mathbf{y}\in{\mathbb{R}}^{p}}f(\mathbf{y})$ and $\mathbf{x}^{+}={\rm argmin}_{\mathbf{y}\in{\mathbb{R}}^{p}}\{f(\mathbf{y})+(\alpha/2)\|\mathbf{y}-\mathbf{x}\|_{2}^{2}\}$ imply that $\mathbf{0}\in\partial f(\mathbf{x}^{*})$ and $-\alpha(\mathbf{x}^{+}-\mathbf{x})\in\partial f(\mathbf{x}^{+})$. Since $f$ is $\rho$-strongly convex in $B(\mathbf{x}^{*},r)$, Lemma F.1 forces

[TABLE]

and thus $\|\mathbf{x}^{+}-\mathbf{x}^{*}\|_{2}\leq\frac{\alpha}{\alpha+\rho}\|\mathbf{x}-\mathbf{x}^{*}\|_{2}$. ∎

Lemma F.4** (Neumann expansion).**

Let $\|\cdot\|$be a submultiplicative matrix norm with $\|\mathbf{I}\|=1$. When $\|\mathbf{M}\|<1$, we have $(\mathbf{I}-\mathbf{M})^{-1}=\sum_{j=0}^{\infty}\mathbf{M}^{j}=\mathbf{I}+\mathbf{M}+\mathbf{M}(\mathbf{I}-\mathbf{M})^{-1}\mathbf{M}$, $\|\mathbf{I}-\mathbf{M}\|\leq 1/(1-\|\mathbf{M}\|)$ and $\|(\mathbf{I}-\mathbf{M})^{-1}-(\mathbf{I}+\mathbf{M})\|\leq\|\mathbf{M}\|^{2}/(1-\|\mathbf{M}\|)$.

Lemma F.5**.**

Let $\mathbf{S}\succ 0$ and $\alpha\geq 0$ be deterministic, $\mathbf{A}=(\mathbf{I}+\alpha\mathbf{S}^{-1})^{-1}$, $\{\mathbf{u}_{i}\}_{i=1}^{n}\subseteq{\mathbb{R}}^{d}$ be i.i.d. sub-gaussian random vectors with zero mean and covariance matrix $\mathbf{A}$, $\bar{\mathbf{u}}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{u}_{i}$ and $\widehat{\mathbf{A}}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{u}_{i}\mathbf{u}_{i}^{T}$. Then the following hold:

There exists some positive constant $C$ that only depends on $\|\mathbf{u}_{1}\|_{\psi_{2}}$, such that $\mathbb{P}(\|\bar{\mathbf{u}}\|_{2}>1/10)\leq e^{-n/C}$ and ${\mathbb{P}}(\|\widehat{\mathbf{A}}-\mathbf{A}\|_{2}>1/10)\leq e^{-n/C}$ hold under any one of the two conditions holds: (i) $n\geq Cd$; (ii) $\alpha\geq C{\rm Tr}(\mathbf{S})/n$. 2. 2.

Suppose that $n\geq cd$ for some constant $c>0$. There exist positive constants $C_{1}^{\prime}$ and $C_{2}^{\prime}$ such that ${\mathbb{P}}(\max\{\lVert\bar{}\mathbf{u}\rVert_{2},\lVert\widehat{}\mathbf{A}-\mathbf{A}\rVert_{2}\}\geq C_{2}^{\prime}\sqrt{d/n})\leq 2e^{-C_{1}^{\prime}d}$.

Proof of Lemma F.5.

Let $\lambda_{1}>\cdots>\lambda_{d}>0$ be the eigenvalues of $\mathbf{S}$, and $C>0$ be some constant to be determined. When $n\geq Cd$, the fact $\|\mathbf{A}\|_{2}=1/(1+\alpha/\lambda_{1})\leq 1$ implies $n\geq C{\rm Tr}(\mathbf{A})$. Also, if $\alpha\geq C{\rm Tr}(\mathbf{S})/n$, then the crude estimate

[TABLE]

also leads to $n\geq C{\rm Tr}(\mathbf{A})$. Hence it suffices to find some proper $C$ and show the desired results given $n\geq C{\rm Tr}(\mathbf{A})$.

Now we prove the first statement. We first study concentration of the sample mean vector $\bar{\mathbf{u}}$. Since $\bar{\mathbf{u}}$ is a sub-gaussian random vector with covariance matrix $n^{-1}\mathbf{A}$, Theorem 2.1 in Hsu et al. (2012) asserts the existence of a constant $c_{1}>0$ such that

[TABLE]

Choose any constant $C_{1}\geq 500c_{1}$. Let $t=n/C_{1}$, and suppose that $n\geq C_{1}{\rm Tr}(\mathbf{A})$. Using $\|\mathbf{A}\|_{2}\leq 1$ and ${\rm Tr}(\mathbf{A}^{2})\leq{\rm Tr}(\mathbf{A})\|\mathbf{A}\|_{2}\leq{\rm Tr}(\mathbf{A})$, we get

[TABLE]

Hence ${\mathbb{P}}\left(\|\bar{}\mathbf{u}\|_{2}>1/10\right)\geq e^{-n/C_{1}}$.

Now we come to concentration of the sample covariance matrix $\widehat{}\mathbf{A}$. Let $r(\mathbf{A})={\rm Tr}(\mathbf{A})/\|\mathbf{A}\|_{2}$. According to Theorem 9 in Koltchinskii et al. (2017), there exists a constant $c_{2}\geq 1$ such that the following holds: for any $t\geq 1$, with probability at least $1-e^{-t}$ we have

[TABLE]

Note that the upper bound above can be rewritten as

[TABLE]

Let $C_{2}=100c_{2}^{2}$. When $n\geq C_{2}{\rm Tr}(\mathbf{A})$, by taking $t=n/C_{2}$ we get $\mathbb{P}(\|\widehat{\mathbf{A}}-\mathbf{A}\|_{2}>1/10)\leq e^{-n/C_{2}}$. The proof of the first statement is then finished by taking $C=\max\{C_{1},C_{2}\}$.

We proceed to prove the second statement. Let $t=C_{1}^{\prime}d$ for some constant $C_{1}^{\prime}$. Note that $d\geq{\rm Tr}(\mathbf{A})$ and $t\geq C_{1}^{\prime}{\rm Tr}(\mathbf{A})$. According to (F.2), (F.3) and (F.4), we obtain that with probability at least $1-2e^{-C_{1}^{\prime}d}$, $\lVert\bar{}\mathbf{u}\rVert_{2}\leq\tilde{C}\sqrt{d/n}$ and $\|\widehat{}\mathbf{A}-\mathbf{A}\|_{2}\leq\tilde{C}\max\left\{\sqrt{d/n},d/n\right\}$ hold with some constant $\tilde{C}$. Since $n\geq cd$, we have

[TABLE]

By combining the inequalities above, we obtain that

[TABLE]

with $C_{2}^{\prime}=\tilde{C}\max\{1,1/\sqrt{c}\}$. ∎