Distributed Sketching for Randomized Optimization: Exact Characterization, Concentration and Lower Bounds

Serverless computing is a relatively new technology that offers computing on the cloud without requiring any server management from end users. Functions in serverless computing can be thought of worker nodes that have very limited resources and a lifetime, but also are very scalable. The algorithms that we study in this paper are particularly suitable for serverless computing platforms, since the algorithms do not require peer-to-peer communication among different worker nodes, and have low memory and compute requirements per node. In numerical simulations, we have evaluated our methods on the serverless computing platform AWS Lambda.

Data privacy is an increasingly important issue in cloud computing, one that has been studied in many recent works including [10, 11, 12, 13]. Distributed sketching comes with the benefit of privacy preservation. To clarify, let us consider a setting where the master node computes sketched data $S_{k}A$, $S_{k}b$, $k=1,\dots,q$ locally where $S_{k}\in\mathbb{R}^{m\times n}$ are the sketching matrices and $A\in\mathbb{R}^{n\times d}$, $b\in\mathbb{R}^{n}$ are the data matrix and the output vector, respectively, for the regression problem $\min_{x}\|Ax-b\|_{2}^{2}$. In the distributed sketching setting, the master node sends only the sketched data to worker nodes for computational efficiency, as well as data privacy preservation. In particular, the mutual information and differential privacy can be controlled when we reveal $S_{k}A$ and keep $A$ hidden. Furthermore, one can trade privacy for accuracy by choosing a suitable sketch dimension $m$.

We use hats $\hat{x}$ to denote the estimator for a single sketch and bars $\bar{x}$ to denote the averaged estimator $\bar{x}:=\frac{1}{q}\sum_{k=1}^{q}\hat{x}_{k}$, where $\hat{x}_{k}$ is the estimator for the $k$’th worker node. Stars are used to denote the optimal solution $x^{*}$. We use $f(\cdot)$ to denote the objective of the optimization problem being considered at that point in the text. The letter $\epsilon$ is used for error while $\varepsilon$ is used as the differential privacy parameter. The notation $\sigma_{\min}(\cdot)$ denotes the smallest nonzero singular value of its argument. $O(\cdot)$ is used for big O notation.

All the expectations in the paper are with respect to the randomness over the sketching matrices, and no randomness assumptions are made for the data. We use $S\in\mathbb{R}^{m\times n}$ to denote random sketching matrices and we often refer to $m$ as the sketch size. For non-iterative distributed algorithms, we use $S_{k}\in\mathbb{R}^{m\times n}$ for the sketching matrix used by worker node $k$. For iterative algorithms, $S_{t,k}\in\mathbb{R}^{m\times n}$ is used to denote the sketching matrix used by worker $k$ in iteration $t$. We assume the sketching matrices are scaled such that $\operatorname{\mathbb{E}}[S^{T}S]=I_{n}$. We omit the subscripts in $S_{k}$ and $S_{t,k}$ for simplicity whenever it does not cause confusion. The notation $m\gtrsim d$ is used to denote that there exists a positive finite constant $c$ such that $m\geq cd$. We denote the thin SVD decomposition of the data matrix as $A=U\Sigma V^{T}$.

For regularized problems, we use $\lambda_{1}$ for the regularization coefficient of the original problem, and $\lambda_{2}$ for the regularization coefficient of the sketched problem. For instance, in the case of the regularized least squares problem, we have

We will derive expressions in the sequel for the optimal selection of the coefficient $\lambda_{2}$ which we denote as $\lambda_{2}^{*}$.

Random projections are a popular way of performing randomized dimensionality reduction, which are widely used in many computational and learning problems [14, 15, 16, 17]. Many works have studied randomized sketching methods for least squares and other optimization problems [18, 19, 2, 20, 21, 3, 22].

Sarlós [23] showed that the relative error of the single Gaussian sketch estimator in (I-B) in the unregularized case $\lambda_{1}=\lambda_{2}=0$ is bounded as $f(\hat{x})/f(x^{*})\leq(1+Cd\log(d)/m)^{2}$ with probability at least $1/3$, where $C$ is a constant. The relative error was improved to $f(\hat{x})/f(x^{*})\leq(1+Cd/m)$ with exponentially small failure probability in subsequent work [7]. In contrast, we derive the expectation of the error exactly and show exponential concentration around this expected value in this paper. As a result, we obtain a tight upper and lower bound for the relative error that holds with exponentially high probability (see Theorem II.8).

