Sketched Iterative Algorithms for Structured Generalized Linear Models

1 Introduction

In this paper, we consider algorithms for efficiently learning generalized linear models (GLMs) of the form:

where xi ∈ Rp is a data point, yi ∈ R is the response, and θ∗ ∈ Rp is the structured parameter that we want to learn. For example, in linear regression yi = xTi θ∗ +wi, where wi is the noise; in logistic regression yi ∈ {-1, +1} and P(yi|xTi θ∗ ) = 1 1+e-yixTi θ∗ . θ∗ can be sparse, low rank, etc [Tibshirani, 1996; Recht et al., 2010]. Estimators for θ∗ can be posed as the following constrained problem [Oymak et al., 2018]

where l : R×R → R is the loss function, and is a function of xTi θ and yi, e.g., in linear regression, l(ˆ y; y) is the square loss 12 (y - ˆ y)2 , in logistic regression, l(ˆ y; y) is the logistic loss log(1 + e-yˆ y ), etc., n is the sample size, and R(.) is a suitable norm inducing structures in θ [Banerjee et al., 2014;

Negahban et al., 2012], e.g., sparsity can be induced by L1 norm or k-support norm [Tibshirani, 1996; Argyriou et al., 2012], low-rank can be induced by nuclear norm [Recht et al., 2010], etc. In this paper, we set λ = R(θ∗ ).

In modern machine learning, problems like (2) are confronted with extremely large datasets, i.e., sample size n is large. Dealing with such large datasets brings about several computational challenges, including large per iteration complexity for methods like gradient descent and memory management challenges for most methods. In this work, inspired by sketching [Woodruff, 2014], a powerful technique for sample size reduction based on random projection, we propose the sketched version of two standard iterative methods, projected gradient descent (PGD) and stochastic variance reduced gradient (SVRG). Since sketching of the data happens once upfront, the iterative algorithms work with a dataset of much smaller effective size.

In this work, we reduce the computation of GLMs while keeping the structure of parameters. For statistical estimation problem (2), at each iteration of PGD and each stage of SVRG, we need to do matrix vector multiplication, and sketching reduces the per iteration time complexity. We show that PGD and SVRG will reach a solution ˆ θ such that ‖ˆ θ - θ∗ ‖2 ≤ Cn ∥∥∥[ ∂ ∂ˆ y l(xTi θ∗ ; yi)]n i=1∥∥∥2 , for constant C > 0. We have two main contributions. First, we present novel sketched iterative algorithms S-PGD and S-SVRG for general loss function of form (2). Our algorithms combine iterative algorithms and sketching. For arbitrary δ > 0, if we choose sketching dimension m≥ c δ2 for constant c > 0 , then the approximate solution ˆ θ of our algorithms satisfies

for positive constants C1 and C2. Term ‖( ∂ ∂ˆ y l(xTi θ∗ ; yi))ni=1 ‖2 is related to the estimation error of model (1). Different choices of δ determine the accuracy of our approximation. In the literature, progress have been made in applying sketching to least squares [Woodruff, 2014]. However, sketching methods for general optimization problems are limited [Tang et al., 2017; Pilanci and Wainwright, 2016]. [Pilanci and Wainwright, 2017] requires approximating the Hessian of program objective function using sketching in every iteration, which can be computationally expensive . Our algorithms approximate the gradient, which effectively alleviate the computational burden for large-scale problems. To the best of our knowledge, it is the first sketched first-order algorithm for GLMs.

Our second main contribution is to establish the theoretical convergence rates of the proposed sketched algorithms, and experiments on synthetic date to illustrate the theoretical results. We characterize the convergence rates of both S-PGD and S-SVRG, with suitable modifications to standard SVRG. Our results hold for any norm R(.). We show that for sketched algorithms, sketching dimension m is a key factor that determines the tradeoff between computational time and approximation accuracy. We analyze both S-PGD and S-SVRG for constrained least squares and S-PGD for general GLMs. We show that as long as a statistical condition called restricted strong convexity holds, geometric convergence rate can be established even when m < p.

