OverSketched Newton: Fast Convex Optimization for Serverless Systems

TL;DR

OverSketched Newton is a randomized optimization algorithm designed for serverless systems, using matrix sketching to efficiently approximate the Hessian, resulting in faster convergence and resilience to system stragglers.

Contribution

It introduces a novel Hessian approximation method using matrix sketching tailored for serverless architectures, with proven convergence guarantees and empirical speedups.

Findings

Achieves ~50% reduction in total running time on AWS Lambda.

Provides convergence guarantees for both strongly convex and non-convex problems.

Demonstrates robustness against stragglers in serverless environments.

Abstract

Motivated by recent developments in serverless systems for large-scale computation as well as improvements in scalable randomized matrix algorithms, we develop OverSketched Newton, a randomized Hessian-based optimization algorithm to solve large-scale convex optimization problems in serverless systems. OverSketched Newton leverages matrix sketching ideas from Randomized Numerical Linear Algebra to compute the Hessian approximately. These sketching methods lead to inbuilt resiliency against stragglers that are a characteristic of serverless architectures. Depending on whether the problem is strongly convex or not, we propose different iteration updates using the approximate Hessian. For both cases, we establish convergence guarantees for OverSketched Newton and empirically validate our results by solving large-scale supervised learning problems on real-world datasets. Experiments demonstrate a reduction of ~50% in total running time on AWS Lambda, compared to state-of-the-art distributed optimization schemes.

Figures

Tables

Dataset	Training Samples	Features	Testing samples
Synthetic	$300,000$	$3000$	$100,000$
EPSILON	$400,000$	$2000$	$100,000$
WEBPAGE	$48,000$	$300$	$15,000$
a9a	$32,000$	$123$	$16,000$
EMNIST	$240,000$	$7840$	$40,000$

Equations

$$f({\mathbf{w}}^{\ast })=\underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }f(\mathbf{w}),$$

$$\displaystyle\mathbf{w}_{t+1}=\arg\min_{\mathbf{w}\in\mathbb{R}^{d}}\Big{\{}f(\mathbf{w}_{t})+\nabla f(\mathbf{w}_{t})^{T}(\mathbf{w}-\mathbf{w}_{t})$$

$$\displaystyle+\frac{1}{2}(\mathbf{w}-\mathbf{w}_{t})^{T}\nabla^{2}f(\mathbf{w}_{t})(\mathbf{w}-\mathbf{w}_{t})\Big{\}}.$$

$$\underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }||\nabla f(\mathbf{w})|{|}^2.$$

$${\mathbf{p}}_t=\arg \underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }||\nabla f({\mathbf{w}}_t)+{\mathbf{H}}_t\mathbf{p}|{|}^2.$$

$${\mathbf{w}}_{t+1}={\mathbf{w}}_t+{\mathbf{p}}_t={\mathbf{w}}_t-[{\nabla }^{2}f({\mathbf{w}}_t){]}^{†}\nabla f({\mathbf{w}}_t).$$

$$\mathbf{S}=\frac{1}{\sqrt{N}}({\mathbf{S}}_1,{\mathbf{S}}_2,\cdots ,{\mathbf{S}}_{N+e}),$$

$${\alpha }_t=\underset{\alpha \le 1}{\max }\alpha \text{ such that }f({\mathbf{w}}_t+\alpha {\mathbf{p}}_t)\le f({\mathbf{w}}_t)+\alpha \beta {\mathbf{p}}_t^T\nabla f({\mathbf{w}}_t),$$

$${\alpha }_t=\underset{\alpha \le 1}{\max }\alpha \text{ such that }||\nabla f({\mathbf{w}}_t+\alpha {\mathbf{p}}_t)|{|}^2\le ||\nabla f({\mathbf{w}}_t)|{|}^2+2\alpha \beta {\mathbf{p}}_t^T{\hat{\mathbf{H}}}_t\nabla f({\mathbf{w}}_t).$$

$${\mathbf{w}}_{t+1}={\mathbf{w}}_t+{\alpha }_t{\mathbf{p}}_t.$$

$${\nabla }^{2}f(\mathbf{w})⪰k\mathbf{I};$$

$${\nabla }^{2}f(\mathbf{w})⪯M\mathbf{I};$$

$$||{\nabla }^{2}f(\mathbf{w}+\mathbf{\Delta })-{\nabla }^{2}f(\mathbf{w})|{|}_2\le L||\mathbf{\Delta }|{|}_2,$$

$$f({\mathbf{w}}_{t+1})-f({\mathbf{w}}^{\ast })\le (1-\rho )(f({\mathbf{w}}_t)-f({\mathbf{w}}^{\ast })),$$

$$||{\mathbf{w}}_{t+1}-{\mathbf{w}}^{\ast }|{|}_2\le \frac{25L}{8\gamma }||{\mathbf{w}}_t-{\mathbf{w}}^{\ast }|{|}_2^2+\frac{5ϵ\beta }{\gamma }||{\mathbf{w}}_t-{\mathbf{w}}^{\ast }|{|}_2\text{for}t=1,2,\cdots ,T,$$

$$||({\nabla }^{2}f(\mathbf{w}){)}^{†}|{|}_2\le 1/\eta .$$

$$||{\mathbf{U}}^{\mathbf{T}}\nabla f(\mathbf{w})|{|}^2\ge \nu ||\nabla f(\mathbf{w})|{|}^2\forall \mathbf{w}\in {\mathbb{R}}^{\mathbf{d}}.$$

$$\displaystyle||\nabla f(\mathbf{w}_{t+1})||^{2}\leq\bigg{(}1-2\beta\alpha\nu\frac{(1-\epsilon)\eta}{M(1+\epsilon)}\bigg{)}||\nabla f(\mathbf{w}_{t})||^{2},$$

$$\min_{\mathbf{w}\in\mathbb{R}^{d}}~{}\bigg{\{}f(\mathbf{w})=\frac{1}{n}\sum_{i=1}^{n}\log(1+e^{-y_{i}\mathbf{w}^{T}\mathbf{x}_{i}})+\frac{\lambda}{2}\|\mathbf{w}\|^{2}_{2}\bigg{\}}.$$

$$\nabla f(\mathbf{w})=\frac{1}{n}\sum _{i=1}^n\frac{-y_i{\mathbf{x}}_i}{1+e^{y_i{\mathbf{w}}_i^T{\mathbf{x}}_i}}+\lambda \mathbf{w}.$$

$$f(\mathbf{W})=\sum _{n=1}^N\left[\sum _{k=1}^K{y}_{kn}{\mathbf{w}}_k^T{\mathbf{x}}_n-\log \sum _{l=1}^K\exp \left({\mathbf{w}}_l^T{\mathbf{x}}_n\right)\right].$$

$$\nabla f_i(\mathbf{W})=\sum _{n=1}^N\left[\frac{\exp \left({\mathbf{w}}_i^T{\mathbf{x}}_n\right)}{{\sum }_{l=1}^K\exp \left({\mathbf{w}}_l^T{\mathbf{x}}_n\right)}-y_{in}\right]{\mathbf{x}}_n\forall i\in [1,k],$$

$${\mathbf{H}}_{ij}(\mathbf{W})=\frac{d}{d{\mathbf{w}}_j}\nabla f_i(\mathbf{W})=\frac{d}{d{\mathbf{w}}_j}\mathbf{X}{\mathbf{\beta }}_{\mathbf{i}}=\mathbf{X}\frac{d}{d{\mathbf{w}}_j}{\mathbf{\beta }}_i=\mathbf{X}{\mathbf{Z}}_{\mathbf{ij}}{\mathbf{X}}^T$$

$$Z_{ij}(n)=\frac{\exp ({\alpha }_{in})}{{\sum }_{l=1}^K\exp ({\alpha }_{ln})}\left(\mathbb{I}(i=j)-\frac{\exp ({\alpha }_{jn})}{{\sum }_{l=1}^K\exp ({\alpha }_{ln})}\right)\forall n\in [1,N],$$

$$\displaystyle\nabla^{2}f(\mathbf{W})=\left[\begin{array}[]{ccc}{\mathbf{H}_{11}}&{\cdots}&{\mathbf{H}_{1K}}\\ {\vdots}&{\ddots}&{\vdots}\\ {\mathbf{H}_{K1}}&{\cdots}&{\mathbf{H}_{KK}}\end{array}\right]=\left[\begin{array}[]{ccc}{\mathbf{X}\mathbf{Z}_{11}}\mathbf{X}^{T}&{\cdots}&\mathbf{X}{\mathbf{Z}_{1K}}\mathbf{X}^{T}\\ {\vdots}&{\ddots}&{\vdots}\\ \mathbf{X}{\mathbf{Z}_{K1}}\mathbf{X}^{T}&{\cdots}&\mathbf{X}{\mathbf{Z}_{KK}\mathbf{X}^{T}}\end{array}\right]=\bar{\mathbf{X}}\bar{\mathbf{Z}}\bar{\mathbf{X}}^{T},$$

$$\underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }\frac{1}{2n}||{\mathbf{X}}^T\mathbf{w}-\mathbf{y}|{|}_2^2+\frac{\lambda }{2}\Vert \mathbf{w}{\Vert }_2^2.$$

$$\underset{\mathbf{Ax}\le \mathbf{b}}{\mathrm{minimize}}{\mathbf{c}}^T\mathbf{x},$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^m}{\min }f(\mathbf{x})=\underset{\mathbf{x}\in {\mathbb{R}}^m}{\min }\left(\tau {\text{c}}^T\mathbf{x}-\sum _{i=1}^n\log (b_i-{\mathbf{a}}_i\mathbf{x})\right),$$

$${\nabla }^{2}f(\mathbf{x})={\mathbf{A}}^T\text{diag}\frac{1}{(b_i-{\alpha }_i{)}^2}\mathbf{A}.$$

$$\underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }\frac{1}{2}||\mathbf{Xw}-\mathbf{y}|{|}_2^2+\lambda ||\mathbf{w}|{|}_1,$$