The work [3] investigates distributed sketching and averaging for regression from optimization and statistical perspectives. The most relevant result in [3], using our notation, can be stated as follows. Setting the sketch dimension $m=O(\mu d(\log d)/\epsilon)$ for uniform sampling, where $\mu$ is row coherence $\mu:=n/\operatorname{rank}(A)\max_{i}\ell_{i}$ and $\ell_{i}$ is the $i$’th row leverage score, and $m=\tilde{O}(d/\epsilon)$ (with $\tilde{O}$ hiding logarithmic factors) for other sketches, the inequality $f(\bar{x})-f(x^{*})\leq(\epsilon/q+\epsilon^{2})f(x^{*})$ holds with high probability. According to this result, for large $q$, the cost at the averaged solution $f(\bar{x})$ will be upper bounded by $\epsilon^{2}f(x^{*})$. In this work we prove that for Gaussian sketch, $f(\bar{x})$ will converge to the optimal cost $f(x^{*})$ as $q$ increases. We also identify the exact expected error for a given number of workers $q$. In addition, our results on regularized least squares regression improve on the results in [3] for averaging multiple sketched solutions. In particular, in [3], the sketched sub-problems use the same regularization parameter as the original problem, which leads to biased solutions. We analyze the bias of the averaged solution, and provide explicit formulas for selecting the regularization parameter of the sketched sub-problems to achieve unbiasedness.

We show that the expected difference between the costs of the averaged solution and the optimal solution has two components, namely variance and squared bias (see Lemma IV.1). This result implies that for the Gaussian sketch, which we prove to be unbiased (see Lemma II.1), the number of workers required for a target error $\epsilon$ scales as $1/\epsilon$. Remarkably, this result does not involve condition numbers. In contrast, for the asynchronous Hogwild algorithm [9], the number of iterations required for error $\epsilon$ scales with $\log(1/\epsilon)/\epsilon$ and also depends on the condition number of the input data, although it addresses a more general class of problems.

One of the crucial results of our work is developing bias correction for averaging sketched solutions. Namely, the setting considered in [22] is based on a distributed second order optimization method that involves averaging approximate update directions. However, in that work, the bias was not taken into account which degrades accuracy and limits applicability to massively parallel computation. We additionally provide results for the unregularized case, which corresponds to the distributed version of the Newton sketch algorithm introduced in [21].

• •

  We characterize the exact expected error of the averaged estimator for Gaussian sketch in closed form as

  $\displaystyle\frac{\operatorname{\mathbb{E}}[f(\bar{x})]-f(x^{*})}{f(x^{*})}=\frac{1}{q}\frac{d}{m-d-1}$

  (2)

  where $\bar{x}=\frac{1}{q}\sum_{k=1}^{q}\hat{x}_{k}$ is the averaged solution and $\hat{x}_{k}=\arg\min_{x}\|S_{k}Ax-S_{k}b\|^{2}_{2}$ is the output of the $k$’th worker node. In addition, we obtain a similar result for the error of the averaged estimator for the least-norm solution in the underdetermined case $n<d$.

• •

  For both the single sketch and averaged estimator, we show that the relative error $\frac{f(\bar{x})-f(x^{*})}{f(x^{*})}$ is concentrated around its expectation given in (2) with exponentially high probability. This provides both an upper and lower bound on the performance of sketching, unlike previous results in the literature. Moreover, it offers a guaranteed recipe to set the sketch size $m$ and number of workers $q$.

• •

  We show that for Gaussian distributed sketch, the expected error of the averaged estimator $\operatorname{\mathbb{E}}[f(\bar{x})]-f(x^{*})$ matches the error lower bound for any unbiased estimator obtained via Fisher information. In addition, we provide a lower bound for general, possibly biased sketching based estimators.