The rest of paper is organized as follows. We first review the relevant background and related works in Section 2. We propose our algorithms in Section 3. In Section 4, we discuss convergence of S-PGD and S-SVRG for constrained least squares. In Section 5, we extend our convergence results to general loss function. In Section 6, we compare our algorithm with SGD. We show our algorithms' experimental results and conclude in Section 7 and 8 respectively.

Notations: Throughout this paper, we denote l′ (ˆ y; y) = ∂ ∂ˆ y l(ˆ y; y) and l′′ (ˆ y; y) = ∂ 2 ∂ˆ y2 l(ˆ y; y), matrix X = [x1, . . . , xn]T the stack up of all data points, vector dn(θ) = [l′ (xT1 θ; y1) . . . l′ (xTn θ; yn)]T the first order derivatives of all loss functions, and vector d(2) n (θ) = [l′′ (xT1 θ; y1) . . . l′′ (xTn θ; yn)]T the second order derivatives of all loss functions. For a given positive integer n, denote [n] = {1, . . . , n}. We use In as the n × n identity matrix. We denote Bn a unit ball and Sn-1 a unit sphere in Rn .

2 Background and RelatedWork

2.1 Projected Gradient Descent (PGD)

One efficient first order algorithm for problem (2) is the projected gradient descent (PGD) [Bertsekas, 2010]. Each step of PGD is given byθt+1 = PK(θt - η∇Ln(θt)), (3)

where η is the learning rate, PK(.) is the projection operator such that PK(z) = argminy∈K ‖y - z‖22 , and K is the constraint set K = {θ : R(θ) ≤ λ}. The computation of PGD includes a gradient update and a projection onto feasible set.

2.2 Stochastic Variance Reduction Gradient (SVRG)

If sample size of problem (2) is large, then a popular modification of PGD is stochastic variance reduction gradient (SVRG) [Johnson and Zhang, 2013; Schmidt et al., 2017; Defazio et al., 2014]. Each outer iteration or stage s of SVRG maintains an estimation of solution to (2) θs , and a full gradient at θs denoted by µs = ∇Ln(θs ). An inner iteration t + 1 is given by the following rule

2.3 Sub-Gaussian Sketch (SGS)

Given matrices A ∈ Rn×p and B ∈ Rn×q , the time complexity of matrix product AT B is O(pqn). Sketching approximates AT B by first sampling a random matrix S ∈ Rm×n with m < n, and then computing SA and SB. Computing the production (SA)T (SB) requires only O(pqm).

We say S is a σ-sub-Gaussian random matrix [Vershynin, 2012] if each row of S is an i.i.d. copy of a random vector x, which satisfies supu∈Sn-1 P(|〈u, x〉| ≥ t) ≤ 2e- nt2 2σ2 . One example of sub-Gaussian random matrix is the Gaussian random matrix, whose each entry is an i.i.d. sample from standard Gaussian distribution. Sub-Gaussian sketch takes O(mn(p + q)) basic operations to compute sketched matrices SA and SB. In this paper, we assume EST S = In.

2.4 RandomOrthogonal Sketch (ROS)

In order to construct sketched matrices faster, we consider random orthogonal sketch [Pilanci and Wainwright, 2015]. Let S = √ nm P · H · D, where H ∈ Rn×n is an orthonormal matrix, P ∈ Rm×n samples m coordinates of Rn uniformly at random, and D ∈ Rn×n is a diagonal re-randomization matrix of Rademacher random variables. One commonly used ROS is random Hadamard sketch (RHS) [Ailon and Liberty, 2008]. H is called a Hadamard matrix if Hi,j ∈ {- 1 √n , 1 √n }. For any x ∈ Rn the mapping Hx can be computed in O(n log n). RHS takes O(n(p + q) log(m)) basic operations to compute sketched matrices SA and SB.

2.5 GaussianWidth

In this paper, both optimization and statistical results are described in terms of a geometric measure called Gaussian width [Banerjee et al., 2014; Vershynin, 2015; Chandrasekaran et al., 2012]. For any set A ∈ Rp , Gaussian width w(A) measures the complexity of A, and it is defined as