$$\underset{\begin{array}{c}||{\mathbf{X}}^T\mathbf{z}|{|}_{\infty }\le \lambda ,\mathbf{z}\in {\mathbb{R}}^n\end{array}}{\min }\frac{1}{2}||\mathbf{y}-\mathbf{z}|{|}_2^2,$$

Full text

1 Introduction

In recent years, there has been tremendous growth in users performing distributed computing operations on the cloud, largely due to extensive and inexpensive commercial offerings like Amazon Web Services (AWS), Google Cloud, Microsoft Azure, etc. Serverless platforms—such as AWS Lambda, Cloud functions and Azure Functions—penetrate a large user base by provisioning and managing the servers on which the computation is performed. These platforms abstract away the need for maintaining servers, since this is done by the cloud provider and is hidden from the user—hence the name serverless. Moreover, allocation of these servers is done expeditiously which provides greater elasticity and easy scalability. For example, up to ten thousand machines can be allocated on AWS Lambda in less than ten seconds [1, 2, 3, 4].

The use of serverless systems is gaining significant research traction, primarily due to its massive scalability and convenience in operation. It is forecasted that the market share of serverless will grow by USD 9.16 billion during 2019-2023 (at a CAGR of 11%) [5]. Indeed, according to the Berkeley view on Serverless Computing [6], serverless systems are expected to dominate the cloud scenario and become the default computing paradigm in the coming years while client-server based computing will witness a considerable decline. For these reasons, using serverless systems for large-scale computation has garnered significant attention from the systems community [3, 4, 7, 8, 9, 10, 11, 12].

Due to several crucial differences between the traditional High Performance Computing (HPC) / serverful and serverless architectures, existing distributed algorithms cannot, in general, be extended to serverless computing. First, unlike serverful computing, the number of inexpensive workers in serverless platforms is flexible, often scaling into the thousands [3, 4]. This heavy gain in the computation power, however, comes with the disadvantage that the commodity workers in serverless architecture are ephemeral and have low memory.111For example, serverless nodes in AWS Lambda, Google Cloud Functions and Microsoft Azure Functions have a maximum memory of 3 GB, 2 GB and 1.5 GB, respectively, and a maximum runtime of 900 seconds, 540 seconds and 300 seconds, respectively (these numbers may change over time). The ephemeral nature of the workers in serverless systems requires that new workers should be invoked every few iterations and data should be communicated to them. Moreover, the workers do not communicate amongst themselves, and instead they read/write data directly from/to a single high-latency data storage entity (e.g., cloud storage like AWS S3 [3]).

Second, unlike HPC/serverful systems, nodes in the serverless systems suffer degradation due to what is known as system noise. This can be a result of limited availability of shared resources, hardware failure, network latency, etc. [13, 14]. This results in job time variability, and hence a subset of much slower nodes, often called stragglers. These stragglers significantly slow the overall computation time, especially in large or iterative jobs. In Fig. 1, we plot the running times for a distributed matrix multiplication job with 3600 workers on AWS Lambda and demonstrate the effect of stragglers on the total job time. In fact, our experiments consistently demonstrate that at least $2\%$ workers take significantly longer than the median job time, severely degrading the overall efficiency of the system.

Due to these issues, first-order methods, e.g., gradient descent and Nesterov Accelerated Gradient (NAG) methods, tend to perform poorly on distributed serverless architectures [15]. Their slower convergence is made worse on serverless platforms due to persistent stragglers. The straggler effect incurs heavy slowdown due to the accumulation of tail times as a result of a subset of slow workers occurring in each iteration.

Compared to first-order optimization algorithms, second-order methods—which use the gradient as well as Hessian information—enjoy superior convergence rates. For instance, Newton’s method enjoys quadratic convergence for strongly convex and smooth problems, compared to the linear convergence of gradient descent [16]. Moreover, second-order methods do not require step-size tuning and unit step-size provably works for most problems. These methods have a long history in optimization and scientific computing (see, e.g., [16]), but they are less common in machine learning and data science. This is partly since stochastic first order methods suffice for downstream problems [17] and partly since naive implementations of second order methods can perform poorly [18]. However, recent theoretical work has addressed many of these issues [19, 20, 21, 22, 23], and recent implementations have shown that high-quality implementations of second order stochastic optimization algorithms can beat state-of-the-art in machine learning applications [24, 25, 26, 27, 28] in traditional systems.

1.1 Main Contributions

In this paper, we argue that second-order methods are highly compatible with serverless systems that provide extensive computing power by invoking thousands of workers but are limited by the communication costs and hence the number of iterations; and, to address the challenges of ephemeral workers and stragglers in serverless systems, we propose and analyze a randomized and distributed second-order optimization algorithm, called OverSketched Newton. OverSketched Newton uses the technique of matrix sketching from Sub-Sampled Newton (SSN) methods [19, 20, 21, 22], which are based on sketching methods from Randomized Numerical Linear Algebra (RandNLA) [29, 30, 31], to obtain a good approximation for the Hessian, instead of calculating the full Hessian.

OverSketched Newton has two key components. For straggler-resilient Hessian calculation in serverless systems, we use the sparse sketching based randomized matrix multiplication method from [32]. For straggler mitigation during gradient calculation, we use the recently proposed technique based on error-correcting codes to create redundant computation [33, 34, 35]. We prove that, for strongly convex functions, the local convergence rate of OverSketched Newton is linear-quadratic, while its global convergence rate is linear. Then, going beyond the usual strong convexity assumption for second-order methods, we adapt OverSketched Newton using ideas from [22]. For such functions, we prove that a linear convergence rate can be guaranteed with OverSketched Newton. To the best of our knowledge, this is the first work to prove convergence guarantees for weakly-convex problems when the Hessian is computed approximately using ideas from RandNLA.

We extensively evaluate OverSketched Newton on AWS Lambda using several real-world datasets obtained from the LIBSVM repository [36], and we compare OverSketched Newton with several first-order (gradient descent, Nesterov’s method, etc.) and second-order (exact Newton’s method [16], GIANT [24], etc.) baselines for distributed optimization. We further evaluate and compare different techniques for straggler mitigation, such as speculative execution, coded computing [33, 34], randomization-based sketching [32] and gradient coding [37]. We demonstrate that OverSketched Newton is at least 9x and 2x faster than state-of-the-art first-order and second-order schemes, respectively, in terms of end-to-end training time on AWS Lambda. Moreover, we show that OverSketched Newton on serverless systems outperforms existing distributed optimization algorithms in serverful systems by at least $30\%$.

1.2 Related Work

Our results tie together three quite different lines of work, each of which we review here briefly.

Existing Straggler Mitigation Schemes: Strategies like speculative execution have been traditionally used to mitigate stragglers in popular distributed computing frameworks like Hadoop MapReduce [38] and Apache Spark [39]. Speculative execution works by detecting workers that are running slower than expected and then allocating their tasks to new workers without shutting down the original straggling task. The worker that finishes first communicates its results. This has several drawbacks, e.g. constant monitoring of tasks is required and late stragglers can still hurt the efficiency.

Recently, many coding-theoretic ideas have been proposed to introduce redundancy into the distributed computation for straggler mitigation (e.g. see [33, 34, 37, 35, 40, 41]). The idea of coded computation is to generate redundant copies of the result of distributed computation by encoding the input data using error-correcting-codes. These redundant copies are then used to decode the output of the missing stragglers. Our algorithm to compute gradients in a distributed straggler-resilient manner uses codes to mitigate stragglers, and we compare our performance with speculative execution.

Approximate Second-order Methods: In many machine learning applications, where the data itself is noisy, using the exact Hessian is not necessary. Indeed, using ideas from RandNLA, one can prove convergence guarantees for SSN methods on a single machine, when the Hessian is computed approximately [23, 19, 20, 21]. To accomplish this, many sketching schemes can be used (sub-Gaussian, Hadamard, random row sampling, sparse Johnson-Lindenstrauss, etc. [29, 30]), but these methods cannot tolerate stragglers, and thus they do not perform well in serverless environments.

This motivates the use of the OverSketch sketch from our recent work in [32]. OverSketch has many nice properties, like subspace embedding, sparsity, input obliviousness, and amenability to distributed implementation. To the best of our knowledge, this is the first work to prove and evaluate convergence guarantees for algorithms based on OverSketch. Our guarantees take into account the amount of communication at each worker and the number of stragglers, both of which are a property of distributed systems.

There has also been a growing research interest in designing and analyzing distributed implementations of stochastic second-order methods  [42, 43, 44, 24, 45, 46]. However, these implementations are tailored for serverful distributed systems. Our focus, on the other hand, is on serverless systems.

Distributed Optimization on Serverless Systems: Optimization over the serverless framework has garnered significant interest from the research community. However, these works either evaluate and benchmark existing algorithms (e.g., see [9, 10, 11]) or focus on designing new systems frameworks for faster optimization (e.g., see [12]) on serverless. To the best of our knowledge, this is the first work that proposes a large-scale distributed optimization algorithm that specifically caters to serverless architectures with provable convergence guarantees. We exploit the advantages offered by serverless systems while mitigating the drawbacks such as stragglers and additional overhead per invocation of workers.

2 Newton’s Method: An Overview

We are interested in solving on serverless systems in a distributed and straggler-resilient manner problems of the form:

[TABLE]

where $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a closed and convex function bounded from below. In the Newton’s method, the update at the $(t+1)$-th iteration is obtained by minimizing the Taylor’s expansion of the objective function $f(\cdot)$ at $\mathbf{w}_{t}$, that is

[TABLE]

For strongly convex $f(\cdot)$, that is, when $\nabla^{2}f(\cdot)$ is invertible, Eq. (2) becomes $\mathbf{w}_{t+1}=\mathbf{w}_{t}-\mathbf{H}_{t}^{-1}\nabla f(\mathbf{w}_{t}),$ where $\mathbf{H}_{t}=\nabla^{2}f(\mathbf{w}_{t})$ is the Hessian matrix at the $t$-th iteration. Given a good initialization and assuming that the Hessian is Lipschitz, the Newton’s method satisfies the update $||\mathbf{w}_{t+1}-\mathbf{w}^{*}||_{2}\leq c||\mathbf{w}_{t}-\mathbf{w}^{*}||^{2}_{2},$ for some constant $c>0$, implying quadratic convergence [16].