• •

  We consider the privacy preserving properties of distributed sketching methods, in which only random projections of the data matrix $A$ need to be shared with worker notes. We derive conditions on the $(\varepsilon,\delta)$-differential privacy when $S$ is the i.i.d. Gaussian sketch matrix. Combined with our results, we show that the approximation error of this distributed privacy preserving regression algorithm scales as $O(1/\varepsilon^{2})$.

• •

  We analyze the convergence rate of a distributed version of the iterative Hessian sketch algorithm which was introduced in [4] and show that the number of iterations required to reach error $\epsilon$ with $q$ workers scales with $\log(1/\epsilon)/\log(q)$.

• •

  We show that $\hat{x}=\arg\min_{x}\|SAx-Sb\|_{2}^{2}+\lambda_{2}\|x\|_{2}^{2}$ is not an unbiased estimator of the optimal solution, i.e. $\operatorname{\mathbb{E}}[A(\hat{x}-x^{*})]\neq 0$ when $\lambda_{2}=\lambda_{1}$ and provide a closed form expression for the regularization coefficient $\lambda_{2}^{*}$ to make the estimator unbiased under certain assumptions.

• •

  In addition to Gaussian sketch, we derive bias bounds for uniform sampling, randomized Hadamard sketch and leverage score sampling. Analysis of the bias term is critical in understanding how close to the optimal solution we can hope to get and establishing the dependence on the sketch dimension $m$. Moreover, we utilize the derived bias bounds and find an upper bound on the error of the averaged estimator for these other sketching methods.

• •

  We discuss averaging for the distributed version of Newton Sketch [21] and show that using the same regularization coefficient as the original problem, i.e., $\lambda_{2}=\lambda_{1}$, as most works in the literature consider, is not optimal. Furthermore, we derive an expression for choosing the regularization coefficient $\lambda_{2}^{*}$ for unbiasedness.

• •

  We provide numerical simulations that illustrate the practicality and scalability of distributed sketching methods on the serverless cloud computing platform AWS Lambda.

We consider various sketching matrices in this work including Gaussian sketch, uniform sampling, randomized Hadamard sketch, Sparse Johnson-Lindenstrauss Transform (SJLT), and hybrid sketch. We now briefly describe each of these sketching methods:

1. 1.

   Gaussian sketch: Entries of $S\in\mathbb{R}^{m\times n}$ are i.i.d. and sampled from the Gaussian distribution. Sketching a matrix $A\in\mathbb{R}^{n\times d}$ using Gaussian sketch requires computing the matrix product $SA$ which has computational complexity equal to $O(mnd)$.

2. 2.

   Randomized Hadamard sketch: The sketch matrix in this case can be represented as $S=PHD$ where $P\in\mathbb{R}^{m\times n}$ is for uniform sampling of $m$ rows out of $n$ rows, $H\in\mathbb{R}^{n\times n}$ is the Hadamard matrix, and $D\in\mathbb{R}^{n\times n}$ is a diagonal matrix with diagonal entries sampled randomly from the Rademacher distribution. Multiplication by $D$ to obtain $DA$ requires $O(nd)$ scalar multiplications. Hadamard transform can be implemented as a fast transform with complexity $O(n\log(n))$ per column, and a total complexity of $O(nd\log(n))$ to sketch all $d$ columns of $DA$.

3. 3.

   Uniform sampling: Uniform sampling randomly selects $m$ rows out of the $n$ rows of $A$ where the probability of any row being selected is the same.

4. 4.

   Leverage score sampling: Row leverage scores of a matrix $A$ are given by $\ell_{i}=\|\tilde{u}_{i}\|_{2}^{2}$ for $i=1,\dots,n$ where $\tilde{u}_{i}\in\mathbb{R}^{d}$ denotes the $i$’th row of $U$. The matrix $U$ is the matrix whose columns are the left singular vectors of $A$, i.e., $A=U\Sigma V^{T}$. There is only one nonzero element in every row $s_{i}\in\mathbb{R}^{n}$ of the sketching matrix $S$ and the probability that the $j$’th entry of $s_{i}$ is nonzero is proportional to the leverage score $\ell_{j}$. More precisely, the rows $s_{1},\dots,s_{m}$ are sampled i.i.d. such that $\mathbb{P}[s_{i}=e_{k}/\sqrt{mp_{k}}]=p_{k},\,\forall i,\,\forall k$ where $p_{k}=\frac{\ell_{k}}{\sum_{j=1}^{n}\ell_{j}}$. Note that we have $\operatorname{\mathbb{E}}[S^{T}S]=\operatorname{\mathbb{E}}[\sum_{i=1}^{m}s_{i}s_{i}^{T}]=m\sum_{k=1}^{n}p_{k}e_{k}e_{k}^{T}/(mp_{k})=I_{n}$.