where g ∼ N(0, In) is a standard Gaussian random vector.

2.6 RelatedWork

Let the minimizer of problem (2) be θ′ . Statistical guarantees and computational methods for (2) has been extensively studied. From statistical perspective, much progress has been made in the regime p ≫ n. [Chandrasekaran et al., 2012] analyzed exact recovery guarantee of atomic norm minimization linear inverse problem in noiseless setting, under restricted strong convexity (RSC) condition [Raskutti et al., 2010]. In noisy setting, [Negahban et al., 2012] established a bound of estimation error ‖θ′ - θ∗ ‖2 for decomposable norm regularized problem; [Banerjee et al., 2014] analyzed general norm regularized estimator using Gaussian width [Vershynin, 2015].

From the computational perspective, two key ideas have been explored in recent works. Firstly, randomized optimization algorithm has been investigated to reduce computation and storage of large scale problems. One line of research is random projection technique including importance sampling [Mahoney, 2011], sub-Gaussian sketch [Pilanci and Wainwright, 2015], fast Johnson-Lindenstrauss transforms [Ailon and Liberty, 2008], and dimensionality reduction [Zhang et al., 2013;

Algorithm 1: Sketched Projected Gradient Descent (S-PGD)

Inputs: X ∈ Rn×p .

Initialize: θ0 = 0.

Generate a random matrix S ∈ Rm×n

for t = 0, 1, . . . , T do

end for

Wang et al., 2017] etc.. In [Pilanci andWainwright, 2017], a second-order method called "Newton-sketch" was proposed, which approximates the Hessian of objective function by sketching. [Tang et al., 2017] extended [Pilanci and Wainwright, 2017] to first-order method for constrained least squares. The other line of research is stochastic optimization algorithms. SGD and SVRG are two popular examples [Johnson and Zhang, 2013; Robbins and Monro, 1951].

The second key idea is to use statistical assumptions to establish fast convergence of PGD. Two kinds of convergence are considered in literature: convergence to statistical optimum θ∗ and convergence to optimization optimum θ′ . [Agarwal et al., 2012] established convergence of PGD by ‖θt -θ′ ‖2 ≤ ρt ‖θ0-θ′ ‖2+‖θ′ -θ∗ ‖2 . [Qu et al., 2017] extended [Agarwal et al., 2012] to SVRG. In statistical machine learning, a more interesting problem is ‖θt -θ∗ ‖2. [Oymak et al., 2018] showed that estimation error converge geometrically to an error bound under RSC by PGD. [Li et al., 2016] showed linear convergence of SVRG to statistical error bound with nonconvex constrain.

3 Sketched Iterative Algorithms

In this section, we propose our sketched projected gradient descent and sketched stochastic variance reduction gradient.

For problem (2), one step of PGD algorithm is

For problems with large sample size, the computation of matrix-vector product XT · dn(θ) is not efficient. Given a sketching matrix S ∈ Rm×n with m < n, we improve the computation by using sketched matrix SX and vector Sdn(θ). The sketched algorithm is given by 1. We can compute and store the sketched data matrix ˆX = SX upfront and use ˆX in every iteration of optimization algorithm, improving the per iteration complexity. The construction of ˆX takes O(mnp) for SGS and O(mp logn) for RHS. For each iteration, we first construct sketched vector Sdn(θ), which requires O(mn) for SGS and O(m log n) for RHS. Then the time complexity of sketched matrix-vector product XT ST Sdn(θ) is O(pm). If we use RHS, then each iteration takes O(n logm+ pm) that is better than O(pn) for direct computation. We call this algorithm sketched projected gradient descent (S-PGD), and it is described as Algorithm 1. For the special case constrained least squares, it is shown [Pilanci and Wainwright, 2015; Tang et al., 2017] that S-PGD 1 finds the solution to the fol-