One shortcoming for the classical Newton’s method is that it works only for strongly convex objective functions. In particular, if $f$ is weakly-convex222For the sake of clarity, we call a convex function weakly-convex if it is not strongly convex., that is, if the Hessian matrix is not positive definite, then the objective function in (2) may be unbounded from below. To address this shortcoming, authors in [22] recently proposed a variant of Newton’s method, called Newton-Minimum-Residual (Newton-MR). Instead of (1), Newton-MR considers the following auxiliary optimization problem:

[TABLE]

Note that the minimizers of this auxiliary problem and (1) are the same when $f(\cdot)$ is convex. Then, the update direction in the $(t+1)$-th iteration is obtained by minimizing the Taylor’s expansion of $||\nabla f(\mathbf{w}_{t}+\mathbf{p})||^{2}$, that is,

[TABLE]

The general solution of the above problem is given by $\mathbf{p}=-[\mathbf{H}_{t}]^{\dagger}\nabla f(\mathbf{w}_{t})+(\mathbf{I}-\mathbf{H}_{t}[\mathbf{H}_{t}]^{\dagger})\textbf{q},~{}\forall~{}\textbf{q}\in\mathbb{R}^{d}$, where $[\cdot]^{\dagger}$ is the Moore-Penrose inverse. Among these, the minimum norm solution is chosen, which gives the update direction in the $t$-th iteration as $\mathbf{p}_{t}=-\mathbf{H}^{\dagger}_{t}\nabla f(\mathbf{w}_{t})$. Thus, the model update is

[TABLE]

OverSketched Newton considers both of these variants.

3 OverSketched Newton

In many applications like machine learning where the training data itself is noisy, using the exact Hessian is not necessary. Indeed, many results in the literature prove convergence guarantees for Newton’s method when the Hessian is computed approximately using ideas from RandNLA for a single machine (e.g. [23, 19, 20, 47]). In particular, these methods perform a form of dimensionality reduction for the Hessian using random matrices, called sketching matrices. Many popular sketching schemes have been proposed in the literature, for example, sub-Gaussian, Hadamard, random row sampling, sparse Johnson-Lindenstrauss, etc. [29, 30]. Inspired from these works, we present OverSketched Newton, a stochastic second order algorithm for solving—on serverless systems, in a distributed, straggler-resilient manner—problems of the form (1).

Distributed straggler-resilient gradient computation: OverSketched Newton computes the full gradient in each iteration by using tools from error-correcting codes [33, 34]. Our key observation is that, for several commonly encountered optimization problems, gradient computation relies on matrix-vector multiplications (see Sec. 4 for examples). We leverage coded matrix multiplication technique from [34] to perform the large-scale matrix-vector multiplication in a distributed straggler-resilient manner. The idea of coded matrix multiplication is explained in Fig. 3; detailed steps are provided in Algorithm 1.

Distributed straggler-resilient approximate Hessian computation: For several commonly encountered optimization problems, Hessian computation involves matrix-matrix multiplication for a pair of large matrices (see Sec. 4 for several examples). For computing the large-scale matrix-matrix multiplication in parallel in serverless systems, we propose to use a straggler-resilient scheme called OverSketch from [32]. OverSketch does blocked partitioning of input matrices where each worker works on square blocks of dimension $b$. Hence, it is more communication efficient than existing coding-based straggler mitigation schemes that do naïve row-column partition of input matrices [48, 40]. We note that it is well known in HPC that blocked partitioning of input matrices can lead to communication-efficient methods for distributed multiplication [32, 49, 50].

OverSketch uses a sparse sketching matrix based on Count-Sketch [29]. It has similar computational efficiency and accuracy guarantees as that of the Count-Sketch, with two additional properties: it is amenable to distributed implementation; and it is resilient to stragglers. More specifically, the OverSketch matrix is constructed as follows [32].

Recall that the Hessian $\nabla^{2}f(\cdot)\in\mathbb{R}^{d\times d}$. First choose the desired sketch dimension $m$ (which depends on $d$), block-size $b$ (which depends on the memory of the workers), and straggler tolerance $\zeta>0$ (which depends on the distributed system). Then, define $N=m/b$ and $e=\zeta N$, for some constant $\zeta>0$. Here $\zeta$ is the fraction of stragglers that we want our algorithm to tolerate. Thus, $e$ is the maximum number of stragglers per $N+e$ workers that can be tolerated. The sketch $\mathbf{S}$ is then given by

[TABLE]

where $\mathbf{S}_{i}\in\mathbb{R}^{n\times b}$, for all $i\in[1,N+e]$, are i.i.d. Count-Sketch matrices333Each of the Count-Sketch matrices $\mathbf{S}_{i}$ is constructed (independently of others) as follows. First, for every row $j$, $j\in[n]$, of $\mathbf{S}_{i}$, independently choose a column $h(j)\in[b]$. Then, select a uniformly random element from $\{-1,+1\}$, denoted as $\sigma(i)$. Finally, set $\mathbf{S}_{i}(j,h(j))=\sigma(i)$ and set $\mathbf{S}_{i}(j,l)=0$ for all $l\neq h(j)$. (See [29, 32] for details.) with sketch dimension $b$. Note that $\mathbf{S}\in\mathbb{R}^{n\times(m+eb)}$, where $m=Nb$ is the required sketch dimension and $e$ is the over-provisioning parameter to provide resiliency against $e$ stragglers per $N+e$ workers. We leverage the straggler resiliency of OverSketch to obtain the sketched Hessian in a distributed straggler-resilient manner. An illustration of OverSketch is provided in Fig. 3; see Algorithm 2 for details.

Model update: Let $\hat{\mathbf{H}}_{t}=\mathbf{A}_{t}^{T}\mathbf{S}_{t}\mathbf{S}_{t}^{T}\mathbf{A}_{t}$, where $\mathbf{A}_{t}$ is the square root of the Hessian $\nabla^{2}f(\mathbf{w}_{t})$, and $\mathbf{S}_{t}$ is an independent realization of (4) at the $t$-th iteration. For strongly-convex functions, the update direction is $\mathbf{p}_{t}=-\hat{\mathbf{H}}_{t}^{-1}\nabla f(\mathbf{w}_{t})$. We use line-search to choose the step-size, that is, find

[TABLE]

for some constant $\beta\in(0,1/2]$. For weakly-convex functions, the update direction (inspired by Newton-MR [22]) is $\mathbf{p}_{t}=-\hat{\mathbf{H}}_{t}^{\dagger}\nabla f(\mathbf{w}_{t})$, where $\hat{\mathbf{H}}^{\dagger}_{t}$ is the Moore-Penrose inverse of $\hat{\mathbf{H}}_{t}$. To find the update $\mathbf{w}_{t+1}$, we find the right step-size $\alpha_{t}$ using line-search in (5), but with $f(\cdot)$ replaced by $||\nabla f(\cdot)||^{2}$ and $\nabla f(\mathbf{w}_{t})$ replaced by $2\hat{\mathbf{H}}_{t}\nabla f(\mathbf{w}_{t})$, according to the objective in $||\nabla f(\cdot)||^{2}$. More specifically, for some constant $\beta\in(0,1/2]$,

[TABLE]

Note that for OverSketched Newton, we use $\hat{\mathbf{H}}_{t}$ in the line-search since the exact Hessian is not available. The update in the $t$-th iteration in both cases is given by

[TABLE]

Note that (5) line-search can be solved approximately in single machine systems using Armijo backtracking line search [51, 16]. OverSketched Newton is concisely described in Algorithm 3. In Section 3.2, we describe how to implement distributed line-search in serverless systems when the data is stored in the cloud. Next, we prove convergence guarantees for OverSketched Newton that uses the sketch matrix in (4) and full gradient for approximate Hessian computation.

3.1 Convergence Guarantees