5. 5.

   Sparse Johnson-Lindenstrauss Transform (SJLT) [24]: The sketching matrix for SJLT is a sparse matrix where each column has exactly $s$ nonzero entries and the columns are independently distributed. The nonzero entries are sampled from the Rademacher distribution. It takes $O(snd/m)$ addition operations to sketch a data matrix using SJLT.

6. 6.

   Hybrid sketch: The method that we refer to as hybrid sketch is a sequential application of two different sketching methods. In particular, it might be computationally feasible for worker nodes to sample as much data as possible, say $m^{\prime}$ rows, and then reduce the dimension of the available data to the final sketch dimension $m$ using another sketch with better error properties than uniform sampling such as Gaussian sketch or SJLT. For instance, a hybrid sketch of uniform sampling followed by Gaussian sketch would have computational complexity $O(mm^{\prime}d)$.

Section II deals with the application of distributed sketching to quadratic problems for Gaussian sketch. Section III deals with the privacy preserving property of distributed sketching. Section IV gives theoretical results for randomized Hadamard sketch, uniform sampling, and leverage score sampling. The tools developed in Section II are applied to second order optimization for unconstrained convex problems in Section V. Section VI provides an overview of applications and example problems where distributed sketching methods could be applied. Section VII presents numerical results and Section VIII concludes the main part of the paper.

We give the proofs for the majority of the lemmas and theorems in the Appendix along with additional numerical results.

We begin with a non-iterative model averaging based algorithm, which we refer to as distributed randomized regression. This algorithm is for the unregularized case $\lambda_{1}=0$. Each of the worker nodes computes an approximate solution $\hat{x}_{k}$ and these are averaged at the master node to compute the final solution $\bar{x}$. This method is outlined in Algorithm 1. The theoretical analysis in this section assumes that the sketch matrix is Gaussian, which is generalized to other sketching matrices in Section IV. Although computing the Gaussian sketch is not as efficient as other fast sketches such as the randomized Hadamard sketch, it has several significant advantages over other sketches: (1) exact relative error expressions can be derived as we show in this section, (2) the solution is unbiased, (3) a differential privacy bound can be provided as we show in Section III, and (4) computation of the sketch can be trivially parallelized.

We first obtain a characterization of the expected error for the single sketch estimator in Lemma II.1.

To the best of our knowledge, this result of exact error characterization is novel in the theory of sketching. Similar tools to those that we used in this proof were used in [26], however, the exact context in which they are used and the end result are different from this work. Furthermore, existing results (see e.g. [6, 3, 16]) characterize a high probability upper bound on the error, whereas the above is a sharp and exact formula for the expected squared norm of the error. Theorem II.2 builds on Lemma II.1 to characterize the expected error for the averaged solution $\bar{x}$.

Theorem II.2 illustrates that the expected error scales as $1/q$, and converges to zero as $q\rightarrow 0$ as long as $m\geq d+2$. This is due to the unbiasedness of the Gaussian sketch. Other sketching methods such as uniform sampling or randomized Hadamard sketch do not have this property, as we further investigate in Section IV.

In this subsection, we show high probability concentration bounds for the error of the Gaussian sketch. Importantly, our result provides both an upper and lower bounds on the relative error with exponentially small failure probability. We begin with the single sketch estimator $\hat{x}$ and extend the results to the averaged estimator $\bar{x}=\frac{1}{q}\sum_{k=1}^{q}\hat{x}_{k}$.

The next theorem states the main concentration result for the single sketch estimator. Both upper and lower bounds are given for the ratio of the cost of the sketched solution to the optimal solution. The full proof is given in the Appendix.