Inputs: X ∈ Rn×p . Initialize: θ0 = 0. Generate a random matrix S ∈ Rm×n Let ˆX = SX for s = 0, 1, . . . do µs = 1n ˆXT Sdn(θs ) for t = 0, 1, . . . , T do Randomly select i from 1, . . . ,m. θs t+1 = PK(θs t + η(ˆxisTi (dn(θt) -dn(θs ) +µs )) end for θs+1 = ∑mt=1 1 2m+1-t θs t end for

lowing approximation problem

where ˆ y = Sy. One iteration of S-PGD reduces per iteration computation from O(pn) to O(pm).

By using ideas similar to S-PGD, we propose sketched stochastic variance reduction gradient (S-SVRG) algorithm 2 . Denote ¯ θs 1 = θs 1 , and we average all the intermediate results by ¯ θs t = 12 ¯ θs t-1 + 12 θs t for all 1 < t ≤ m. At last we output θs+1 = ¯ θsm. For large scale problems, S-SVRG applies an one shot sketching, therefore less memory is required in computation. The S-SVRG is given by Algorithm 2 .

4 Constrained Least Squares

In this section, we show the convergence of S-PGD and S-SVRG for constrained least squares. By feasibility, the solution of problem (2) lies in error cone C(θ∗ ), which is the smallest closed cone contains set K - θ denoted as C(θ∗ ) = cone(K - θ∗ ), In this paper, we use C short for C(θ∗ ). We present both optimization and statistical results in terms of Gaussian width w(XC ∩ Sn-1 ). We make the following assumptions:

Assumption 1 (Restricted Eigenvelue (RE) Condition) There is a κn(C) > 0 such that for all v ∈ C

Assumption 2 (Boundedness (BD) Condition) There is a Ln > 0 such that for all θ1, θ2 ∈ K,

For constrained least squares, RE and BD conditions corresponds to the RSC and smoothness condition for general loss function. In literature, the RE and BD conditions have been comprehensively studied. [Negahban et al., 2012; Banerjee et al., 2014] have shown that when X is a 1-subGaussian random matrix, then κn(C) ≥ 1 - w2 (C∩Sp-1 ) n with high probability. [Vershynin, 2012] has shown that when X is sampled from a general distribution, then Ln = O(1 + c np ). We assume X satisfies both RE and BD conditions.

4.1 Convergence of S-PGD

We give the main theoretical result of S-PGD under our assumptions above. Our results precisely characterize the convergence rate and statistical guarantee of S-PGD for constrained least squares estimator .

Theorem 1 Let C = C(θ∗ ) and S ∈ Rm×n be a σ-sub Gaussian random matrix. Let θ∗ ∈ Rp be an arbitrary vector, and y = Xθ∗ + w. With θ0 = 0, we apply sketched PGD update θt+1 = PK(θt + ηn XT ST S(y - Xθt)). Let us assume the RE and BD conditions hold. Set learning rate η = cmn Ln(√m+√n)2 , if m > w2 (XC∩Sm-1 ) δ2 , then with proba-

where ξ(C) = supu∈C∩Bn uT XT w ‖w‖2 .

Theorem 1 shows the precise geometric convergence rate of S-PGD under RSC condition. For sub-Gaussian noise w, we know from literature [Banerjee et al., 2014] that ‖w‖2 = Θ(√n) with high probability. If we choose δ = w(XC∩Sn-1 ) √m , then the second term on the right hand of (10) goes to 1 √nκn ξ(C) with the growth of m. If we also assume X is a sub-Gaussian random matrix, then statistical error bound ξ(C) = O(w(C ∩ Bn)) and RE condition κn(C) = Ω(1 - w(C∩Sp-1 ) n ) with high probability[Banerjee et al., 2014]. Therefore, as n increases, 1 √nκn ξ(C) goes down to zero. The result of [Oymak et al., 2018] shows that the convergence rate ρ(C) of PGD is 1 - κn(C) Ln , thus S-PGD have slower convergence rate ρ(C) due to the effect of sketching. But each update of S-PGD can be computed more efficiently.

4.2 Convergence of S-SVRG

We characterize convergence rate for constrained least squares using S-SVRG. Our conclusions precisely describe how sketching impacts the convergence rate of S-SVRG and the statistical error bound of our estimator. The convergence of S-SVRG is presented in the following theorem:

Theorem 2 Let C = C(θ∗ ), S ∈ Rm×n be a σ-sub Gaussian random matrix. Let θ∗ ∈ Rp be an arbitrary vector, and y = Xθ∗ + w. With θ0 = 0, we apply sketched SVRG update θs t+1 = PK(θs t + η(ˆxiˆxTi (θs - θs t ) - ˆµs )). Denote the i-th row of S as si, ˆxi is given by ˆxi = 1 √n XT si. Let us assume the RE and BD conditions, and Ln (1-δ)κn(C) ≤ 32 hold. Set learning rate η = 1 2Ln , if m > w2 (XC∩Sm-1 ) δ2 , then starting from θ0 = 0, with probability at least 1 - 7 exp (-c1 mδ2 σ4 ),

Theorem 2 indicates the linear convergence of S-SVRG. Similar to PGD, sketching brings approximation error √Lnδ. To get better approximation solution in other words smaller δ, we should increase sketching dimension m such that m > w2 (XC∩Sm-1 ) δ2 . We also need to choose δ ≤ 1 - 2Ln 3κn(C) such that bound (11) is valid.

Example (Lasso): If the structured parameter θ∗ is s-sparse, we can use l1 norm ball as constraint. l1 constrained least squares is also known as Lasso [Tibshirani, 1996].

Let S ∈ Rm×n be a σ-sub Gaussian random matrix. Let θ∗ ∈ Rp be an s-sparse vector, and y = Xθ∗ + w. With θ0 = 0, we apply sketched PGD update. Let us assume the RSC and BD conditions hold. Take η = cmn Ln(√m+√n)2 , if m > c · s log p δ2 for some c > 0 [Banerjee et al., 2014; Pilanci and Wainwright, 2015], then starting from θ0 = 0, with high probability

Example (Low-Rank Matrix Regression): In low-rank matrix regression, the true parameter Θ∗ ∈ Rd×p is a rank r matrix with r ≪ min{d, p}. n observations can be accessed by the noisy linear model yi = 〈Xi,Θ∗ 〉 + wi, where 〈., .〉 denotes the entry wise inner product of two matrices, and Xi is a data point. One commonly used constraint is nuclear norm ball [Recht et al., 2010]. Denote ‖.‖∗ the nuclear norm, we set R(θ) = ‖θ‖∗. Let us assume the RE and BD conditions hold. Take η = cmn Ln(√m+√n)2 , if m > c · r(d+p) δ2 for some c > 0 [Chandrasekaran et al., 2012], then starting from Θ0 = 0, with high probability

5 General Loss Functions

In previous sections, we have discussed sketched iterative algorithms for constrained least squares. Unlike sketched least squares, Algorithm 1 does not solve a specific optimization problem. Instead, we find that if Algorithm 1 converges, the convergent point will be a solution to a system of nonlinear equations. The conclusion follows Lemma 3.

Lemma 3 Assume there is an η > 0 and a γ < 1 such that

for all θ1, θ2 ∈ Rp , then starting from arbitrary point, the iteration θt+1 = θt - ηXT ST S · dn(θt), has a limitation ˆ θ which satisfies XT ST Sdn(ˆ θ) = 0.

From Lemma 3, S-PGD approximates first order optimality condition XT dn(θ) = 0, which is the same as least squares.

If sketching dimension m is large enough, then ST S will be approximately identity, and we would expect ˆ θ and (2)'s optimum to be very close.

Our results requires the following assumptions:

Assumption 3 (Strong Convexity and Smoothness of Loss Function) All l(.; yi) are strongly convex and smooth. That is there is an α > 0 and a β > 0 such that α ≤ l′′ (xTi θ; yi) ≤ β, for any i ∈ [n], θ ∈ K.

Assumption 4 There is a 0 < δ1 < 1 such that for any u ∈ XC, and diagonal matrix D = diag(d) with d ∈ [ αβ , 1]n

Assumption 5 There is a δ2 > 0 such that for any u ∈ XC ∩ Bn, v ∈ span(X) ∩ Bn, and diagonal matrix D = diag(d) with d ∈ [ αβ , 1]n , where span(X) is the subspace spanned by the columns of X

Assumptions 4 and 5 implies that matrix ST S is close to identity in some scale cones. From literature [Pilanci and Wainwright, 2017; Pilanci and Wainwright, 2015; Woodruff, 2014], sketch dimension m determines the size of δ1 and δ2 in assumptions 5, 4. Assumption 4 also indicates that m should be in proportional to loss function condition number βα .

We then characterize the convergence of S-PGD.

Theorem 4 Let θ∗ be an arbitrary vector, and y is generated according to (1). With θ0 = 0, we apply S-PGD. Let us assume the RE and BD conditions hold for X. Set learning rate η = α(1-δ2)κn(C) β2 (1+δ2)2L2n , if assumptions 3, 5 , and 4 hold, we have at step t + 1

In Theorem 4 we use the same RE and BD conditions as least squares (Theorem 1). Note that in general case S-PGD convergence rate ρ(C; S) = √1- κ2n (C) L2n is slightly worse than ρ(C) = 1 - κn(C) Ln of least squares (Theorem 1) .

Example (Structured Logistic Regression): In structured logistic regression, its loss function is logistic loss

where label yi ∈ {+1,-1}. If the structured parameter θ∗ can be captured by norm R(.), then the optimization problem is posed as following

The convexity and smoothness constant in Assumption 3 can be bounded by the following conclusion

Lemma 5 When l(xTi θ; yi) is the logistic loss (16), we have α ≤ l′′ (xTi θ; yi) ≤ 14 , where α = ((1 + e-MR(θ∗ ) )(1 + eMR(θ∗ ) ))-1 for all i and θ ∈ K, where M = maxi R∗ (xi), and R∗ (.) is the dual norm of R(.).

Let S ∈ Rm×n be a sketching matrix, the main step of S-PGD for (17) is θt+1 = PR(θ)≤R(θ∗ )(θt - ηn (SX)T S · dn(θt)), where dn,i(θt) = -yi 1+eyixTi θt . Let us assume there exists a RSC constant κn(C) > 0 and smoothness constant Ln > 0, where error cone C = cone{∆ : R(θ∗ + ∆) ≤ R(θ∗ )}. Take η = α(1-δ1)κn(C) β2 (1+δ2)2L2n , then starting from θ0 = 0,

6 Comparison with SGD

In this section, we compare S-PGD with SGD. SGD can be treated a special case of S-PGD. We introduce leverage score and discuss how different types of sketching affects the leverage score. Finally, we compare the performance of SGD and Hadamard sketching using leverage score.

Let S = √ nm P · In, where In ∈ Rn is the identity matrix, then one iteration of SGD corresponds to S-PGD using the S defined above. In each iteration, SGD sample a new sketching matrix S in every iteration, while S-PGD only sample one S.

Let us consider ROS without re-randomization: S = √ nm P · H, we show that the performance of different S can be measured by leverage score. For S-PGD (1), the ith leverage score [Mahoney, 2011] is defined as pi(θ) = ‖hTi X‖22 (hTi d(θ))2 ∑nj=1 ‖hTj X‖22 (hTj d(θ))2 , where hi is the i-th row of H. Let µ(θ) = n · maxi pi(θ), we have the following approximation bound from [Mahoney, 2011; Drineas et al., 2006]

Theorem 6 Let δ ∈ (0, 1), with probability at least 1-δ, the difference between gradient and sketching gradient can be bounded by

From Theorem 6, the performance of SGD depends on µ(θ). When d(θ) is close to sparse, µ(θ) is large and sketching update XT ST Sd(θ) will not work. We set H to be the Hadamard matrix and re-randomize S using D, then XT ST Sd(θ) will provide a desired stochastic gradient,

7 Experimental Results

In this section, we show the experimental results of our algorithms on synthetic dataset, and how the choice of m affects computational efficiency and statistical guarantee.

In our first set of experiments, we show how m affects the performance of our algorithm using least squares. We draw design matrix X ∈ Rn×p randomly from Gaussian distribution. We choose parameter θ∗ to be an s-sparse vector, non-zero entries of θ∗ are drawn from standard Gaussian distribution. Response y is given by y = Xθ∗ + σ · w, where σ > 0 is a constant and w is drawn from standard Gaussian N(0, 1). We sample sketching matrix S fromN(0, 1 √m ) elementwisely.

In the first experiment, we choose n = 131, 072(217 ), p = 2000, m = 1000, and s = 10. We plot mean running time for each iteration versus mean estimation error ‖θt - θ∗ ‖2 and error bar of estimation error in Figure 1. We can see that both algorithms converge linearly. In the second experiment we choose n = 20000, p = 5000, s = 10, and m ∈ {20000, 18000, 16000, 14000, 12000, 10000}. For each m we run S-PGD 10 times. We stop our algorithm when relative error ‖θt+1-θt‖2 ‖θt‖2 is smaller than 0.01. We plot ratio mn versus estimation error ‖ˆ θ - θ∗ ‖2 and running time for each iteration, estimation error in Figure 2 left. We can see that as m increases, S-PGD converges slower but the output becomes more accurate.

In our second set of experiments, we compare RHS based S-PGD with PGD for binary logistic regression. Each data point (yi, xi) is generated by xi ∼ N(0, Ip), P(yi|xTi θ∗ ) = 1 1+e-yixTi θ∗ . Let p = 1000, n = 1, 048, 576(220 ), and m = 2000. We run 900 iterations for both algorithms. We collect the running time and estimation error every 10 iterations. The experimental results are shown in Figure 3. We can see that PGD takes less iterations to reach a smaller error, but in about the first 600 seconds S-PGD outperforms PGD. In summary, PGD has faster convergence rate but slower per iteration computation than S-PGD. S-PGD do more iterations in less time than PGD.

In our third set of experiments, we compare RHS based S-PGD with SGD for binary logistic regression. We choose n = 131, 072(217 ), p = 1000, and m = 1000. Firstly, we sample uniformly from the unit sphere Sp-1 ⊂ Rp a vector θ0. Denote S+p-1 a sphere gap such that for all u ∈ S+p-1 the angle between θ0 and u is in [0, π3 ]. Let y ∈ {1,-1}n be the label vector and X ∈ Rn×p be the data matrix. We set yi = 1 and sample xi uniformly from S+p-1 if i = 2, . . . , 65, 536(216 ); We set yi = -1 if i = 65, 538, . . . , 131, 072 and sample xi uniformly from -S+p-1 = {v : -v ∈ S+p-1 }. Then we set y1 = 1, y65537 = -1 and randomly sample x1 and x65537 such that xi = 1000 · gi, i = 1, 65, 537 where gi ∈ Rp is a random vector that satisfies gT i θ0 = 0. Finally, we constructed a dataset with two outliers. We set sketching dimension m = 1 and initialize both algorithms by w0. We compare the performance of S-PGD and SGD in terms of Ln(θ) and iterations. We run both S-PGD and SGD 100 times, and plot the average of training loss versus number of iterations.

From Figure 2 right, the loss function trained by S-PGD keeps going down while SGD stops, which implies that SGD failed to find an effective sample in fixed iterations. Vector d(θ0) is close to sparse since only two samples are misclassified, and µ(θ0) is large because the misclassified samples have greater ‖x‖2. RHS based S-PGD finds descent directions and move towards the real optimum.

8 Conclusion

In this paper we discussed sketching based iterative algorithms for generalized linear models. We show that under proper assumptions, these algorithms converge linearly to an error bound. The rate of convergence and error bound is determined by size of random projection and property of design matrix. In future works we want to explore if there exists similar results for nonlinear problems.

Acknowledgments

The research was also supported by NSF grants IIS-1563950, IIS-1447566, IIS-1447574, IIS-1422557, CCF-1451986, CNS- 1314560, IIS-0953274, IIS-1029711, NASA grant NNX12AQ39A, and gifts from Adobe, IBM, and Yahoo.