First, we focus our attention to strongly convex functions. We consider the following assumptions. We note that these assumptions are standard for analyzing approximate Newton methods, (e.g., see  [23, 19, 20].

Assumptions:

$f$ is twice-differentiable; 2. 2.

$f$ is $k$-strongly convex ($k>0$), that is,

[TABLE] 3. 3.

$f$ is $M$-smooth ($k\leq M<\infty$), that is,

[TABLE] 4. 4.

the Hessian is $L$-Lipschitz continuous, that is, for any $\boldsymbol{\Delta}\in\mathbb{R}^{d}$

[TABLE]

where $||\cdot||_{2}$ is the spectral norm for matrices.

We first prove the following “global” convergence guarantee which shows that OverSketched Newton would converge from any random initialization of $\mathbf{w}_{0}\in\mathbb{R}^{d}$ with high probability.

Theorem 3.1** (Global convergence for strongly-convex $f$).**

Consider Assumptions 1, 2, and 3 and step-size $\alpha_{t}$ given by Eq. (5). Let $\mathbf{w}^{*}$ be the optimal solution of (1). Let $\epsilon$ and $\mu$ be positive constants. Then, using the sketch in (4) with a sketch dimension $Nb+eb=\Omega(\frac{d^{1+\mu}}{\epsilon^{2}})$ and the number of column-blocks $N+e=\Theta_{\mu}(1/\epsilon)$, the updates for OverSketched Newton, for any $\mathbf{w}_{t}\in\mathbb{R}^{d}$, satisfy

[TABLE]

with probability at least $1-1/d^{\tau}$, where $\rho=\frac{2\alpha_{t}\beta k}{M(1+\epsilon)}$ and $\tau>0$ is a constant depending on $\mu$ and constants in $\Omega(\cdot)$ and $\Theta(\cdot)$. Moreover, $\alpha_{t}$ satisfies $\alpha_{t}\geq\frac{2(1-\beta)(1-\epsilon)k}{M}$.

Proof.

See Section 6.1. ∎

Theorem 3.1 guarantees the global convergence of OverSketched Newton starting with any initial estimate $\mathbf{w}_{0}\in\mathbb{R}^{d}$ to the optimal solution $\mathbf{w}^{*}$ with at least a linear rate.

Next, we can also prove an additional “local” convergence guarantee for OverSketched Newton, under the assumption that $\mathbf{w}_{0}$ is sufficiently close to $\mathbf{w}^{*}$.

Theorem 3.2** (Local convergence for strongly-convex $f$).**

Consider Assumptions 1, 2, and 4 and step-size $\alpha_{t}=1$. Let $\mathbf{w}^{*}$ be the optimal solution of (1) and $\gamma$ and $\beta$ be the minimum and maximum eigenvalues of $\nabla^{2}f(\mathbf{w}^{*})$, respectively. Let $\epsilon\in(0,\gamma/(8\beta)]$ and $\mu>0$. Then, using the sketch in (4) with a sketch dimension $Nb+eb=\Omega(\frac{d^{1+\mu}}{\epsilon^{2}})$ and the number of column-blocks $N+e=\Theta_{\mu}(1/\epsilon)$, the updates for OverSketched Newton, with initialization $\mathbf{w}_{0}$ such that $||\mathbf{w}_{0}-\mathbf{w}^{*}||_{2}\leq\frac{\gamma}{8L}$, follow

[TABLE]

with probability at least $1-T/d^{\tau}$, where $\tau>0$ is a constant depending on $\mu$ and constants in $\Omega(\cdot)$ and $\Theta(\cdot)$.

Proof.

See Section 6.2. ∎

Theorem 3.2 implies that the convergence is linear-quadratic in error $\Delta_{t}=\mathbf{w}_{t}-\mathbf{w}^{*}$. Initially, when $||\Delta_{t}||_{2}$ is large, the first term of the RHS will dominate and the convergence will be quadratic, that is, $||\Delta_{t+1}||_{2}\lesssim\frac{25L}{8\gamma}||\Delta_{t}||^{2}_{2}$. In later stages, when $||\mathbf{w}_{t}-\mathbf{w}^{*}||_{2}$ becomes sufficiently small, the second term of RHS will start to dominate and the convergence will be linear, that is, $||\Delta_{t+1}||_{2}\lesssim\frac{5\epsilon\beta}{\gamma}||\Delta_{t}||_{2}$. At this stage, the sketch dimension can be increased to reduce $\epsilon$ to diminish the effect of the linear term and improve the convergence rate in practice. Note that, for second order methods, the number of iterations $T$ is in the order of tens in general while the number of features $d$ is typically in thousands. Hence, the probability of failure is generally small (and can be made negligible by choosing $\tau$ appropriately).

Though the works [23, 19, 20, 47, 52] also prove convergence guarantees for approximate Hessian-based optimization, no convergence results exist for the OverSketch matrix in Eq. (4) to the best of our knowledge. OverSketch has many nice properties like sparsity, input obliviousness, and amenability to distributed implementation, and our convergence guarantees take into account the block-size $b$ (that captures the amount of communication at each worker) and the number of stragglers $e$, both of which are a property of the distributed system. On the other hand, algorithms in [23, 19, 20, 47, 52] are tailored to run on a single machine.

Next, we consider the case of weakly-convex functions. For this case, we consider two more assumptions on the Hessian matrix, similar to [22]. These assumptions are a relaxation of the strongly-convex case.

Assumptions:

There exists some $\eta>0$ such that, $\forall~{}\mathbf{w}\in\mathbb{R}^{d}$,

[TABLE]

This assumption establishes regularity on the pseudo-inverse of $\nabla^{2}f(\mathbf{x})$. It also implies that $||\nabla^{2}f(\mathbf{w})\mathbf{p}||\geq\eta||\mathbf{p}||~{}\forall~{}p\in\text{ Range}(\nabla^{2}f(\mathbf{w}))$, that is, the minimum ‘non-zero’ eigenvalue of $\nabla^{2}f(\mathbf{w})$ is lower bounded by $\eta$; just as in the $k$-strongly convex case, the smallest eigenvalue is greater than $k$. 2. 6.

Let $\mathbf{U}\in\mathbb{R}^{d\times d}$ be any arbitrary orthogonal basis for $\text{Range}(\nabla^{2}f(\mathbf{w}))$, there exists $0<\nu\leq 1$, such that,

[TABLE]

This assumption ensures that there is always a non-zero component of the gradient in the subspace spanned by the Hessian, and, thus, ensures that the model update $-\hat{\mathbf{H}}_{t}^{\dagger}\nabla f(\mathbf{w}_{t})$ will not be zero.

Note that the above assumptions are always satisfied by strongly-convex functions. Next, we prove global convergence of OverSketched Newton when the objective is weakly-convex.

Theorem 3.3** (Global convergence for weakly-convex $f$).**

Consider Assumptions 1,3,4,5 and 6 and step-size $\alpha_{t}$ given by Eq. (6). Let $\epsilon\in\left(0,\frac{(1-\beta)\nu\eta}{2M}\right]$ and $\mu>0$. Then, using an OverSketch matrix with a sketch dimension $Nb+eb=\Omega(\frac{d^{1+\mu}}{\epsilon^{2}})$ and the number of column-blocks $N+e=\Theta_{\mu}(1/\epsilon)$, the updates for OverSketched Newton, for any $\mathbf{w}_{t}\in\mathbb{R}^{d}$, satisfy

[TABLE]

with probability at least $1-1/d^{\tau}$, where \alpha=\frac{\eta}{2Q}\big{[}(1-\beta)\nu\eta-2\epsilon M\big{]},Q=(L||\nabla f(\mathbf{w}_{0})||+M^{2}), $\mathbf{w}_{0}$ is the initial iterate of the algorithm and $\tau>0$ is a constant depending on $\mu$ and constants in $\Omega(\cdot)$ and $\Theta(\cdot)$.

Proof.

See Section 6.3. ∎

Even though we present the above guarantees for the sketch matrix in Eq. (4), our analysis is valid for any sketch that satisfies the subspace embedding property (Lemma 6.1; see [29] for details on subspace embedding property of sketches). To the best of our knowledge, this is the first work to prove the convergence guarantees for weakly-convex functions when the Hessian is calculated approximately using sketching techniques. Later, authors in [53] extended the analysis to the case of general Hessian perturbations with additional assumptions on the type of perturbation.

3.2 Distributed Line Search

Here, we describe a line-search procedure for distributed serverless optimization, which is inspired by the line-search method from [24] for serverful systems. To solve for the step-size $\alpha_{t}$ as described in the optimization problem in (5), we set $\beta=0.1$ and choose a candidate set $\mathcal{S}=\{4^{0},4^{1},\cdots,4^{-5}\}$. After the master calculates the descent direction $\mathbf{p}_{t}$ in the $t$-th iteration, the $i$-th worker calculates $f_{i}(\mathbf{w}_{t}+\alpha\mathbf{p}_{t})$ for all values of $\alpha$ in the candidate set $\mathcal{S}$, where $f_{i}(\cdot)$ depends on the local data available at the $i$-th worker and $f(\cdot)=\sum_{i}f_{i}(\cdot)$444For the weakly-convex case, the workers calculate $\nabla f_{i}(\cdot)$ instead of $f_{i}(\cdot)$, and the master calculates $||\nabla f(\cdot)||^{2}$ instead of $f(\cdot)$..

The master then sums the results from workers to obtain $f(\mathbf{w}_{t}+\alpha\mathbf{p}_{t})$ for all values of $\alpha$ in $\mathcal{S}$ and finds the largest $\alpha$ that satisfies the Armijo condition in (5)555Note that codes can be used to mitigate stragglers during distributed line-search in a manner similar to the gradient computation phase.. Note that line search requires an additional round of communication where the master communicates $\mathbf{p}_{t}$ to the workers through cloud and the workers send back the function values $f_{i}(\cdot)$. Finally, the master finds the best step-size from set $\mathcal{S}$ and finds the model estimate $\mathbf{w}_{t+1}$.

4 OverSketched Newton on Serverless Systems: Examples

Here, we describe several examples where our general approach can be applied.

4.1 Logistic Regression using OverSketched Newton

The optimization problem for supervised learning using Logistic Regression takes the form

[TABLE]

Here, $\mathbf{x}_{1},\cdots,\mathbf{x}_{n}\in\mathbb{R}^{d\times 1}$ and $y_{1},\cdots,y_{n}\in\mathbb{R}$ are training sample vectors and labels, respectively. The goal is to learn the feature vector $\mathbf{w}^{*}\in\mathbb{R}^{d\times 1}$. Let $\mathbf{X}=[\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n}]\in\mathbb{R}^{d\times n}$ and $\mathbf{y}=[y_{1},\cdots,y_{n}]\in\mathbb{R}^{n\times 1}$ be the example and label matrices, respectively. The gradient for the problem in (7) is given by

[TABLE]

Calculation of $\nabla f(\mathbf{w})$ involves two matrix-vector products, $\boldsymbol{\alpha}=\mathbf{X}^{T}\mathbf{w}$ and $\nabla f(\mathbf{w})=\frac{1}{n}\mathbf{X}\boldsymbol{\beta}+\lambda\mathbf{w}$, where $\beta_{i}=\frac{-y_{i}}{1+e^{y_{i}\alpha_{i}}}~{}\forall~{}i\in[1,\cdots,n]$. When the example matrix is large, these matrix-vector products are performed distributedly using codes. Faster convergence is obtained by second-order methods which will additionally compute the Hessian $\mathbf{H}=\frac{1}{n}\mathbf{X}\boldsymbol{\Lambda}\mathbf{X}^{T}+\lambda\mathbf{I}_{d}$, where $\boldsymbol{\Lambda}$ is a diagonal matrix with entries given by $\Lambda(i,i)=\frac{e^{y_{i}\alpha_{i}}}{(1+e^{y_{i}\alpha_{i}})^{2}}$. The product $\mathbf{X}\boldsymbol{\Lambda}\mathbf{X}^{T}$ is computed approximately in a distributed straggler-resilient manner using the sketch matrix in (4). Using the result of distributed multiplication, the Hessian matrix H is calculated at the master and the model is updated as $\mathbf{w}_{t+1}=\mathbf{w}_{t}-\textbf{H}^{-1}\nabla f(\mathbf{w}_{t})$. In practice, efficient algorithm like conjugate gradient, that provide a good estimate in a small number of iterations, can be used locally at the master to solve for $\mathbf{w}_{t+1}$ [54].666Note that here we have assumed that the number of features is small enough to perform the model update locally at the master. This is not necessary, and straggler resilient schemes, such as in [35], can be used to perform distributed conjugate gradient in serverless systems.