Then we note that after applying Lemma II.5, the conditional expectation of the error is given by $\operatorname{\mathbb{E}}[\|A(\hat{x}-x^{*})\|_{2}^{2}\,|SA]=\frac{f(x^{*})}{m}\operatorname{tr}((U^{T}S^{T}SU)^{-1})$. Next, we focus on the trace term $\operatorname{tr}((U^{T}S^{T}SU)^{-1})$ and relate it to the Stieltjes transform.

We derive a high probability bound for the trace term around its expectation by leveraging the concentration of the empirical Stieltjes transform to its expectation, which is given in Lemma II.7 below.

The proof of Lemma II.7 is provided in the Appendix. Combining the concentration of the Stieltjes transform with the concentration of Gaussian quadratic forms completes the proof of Theorem II.4.

Next, we show that the relative error of the averaged estimator is also concentrated with exponentially small failure probability.

The above result shows that the relative error of the averaged estimator using distributed sketching concentrates considerably faster compared to a single sketch. Moreover, the above bound offers a method to choose the values of the sketch size $m$ and number of workers $q$ to achieve a desired relative error with exponentially small error probability.

Figure 1 shows the histograms for the single sketch ($q=1$) and averaged estimator ($q=4$ and $q=10)$ for Gaussian and randomized Hadamard sketches. This experiment demonstrates that the optimality ratio $f(\bar{x})/f(x^{*})$ is concentrated around its mean and its mean and concentration improve as we increase the number of workers $q$.

We now present an error lower bound result for the Gaussian sketch. We first consider single sketch estimators and then discuss the distributed sketching case. Two different lower bounds can be obtained depending on whether the estimator is restricted to be unbiased or not. Lemma II.9 and II.10 provide error lower bounds for all unbiased and general (i.e., possibly biased) estimators, respectively. The proof of Lemma II.9 is based on Fisher information and Cramér-Rao lower bound while Lemma II.10 additionally employs the van Trees inequality. We refer the reader to [28] for details on these single sketch lower bounds.

We note that the expected error of the single sketch estimator that we provide in Lemma II.1 exactly matches the error lower bound in Lemma II.9. Therefore, we conclude that no other unbiased estimator based on a single sketch $SA$, $Sb$ can achieve a better expected relative error.

We now provide a novel generalized lower bound that applies to averaged estimators $\bar{x}$ obtained via distributed sketching. Theorem II.11 states the main lower bound result for the averaged estimator.

We note that the expected error of the averaged estimator given in Theorem II.2 matches exactly the lower bound in Theorem II.11 for unbiased estimators. For general estimators, we observe that the lower bound matches the upper bound for large $m$, e.g., when $m-d-1=O(m)$.

In this section, we consider an iterative algorithm for solving the unregularized least squares problem in (3), where $\lambda_{1}=0$ with higher accuracy. Note that, Newton’s method terminates in one step when applied to this problem since the Hessian is $A^{T}A$ and

reduces to directly solving the normal equations $(A^{T}A)x=A^{T}b$. However, the computational cost of this direct solution is often prohibitive for large scale problems. To remedy this, the method of iterative Hessian sketch introduced in [4] employs a randomly sketched Hessian $A^{T}S_{t}^{T}S_{t}^{T}A$ as follows

where $S_{t}$ corresponds to the sketching matrix at iteration $t$. Sketching reduces the row dimension of the data from $n$ to $m$ and hence computing an approximate Hessian $A^{T}S_{t}^{T}S_{t}A$ is computationally cheaper than the exact Hessian $A^{T}A$. Moreover, for regularized problems one can choose $m$ to be smaller than $d$ as we investigate in Section II-F.

In a distributed computing setting, one can obtain more accurate update directions by averaging multiple trials, where each worker node computes an independent estimate of the update direction. These approximate update directions can be averaged at the master node and the following update takes place

Here $S_{t,k}$ is the sketching matrix for the $k$’th node at iteration $t$. The details of the distributed IHS algorithm are given in Algorithm 2. We note that although the update equation involves a matrix inverse term, in practice, this can be replaced with an approximate linear system solver. In particular, it might be computationally more efficient for worker nodes to compute their approximate update directions using indirect methods such as conjugate gradient.

We note that in Algorithm 2, worker nodes communicate their approximate update directions and not the approximate Hessian matrix, which reduces the communication complexity from $O(d^{2})$ to $O(d)$ for each worker per iteration.

