--- title: "Sketched Iterative Algorithms for Structured Generalized Linear Models" authors: - "Qilong Gu" - "Arindam Banerjee" doi: "10.24963/ijcai.2019/332" url: "https://ijcai.org/proceedings/2019/332" pdf: "https://ijcai.org/proceedings/2019/0332.pdf" venue: "Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence" year: 2019 pages: "2392-2398" track: "Main track" keywords: - "Machine Learning: Dimensionality Reduction and Manifold Learning" - "Machine Learning: Online Learning" --- # Sketched Iterative Algorithms for Structured Generalized Linear Models Qilong Gu, Arindam Banerjee *Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence*, pages 2392-2398. DOI: [10.24963/ijcai.2019/332](https://doi.org/10.24963/ijcai.2019/332) ## Abstract Recent years have seen advances in optimizing large scale statistical estimation problems. In statistical learning settings iterative optimization algorithms have been shown to enjoy geometric convergence. While powerful, such results only hold for the original dataset, and may face computational challenges when the sample size is large. In this paper, we study sketched iterative algorithms, in particular sketched-PGD (projected gradient descent) and sketched-SVRG (stochastic variance reduced gradient) for structured generalized linear model, and illustrate that these methods continue to have geometric convergence to the statistical error under suitable assumptions. Moreover, the sketching dimension is allowed to be even smaller than the ambient dimension, thus can lead to significant speed-ups. The sketched iterative algorithms introduced provide an additional dimension to study the trade-offs between statistical accuracy and time. ## Full text ## 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 ``` Algorithm 2: Sketched Stochastic Variance Reduction Gradient (S-SVRG) 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. Figure 1: Plot of running time versus estimation error ‖θt - θ∗ ‖2 for PGD (left) and SVRG (right). We set sample size n to be 217 and dimension p to be 2000. We set sketching dimension m to be 1000. Estimation error is presented in log scale. We can see that both algorithms converge linearly. ![Image](https://ijcai.org/proceedings/2019/fulltext/332/image_000000_ae19430ae21e3d6e4d045f719da13ce48d8387cdae9c90a398ed415a184ef719.png) Figure 2: (left) Dimension reduction ratio mn versus running time (left, blue) and statistical error ‖θt -θ∗ ‖2 of least squares (right, red) using S-PGD. As we increase the size of m, S-PGD uses more time to converge but the output becomes more accurate. (right) Comparison of S-PGD and SGD on our designed synthetic dataset start from special initial point. S-PGD converges while SGD increases the training loss. ![Image](https://ijcai.org/proceedings/2019/fulltext/332/image_000001_feb2ae0d7a3e4f5be902181457e262ae4eb3439014a9870c40eb95fbfd3bf036.png) 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. Figure 3: Convergence comparison of PGD and S-PGD. For both algorithms we plot running time versus estimation error ‖θt - θ∗ ‖2 (left) and number of iterations versus estimation error (right). We set sample size n to be 220 and dimension p to be 1000. We set sketching dimension m to be 2000. We can see from the left that in about the first 600 seconds S-PGD achieves lower estimation error. By comparing the result of both left and right we can see that S-PGD has faster per iteration computation. ![Image](https://ijcai.org/proceedings/2019/fulltext/332/image_000002_d61f2a64153c0a289d784bb95cc4184a6e76f9c9d877598c5f211598a3918ace.png) 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. ## References 1. [Agarwal et al., 2012] Alekh Agarwal, Sahand Negahban, and Martin J. Wainwright. Fast global convergence of gradient methods for high-dimensional statistical recovery. The Annals of Statistics, 40(5):2452-2482, 10 2012. 2. [Ailon and Liberty, 2008] Nir Ailon and Edo Liberty. Fast dimension reduction using Rademacher series on dual BCH codes. Discrete & Computational Geometry, 42(4):615, Sep 2008. 3. [Argyriou et al., 2012] Andreas Argyriou, Rina Foygel, and Nathan Srebro. Sparse Prediction with the k-Support Norm. In NIPS, pages 1457-1465, 2012. 4. [Banerjee et al., 2014] Arindam Banerjee, Sheng Chen, Farideh Fazayeli, and Vidyashankar Sivakumar. Estimation with Norm Regularization. In NIPS, pages 1556-1564, 2014. 5. [Bertsekas, 2010] D. P. Bertsekas. Nonlinear Programming. Athena Scientific, 2010. 6. [Chandrasekaran et al., 2012] Venkat Chandrasekaran, Ben- jamin Recht, Pablo A. Parrilo, and Alan S. Willsky. The Convex Geometry of Linear Inverse Problems. Foundations of Computational Mathematics, 12(6):805-849, 2012. 7. [Defazio et al., 2014] Aaron Defazio, Francis Bach, and Si- mon Lacoste-julien. Saga: A fast incremental gradient method with support for non-strongly convex composite objectives. In NIPS, pages 1646-1654, 2014. 8. [Drineas et al., 2006] P. Drineas, R. Kannan, and M. Ma- honey. Fast monte carlo algorithms for matrices i: Approxi- mating matrix multiplication. SIAM Journal on Computing, 36(1):132-157, 2006. 9. [Johnson and Zhang, 2013] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In NIPS, pages 315-323, 2013. 10. [Li et al., 2016] Xingguo Li, Tuo Zhao, Raman Arora, Han Liu, and Jarvis Haupt. Stochastic variance reduced opti- mization for nonconvex sparse learning. In ICML, pages 917-925, 2016. 11. [Mahoney, 2011] Michael W. Mahoney. Randomized algo- rithms for matrices and data. Foundations and Trends® in Machine Learning, 3(2):123-224, 2011. 12. [Negahban et al., 2012] Sahand N. Negahban, Pradeep Ravikumar, Martin J. Wainwright, and Bin Yu. A Uni- fied Framework for High-Dimensional Analysis of M - Estimators with Decomposable Regularizers. Statistical Science, 27(4):538-557, 2012. 13. [Oymak et al., 2018] Samet Oymak, Benjamin Recht, and Mahdi Soltanolkotabi. Sharp time-data tradeoffs for lin- ear inverse problems. IEEE Transactions on Information Theory, pages 4129-4158, 2018. 14. [Pilanci and Wainwright, 2015] Mert Pilanci and Martin J. Wainwright. Randomized sketches of convex programs with sharp guarantees. IEEE Transactions on Information Theory, pages 5096-5115, 2015. 15. [Pilanci and Wainwright, 2016] Mert Pilanci and Martin J. Wainwright. Iterative Hessian sketch: Fast and accurate so- lution approximation for constrained least-squares. Journal of Machine Learning Research (JMLR), 17(1):1842-1879, January 2016. 16. [Pilanci and Wainwright, 2017] Mert Pilanci and Martin J. Wainwright. Newton sketch: A linear-time optimization algorithm with linear-quadratic convergence. SIAM Journal of Optimization, 27(1):205-245, 2017. 17. [Qu et al., 2017] Chao Qu, Yan Li, and Huan Xu. Lin- ear Convergence of SVRG in Statistical Estimation. arXiv:1611.01957, 2017. 18. [Raskutti et al., 2010] Garvesh Raskutti, Martin J. Wain- wright, and Bin Yu. Restricted Eigenvalue Properties for Correlated Gaussian Designs. Journal of Machine Learning Research (JMLR), 11:2241-2259, 2010. 19. [Recht et al., 2010] Benjamin Recht, Maryam Fazel, and Pablo A. Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM Review, 52(3):471-501, 2010. 20. [Robbins and Monro, 1951] Herbert Robbins and Sutton Monro. A stochastic approximation method. The Annals of Mathematical Statistics, 22(3):400-407, 09 1951. 21. [Schmidt et al., 2017] Mark Schmidt, Nicolas Le Roux, and Francis Bach. Minimizing finite sums with the stochastic average gradient. Mathematical Programming, 162(1):83- 112, Mar 2017. 22. [Tang et al., 2017] Junqi Tang, Mohammad Golbabaee, and Mike E. Davies. Gradient projection iterative sketch for large-scale constrained least-squares. In ICML, pages 3377- 3386, 2017. 23. [Tibshirani, 1996] Robert Tibshirani. Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society, 58(1):267-288, 1996. 24. [Vershynin, 2012] Roman Vershynin. Introduction to the non- asymptotic analysis of random matrices. In Compressed Sensing, pages 210-268. Cambridge University Press, Cam- bridge, nov 2012. 25. [Vershynin, 2015] Roman Vershynin. Estimation in High Di- mensions: A Geometric Perspective, pages 3-66. Springer International Publishing, Cham, 2015. 26. [Wang et al., 2017] Jialei Wang, Jason Lee, Mehrdad Mah- davi, Mladen Kolar, and Nati Srebro. Sketching Meets Random Projection in the Dual: A Provable Recovery Al- gorithm for Big and High-dimensional Data. In AISTATS, pages 1150-1158, 2017. 27. [Woodruff, 2014] David P. Woodruff. Sketching as a tool for numerical linear algebra. Foundations and Trends® in Theoretical Computer Science, 10(1-2):1-157, 2014. 28. [Zhang et al., 2013] Lijun Zhang, Mehrdad Mahdavi, Rong Jin, Tianbao Yang, and Shenghuo Zhu. Recovering the Optimal Solution by Dual Random Projection. In COLT, pages 135-157. 2013.