We provide a detailed description of OverSketched Newton for large-scale logistic regression for serverless systems in Algorithm 4. Steps 4, 8, and 14 of the algorithm are computed in parallel on AWS Lambda. All other steps are simple vector operations that can be performed locally at the master, for instance, the user’s laptop. Steps 4 and 8 are executed in a straggler-resilient fashion using the coding scheme in [34], as illustrated in Fig. 1 and described in detail in Algorithm 1.

We use the coding scheme in [34] since the encoding can be implemented in parallel and requires less communication per worker compared to the other schemes, for example schemes in [33, 40], that use Maximum Distance Separable (MDS) codes. Moreover, the decoding scheme takes linear time and is applicable on real-valued matrices. Note that since the example matrix $\mathbf{X}$ is constant in this example, the encoding of $\mathbf{X}$ is done only once before starting the optimization algorithm. Thus, the encoding cost can be amortized over iterations. Moreover, decoding over the resultant product vector requires negligible time and space, even when $n$ is scaling into the millions.

The same is, however, not true for the matrix multiplication for Hessian calculation (step 14 of Algorithm 4), as the matrix A changes in each iteration, thus encoding costs will be incurred in every iteration if error-correcting codes are used. Moreover, encoding and decoding a huge matrix stored in the cloud incurs heavy communication cost and becomes prohibitive. Motivated by this, we use OverSketch in step 14, as described in Algorithm 2, to calculate an approximate matrix multiplication, and hence the Hessian, efficiently in serverless systems with inbuilt straggler resiliency.777We also evaluate the exact Hessian-based algorithm with speculative execution, i.e., recomputing the straggling jobs, and compare it with OverSketched Newton in Sec. 5.

4.2 Softmax Regression using OverSketched Newton

We take unregulairzed softmax regression as an illustrative example for the weakly convex case. The goal is to find the weight matrix $\mathbf{W}=\left[\mathbf{w}_{1},\cdots,\mathbf{w}_{K}\right]$ that fit the training data $\mathbf{X}\in\mathbb{R}^{d\times N}$ and $\mathbf{y}\in\mathbb{R}^{K\times N}$. Here $\mathbf{w}_{i}\in\mathbb{R}^{d}$ represesents the weight vector for the $k$-th class for all $i\in[1,K]$ and $K$ is the total number of classes. Hence, the resultant feature dimension for softmax regression is $dK$. The optimization problem is of the form

[TABLE]

The gradient vector for the $i$-th class is given by

[TABLE]

which can be written as matrix products $\boldsymbol{\alpha_{i}}=\mathbf{X}^{T}\mathbf{w}_{i}$ and $\nabla f_{i}(\mathbf{W})=\mathbf{X}\boldsymbol{\beta}_{i}$, where the entries of $\boldsymbol{\beta}_{i}\in\mathbb{R}^{N}$ are given by $\beta_{in}=\left(\frac{\exp(\alpha_{in})}{\sum_{l=1}^{K}\exp(\alpha_{ln})}-y_{in}\right)$. Thus, the full gradient matrix is given by $\nabla f(\mathbf{W})=\mathbf{X}\boldsymbol{\beta}$ where the entries of $\boldsymbol{\beta}\in\mathbb{R}^{N\times K}$ are dependent on $\boldsymbol{\alpha}\in\mathbb{R}^{N\times K}$ as above and the matrix $\boldsymbol{\alpha}$ is given by $\boldsymbol{\alpha}=\mathbf{X}^{T}\mathbf{W}$. We assume that the number of classes $K$ is small enough such that tall matrices $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are small enough for the master to do local calculations on them.

Since the effective number of features is $d\times K$, the Hessian matrix is of dimension $dK\times dK$. The $(i,j)$-th component of the Hessian, say $\mathbf{H}_{ij}$, is

[TABLE]

where $\mathbf{Z}_{ij}\in\mathbb{R}^{N\times N}$ is a diagonal matrix whose $n$-th diagonal entry is

[TABLE]

where $\mathbb{I}(\cdot)$ is the indicator function and $\boldsymbol{\alpha}=\mathbf{X}\mathbf{W}$ was defined above. The full Hessian matrix is obtained by putting together all such $\mathbf{H}_{ij}$’s in a $dK\times dK$ matrix and can be expressed in a matrix-matrix multiplication form as

[TABLE]

where $\bar{\mathbf{X}}\in\mathbb{R}^{dK\times NK}$ is a block diagonal matrix that contains $\mathbf{X}$ in the diagonal blocks and $\bar{\mathbf{Z}}\in\mathbb{R}^{NK\times NK}$ is formed by stacking all the $\mathbf{Z}_{ij}$’s for $i,j\in[1,K]$. In OverSketched Newton, we compute this multiplication using sketching in serverless systems for efficiency and resiliency to stragglers. Assuming $d\times K$ is small enough, the master can then calculate the update $\mathbf{p}_{t}$ using efficient algorithms such the minimum-residual method [55, 22].

4.3 Other Example Problems

In this section, we describe several other commonly encountered optimization problems that can be solved using OverSketched Newton.

Ridge Regularized Linear Regression: The optimization problem is

[TABLE]

The gradient in this case can be written as $\frac{1}{n}\mathbf{X}(\boldsymbol{\beta}-\mathbf{y})+\lambda\mathbf{w}$, where $\boldsymbol{\beta}=\mathbf{X}^{T}\mathbf{w}$, where the training matrix $\mathbf{X}$ and label vector $\mathbf{y}$ were defined previously. The Hessian is given by $\nabla^{2}f(\mathbf{w})=\mathbf{X}\mathbf{X}^{T}+\lambda\mathbf{I}$. For $n\gg d$, this can be computed approximately using the sketch matrix in (4).

Linear programming via interior point methods: The following linear program can be solved using OverSketched Newton

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{m\times 1},\mathbf{c}\in\mathbb{R}^{m\times 1},\mathbf{b}\in\mathbb{R}^{n\times 1}$ and $\mathbf{A}\in\mathbb{R}^{n\times m}$ is the constraint matrix with $n>m$. In algorithms based on interior point methods, the following sequence of problems using Newton’s method

[TABLE]

where $\mathbf{a}_{i}$ is the $i$-th row of $\mathbf{A}$, $\tau$ is increased geometrically such that when $\tau$ is very large, the logarithmic term does not affect the objective value and serves its purpose of keeping all intermediates solution inside the constraint region. The update in the $t$-th iteration is given by $\mathbf{x}_{t+1}=\mathbf{x}_{t}-(\nabla^{2}f(\mathbf{x}_{t}))^{-1}\nabla f(\mathbf{x}_{t})$, where $\mathbf{x}_{t}$ is the estimate of the solution in the $t$-th iteration. The gradient can be written as $\nabla f(\mathbf{x})=\tau\textbf{c}+\mathbf{A}^{T}\boldsymbol{\beta}$ where $\beta_{i}=1/(b_{i}-\alpha_{i})$ and $\boldsymbol{\alpha}=\mathbf{A}\mathbf{x}$.

The Hessian for the objective in (21) is given by

[TABLE]

The square root of the Hessian is given by $\nabla^{2}f(\mathbf{x})^{1/2}=\text{diag}\frac{1}{|bi-\alpha_{i}|}\mathbf{A}$. The computation of Hessian requires $O(nm^{2})$ time and is the bottleneck in each iteration. Thus, we can use sketching to mitigate stragglers while evaluating the Hessian efficiently, i.e. $\nabla^{2}f(\mathbf{x})\approx(\mathbf{S}\nabla^{2}f(\mathbf{x})^{1/2})^{T}\times(\mathbf{S}\nabla^{2}f(\mathbf{x})^{1/2})$, where $\mathbf{S}$ is the OverSketch matrix defined in (4).

Lasso Regularized Linear Regression: The optimization problem takes the following form

[TABLE]

where $\mathbf{X}\in\mathbb{R}^{n\times d}$ is the measurement matrix, the vector $\mathbf{y}\in\mathbb{R}^{n}$ contains the measurements, $\lambda\geq 0$ and $d\gg n$. To solve (23), we consider its dual variation

[TABLE]

which is amenable to interior point methods and can be solved by optimizing the following sequence of problems where $\tau$ is increased geometrically

[TABLE]

where $\mathbf{x}_{j}$ is the $j$-th column of $\mathbf{X}$. The gradient can be expressed in few matrix-vector multiplications as $\nabla f(\mathbf{z})=\tau(\mathbf{z}-\mathbf{y})+\mathbf{X}(\boldsymbol{\beta}-\boldsymbol{\gamma}),$ where $\beta_{i}=1/(\lambda-\alpha_{i})$, $\gamma_{i}=1/(\lambda+\alpha_{i})$, and $\boldsymbol{\alpha}=\mathbf{X}^{T}\mathbf{z}$. Similarly, the Hessian can be written as $\nabla^{2}f(\mathbf{z})=\tau\mathbf{I}+\mathbf{X}\boldsymbol{\Lambda}\mathbf{X}^{T}$, where $\boldsymbol{\Lambda}$ is a diagonal matrix whose entries are given by $\Lambda_{ii}=1/(\lambda-\alpha_{i})^{2}+1/(\lambda+\alpha_{i})^{2}~{}\forall~{}i\in[1,n]$.

Other common problems where OverSketched Newton is applicable include Linear Regression, Support Vector Machines (SVMs), Semidefinite programs, etc.

5 Experimental Results