We establish the convergence rate for Gaussian sketch in Theorem II.13, which provides an exact characterization of the expected error. First, we give the definition of the error in Definition II.12.

To state the result, we first introduce the following moments of the inverse Wishart distribution (see Appendix).

Next, Corollary II.14 builds on Theorem II.13 to characterize the number of iterations for Algorithm 2 to achieve an error of $\epsilon$. The number of iterations required for error $\epsilon$ scales with $\log(1/\epsilon)/\log(q)$.

In this section, we consider the underdetermined case where $n<d$ and applying the sketching matrix from the right, i.e., on the features. We will refer to this method as right sketch as well. Let us define the minimum norm solution

The above problem has a closed-form solution given by $x^{*}=A^{T}(AA^{T})^{-1}b$ when the matrix $A$ is full row rank. We will assume that the full row rank condition holds in the sequel. Let us denote the optimal value of the minimum norm objective as $f(x^{*})=\|x^{*}\|_{2}^{2}=b^{T}(AA^{T})^{-1}b$. The $k$’th worker node will compute the approximate solution via $\hat{x}_{k}=S_{k}^{T}\hat{z}_{k}$ where

where $S_{k}\in\mathbb{R}^{m\times d}$ and $z\in\mathbb{R}^{m}$. The averaged solution is then computed at the master node as $\bar{x}=\frac{1}{q}\sum_{k=1}^{q}\hat{x}_{k}$. We will assume that the sketch matrices $S_{k}$ are i.i.d. Gaussian in deriving the error expressions. Lemma II.16 establishes the approximation error for a single right sketch estimator.

An exact formula for averaging multiple outputs in right sketch that parallels Theorem II.2 can be obtained in a similar fashion:

Hence, for the distributed Gaussian right sketch, we establish the approximation error as stated in Theorem II.17.

We have previously studied distributed randomized regression for unregularized problems and showed that Gaussian sketch leads to unbiased estimators. In this section, we focus on the regularized case and show that using the original regularization coefficient for the sketched problems causes the estimators to be biased. In addition, we provide a bias correction procedure. More precisely, the method described in this section is based on non-iterative averaging for solving the linear least squares problem with $\ell_{2}$ regularization, i.e., ridge regression.

Note that we have $\lambda_{1}$ as the regularization coefficient of the original problem and $\lambda_{2}$ for the sketched sub-problems. If $\lambda_{2}$ is chosen to be equal to $\lambda_{1}$, then this scheme reduces to the framework given in the work of [3] and we show in Theorem II.20 that $\lambda_{2}=\lambda_{1}$ leads to a biased estimator, which does not converge to the optimal solution.

We first introduce the following results on traces involving random Gaussian matrices which are instrumental in our result.

We list the distributed randomized ridge regression method in Algorithm 3. This algorithm assumes that our goal is to solve a large scale regression problem for a given value of regularization coefficient $\lambda_{1}$. Theorem II.20 states the main result of this subsection. In short, for the averaged result to converge to the optimal solution $x^{*}$, the regularization coefficient needs to be modified to $\lambda_{2}^{*}$.

Figure 2 illustrates the practical implications of Theorem II.20. If $\lambda_{2}$ is chosen according to the formula in (25), then the averaged solution $\bar{x}$ gives a significantly better approximation to $x^{*}$ than if we had used $\lambda_{2}=\lambda_{1}$. The data matrix $A$ in Figure 2(a) has identical singular values, and 2(b) shows the case where the singular values of $A$ are not identical. When the singular values of $A$ are not all equal to each other, we set $\sigma$ to the mean of the singular values of $A$ as a heuristic, which works extremely well as shown in Figure 2(b). According to the formula in (25), the value of $\lambda_{2}$ that we need to use to achieve zero bias is found to be $\lambda_{2}^{*}=0.833$ whereas $\lambda_{1}=5$. The plot in Figure 2(b) illustrates that even if the assumption that $\Sigma=\sigma I_{d}$ in Theorem II.20 is violated, the proposed bias corrected averaging method outperforms vanilla averaging in [3] where $\lambda_{2}=\lambda_{1}$.