In this section, we evaluate OverSketched Newton on AWS Lambda using real-world and synthetic datasets, and we compare it with state-of-the-art distributed optimization algorithms888A working implementation of OverSketched Newton is available at https://github.com/vvipgupta/OverSketchedNewton. We use the serverless computing framework, Pywren [3]. Our experiments are focused on logistic and softmax regression, which are popular supervised learning problems, but they can be reproduced for other problems described in Section 4. We present experiments on the following datasets:

For comparison of OverSketched Newton with existing distributed optimization schemes, we choose recently-proposed Globally Improved Approximate Newton Direction (GIANT) [24]. The reason is that GIANT boasts a better convergence rate than many existing distributed second-order methods for linear and logistic regression, when $n\gg d$. In GIANT, and other similar distributed second-order algorithms, the training data is evenly divided among workers, and the algorithms proceed in two stages. First, the workers compute partial gradients using local training data, which is then aggregated by the master to compute the exact gradient. Second, the workers receive the full gradient to calculate their local second-order estimate, which is then averaged by the master. An illustration is shown in Fig. 4.

For straggler mitigation in such serverful systems based algorithms, [37] proposes a scheme for coding gradient updates called gradient coding, where the data at each worker is repeated multiple times to compute redundant copies of the gradient. See Figure 5(b) for illustration. Figure 5(a) illustrates the scheme that waits for all workers and Figure 5(c) illustrates the ignoring stragglers approach. We use the three schemes for dealing with stragglers illustrated in Figure 5 during the two stages of GIANT, and we compare their convergence with OverSketched Newton. We further evaluate and compare the convergence exact Newton’s method (employed with speculative execution, that is, reassigning and recomputing the work for straggling workers).

5.1 Comparisons with Existing Second-Order Methods on AWS Lambda

In Figure 6, we present our results on a synthetic dataset with $n=300,000$ and $d=3000$ for logistic regression on AWS Lambda. Each column $\mathbf{x}_{i}\in\mathbb{R}^{d}$, for all $i\in[1,n]$, is sampled uniformly randomly from the cube $[-1,1]^{d}$. The labels $y_{i}$ are sampled from the logistic model, that is, $\mathbb{P}[y_{i}=1]=1/(1+exp(\mathbf{x}_{i}\mathbf{w}+b))$, where the weight vector $\mathbf{w}$ and bias $b$ are generated randomly from the normal distribution. ht vector $\mathbf{w}$ and bias $b$ are generated randomly from the normal distribution.

The orange, blue and red curves demonstrate the convergence for GIANT with the full gradient (that waits for all the workers), gradient coding and mini-batch gradient (that ignores the stragglers while calculating gradient and second-order updates) schemes, respectively. The purple and green curves depict the convergence for the exact Newton’s method and OverSketched Newton, respectively. The gradient coding scheme is applied for one straggler, that is the data is repeated twice at each worker. We use 60 Lambda workers for executing GIANT in parallel. Similarly, for Newton’s method, we use 60 workers for matrix-vector multiplication in steps 4 and 8 of Algorithm 4, $3600$ workers for exact Hessian computation and $600$ workers for sketched Hessian computation with a sketch dimension of $10d=30,000$ in step 14 of Algorithm 4. In all cases, unit step-size was used to update the model999Line-search in Section 3 was mainly introduced to prove theoretical guarantees. In our experiments, we observe that constant step-size works well for OverSketched Newton.

Remark 1**.**

In our experiments, we choose the number of workers in such a way that each worker receives approximately the same amount of data to work with, regardless of the algorithm. This is motivated by the fact that the memory at each worker is the bottleneck in serverless systems (e.g., in AWS Lambda, the memory at each worker can be as low as 128 MB). Note that this is unlike serverful/HPC systems, where the number of workers is the bottleneck.

An important point to note from Fig. 6 is that the uncoded scheme (that is, the one that waits for all stragglers) has the worst performance. The implication is that good straggler/fault mitigation algorithms are essential for computing in the serverless setting. Secondly, the mini-batch scheme outperforms the gradient coding scheme by $25\%$. This is because gradient coding requires additional communication of data to serverless workers (twice when coding for one straggler, see [37] for details) at each invocation to AWS Lambda. On the other hand, the exact Newton’s method converges much faster than GIANT, even though it requires more time per iteration.

The number of iterations needed for convergence for OverSketched Newton and exact Newton (that exactly computes the Hessian) is similar, but OverSketched Newton converges in almost half the time due to an efficient computation of (approximate) Hessian (which is the computational bottleneck and thus reduces time per iteration).

5.1.1 Logistic Regression on EPSILON, WEBPAGE and a9a Datasets

In Figure 7, we repeat the above experiment with EPSILON classification dataset obtained from [36], with $n=0.4~{}million$ and $d=2000$. We plot training and testing errors for logistic regression for the schemes described in the previous section. Here, we use $100$ workers for GIANT, and $100$ workers for matrix-vector multiplications for gradient calculation in OverSketched Newton. We use gradient coding designed for three stragglers in GIANT. This scheme performs worse than uncoded GIANT that waits for all the stragglers due to the repetition of training data at workers. Hence, one can conclude that the communication costs dominate the straggling costs. In fact, it can be observed that the mini-batch gradient scheme that ignores the stragglers outperforms the gradient coding and uncoded schemes for GIANT.

During exact Hessian computation, we use $10,000$ serverless workers with speculative execution to mitigate stragglers (i.e., recomputing the straggling jobs) compared to OverSketched Newton that uses $1500$ workers with a sketch dimension of $15d=30,000$. OverSketched Newton requires a significantly smaller number of workers, as once the square root of Hessian is sketched in a distributed fashion, it can be copied into local memory of the master due to dimension reduction, and the Hessian can be calculated locally. Testing error follows training error closely, and important conclusions remain the same as in Figure 6. OverSketched Newton outperforms GIANT and exact Newton-based optimization by at least $46\%$ in terms of running time.

We repeated the above experiments for classification on the WEBPAGE ($n=49,749$ and $d=300$) and a9a ($n=32,561$ and $d=123$) datasets [36]. For both datasets, we used 30 workers for each iteration in GIANT and any matrix-vector multiplications. Exact hessian calculation invokes $900$ workers as opposed to $300$ workers for OverSketched Newton, where the sketch dimension was $10d=3000$. The results for training loss on logistic regression are shown in Figure 8. Testing error closely follows the training error in both cases. OverSketched Newton outperforms exact Newton and GIANT by at least $\sim 25\%$ and $\sim 75\%$, respectively, which is similar to the trends witnessed heretofore.

Remark 2**.**

Note that conventional distributed second-order methods for serverful systems—which distribute training examples evenly across workers (such as [42, 44, 43, 24, 45, 46])—typically find a “localized approximation” (localized to each machine) of second-order update at each worker and then aggregate it. OverSketched Newton, on the other hand, uses the massive storage and compute power in serverless systems to find a more “globalized approximation” (globalized in the sense of across machine). Thus, it performs better in practice.

5.2 Softmax Regression on EMIST

In Fig. 12, we solve unregularized softmax regression, which is weakly convex (see Sec. 4.2 for details). We use the Extended MNIST (EMNIST) dataset [56] with $N=240,000$ training examples, $d=784$ features and $K=10$ classes. Note that GIANT cannot be applied here as the objective function is not strongly convex. We compare the convergence rate of OverSketched Newton, exact Hessian and gradient descent based schemes.

For gradient computation in all three schemes, we use 60 workers. However, exact Newton scheme requires 3600 workers to calculate the $dK\times dK$ Hessian and recomputes the straggling jobs, while OverSketched Newton requires only 360 workers to calculate the sketch in parallel with sketch dimension $6dK=47,040$. The approximate Hessian is then computed locally at the master using its sketched square root, where the sketch dimension is $6dK=47,040$. The step-size is fixed and is determined by hyperparamter tuning before the start of the algorithm. Even for the weakly-convex case, second-order methods tend to perform better. Moreover, the runtime of OverSketched Newton outperforms both gradient descent and Exact Newton based methods by $\sim 75\%$ and $\sim 50\%$, respectively.

5.3 Coded computing versus Speculative Execution

In Figure 12, we compare the effect of straggler mitigation schemes, namely speculative execution, that is, restarting the jobs with straggling workers, and coded computing on the convergence rate during training and testing. We regard OverSketch based matrix multiplication as a coding scheme in which some redundancy is introduced during “over” sketching for matrix multiplication. There are four different cases, corresponding to gradient and hessian calculation using either speculative execution or coded computing. For speculative execution, we wait for at least $90\%$ of the workers to return (this works well as the number of stragglers is generally less than $10\%$) and restart the jobs that did not return till this point.

For both exact Hessian and OverSketched Newton, using codes for distributed gradient computation outperforms speculative execution based straggler mitigation. Moreover, computing the Hessian using OverSketch is significantly better than exact computation in terms of running time as calculating the Hessian is the computational bottleneck in each iteration.

5.4 Comparison with First-Order Methods on AWS Lambda

In Figure 12, we compare gradient descent and Nesterov Accelerated Gradient (NAG) (while ignoring the stragglers) with OverSketched Newton for logistic regression on EPSILON dataset. We observed that for first-order methods, there is only a slight difference in convergence for a mini-batch gradient when the batch size is $95\%$. Hence, for gradient descent and NAG, we use 100 workers in each iteration while ignoring the stragglers.101010We note that stochastic methods such as SGD perform worse that gradient descent since their update quality is poor, requiring more iterations (hence, more communication) to converge while not using the massive compute power of serverless. For example, 20% minibatch SGD in the setup of Fig. 12 requires $1.9\times$ more time than gradient descent with same number of workers. These first-order methods were given the additional advantage of backtracking line-search, which determined the optimal amount to move in given a descent direction.111111We remark that backtracking line-search required $\sim 13\%$ of the total time for NAG. Hence, as can be seen from Fig. 12, any well-tuned step-size method would still be significantly slower than OverSketched Newton. Overall, OverSketched Newton with unit step-size significantly outperforms gradient descent and NAG with backtracking line-search.

5.5 Comparison with Serverful Optimization

In Fig. 12, we compare OverSketched Newton on AWS Lambda with existing distributed optimization algorithm GIANT in serverful systems (AWS EC2). The results are plotted on synthetically generated data for logistic regression. For serverful programming, we use Message Passing Interface (MPI) with one c3.8xlarge master and $60$ t2.medium workers in AWS EC2. In [4], the authors observed that many large-scale linear algebra operations on serverless systems take at least $30\%$ more time compared to MPI-based computation on serverful systems. However, as shown in Fig. 12, we observe a slightly surprising trend that OverSketched Newton outperforms MPI-based optimization (that uses existing state-of-the-art optimization algorithm). This is because OverSketched Newton exploits the flexibility and massive scale at disposal in serverless, and thus produces a better approximation of the second-order update than GIANT.121212We do not compare with exact Newton in serverful sytems since the data is large and stored in the cloud. Computing the exact Hessian would require a large number of workers (e.g., we use 10,000 workers for exact Newton in EPSILON dataset) which is infeasible in existing serverful systems.

6 Proofs

To complete the proofs in this section, we will need the following lemma.

Lemma 6.1**.**

Let $\hat{\mathbf{H}}_{t}=\mathbf{A}_{t}^{T}\mathbf{S}_{t}\mathbf{S}_{t}^{T}\mathbf{A}_{t}$ where $\mathbf{S}_{t}$ is the sparse sketch matrix in (4) with sketch dimension $m=\Omega(d^{1+\mu}/\epsilon^{2})$ and $N=\Theta_{\mu}(1/\epsilon)$. Then, the following holds

[TABLE]

with probability at least $1-\frac{1}{d^{\tau}}$, where $\tau>0$ is a constant depending on $\mu$ and the constants in $\Theta(\cdot)$ and $\Omega(\cdot)$, and $\lambda_{\max}(\cdot)$ and $\lambda_{\min}(\cdot)$ denote the maximum and minimum eigenvalues, respectively. In general,

[TABLE]

where $\lambda_{i}(\cdot)$ is the $i$-th eigenvalue.

Proof.

We note than $N$ is the number of non-zero elements per row in the sketch $\mathbf{S}_{t}$ in (4) after ignoring stragglers. We use Theorem 8 in [57] to bound the singular values for the sparse sketch $\mathbf{S}_{t}$ in (4) with sketch dimension $m=\Omega(d^{1+\mu}/\epsilon^{2})$ and $N=\Theta(1/\epsilon)$. It says that $\mathbb{P}(\forall~{}\mathbf{x}~{}\in~{}\mathbb{R}^{n},||\mathbf{S}_{t}\mathbf{x}||_{2}\in(1\pm\epsilon/3)||\mathbf{x}||_{2})>1-1/d^{\tau}$, where $\tau>0$ depends on $\mu$ and the constants in $\Theta(\cdot)$ and $\Omega(\cdot)$. Thus, $||\mathbf{S}_{t}\mathbf{x}||_{2}\in(1\pm\epsilon/3)||\mathbf{x}||_{2}$, which implies that

[TABLE]

with probability at least $1-1/d^{\tau}$. For $\epsilon\leq 1/2$, this leads to the following inequality

[TABLE]

with probability at least $1-1/d^{\tau}$. Also, since $(1-\epsilon)\mathbf{x}^{T}\mathbf{x}\leq\mathbf{x}^{T}\mathbf{S}_{t}\mathbf{S}_{t}^{T}\mathbf{x}\leq(1+\epsilon)\mathbf{x}^{T}\mathbf{x}~{}\forall~{}x\in\mathbb{R}^{n}$ by the inequality above, replacing $\mathbf{x}$ by $\mathbf{A}\mathbf{y}$, we get

[TABLE]

with probability at least $1-1/d^{\tau}$. Let $\mathbf{y}_{1}$ be the unit norm eigenvector corresponding to the minimum eigenvalue of $\hat{\mathbf{H}}_{t}=\mathbf{A}^{T}_{t}\mathbf{S}_{t}\mathbf{S}_{t}^{T}\mathbf{A}_{t}$. Since the above inequality is true for all $\mathbf{y}$, we have

[TABLE]

with probability at least $1-1/d^{\tau}$. Along similar lines, we can prove that $\lambda_{\max}(\hat{\mathbf{H}}_{t})\leq(1-\epsilon)\lambda_{\max}(\nabla^{2}f(\mathbf{w}_{t}))$ with probability at least $1-1/d^{\tau}$ using the right hand inequality in (27). Together, these prove the first result.

In general, Eq. (26) implies that the eigenvalues of $(\mathbf{S}_{t}\mathbf{S}_{t}^{T}-\mathbf{I})$ are in the set $[-\epsilon,\epsilon]$. Thus, all the eigenvalues of $\mathbf{A}^{T}_{t}(\mathbf{S}_{t}\mathbf{S}_{t}^{T}-\mathbf{I})\mathbf{A}_{t}$ are in the set $[-\epsilon\lambda_{\max}(\nabla^{2}f(\mathbf{w}_{t})),\epsilon\lambda_{\max}(\nabla^{2}f(\mathbf{w}_{t}))]$ Also, we can write

[TABLE]

Now, applying Weyl’s inequality (see [58], Section 1.3) on symmetric matrices $\hat{\mathbf{H}}_{t}=\mathbf{A}_{t}^{T}\mathbf{S}_{t}\mathbf{S}_{t}^{T}\mathbf{A}_{t}$, $\nabla^{2}f(\mathbf{w}_{t})=\mathbf{A}_{t}^{T}\mathbf{A}_{t}$ and $\mathbf{A}_{t}^{T}(\mathbf{S}_{t}\mathbf{S}_{t}^{T}-\mathbf{I})\mathbf{A}_{t}$, we get

[TABLE]

which proves the second result.

∎

6.1 Proof of Theorem 3.1

Let’s define $\mathbf{w}_{\tau}=\mathbf{w}_{t}+\tau\mathbf{p}_{t}$, where the descent direction $\mathbf{p}_{t}$ is given by $\mathbf{p}_{t}=-\hat{\mathbf{H}}_{t}^{-1}\nabla f(\mathbf{w}_{t})$. Also, from Lemma 6.1, we have

[TABLE]

with probability at least $1-1/d^{\tau}.$ Using the above inequalities and the fact that $f(\cdot)$ is $k$-strongly convex and $M$-smooth, we get

[TABLE]

with probability at least $1-1/d^{\tau}.$

Next, we show that there exists an $\alpha>0$ such that the Armijo line search condition in (5) is satisfied. From the smoothness of $f(\cdot)$, we get (see [59], Theorem 2.1.5)

[TABLE]

Now, for $\mathbf{w}_{\alpha}$ to satisfy the Armijo rule, $\alpha$ should satisfy

[TABLE]

where the last inequality follows from the definition of $\mathbf{p}_{t}$. Now, using the lower bound from (28), $\mathbf{w}_{\alpha}$ satisfies Armijo rule for all

[TABLE]

Hence, we can always find an $\alpha_{t}\geq\frac{2(1-\beta)(1-\epsilon)k}{M}$ using backtracking line search such that $\mathbf{w}_{t+1}$ satisfies the Armijo condition, that is,

[TABLE]

which in turn implies

[TABLE]

with probability at least $1-1/d^{\tau}$. Here the last inequality follows from the bound in (28). Moreover, $k$-strong convexity of $f(\cdot)$ implies (see [59], Theorem 2.1.10)

[TABLE]

Using the inequality from (29) in the above inequality, we get

[TABLE]

where $\rho=\frac{2\alpha_{t}\beta k}{M(1+\epsilon)}$. Rearranging, we get

[TABLE]

with probability at least $1-1/d^{\tau}$, which proves the desired result.

6.2 Proof of Theorem 3.2

According to OverSketched Newton update, $\mathbf{w}_{t+1}$ is obtained by solving

[TABLE]

Thus, we have, for any $\mathbf{w}\in\mathbb{R}^{d}$,

[TABLE]

Substituting $\mathbf{w}$ by $\mathbf{w}^{*}$ in the above expression and calling $\Delta_{t}=\mathbf{w}^{*}-\mathbf{w}_{t}$, we get

[TABLE]

Now, due to the optimality of $\mathbf{w}^{*}$, we have $\nabla f(\mathbf{w}^{*})^{T}\Delta_{t+1}\geq 0$. Hence, we can write

[TABLE]

Next, substituting \nabla f(\mathbf{w}_{t})-\nabla f(\mathbf{w}^{*})=\Big{(}\int_{0}^{1}\nabla^{2}f(\mathbf{w}^{*}+p(\mathbf{w}_{t}-\mathbf{w}^{*}))dp\Big{)}(\mathbf{w}_{t}-\mathbf{w}^{*}) in the above inequality, we get

[TABLE]

Using Cauchy-Schwartz inequality in the LHS above, we get

[TABLE]

Now, using the $L$-Lipschitzness of $\nabla^{2}f(\cdot)$ in the inequality above, we get

[TABLE]

Note that for the positive definite matrix $\nabla^{2}f(\mathbf{w}_{t})=\mathbf{A}_{t}^{T}\mathbf{A}_{t}$, we have $||\mathbf{A}_{t}||^{2}_{2}=||\nabla^{2}f(\mathbf{w}_{t})||_{2}$. Moreover,

[TABLE]

Now, using Equation 26 from the proof of Lemma 6.1, we get $||\mathbf{S}_{t}\mathbf{S}_{t}^{T}-\mathbf{I}||_{2}=\lambda_{\max}(\mathbf{S}_{t}\mathbf{S}_{t}^{T}-\mathbf{I})\leq\epsilon$. Using this to bound the RHS of (6.2), we have, with probability at least $1-1/d^{\tau}$,

[TABLE]

where the last inequality follows from $\hat{\mathbf{H}}_{t}=\mathbf{A}_{t}^{T}\mathbf{S}_{t}^{T}\mathbf{S}_{t}\mathbf{A}_{t}$. Now, since the sketch dimension $m=\Omega(d^{1+\mu}/\epsilon^{2})$, using Eq. (26) from the proof of Lemma 1 in above inequality, we get, with probability at least $1-1/d^{\tau}$,

[TABLE]

Now, since $\gamma$ and $\beta$ are the minimum and maximum eigenvalues of $\nabla^{2}f(\mathbf{w}^{*})$, we get

[TABLE]

by the Lipschitzness of $\nabla^{2}f(\mathbf{w})$, that is, $|\Delta_{t+1}^{T}(\nabla^{2}f(\mathbf{w}_{t})-\nabla^{2}f(\mathbf{w}^{*}))\Delta_{t+1}|\leq L||\Delta_{t}||_{2}||\Delta||_{t+1}^{2}$. Rearranging for $\epsilon\leq\gamma/(8\beta)<1/2$, we get

[TABLE]

with probability at least $1-1/d^{\tau}$.

Let $\xi_{T}$ be the event that the above inequality (in (31)) is true for $t=0,1,\cdots,T$. Thus,

[TABLE]

where the second inequality follows from Bernoulli’s inequality. Next, assuming that the event $\xi_{T}$ holds, we prove that $||\Delta_{t}||_{2}\leq\gamma/5L$ using induction. We can verify the base case using the initialization condition, i.e. $||\Delta_{0}||_{2}\leq\gamma/8L$. Now, assuming that $||\Delta_{t-1}||_{2}\leq\gamma/5L$ and using it in the inequality (31), we get

[TABLE]

where the last inequality uses the fact that $\epsilon\leq\gamma/(8\beta)$. Thus, by induction,

[TABLE]

Using this in (31), we get the desired result, that is,

[TABLE]

with probability at least $1-T/d^{\tau}$.

6.3 Proof of Theorem 3.3

Let us define a few short notations for convenience. Say $\mathbf{g}_{t}=\nabla f(\mathbf{w}_{t})$ and $\mathbf{H}_{t}=\nabla^{2}f(\mathbf{w}_{t})=\mathbf{A}_{t}^{T}\mathbf{A}_{t}$, and we know that $\hat{\mathbf{H}}_{t}=\mathbf{A}^{T}_{t}\mathbf{S}_{t}\mathbf{S}^{T}_{t}\mathbf{A}_{t}$. Moreover, all the results with approximate Hessian $\hat{\mathbf{H}}_{t}$ hold with probability $1-1/d^{\tau}$. We skip its mention in most of the proof for brevity. The following lemmas will assist us in the proof.

Lemma 6.2**.**

$M$-smoothness of $f(\cdot)$ and $L$-Lipchitzness of $\nabla^{2}f(\cdot)$ imply

[TABLE]

for all $\mathbf{x}\in\mathbb{R}^{d}$, $Q=(L\delta+M^{2})$, where $\mathbf{y}\in\mathcal{Y}$, where $\mathcal{Y}=\{\mathbf{y}\in\mathbb{R}^{d}|~{}||\nabla f(\mathbf{y})||\leq\delta\}$ and $\delta>0$ is some constant.

Proof.

We have

[TABLE]

By applying triangle inequality and Cauchy-Schwarz to above equation, we get

[TABLE]

From the smoothness of $f(\cdot)$, that is, Lipshitzness of gradient, we get $||\nabla^{2}f(\mathbf{x})||_{2}\leq M~{}\forall~{}x\in\mathbb{R}^{d}$. Additionally, using Lipshitzness of Hessian, we get

[TABLE]

for $\mathbf{y}\in\mathcal{Y}$. This proves the desired result. ∎

Lemma 6.3**.**

Let $\mathbf{A}^{T}=\mathbf{U}\sqrt{\Sigma}\mathbf{V}^{T}$ and $\mathbf{A}^{T}\mathbf{S}_{t}=\hat{\mathbf{U}}\sqrt{\hat{\Sigma}}\hat{\mathbf{V}}^{T}$ be the truncated Singular Value Decompositions (SVD) of $\mathbf{A}^{T}$ and $\mathbf{A}^{T}\mathbf{S}_{t}$, respectively. Thus, $\mathbf{H}_{t}=\mathbf{U}\Sigma\mathbf{U}^{T}$ and $\hat{\mathbf{H}}_{t}=\hat{\mathbf{U}}\hat{\Sigma}\hat{\mathbf{U}}^{T}$. Then, for all $\mathbf{g}\in\mathbb{R}^{d}$, we have

[TABLE]

where $\eta$ is defined in Assumption (5).

Proof.

For all $\mathbf{g}\in\mathbb{R}^{d}$, using the fact that $\mathbf{A}=\bf{V}\sqrt{\Sigma}\mathbf{U}^{T}$, we get

[TABLE]

where the last inequality uses Assumption (5). In a similar fashion, we can obtain

[TABLE]

where the last inequality uses $M$-smoothness of $f(\cdot)$ and Lemma 6.1. Also, from the subspace embedding property of $\mathbf{S}_{t}$ (see Lemma 6.1), we have

[TABLE]

Now, using the above inequality and Eqs. (6.3) and (6.3), we get

[TABLE]

which is the desired result. ∎

Now we are ready to prove Theorem 3.3. Let $\mathbf{H}_{t}=\mathbf{U}\Sigma\mathbf{U}^{T}$ and $\hat{\mathbf{H}}_{t}=\hat{\mathbf{U}}\hat{\Sigma}\hat{\mathbf{U}}^{T}$ be the truncated SVDs of $\mathbf{H}_{t}$ and $\hat{\mathbf{H}}_{t}$, respectively. Also, let $\alpha_{t}$ be the step-size obtained using line-search in (6) in the $t$-th iteration. Thus, Eq. (6) with the update direction $\mathbf{p}_{t}=-\hat{\mathbf{H}}^{\dagger}_{t}\mathbf{g}_{t}$ implies

[TABLE]

where the last equality uses the fact that $\hat{\mathbf{H}}^{\dagger}_{t}$ can be expressed as $\hat{\mathbf{H}}^{\dagger}_{t}=\hat{\mathbf{U}}\hat{\Sigma}^{-1}\hat{\mathbf{U}}^{T}$. Note that Lemma 6.2 implies that the function $||\nabla f(\mathbf{y})||^{2}/2$ is smooth for all $y\in\mathcal{Y}$, where $\mathcal{Y}=\{\mathbf{y}\in\mathbb{R}^{d}|~{}||\nabla f(\mathbf{y})||\leq\delta\}$. Smoothness in turn implies the following property (see [59], Theorem 2.1.10)

[TABLE]

where $Q=L\delta+M^{2}$. We take $\delta=||\nabla f(\mathbf{w}_{0})||$ where $\mathbf{w}_{0}$ is the initial point of our algorithm. Due to line-search condition in (6), it holds that $||\nabla f(\mathbf{w}_{t})||\leq||\nabla f(\mathbf{w}_{0})||~{}\forall~{}t>0$. Thus, substituting $\mathbf{x}=\mathbf{w}_{t}$ and $\mathbf{y}=\mathbf{w}_{t+1}=\mathbf{w}_{t}+\alpha_{t}\mathbf{p}_{t}$, we get

[TABLE]

where

[TABLE]

Also, since the minimum non-zero eigenvalue of $\mathbf{H}_{t}\geq\eta$ from Assumption (5), the minimum non-zero eigenvalue of $\hat{\mathbf{H}}_{t}$ is at least $\eta-\epsilon M$ from Lemma 6.1. Thus,

[TABLE]

|Moreover,

[TABLE]

Using this in (6.3), we get

[TABLE]

Now,

[TABLE]

where the last inequality is obtained by applying the triangle inequality and Cauchy-Schwartz inequality. This can be further simplified using Lemma 6.1 and Eq. (40) as

[TABLE]

Using the above in Eq. (42), we get

[TABLE]

Note that the upper bound in Eq. (43) always holds. Also, we want the inequality in (6.3) to hold for some $\alpha_{t}>0$. Therefore, we want $\alpha_{t}$ to satisfy the following (and hope that it is always satisfied for some $\alpha_{t}>0$)

[TABLE]

Thus, any $\alpha_{t}$ satisfying the above inequality would definitely satisfy the line-search termination condition in

Now, using Lemma 6.3 and Assumption (6), we have

[TABLE]

Using the above in Eq. (6.3) to find an iteration independent bound on $\alpha_{t}$, we get

[TABLE]

Hence, line-search will always terminate for all $\alpha_{t}$ that satisfy the above inequality. This can be further simplified by assuming that $\epsilon$ is small enough such that $\epsilon<\eta/2M$. Thus, $\eta-M\epsilon>\eta/2$, and the sufficient condition on $\alpha_{t}$ in (46) becomes

[TABLE]

For a positive $\alpha_{t}$ to always exist, we require $\epsilon$ to further satisfy

[TABLE]

which is tighter than the initial upper bound on $\epsilon$. Now, Eqs. (6.3) and (45) proves the desired result, that is

[TABLE]

Thus, OverSketched Newton for the weakly-convex case enjoys a uniform linear convergence rate of decrease in $||\nabla f(\mathbf{w})||^{2}$.

7 Conclusions

We proposed OverSketched Newton, a straggler-resilient distributed optimization algorithm for serverless systems. It uses the idea of matrix sketching from RandNLA to find an approximate second-order update in each iteration. We proved that OverSketched Newton has a local linear-quadratic convergence rate for the strongly-convex case, where the dependence on the linear term can be made to diminish by increasing the sketch dimension. Moreover, it has a linear global convergence rate for weakly-convex functions. By exploiting the massive scalability of serverless systems, OverSketched Newton produces a global approximation of the second-order update. Empirically, this translates into faster convergence than state-of-the-art distributed optimization algorithms on AWS Lambda.

Acknowledgments

This work was partially supported by NSF grants CCF-1748585 and CNS-1748692 to SK, and NSF grants CCF-1704967 and CCF- 0939370 (Center for Science of Information) to TC, and ARO, DARPA, NSF, and ONR grants to MWM, and NSF Grant CCF-1703678 to KR. The authors would like to additionally thank Fred-Roosta and Yang Liu for helpful discussions regarding our proof techniques and AWS for providing promotional cloud credits for research.