High-Dimensional Optimization in Adaptive Random Subspaces
Jonathan Lacotte Mert Pilanci Marco Pavone











Quiescent galaxies at z ≳ 2.5: observations vs. models
    Giacomo Girelli
    December 30, 2023
======================================================























    We propose a new randomized optimization method for high-dimensional problems which can be seen as a generalization of coordinate descent to random subspaces. We show that an adaptive sampling strategy for the random subspace significantly outperforms the oblivious sampling method, which is the common choice in the recent literature. The adaptive subspace can be efficiently generated by a correlated random matrix ensemble whose statistics mimic the input data. We prove that the improvement in the relative error of the solution can be tightly characterized in terms of the spectrum of the data matrix, and provide probabilistic upper-bounds. We then illustrate the consequences of our theory with data matrices of different spectral decay. Extensive experimental results show that the proposed approach offers significant speed ups in machine learning problems including logistic regression, kernel classification with random convolution layers and shallow neural networks with rectified linear units. Our analysis is based on convex analysis and Fenchel duality, and establishes connections to sketching and randomized matrix decomposition. 






§ INTRODUCTION

Random Fourier features, Nystrom method and sketching techniques have been successful in large scale machine learning problems. The common practice is to employ oblivious sampling or sketching matrices, which are typically randomized and fixed ahead of the time. However, it is not clear whether one can do better by adapting the sketching matrices to data. In this paper, we show that adaptive sketching matrices can significantly improve the approximation quality. We characterize the approximation error on the optimal solution in terms of the smoothness of the function, and spectral properties of the data matrix. 

Many machine learning problems end up being high dimensional optimization problems, which typically follow from forming the kernel matrix of a large dataset or mapping the data trough a high dimensional feature map, such as random Fourier features <cit.> or convolutional neural networks <cit.>. Such high dimensional representations induce higher computational and memory complexities, and result in slower training of the models. Random projections are a classical way of performing dimensionality reduction, and are widely used in many algorithmic contexts <cit.>. Nevertheless, only recently these methods have captured great attention as an effective way of performing dimensionality reduction in convex optimization. In the context of solving a linear system Ax=b and least-squares optimization, the authors of <cit.> propose a randomized iterative method with linear convergence rate, which, at each iteration, performs a proximal update x^(k+1)=_x ∈ Tx - x^(k)_2^2, where the next iterate x^(k+1) is restricted to lie within an affine subspace T = x^(k) + range(A^⊤ S), and S is a n × m dimension-reduction matrix with m min{n, d}. In the context of kernel ridge regression, the authors of <cit.> propose to approximate the n-dimensional kernel matrix by sketching its columns to a lower m-dimensional subspace, chosen uniformly at random. From the low dimensional kernel ridge solution α^* ∈ℝ^m, they show how to reconstruct an approximation x∈ℝ^n of the high dimensional solution x^* ∈ℝ^n. Provided that the sketching dimension m is large enough – as measured by the spectral properties of the kernel matrix K –, the estimate x retains some statistical properties of x^*, e.g., minimaxity. Similarly, in the broader context of classification through convex loss functions, the authors of <cit.> propose to project the d-dimensional features of a given data matrix A to a lower m-dimensional subspace, chosen independently of the data. After computing the optimal low-dimensional classifier α^* ∈ℝ^m, their algorithm returns an estimate x∈ℝ^d of the optimal classifier x^* ∈ℝ^d. Even though they provide formal guarantees on the estimation error x-x^*_2, their results rely on several restrictive assumptions, that is, the data matrix A must be low rank, or, the classifier x^* must lie in the span of the top few left singular vectors of A. Further, random subspace optimization has also been explored for large-scale trust region problems <cit.>, also using a subspace chosen uniformly at random. Our proposed approach draws connections with the Gaussian Kaczmarz method proposed in <cit.> and the kernel sketching method in <cit.>. Differently, we are interested in smooth, convex optimization problems with ridge regularization. In contrast to <cit.>, we do not make any assumption on the optimal solution x^*.

Our work relates to the considerable amount of literature on randomized approximations of high dimensional kernel matrices K. The typical approach consists of building a low-rank factorization of the matrix K, using a random subset of its columns <cit.>. The so-called Nystrom method has proven to be effective empirically <cit.>, and many research efforts have been devoted to improving and analyzing the performance of its many variants (e.g., uniform column sub-sampling, leverage-score based sampling), especially in the context of kernel ridge regression <cit.>. In a related vein, sketching methods have been proposed to reduce the space complexity of storing high-dimensional data matrices <cit.>, by projecting their rows to a randomly chosen lower dimensional subspace. Our theoretical findings build on known results for low-rank factorization of positive semi-definite (p.s.d.) matrices <cit.>, and show intimate connections with kernel matrices sketching <cit.>. Lastly, our problem setting also draws connections with compressed sensing <cit.> where the goal is to recover a high dimensional structured signal from a small number of randomized, and usually oblivious, measurements. 















 §.§ Contributions

In this work, we propose a novel randomized subspace optimization method with strong solution approximation guarantees which outperform oblivious sampling methods. We derive probabilistic bounds on the error of approximation for general convex functions. We show that our method provides a significant improvement over the oblivious version, and theoretically quantify this observation as function of the spectral properties of the data matrix. We also introduce an iterative version of our method, which converges to the optimal solution by iterative refinement.









 §.§ An overview of our results








Let f : ℝ^n →ℝ be a convex and μ-strongly smooth function, i.e., ∇^2 f(w) ≼μ I_n for all w ∈ℝ^n, and A ∈ℝ^n × d a high-dimensional matrix. We are interested in solving the primal problem


    x^* = _x ∈ℝ^d f(Ax) + λ/2x_2^2 ,


Given a random matrix S ∈ℝ^d × m with m ≪ d, we consider instead the sketched primal problem


    α^* ∈_α∈ℝ^m f(ASα) + λ/2α^⊤ S^⊤ S α ,


where we effectively restrict the optimization domain to a lower m-dimensional subspace. In this work, we explore the following questions: How can we estimate the original solution x^* given the sketched solution α^*? Is a uniformly random subspace the optimal choice, e.g., S∼ Gaussian i.i.d.? Or, can we come up with an adaptive sampling distribution that is related to the matrix A, which yields stronger guarantees? 

By Fenchel duality analysis, we exhibit a natural candidate for an approximate solution to x^*, given by x = -λ^-1 A^⊤∇ f(ASα^*). Our main result (Section <ref>) establishes that, for an adaptive sketching matrix of the form S = A^⊤S where S is typically Gaussian i.i.d., the relative error satisfies a high-probability guarantee of the form x-x^*_2/x^*_2 ε, with ε < 1. Our error bound ε depends on the smoothness parameter μ, the regularization parameter λ, the shape of the domain of the Fenchel conjugate f^* and the spectral decay of the matrix A. Further, we show that this error can be explicitly controlled in terms of the singular values of A, and we derive concrete bounds for several standard spectral profiles, which arise in data analysis and machine learning. In particular, we show that using the adaptive matrix S=A^⊤S provides much stronger guarantees than oblivious sketching, where S is independent of A. Then, we take advantage of the error contraction (i.e., ε < 1), and extend our adaptive sketching scheme to an iterative version (Section <ref>), which, after T iterations, returns a higher precision estimate x^(T) that satisfies x^(T)-x^*_2 / x^*_2 ε^T. Throughout this work, we specialize our formal guarantees and empirical evaluations (Section <ref>) to Gaussian matrices S, which is a standard choice and yields the tightest error bounds. However, our approach extends to a broader class of matrices S, such as Rademacher matrices, sub-sampled randomized Fourier (SRFT) or Hadamard (SRHT) transforms, and column sub-sampling matrices. Thus, it provides a general framework for random subspace optimization with strong solution guarantees.










§ CONVEX OPTIMIZATION IN ADAPTIVE RANDOM SUBSPACES











We introduce the Fenchel conjugate of f, defined as f^*(z) := sup_w ∈ℝ^n{ w^⊤ z - f(w) }, which is convex and its domain domf^* := {z ∈ℝ^n | f^*(z) < +∞} is a closed, convex set. Our control of the relative error x - x^*_2 / x^*_2 is closely tied to controlling a distance between the respective solutions of the dual problems of (<ref>) and (<ref>). The proof of the next two Propositions follow from standard convex analysis arguments <cit.>, and are deferred to Appendix <ref>.


	
	Under the previous assumptions on f, it holds that
	
	
    min_x f(Ax) + λ/2x_2^2 = max_z - f^*(z) - 1/2λA^⊤ z_2^2 .

	
	There exist an unique primal solution x^* and an unique dual solution z^*. Further, we have Ax^* ∈∂ f^*(z^*), z^* = ∇ f(Ax^*) and x^* = -1/λ A^⊤ z^*.
	



	
	Strong duality holds for the sketched program
	
	
    min_α f(ASα) + λ/2Sα_2^2 = max_y - f^*(y) - 1/2 λP_S A^⊤ y_2^2 ,

	
	where P_S = S(S^⊤ S)^† S^⊤ is the orthogonal projector onto the range of S. There exist a sketched primal solution α^* and an unique sketched dual solution y^*. Further, for any solution α^*, it holds that ASα^* ∈∂ f^*(y^*) and y^* = ∇ f(AS α^*).
	


We define the following deterministic functional Z_f which depends on f^*, the data matrix A and the sketching matrix S, and plays an important role in controlling the approximation error,


    Z_f ≡ Z_f(A, S) = sup_Δ∈ (domf^* - z^*)(Δ^⊤ A P_S^⊥ A^⊤Δ/Δ_2^2)^1/2,


where P_S^⊥ = I - P_S is the orthogonal projector onto range(S)^⊥. The relationship x^* = -λ^-1 A^⊤∇ f(Ax^*) suggests the point x = -λ^-1 A^⊤∇ f(ASα^*) as a candidate for approximating x^*. The Fenchel dual programs of (<ref>) and (<ref>) only differ in their quadratic regularization term, A^⊤ z_2^2 and P_S A^⊤ y_2^2, which difference is tied to the quantity P_S^⊥ A^⊤ (z-y)_2. As it holds that x-x^*_2 = λ^-1A^⊤(z^*-y^*)_2, we show that the error x - x^*_2 can be controlled in terms of the spectral norm P_S^⊥ A^⊤_2, or more sharply, in terms of the quantity Z_f, which satisfies Z_f P^⊥_S A^⊤_2. We formalize this statement in our next result, which proof is deferred to Appendix <ref>.


	
	Let α^* be any minimizer of the sketched program (<ref>). Then, under the condition λ 2 μ Z_f^2, we have
	
	
    x-x^*_2 √(μ/2λ) Z_f x^*_2 ,

	
	which further implies
	
	
    x-x^*_2 √(μ/2λ)P_S^⊥ A^⊤_2 x^*_2 .

	


For an adaptive sketching matrix S=A^⊤S, we rewrite P_S^⊥ A^⊤^2_2=K - K S (S^⊤ K S)^†S^⊤ K_2, where K = AA^⊤ is p.s.d. Combining our deterministic bound (<ref>) with known results <cit.> for randomized low-rank matrix factorization in the form KS (S^⊤ K S)^†S^⊤ K of p.s.d. matrices K, we can give guarantees with high probability (w.h.p.) on the relative error for various types of matrices S. For conciseness, we specialize our next result to adaptive Gaussian sketching, i.e., S Gaussian i.i.d. Given a target rank k  2, we introduce a measure of the spectral tail of A as R_k(A) = (σ^2_k + 1/k∑_j = k+1^ρσ^2_j)^1/2, where ρ is the rank of the matrix A and σ_1 σ_2 σ_ρ its singular values. The proof of the next result follows from a combination of Theorem <ref> and Corollary 10.9 in <cit.>, and is deferred to Appendix <ref>.


	
	Given k min(n,d)/2 and a sketching dimension m=2k, let S = A^⊤S, with S∈ℝ^n × m Gaussian i.i.d. Then, for some universal constant c_0  36, provided λ 2 μ c^2_0 R_k^2(A), it holds with probability at least 1-12e^-k that
	
	
    x - x^*_2  c_0 √(μ/2λ) R_k(A) x^*_2 .

	
	



	The quantity Z_f := sup_Δ∈ (domf^* - z^*)(Δ^⊤ A P_S^⊥ A^⊤Δ/Δ_2^2)^1/2 is the eigenvalue of the matrix P_S^⊥ A^⊤, restricted to the spherical cap 𝒦 := (dom f^* - z^*) ∩𝒮^n-1, where 𝒮^n-1 is the unit sphere in dimension n. Thus, depending on the geometry of 𝒦, the deterministic bound (<ref>) might be much tighter than (<ref>), and yield a probabilistic bound better than (<ref>). The investigation of such a result is left for future work.












 §.§ Theoretical predictions as a function of spectral decay











We study the theoretical predictions given by (<ref>) on the relative error, for different spectral decays of A and sketching methods, in particular, adaptive Gaussian sketching versus oblivious Gaussian sketching and leverage score column sub-sampling <cit.>. We denote ν_k = σ_k^2 the eigenvalues of A A^⊤. For conciseness, we absorb μ into the eigenvalues by setting ν_k ≡μν_k and μ≡ 1. This re-scaling leaves the right-hand side of the bound (<ref>) unchanged, and does not affect the analysis below. Then, we assume that ν_1= 𝒪(1), λ∈ (ν_ρ, ν_1), and λ→ 0 as n → +∞. These assumptions are standard in empirical risk minimization and kernel regression methods <cit.>, which we focus on in Sections <ref> and <ref>. We consider three decaying schemes of practical interest. The matrix A has either a finite-rank ρ, a κ-exponential decay where ν_j ∼ e^-κ j and κ > 0, or, a β-polynomial decay where ν_k ∼ j^-2β and β > 1/2. Among other examples, these decays are characteristic of various standard kernel functions, such as the polynomial, Gaussian and first-order Sobolev kernels <cit.>. Given a precision ε > 0 and a confidence level η∈ (0,1), we denote by m_A (resp. m_O, m_S) a sufficient dimension for which adaptive (resp. oblivious, leverage score) sketching yields the following (ε, η)-guarantee on the relative error. That is, with probability at least 1-η, it holds that x-x^*_2/x^*_2 ε. 

We determine m_A from our probabilistic regret bound (<ref>). For m_S, using our deterministic regret bound (<ref>), it then suffices to bound the spectral norm P_A^⊤S^⊥ A^⊤_2 in terms of the eigenvalues ν_k, when S is a leverage score column sub-sampling matrix. To the best of our knowledge, the tightest bound has been given by <cit.> (see Lemma 5). For m_O, we leverage results from <cit.>. The authors provide an upper bound on the relative error x-x^*_2/x^*_2, when S is Gaussian i.i.d. with variance 1/d. It should be noted that their sketched solution α^* is slightly different from ours. They solve α^* =  f(ASα) + (2λ)^-1α_2^2, whereas we do include the matrix S in the regularization term. One might wonder which regularizer works best when S is Gaussian i.i.d. Through extensive numerical simulations, we observed a strongly similar performance. Further, standard Gaussian concentration results yields that Sα^*_2^2 ≈α^*_2^2.

Our theoretical findings are summarized in Table <ref>, and we give the mathematical details of our derivations in Appendix <ref>. For the sake of clarity, we provide in Table <ref> lower bounds on the predicted values m_O and m_S, and, thus, lower bounds on the ratios m_O / m_A and m_S/m_A. Overall, adaptive Gaussian sketching provides stronger guarantees on the relative error x-x^*_2/x^*_2.



We illustrate numerically our predictions for adaptive Gaussian sketching versus oblivious Gaussian sketching. With n=1000 and d=2000, we generate matrices A^exp and A^poly, with spectral decay satisfying respectively ν_j ∼ n e^-0.1 j and ν_j ∼ n j^-2. First, we perform binary logistic regression, with f(Ax) = n^-1∑_i=1^n ℓ_y_i(a_i^⊤ x) where ℓ_y_i(z) = y_i log(1+e^-z) + (1-y_i) log(1+e^z), y ∈{0,1}^n and a_i is the i-th row of A. For the polynomial (resp. exponential) decay, we expect the relative error x-x^*_2 / x^*_2 to follow w.h.p. a decay proportional to m^-1 (resp. e^-0.05 m). Figure <ref> confirms those predictions. We repeat the same experiments with a second loss function, f(Ax) = (2n)^-1∑_i=1^n (a_i^⊤ x)^2_+ - 2(a_i^⊤ x) y_i. The latter is a convex relaxation of the penalty 1/2(Ax)_+ - y_2^2 for fitting a shallow neural network with a ReLU non-linearity. Again, Figure <ref> confirms our predictions, and we observe that the adaptive method performs much better than the oblivious sketch.








§ ALGORITHMS FOR ADAPTIVE SUBSPACE SKETCHING









 §.§ Numerical conditioning and generic algorithm






A standard quantity to characterize the capacity of a convex program to be solved efficiently is its condition number <cit.>, which, for the primal (<ref>) and (adaptive) sketched program (<ref>), are given by


    κ = λ + sup_xσ_1( A^⊤∇^2 f(Ax) A)/λ + inf_xσ_d( A^⊤∇^2 f(Ax) A),   κ_S = sup_ασ_1( S^⊤ A (λ I + A^⊤∇^2 f(A A^⊤Sα) A) A^⊤S)/inf_ασ_d( S^⊤ A ( λ I + A^⊤∇^2 f(A A^⊤Sα) A) A^⊤S).


The latter can be significantly larger than κ, up to κ_S ≈κσ_1(S^⊤ A A^⊤S)/σ_m(S^⊤ A A^⊤S)≫κ. A simple change of variable overcomes this issue. With A_S,† = AS (S^⊤ S)^-1/2, we solve instead the optimization problem


    α_†^* = _α_†∈ℝ^m f(A_S,†α_†) + λ/2α_†_2^2.


It holds that x = -λ^-1 A^⊤∇ f(A_S,†α_†^*). The additional complexity induced by this change of variables comes from computing the (square-root) pseudo-inverse of S^⊤ S, which requires 𝒪(m^3) flops via a singular value decomposition. When m is small, this additional computation is negligible and numerically stable, and the re-scaled sketched program (<ref>) is actually better conditioned that the original primal program (<ref>), as stated in the next result that we prove in Appendix <ref>.


	
	Under adaptive sketching, the condition number κ_† of the re-scaled sketched program (<ref>) satisfies κ_†κ with probability 1.




We observed a drastic practical performance improvement between solving the sketched program as formulated in (<ref>) and its well-conditioned version (<ref>).


If the chosen sketch dimension m is itself prohibitively large for computing the matrix Q = (S^⊤ S)^-1/2, one might consider a pre-conditioning matrix Q, which is faster to compute, and such that the matrix SQ is well-conditioned. Typically, one might compute a matrix Q based on an approximate singular value decomposition of the matrix S^⊤ S. Then, one solves the optimization problem α^*_Q = _α∈ℝ^m f(AS Qα) + λ/2S Q α_2^2. Provided that Q is invertible, it holds that x satisfies x = -λ^-1 A^⊤∇ f(ASQ α^*_Q). 












 §.§ Error contraction and almost exact recovery of the optimal solution






The estimate x satisfies a guarantee of the form x-x^*_2 εx^*_2 w.h.p., and, with ε < 1 provided that λ is large enough. Here, we extend Algorithm <ref> to an iterative version which takes advantage of this error contraction, and which is relevant when a high-precision estimate x is needed.




A key advantage is that, at each iteration, the same sketching matrix S is used. Thus, the sketched matrix A_S,† has to be computed only once, at the beginning of the procedure.  The output x^(T) satisfies the following recovery property, which empirical benefits are illustrated in Figure <ref>.


	
	After T iterations of Algorithm <ref>, provided that λ 2 μ Z_f^2, it holds that
	
	
    x^(T) - x^*_2 (μ Z_f^2/2λ)^T/2x^*_2 .

	
	Further, if S = A^⊤S where S∈ℝ^n × m with i.i.d. Gaussian entries and m = 2k for some target rank k  2, then, for some universal constant c_0  36, after T iterations of Algorithm <ref>, provided that λ 2 c_0^2 μ R_k^2(A), the approximate solution x^(T) satisfies with probability at least 1-12e^-k,
	
	
    x^(T) - x^*_2 (c_0^2 μ R_k^2(A)/2λ)^T/2x^*_2 .

	













	
	An immediate extension of Algorithms <ref> and <ref> consists in using the power method <cit.>. Given q ∈ℕ, one uses the sketching matrix S=(A^⊤ A)^q A^⊤S. The larger q, the smaller the approximation error AA^⊤ - AS(S^⊤ S)^† S^⊤ A^⊤_2 (see Corollary 10.10 in <cit.>). Of pratical interest are data matrices A with a spectral profile starting with a fast decay, and then becoming flat. This happens typically for A of the form A = A + W, where A has a fast decay and W is a noise matrix with, for instance, independent subgaussian rows <cit.>. Our results easily extend to this setting and we illustrate its empirical benefits in Figure <ref>.













§ APPLICATION TO EMPIRICAL RISK MINIMIZATION AND KERNEL METHODS












By the representer theorem, the primal program (<ref>) can be re-formulated as 


    w^* ∈_w ∈ℝ^n f(K w) + λ/2 w^⊤ K w ,


where K=AA^⊤. Clearly, it holds that x^* = A^⊤ w^*. Given a matrix S with i.i.d. Gaussian entries, we consider the sketched version of the kernelized primal program (<ref>), 


    α^* ∈_α∈ℝ^m f(K Sα) + λ/2α^⊤S^⊤ K Sα .


The sketched program (<ref>) is exactly our adaptive Gaussian sketched program (<ref>). Thus, setting w = -λ^-1∇ f(K Sα^*), it holds that x = A^⊤w. Since the relative error x-x^*_2/x^*_2 is controlled by the decay of the eigenvalues of K, so does the relative error A^⊤ (w-w^*)_2 / A^⊤ w^*_2. More generally, the latter statements are still true if K is any positive semi-definite matrix, and, if we replace A^⊤ by any square-root matrix of K. Here, we denote Z_f ≡ Z_f(K^1/2, K^1/2S) (see Eq. (<ref>)).


	
	Let K ∈ℝ^n × n be any positive semi-definite matrix. Let w^* be any minimizer of the kernel program (<ref>) and α^* be any minimizer of its sketched version (<ref>). Define the approximate solution w = -1/λ∇ f(K Sα^*). If λ 2 μ Z_f^2, then it holds that
	
	
    K^1/2(w-w^*)_2 √(μ/2 λ) Z_f K^1/2w^*_2 .

	


For a positive definite kernel k : ℝ^d ×ℝ^d →ℝ and a data matrix A = [a_1,,a_n]^⊤∈ℝ^n × d, let K be the empirical kernel matrix, with K_ij=k(a_i,a_j). Let φ(·) ∈ℝ^D be a random feature map <cit.>, such as random Fourier features or a random convolutional neural net. We are interested in the computational complexities of forming the sketched versions of the primal (<ref>), the kernel primal (<ref>) and the primal (<ref>) with φ(A) instead of A. We compare the complexities of adaptive and oblivious sketching and uniform column sub-sampling. Table <ref> shows that all three methods have similar complexities for computing AS and φ(A)S. Adaptive sketching exhibits an additional factor 2 that comes from computing the correlated sketching matrices S=A^⊤S and S=φ(A)^⊤S. In practice, the latter is negligible compared to the cost of forming φ(A) which, for instance, corresponds to a forward pass over the whole dataset in the case of a convolutional neural network. On the other hand, uniform column sub-sampling is significantly faster in order to form the sketched kernel matrix K S, which relates to the well-known computational advantages of kernel Nystrom methods <cit.>.














§ NUMERICAL EVALUATION OF ADAPTIVE GAUSSIAN SKETCHING












We evaluate Algorithm <ref> on MNIST and CIFAR10. First, we aim to show that the sketching dimension can be considerably smaller than the original dimension while retaining (almost) the same test classification accuracy. Second, we aim to get significant speed-ups in achieving a high-accuracy classifier. To solve the primal program (<ref>), we use two standard algorithms, stochastic gradient descent (SGD) with (best) fixed step size and stochastic variance reduction gradient (SVRG) <cit.> with (best) fixed step size and frequency update of the gradient correction. To solve the adaptive sketched program (<ref>), we use SGD, SVRG and the sub-sampled Newton method <cit.> – which we refer to as Sketch-SGD, Sketch-SVRG and Sketch-Newton. The latter is well-suited to the sketched program, as the low-dimensional Hessian matrix can be quickly inverted at each iteration. For both datasets, we use 50000 training and 10000 testing images. We transform each image using a random Fourier feature map φ(·) ∈ℝ^D, i.e., ⟨φ(a), φ(a^') ⟩≈exp(-γa-a^'^2_2 ) <cit.>. For MNIST and CIFAR10, we choose respectively D=10000 and γ=0.02, and, D=60000 and γ=0.002, so that the primal is respectively 10000-dimensional and 60000-dimensional. Then, we train a classifier via a sequence of binary logistic regressions – which allow for efficient computation of the Hessian and implementation of the Sketch-Newton algorithm –, using a one-vs-all procedure. 

First, we evaluate the test classification error of x. We solve to optimality the primal and sketched programs for values of λ∈{10^-4, 5 · 10^-5, 10^-5, 5 · 10^-6} and sketching dimensions m ∈{64, 128, 256, 512, 1024}. In Table <ref> are reported the results, which are averaged over 20 trials for MNIST and 10 trials for CIFAR10, and, empirical variances are reported in Appendix <ref>. Overall, the adaptive sketched program yields a high-accuracy classifier for most couples (λ, m). Further, we match the best primal classifier with values of m as small as 256 for MNIST and 512 for CIFAR10, which respectively corresponds to a dimension reduction by a factor ≈ 40 and ≈ 120. These results additionally suggest that adaptive Gaussian sketching introduces an implicit regularization effect, which might be related to the benefits of spectral cutoff estimators. For instance, on CIFAR10, using λ=10^-5 and m=512, we obtain an improvement in test accuracy by more than 2% compared to x^*. Further, over some sketching dimension threshold under which the performance is bad, as the value of m increases, the test classification error of x increases to that of x^*, until matching it.



Further, we evaluate the test classification error of two sketching baselines, that is, oblivious Gaussian sketching for which the matrix S has i.i.d. Gaussian entries, and, adaptive column sub-sampling (Nystrom method) for which S=A^⊤S with S a column sub-sampling matrix. As reported in Table <ref>, adaptive Gaussian sketching performs better for a wide range of values of sketching size m and regularization parameter λ.



Then, we compare the test classification error versus wall-clock time of the optimization algorithms mentioned above. Figure <ref> shows results for some values of m and λ. We observe some speed-ups on the 10000-dimensional MNIST problem, in particular for Sketch-SGD and for Sketch-SVRG, for which computing the gradient correction is relatively fast. Such speed-ups are even more significant on the 60000-dimensional CIFAR10 problem, especially for Sketch-Newton. A few iterations of Sketch-Newton suffice to almost reach the minimum x, with a per-iteration time which is relatively small thanks to dimensionality reduction. Hence, it is more than 10 times faster to reach the best test accuracy using the sketched program. In addition to random Fourier features mapping, we carry out another set of experiments with the CIFAR10 dataset, in which we pre-process the images. That is, similarly to <cit.>, we map each image through a random convolutional layer. Then, we kernelize these processed images using a Gaussian kernel with γ=2 · 10^-5. Using our implementation, the best test accuracy of the kernel primal program (<ref>) we obtained is 73.1%. Sketch-SGD, Sketch-SVRG and Sketch-Newton – applied to the sketched kernel program (<ref>) – match this test accuracy, with significant speed-ups, as reported in Figure <ref>.













§ ACKNOWLEDGEMENTS

This work was partially supported by the National Science Foundation under grant IIS-1838179 and Office of Naval Research, ONR YIP Program, under contract N00014-17-1-2433.









10
	
	rahimi2008random
	Ali Rahimi and Benjamin Recht.
	Random features for large-scale kernel machines.
	In Advances in neural information processing systems, pages
	1177–1184, 2008.
	
	goodfellow2016deep
	Ian Goodfellow, Yoshua Bengio, and Aaron Courville.
	Deep learning.
	MIT press, 2016.
	
	vempala2005random
	Santosh S Vempala.
	The random projection method, volume 65.
	American Mathematical Soc., 2005.
	
	yang2017randomized
	Yun Yang, Mert Pilanci, Martin J Wainwright, et al.
	Randomized sketches for kernels: Fast and optimal nonparametric
	regression.
	The Annals of Statistics, 45(3):991–1023, 2017.
	
	zhang2013recovering
	Lijun Zhang, Mehrdad Mahdavi, Rong Jin, Tianbao Yang, and Shenghuo Zhu.
	Recovering the optimal solution by dual random projection.
	In Conference on Learning Theory, pages 135–157, 2013.
	
	zhang2014random
	Lijun Zhang, Mehrdad Mahdavi, Rong Jin, Tianbao Yang, and Shenghuo Zhu.
	Random projections for classification: A recovery approach.
	IEEE Transactions on Information Theory, 60(11):7300–7316,
	2014.
	
	vu2017random
	Ky Vu, Pierre-Louis Poirion, Claudia D'Ambrosio, and Leo Liberti.
	Random projections for trust region subproblems.
	arXiv preprint arXiv:1706.02730, 2017.
	
	williams2001using
	Christopher KI Williams and Matthias Seeger.
	Using the nyström method to speed up kernel machines.
	In Advances in neural information processing systems, pages
	682–688, 2001.
	
	smola2004tutorial
	Alex J Smola and Bernhard Schölkopf.
	A tutorial on support vector regression.
	Statistics and computing, 14(3):199–222, 2004.
	
	drineas2005nystrom
	Petros Drineas and Michael W Mahoney.
	On the nyström method for approximating a gram matrix for
	improved kernel-based learning.
	journal of machine learning research, 6(Dec):2153–2175, 2005.
	
	kumar2012sampling
	Sanjiv Kumar, Mehryar Mohri, and Ameet Talwalkar.
	Sampling methods for the nyström method.
	Journal of Machine Learning Research, 13(Apr):981–1006, 2012.
	
	yang2012nystrom
	Tianbao Yang, Yu-Feng Li, Mehrdad Mahdavi, Rong Jin, and Zhi-Hua Zhou.
	Nyström method vs random fourier features: A theoretical and
	empirical comparison.
	In Advances in neural information processing systems, pages
	476–484, 2012.
	
	alaoui2015fast
	Ahmed Alaoui and Michael W Mahoney.
	Fast randomized kernel ridge regression with statistical guarantees.
	In Advances in Neural Information Processing Systems, pages
	775–783, 2015.
	
	bach2013sharp
	Francis Bach.
	Sharp analysis of low-rank kernel matrix approximations.
	In Conference on Learning Theory, pages 185–209, 2013.
	
	pilanci2015randomized
	Mert Pilanci and Martin J Wainwright.
	Randomized sketches of convex programs with sharp guarantees.
	IEEE Transactions on Information Theory, 61(9):5096–5115,
	2015.
	
	halko2011
	Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp.
	Finding structure with randomness: Probabilistic algorithms for
	constructing approximate matrix decompositions.
	SIAM review, 53(2):217–288, 2011.
	
	boutsidis2013
	Christos Boutsidis and Alex Gittens.
	Improved matrix algorithms via the subsampled randomized hadamard
	transform.
	SIAM Journal on Matrix Analysis and Applications,
	34(3):1301–1340, 2013.
	
	gittens2016revisiting
	Alex Gittens and Michael W Mahoney.
	Revisiting the nyström method for improved large-scale machine
	learning.
	The Journal of Machine Learning Research, 17(1):3977–4041,
	2016.
	
	witten2015randomized
	Rafi Witten and Emmanuel Candes.
	Randomized algorithms for low-rank matrix factorizations: sharp
	performance bounds.
	Algorithmica, 72(1):264–281, 2015.
	
	candes2006stable
	Emmanuel J Candes, Justin K Romberg, and Terence Tao.
	Stable signal recovery from incomplete and inaccurate measurements.
	Communications on Pure and Applied Mathematics: A Journal Issued
		by the Courant Institute of Mathematical Sciences, 59(8):1207–1223, 2006.
	
	rockafellar2015convex
	Ralph Tyrell Rockafellar.
	Convex analysis.
	Princeton university press, 2015.
	
	friedman2001elements
	Jerome Friedman, Trevor Hastie, and Robert Tibshirani.
	The elements of statistical learning, volume 1.
	Springer, 2001.
	
	kernelprobastats
	A. Berlinet and C. Thomas-Agnan.
	Reproducing kernel Hilbert spaces in probability and
		statistics.
	Kluwer Adademic, 2004.
	
	boyd2004convex
	Stephen Boyd and Lieven Vandenberghe.
	Convex optimization.
	Cambridge university press, 2004.
	
	coates2012learning
	Adam Coates and Andrew Y Ng.
	Learning feature representations with k-means.
	In Neural networks: Tricks of the trade, pages 561–580.
	Springer, 2012.
	
	johnson2013accelerating
	Rie Johnson and Tong Zhang.
	Accelerating stochastic gradient descent using predictive variance
	reduction.
	In Advances in neural information processing systems, pages
	315–323, 2013.
	
	bollapragada2018exact
	Raghu Bollapragada, Richard H Byrd, and Jorge Nocedal.
	Exact and inexact subsampled newton methods for optimization.
	IMA Journal of Numerical Analysis, 39(2):545–578, 2018.
	
	erdogdu2015convergence
	Murat A Erdogdu and Andrea Montanari.
	Convergence rates of sub-sampled newton methods.
	arXiv preprint arXiv:1508.02810, 2015.
	
	scikit-learn
	F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel,
	M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos,
	D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.
	Scikit-learn: Machine learning in Python.
	Journal of Machine Learning Research, 12:2825–2830, 2011.
	
	zhang2016understanding
	Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.
	Understanding deep learning requires rethinking generalization.
	arXiv preprint arXiv:1611.03530, 2016.
	
	tu2016large
	Stephen Tu, Rebecca Roelofs, Shivaram Venkataraman, and Benjamin Recht.
	Large scale kernel learning using block coordinate descent.
	arXiv preprint arXiv:1602.05310, 2016.
	
	recht2018cifar
	Benjamin Recht, Rebecca Roelofs, Ludwig Schmidt, and Vaishaal Shankar.
	Do cifar-10 classifiers generalize to cifar-10?
	arXiv preprint arXiv:1806.00451, 2018.

	vershynin2018high
	Roman Vershynin.
	High-dimensional probability: An introduction with applications in data science, volume 47.
	Cambridge University Press, 2018.
	
	gower2015randomized
	Robert Gower and Peter Richtárik.
	Randomized iterative methods for linear systems.
	SIAM Journal on Matrix Analysis and Applications, 36(4):1660-1690, 2015.
	
















§ ADDITIONAL EXPERIMENTAL RESULTS AND IMPLEMENTATION DETAILS














 §.§ Synthetic examples (Figure <ref>)






With n=1000 and d=2000, we sample two matrices with orthonormal columns U ∈ℝ^n × n and V ∈ℝ^d × n, uniformly at random from their respective spaces. We construct two diagonal matrices Σ^poly, Σ^exp∈ℝ^n × n, such that their respective diagonal elements are Σ^poly_jj = √(n) j^-1 and Σ^exp_jj = √(n) e^-0.05 j. We set A^exp = U Σ^exp V^⊤ and A^poly = U Σ^poly V^⊤, and we sample a planted vector x_gd∈ℝ^d with iid entries 𝒩(0,1). 

In the case of binary logistic regression, for each A ∈{A^exp, A^poly}, we set y_i = 0.5 (sign(⟨ a_i, x_gd⟩)+1), for i=1,,n. 

For the ReLU model, we set y_i = (⟨ a_i, x_gd⟩)_+, where (z)_+ = max(0, z). Hence, each observation y_i is the result of a linear operation z_i = ⟨ x, a_i ⟩ and a non-linear operation y_i = (z_i)_+. Additionally, it can be shown that the global minimum of the optimization problem 



    min_x 1/2∑_i=1^n (a_i^⊤ x)_+^2 - 2 (a_i^⊤ x) y_i,


is equal to x_gd, which motivates using such a convex relaxation.












 §.§ Numerical illustration of the iterative and power methods (Figure <ref>)






We use the MNIST dataset with 50000 training images and 10000 testing images. We rescale the pixel values between [0,1]. Each image is mapped through random cosines φ(·) ∈ℝ^D which approximate the Gaussian kernel, i.e., ⟨φ(a), φ(a^') ⟩≈exp(-γa-a^2_2^2). We choose D=10000 and γ = 0.02. 

We perform binary logistic regression for even-vs-odd classification of the digits.

For the iterative method, we use the sketching matrix S = (A^⊤ A)^2 A^⊤S, where S is Gaussian iid. That is, we run the iterative method on top of the power method, with q=2.








 §.§ Adaptive Gaussian sketching on MNIST and CIFAR10 datasets (Table <ref> and Figure <ref>)








Experiments were run in Python on a workstation with 20 cores and 256 GB of RAM. The MNIST and CIFAR10 datasets were downloaded through the PyTorch.torchvision module and converted to NumPy arrays. We use the Sklearn.kernel_approximation.RBFSampler module to generate random cosines. We use our own implementation of each algorithm for a fair comparison.

For SGD, we use a batch size equal to 128. For SVRG, we use a batch size equal to 128 and update the gradient correction every 400 iterations. For Sketch-SGD, we use a batch size equal to 1024. For Sketch-SVRG, we use a batch size equal to 64 and update the gradient correction every 200 iterations. Each iteration of the sub-sampled Newton method (Sketch-Newton) computes a full-batch gradient, and, the Hessian with respect to a batch of size 1500. 

For SGD and SVRG, we considered step sizes η between 10^-2 and 10^2. We obtained best performance for η = 10^1. For the sub-sampled Newton method, we use a step size η=1, except for the first 5 iterations, for which we use η=0.2.

In Figure <ref>, we did not report results for SVRG for solving the primal (<ref>) on CIFAR, as the computation time for reaching a satisfying performance was significantly larger than for the other algorithms.

In Table <ref>, we did not investigate results for CIFAR with λ=10^-4, as the primal classifier had a test error significantly larger than smaller values of λ.
















§ PROOF OF MAIN RESULTS






Here, we establish our main technical results, that is, the deterministic regret bounds (<ref>) and (<ref>) stated in Theorem <ref> and its high-probability version stated in Corollary <ref>, along with its extension to the iterative Algorithm <ref> as given in Theorem <ref>, and its variant for kernel methods, given in Theorem <ref>. Our analysis is based on convex analysis and Fenchel duality arguments.












 §.§ Proof of Theorem <ref>










We introduce the Fenchel dual program of (<ref>), 


    min_z f^*(z) + 1/2λA^⊤ z_2^2.


For a sketching matrix S ∈ℝ^d × m, the Fenchel dual program of (<ref>) is 


    min_y f^*(y) + 1/2λP_S A^⊤ y_2^2.


Let α^* be any minimizer of the sketched program (<ref>). Then, according to Proposition <ref>, the unique solution of the dual sketched program (<ref>) is


    y^* = ∇ f(ASα^*)


and the subgradient set ∂ f^*(y^*) is non-empty. We fix g_y^*∈∂ f^*(y^*). 

According to Proposition <ref>, the dual program (<ref>) admits a unique solution z^*, which satisfies


    z^* = ∇ f(Ax^*),


and which subgradient set ∂ f^*(z^*) is non-empty. We fix g_z^*∈∂ f^*(z^*).

We denote the error between the two dual solutions by Δ = y^* - z^*. By optimality of y^* with respect to the sketched dual (<ref>) and by feasibility of z^*, first-order optimality conditions imply that


    ⟨1/λ AP_SA^⊤ y^* + g_y^*, Δ⟩ 0.


Similarly, by optimality of z^* with respect to the dual (<ref>) and by feasibility of y^*, we get by first-order optimality conditions that


    ⟨1/λ AA^⊤ z^* + g_z^*, Δ⟩ 0.


It follows that


    ⟨1/λ AP_S A^⊤Δ, Δ⟩   = ⟨1/λ AP_SA^⊤ y^*, Δ⟩ - ⟨1/λ AP_SA^⊤ z^*, Δ⟩
       = ⟨1/λ AP_SA^⊤ y^* + g_y^*, Δ⟩_ 0 + ⟨ g_z^* - g_y^*, Δ⟩ - ⟨1/λ AP_SA^⊤ z^* + g_z^*, Δ⟩
       ⟨ g_z^* - g_y^*, Δ⟩ - ⟨1/λ AP_SA^⊤ z^* + g_z^*, Δ⟩
        = ⟨ g_z^* - g_y^*, Δ⟩ + ⟨1/λ AP_S^⊥ A^⊤ z^*, Δ⟩ - ⟨1/λ A A^⊤ z^* + g_z^*, Δ⟩_ 0
       ⟨ g_z^* - g_y^*, Δ⟩ + ⟨1/λ AP_S^⊥ A^⊤ z^*, Δ⟩.


Strong μ-smoothness of f implies that the function f^* is 1/μ-strongly convex. Hence, it follows that 


    ⟨ g_z^* - g_y^*, Δ⟩ + 1/μΔ^2_2  0.


Therefore, combining (<ref>) with the previous set of inequalities (<ref>), we get


    ⟨1/λ AP_S A^⊤Δ, Δ⟩ + 1/μΔ^2_2 ⟨1/λ AP_S^⊥ A^⊤ z^*, Δ⟩,


and, multiplying both sides by λ,


    ⟨ AP_S A^⊤Δ, Δ⟩ + λ/μΔ^2_2 ⟨ AP_S^⊥ A^⊤ z^*, Δ⟩.


By definition of Z_f and since Δ∈dom f^* - z^*, it holds that


    Δ^⊤ A P_S^⊥ A^⊤Δ/Δ_2^2 Z_f^2,


which we can rewrite as


    ⟨ AP_S A^⊤Δ, Δ⟩⟨ A A^⊤Δ, Δ⟩ - Z_f^2 Δ_2^2.


Hence, combining (<ref>) and (<ref>), we obtain


    (λ/μ - Z_f^2) Δ_2^2 + A^⊤Δ_2^2 ⟨ A P_S^⊥ A^⊤ z^*, Δ⟩,


Under the assumption that λ 2 μ Z_f^2, it holds that λ / μ - Z_f^2 λ/(2μ). Thus,


    (λ/μ - Z_f^2) Δ_2^2 + A^⊤Δ_2^2    λ/2μΔ_2^2 + A^⊤Δ_2^2
       √(2 λ/μ)Δ_2 A^⊤Δ_2.


where we used the fact that for any a, b  0, a + b  2 √(ab), with a = λ/2μΔ_2^2 and b = A^⊤Δ_2^2. Combining the former inequality with inequality (<ref>), we obtain


    √(2 λ/μ)Δ_2 A^⊤Δ_2 ⟨ A P_S^⊥ A^⊤ z^*, Δ⟩.


The right-hand side of the latter inequality can be bounded as


    ⟨ A P_S^⊥ A^⊤ z^*, Δ⟩   = ⟨ A^⊤ z^*, P_S^⊥ A^⊤Δ⟩
       (i)A^⊤ z^*_2 P_S^⊥ A^⊤Δ_2 
       A^⊤ z^*_2 Δ_2 sup_Δ^'∈domf^* - z^*(P_S^⊥ A^⊤Δ^'_2/|Δ^'_2)
       (ii)=A^⊤ z^*_2 Δ_2 Z_f.


where (i) follows from Cauchy-Schwarz inequality, and (ii) holds by definition of Z_f. Thus, inequality (<ref>) becomes


    √(2 λ/μ)Δ_2 A^⊤Δ_2 A^⊤ z^*_2 Δ_2 Z_f.


From Propositions <ref> and <ref>, we have that A^⊤ z^* = - λ x^* and y^* = ∇ f(ASα^*). By definition of x, it follows that A^⊤Δ = -λ(x-x^*). Then, rearranging inequality (<ref>), 


    x-x^*_2 √(μ/2 λ) Z_f x^*_2,


which is exactly the desired regret bound (<ref>). The regret bound (<ref>) immediately follows from the fact that Z_f P_S^⊥ A^⊤_2.












 §.§ Proof of Corollary <ref>












The proof combines our deterministic regret bound (<ref>), along with the following result, which is a re-writing of Corollary 10.9, in <cit.>.



Let k  2 be a target rank and m  1 a sketching dimension such that k < m min(n,d). Let S be an n × m random matrix with iid Gaussian entries. Define the oversampling ratio r = (m-k)/k, and the sketching matrix S=A^⊤S. Then, provided r k  4, it holds with probability at least 1 - 6 e^-r k that


    P_S^⊥ A^⊤_2 c_0/√(2)√(r+1/r)( σ_k+1^2 + 1/rk∑_j = k+1^min(n,d)σ_j^2 )^1/2.


where σ_1 σ_2 are the singular values of A, and the universal constant c_0 satisfies c_0  36. In particular, if m=2k, then it holds with probability at least 1-6e^-k that


    P_S^⊥ A^⊤_2     c_0 ( σ_k+1^2 + 1/k∑_j = k+1^min(n,d)σ_j^2 )^1/2
       = c_0 R_k(A).




From Theorem <ref>, if λ 2μ Z_f^2, then


    x-x^*_2 √(μ/2λ)P_S^⊥ A^⊤_2 x^*_2.


Hence, combining the latter inequality with Lemma <ref>, provided 2  k 1/2min(n,d), m=2k and λ 2μ Z_f^2, it holds with probability at least 1-6e^-k that


    x-x^*_2  c_0 √(μ/2λ) R_k(A) x^*_2.


We want to establish the latter inequality, but under the condition λ 2 c^2_0 μ R_k^2(A). But, by Lemma <ref>, the condition λ 2 c^2_0 μ R_k^2(A) implies that λ 2 μ Z_f^2 with probability at least 1-6e^-k. By union bound, it follows that if λ 2 c^2_0 μ R_k^2(A), then 


    x-x^*_2  c_0 √(μ/2λ) R_k(A) x^*_2,


with probability at least 1-12e^-k.












 §.§ Proof of Theorem <ref>












First, we show that for any t  0, provided λ 2μ Z_f^2,


    x^(t+1) - x^*_2 √(μ/2λ)Z_f x^(t) - x^*_2.


It should be noted that for t=0, the latter inequality is exactly the regret bound (<ref>). The proof for t > 0 follows similar steps.

Fix t  0. Consider the optimization problem


    min_δ∈ℝ^d f(Aδ + Ax^(t)) + λ/2δ + x^(t)_2^2.


which is equivalent to the primal program (<ref>), up to a translation of the optimization variable. Thus, the unique optimal solution of (<ref>) – which exists by strong convexity of the objective – is given by δ^* = x^* - x^(t). By Fenchel duality (Corollary 31.2.1,  <cit.>), it holds that


    min_δ f(Aδ + Ax^(t)) + λ/2δ + x^(t)_2^2 = max_z - f^*(z) - 1/2λz^⊤ A A^⊤ z,


and the optimal dual solution z^* exists and is unique (by strong concavity of the dual objective). Further, by the Karush-Kuhn-Tucker conditions (Theorem 31.3, <cit.>), we have 


    δ^* = - x^(t) - 1/λ A^⊤ z^*
    
        z^* = ∇ f(Aδ^* + Ax^(t)).


Observe that, by using the change of variables α = (S^⊤ S)^-1/2α_†, the optimization problem (<ref>) can be rewritten as


    min_α∈ℝ^m f(A S α + A x^(t)) + λα^⊤ S^⊤x^(t) + λ/2S α_2^2
    ≡   min_α∈ℝ^m f(A S α + A x^(t)) + λ/2S α + x^(t)_2^2.
 

Let α_†^(t+1) be the unique solution of (<ref>). Then, setting α^(t+1) = (S^⊤ S)^-1/2α_†^(t+1), we have


    α^(t+1)∈_α∈ℝ^m f(AS α + A x^(t)) + λ/2S α  + x^(t)_2^2.


By Fenchel duality, we get


    min_α f(AS α + A x^(t)) + λ/2S α  + x^(t)_2^2 = max_y - f^*(y) + y^⊤ A x^(t) - 1/2λy^⊤ A P_S A^⊤ y + λ/2x^(t) P_S^⊥x^(t).


By strong concavity of the dual objective, there exists a unique maximizer y^*. Further, by the Karush-Kuhn-Tucker conditions (Theorem 31.3,  <cit.>), we have 


    ASα^(t+1) + Ax^(t)∈∂ f^*(y^*)


and, thus, y^* = ∇ f(ASα^(t+1) + Ax^(t)). 

We define Δ = y^* - z^*. Following similar steps as in the proof of Theorem <ref>, we obtain


    A^⊤Δ_2^2 + (λ/μ - Z_f^2) Δ_2^2    ⟨ A P_S^⊥ (λx^(t) + A^⊤ z^*), Δ⟩
       = -λ⟨ AP_S^⊥δ^*, Δ⟩.


Since λ 2 μ Z_f^2, it follows that λ / μ - Z_f^2 λ / (2μ). Using the fact that for any a, b  0, we have 2 √(ab) a + b, we obtain the inequality


    √(2λ/μ)A^⊤Δ_2 Δ_2     -λ⟨ AP_S^⊥δ^*, Δ⟩
       λδ^*_2 Z_f Δ_2.


Dividing both sides by λΔ_2 and using the identities δ^* = x^* - x^(t) and A^⊤Δ / λ = x^* - x^(t+1), we obtain the desired contraction inequality


    x^(t+1) - x^*_2 √(μ/2λ)Z_f x^(t) - x^*_2.


By induction, it immediately follows that for any number of iterations T  1, 


    x^(T) - x^*_2 (μ/2λ Z_f^2)^T/2x^*_2.


The high-probability version follows by immediate application of Lemma <ref> to the previous inequality.












 §.§ Proof of Theorem <ref>











Define


    x^* = _x f(K^1/2x) + λ/2x_2^2.


Set x = - 1/λK^1/2∇ f(K Sα^*), and


    Z_f = Z_f(K^1/2, K^1/2S).


Then, by application of Theorem <ref> with S=K^1/2S, it holds that x-x^*_2 √(μ/(2 λ)) Z_f x^*_2, provided that λ 2μ Z_f^2. We conclude by using the facts that x = K^1/2w and x^* = K^1/2w^*. 











§ PROOFS OF INTERMEDIATE RESULTS















 §.§ Proposition <ref> – Strong duality and Karush-Kuhn-Tucker conditions of the primal objective (<ref>)






Denote g(x) = λ/2x_2^2. The functions f and g are proper, closed, convex and their domains are respectively equal to ℝ^n and ℝ^d. It is then trivial that for any x ∈ℝ^d, we have x ∈dom(g) and Ax ∈dom(f). Hence, all conditions to apply strong Fenchel duality results hold (Theorem 31.2, <cit.>). Using the fact that g^*(z) = 1/2 λz_2^2, we get


    inf_x f(Ax) + λ/2x_2^2 = sup_z - f^*(z) - 1/2λA^⊤ z_2^2,


and the supremum is attained for some z^* ∈dom(f^*). The uniqueness of z^* follows from strong concavity of the dual objective, which comes from the 1/μ-strong convexity of f^*. 

Further, the primal objective is also strongly convex over ℝ^d, which implies the existence and uniqueness of a minimizer x^*. 

The Karush-Kuhn-Tucker conditions (Theorem 31.3, <cit.>) imply that Ax^* ∈∂ f^*(z^*). Since ∂ f^* = (∇ f)^-1 (Theorem 23.5, <cit.>), it follows that z^* = ∇ f(Ax^*). Finally, by first-order optimality conditions of x^*, we have that A^⊤∇ f(Ax^*) + λ x^* = 0, i.e., x^* = -λ^-1 A^⊤ z^*.







 §.§ Proposition <ref> – Strong duality and Karush-Kuhn-Tucker conditions of the sketched primal objective (<ref>)






Denote g(α) = λ/2Sα_2^2. The functions f and g are proper, closed, convex and their domains are respectively equal to ℝ^n and ℝ^m. It is then trivial that for any α∈ℝ^m, we have α∈dom(g) and ASα∈dom(f). Hence, all conditions to apply strong Fenchel duality results hold (Theorem 31.2, <cit.>). Using the fact that g^*(y) = 1/2 λ y^⊤ (S^⊤ S)^† y, we get


    inf_α f(ASα) + λ/2Sα_2^2 = sup_y - f^*(y) - 1/2λP_S A^⊤ y_2^2,


and the supremum is attained for some y^* ∈dom(f^*). The uniqueness of y^* follows from strong concavity of the dual objective, which comes from the 1/μ-strong convexity of f^*. 

We establish the existence of a minimizer α of α↦ f(ASα) + λ/2Sα_2^2. The latter function is strongly convex over the subspace (KerS)^⊥. Thus, there exists a unique minimizer α over (KerS)^⊥. Then, for any α∈ℝ^m, writing α = α_⊥ + α_∥ where α_⊥∈Ker(S)^⊥ and α_∥∈Ker(S), we have


    f(ASα) + λ/2Sα_2^2    = f(ASα^⊥)+λ/2Sα^⊥_2^2
        f(ASα) + λ/2Sα^2_2.


Thus, the point α is a minimizer.

Let α^* be any minimizer. The Karush-Kuhn-Tucker conditions (Theorem 31.3, <cit.>) imply that ASα^* ∈∂ f^*(y^*). Since ∂ f^* = (∇ f)^-1 (Theorem 23.5, <cit.>), it follows that y^* = ∇ f(ASα^*).












 §.§ Proof of Proposition <ref> – Numerical conditioning of the re-scaled sketched program







The condition number κ_† of the re-scaled sketched program is equal to


    κ_† = sup_αλ + σ_1( ℐ_S^⊤ A^⊤∇^2 f(A S α) ℐ_S )/inf_αλ + σ_m(ℐ_S^⊤ A^⊤∇^2 f(A S α) ℐ_S ),


where ℐ_S = S (S^⊤ S)^-1/2. 

In order to show that κ_†κ, it suffices to upper bound the numerator in the definition of κ_† by the numerator of κ and to lower bound the denominator of κ_† by the denominator of κ, i.e., it suffices to show that 


    sup_ασ_1( ℐ_S^⊤ A^⊤∇^2 f(A S α) A ℐ_S ) sup_x σ_1(A^⊤∇^2 f(Ax) A),
       inf_ασ_m( ℐ_S^⊤ A^⊤∇^2 f(A S α)A ℐ_S ) inf_x σ_d(A^⊤∇^2 f(Ax) A).


By the trivial inclusion { Sα|α∈ℝ^m }⊆ℝ^d, it holds that 


    sup_ασ_1( ℐ_S^⊤ A^⊤∇^2 f(A S α) A ℐ_S ) sup_x σ_1( ℐ_S^⊤ A^⊤∇^2 f(Ax) A ℐ_S ),
       inf_ασ_m(ℐ_S^⊤ A^⊤∇^2 f(A S α)A ℐ_S ) inf_x σ_m(ℐ_S^⊤ A^⊤∇^2 f(Ax)A ℐ_S ).


Therefore, to establish that κ_†κ, it is sufficient to show that for any x ∈ℝ^d,


    σ_1( ℐ_S^⊤ A^⊤∇^2 f(A x) Aℐ_S ) σ_1(A^⊤∇^2 f(Ax) A),
       σ_m( ℐ_S^⊤ A^⊤∇^2 f(Ax) A ℐ_S ) σ_d(A^⊤∇^2 f(Ax) A).


The first inequality follows from the fact that ℐ_S_2  1. Hence, we have


    σ_1( ℐ_S^⊤ A^⊤∇^2 f(A x) Aℐ_S )    = sup_w ≠ 0w^⊤ℐ_S^⊤ A^⊤∇^2 f(A x) A ℐ_S w/w_2
       = sup_w ≠ 0(ℐ_S w)^⊤ A^⊤∇^2 f(A x) A (ℐ_S w)/ℐ_S w_2ℐ_S w_2/w_2_ 1
       sup_z ≠ 0z^⊤ A^⊤∇^2 f(A x) A z/z_2
       = σ_1(A^⊤∇^2 f(Ax) A).


For the second inequality, we distinguish two cases. 

If the sketching matrix S ∈ℝ^d × m is full-column rank, then, the matrix ℐ_S is actually an isometry, i.e., for any w ∈ℝ^m, we have S (S^⊤ S)^-1/2 w_2 = w_2, which implies that 


    σ_m( ℐ_S^⊤ A^⊤∇^2 f(A x) Aℐ_S )    = inf_w ≠ 0w^⊤ℐ_S^⊤ A^⊤∇^2 f(A x) A ℐ_S w/w_2
       = inf_w ≠ 0(ℐ_S w)^⊤ A^⊤∇^2 f(A x) A (ℐ_S w)/ℐ_S w_2ℐ_S w_2/w_2_= 1
       inf_z ≠ 0z^⊤ A^⊤∇^2 f(A x) A z/z_2
       = σ_d(A^⊤∇^2 f(Ax) A).


Suppose now that the sketching matrix S ∈ℝ^d × m is not full column-rank. By assumption, S = A^⊤S where S∈ℝ^n × m is Gaussian iid, hence, full-column rank almost surely. It implies that there exists a vector v ≠ 0 such that A v =0. Indeed, let v ≠ 0 be a vector such that S^⊤ v = 0, which exists since m < d. The equation S^⊤ v = 0 can be rewritten as S^⊤ A v = 0. Since S^⊤ is full row-rank, we get that Av=0, i.e., Ker A ≠{0}.

From Ker A ≠{0}, we get σ_d(A^⊤∇^2 f(Ax) A) = 0 and σ_m( ℐ_S^⊤ A^⊤∇^2 f(Ax) A ℐ_S )=0, which concludes the proof.











§ PROOF OF BOUNDS IN TABLE <REF>














 §.§ Adaptive Gaussian sketching











From Corollary <ref>, we have that for a target rank k and a sketching dimension m_A=2k, with probability at least 1-12e^-k,


    x-x^*_2/x^*_2≲λ^-1/2(ν_k + 1/k∑_j=k+1^ρν_j )^1/2.


For a matrix A with rank ρ≪min(n,d), with k ρ + 1, the right hand side of the latter equation is equal to 0. In order to achieve this with probability at least 1-η, it is sufficient to oversample by an amount log(12/η), that is, m_A=ρ+1+log(12/η) is sufficient to achieve a (ε=0, η)-guarantee.

For a κ-exponential decay with κ > 0, we have ν_j ∼ e^-κ j, and


    (ν_k + 1/k∑_j=k+1^ρν_j)^1/2∼ e^-κ (k+1)/2,


and it is sufficient for the sketching dimension m_A to satisfy


    m_A ≳κ^-1log( 1/λε) + log(12/η).


For a β-polynomial decay with β > 1/2, we have ν_j ∼ j^-2β and 


    (ν_k + 1/k∑_j=k+1^ρν_j)^1/2∼ k^-β,


and it is sufficient to have


    m_A ≳λ^-1/2βε^-1/β + log(12/η).













 §.§ Oblivious Gaussian sketching
 










For a ρ-rank matrix A, it has been shown in <cit.> that, provided the sketching dimension m_O satisfies


    m_O ≳(ρ+1)log(2ρ/η)/ε^2,


then,


    x-x^*_2/x^*_2≲ε,


with probability at least 1-η, for any ε∈ (0,1/2).

We now justify the bounds for the κ-exponential and β-polynomial decays. Let ρ be the effective rank of the matrix AA^⊤, defined as 


    ρ = ∑_i=1^ρν_i/λ + ν_i.


In <cit.>, the authors have shown that, provided 


    m_O ≳ρ/ε^2 (λ+1)log(2d/η),


then the relative error satisfies


    x-x^*_2/x^*_2≲ε(1+√(λ/ν_k)),


with probability at least 1-η, under the additional condition that the minimizer x^* lies in the subspace spanned by the top k-left singular vectors of A. For simplicity of comparison, we neglect the latter (restrictive) requirement on x^*, and the term √(λ/ν_k) in the latter upper bound, which yields a smaller lower bound on a sufficient sketching size m_O to achieve a (ε, η)-guarantee. Based on those simplifications, oblivious Gaussian sketching yields a relative error such that


    x-x^*_2/x^*_2ε,


with probability at least 1-η, provided that m_O ρε^-2log(2d/η).

For a κ-exponential decay, it holds that


    ρ = ∑_i=1^ρe^-κ i/e^-κ i+λ∫_1^ρe^-κ t/e^-κ t+λ dt
        = ∫_e^-κρ^e^-κ1/κ1/u+λ du
        = 1/κlog(e^-κ+λ/e^-κρ + λ).


Since λ∈ (ν_ρ, 1) = (e^-κρ, 1), it follows that 


    ρ≳κ^-1log1/λ.


Hence, their theoretical predictions state that the sketching dimension m_O must be greater than


    m_O ≳κ^-1ε^-2log(1/λ) log(2d/η)


in order to achieve a (ε, η)-guarantee.

For a β-polynomial decay, it holds that


    ρ = ∑_i=0^ρ1/1+λ i^2β≥ -1 + ∫_0^ρ1/1+λ t^2β dt
            = -1 + λ^-1/2β/2 β∫_0^λρ^2βu^1/2β-1/1+u du
       ≥ -1 + λ^-1/2β/2 β∫_0^1u^1/2β-1/1+u du,


where the last inequality is justified by the fact that the integrand is non-negative, and the fact that λρ^-2β. Since the integral is finite and independent of λ, it follows that


    ρ≳λ^-1/2β,


and the sketching dimension m_O must satisfy 


    m_O ≳λ^-1/2βε^-2log(2d/η)


in order to achieve a (ε, η)-guarantee, according to their theoretical predictions.












 §.§ Leverage score column sampling.
 










Let A = U Σ V^⊤ be a singular value decomposition of the matrix A, where Σ = diag(σ_1, σ_2, , σ_ρ), and σ_1 σ_2 σ_ρ. For a given target rank k, let u_1, , u_k be the first k columns of the matrix U, and denote U_1 = [u_1, , u_k] ∈ℝ^n × k. For j=1,,n, define p_j = k^-1U_1,j_2^2, where U_1, j is the j-th row of the matrix U_1. By orthonormality of the family (u_1, , u_k), it holds that ∑_j=1^n p_j = 1, and p_j  0. The family {p_j}_j=1^n is called the leverage score probability distribution of the Gram matrix A A^⊤. 

Leverage based column sampling consists in, first, computing the exact or approximated leverage score distribution of the matrix AA^⊤, and, second, sampling m columns of AA^⊤ from the latter probability distribution, with replacement. Precisely, the sketching matrix S ∈ℝ^d × m is given as


    S = A^⊤ R D,


where R ∈ℝ^n × m is a column selecting matrix drawn according to the leverage scores, and D ∈ℝ^m × m is a diagonal rescaling matrix, with D_jj = (m p_i)^-1/2, if R_ij = 1.

In order to compare the theoretical guarantees of adaptive Gaussian sketching and leverage-based column sampling, we assume that the leverage scores are computed exactly. Note that if this is not the case, then the sketching size increases as the quality of approximation of the leverage score distribution decreases. As our primary goal is to lower bound the ratio m_S / m_A, our qualitative comparison is not affected (at least not in the favor of adaptive Gaussian sketching) by this assumption.

The authors of <cit.> showed that given δ∈ (0,1], provided m_S satisfies


    m_S ≳δ^-2 k log(k/η),


then, with probability at least 1-η,


    P^⊥_S A^⊤_2 ν_k^1/2 + δ^2 ∑_j=k+1^ρν_j^1/2.


Using the combination of the latter concentration bound with our deterministic regret bound (<ref>) on the relative error, it follows that, under the latter condition on m_S, with probability at least 1-η,


    x-x^*_2/x^*_2λ^-1/2(ν_k^1/2 + δ^2 ∑_j=k+1^ρν_j^1/2).


For a matrix A with rank ρ≪min(n,d), if the sketching size m_S is greater than (ρ + 1) log(ρ + 1/η), then the relative error satisfies an (ε=0, η)-guarantee.

For a κ-exponential decay, we have


    ν_k^1/2 + δ^2 ∑_j=k+1^ρν_j^1/2∼(1 + 2δ^2/κ) e^-κ (k+1)/2.


Taking δ = 1/2, it follows that the sketching size m_S must be greater than κ^-1log(1/λε) log(1/η) to satisfy


    x-x^*_2/x^*_2≲ε


with probability at least 1-η.

For a β-polynomial decay (with β > 1), we have 


    ν_k+1 + δ^2 ∑_j=k+1^ρν_j ∼ k^-β + δ^2 β^-1 k^1-β,


and, provided m_S ≳δ^2 k log(k/η),


    x-x^*_2/x^*_2λ^-1/2(k^-β + δ^2 β^-1 k^1-β).


To achieve a precision ε, it is sufficient to have


    ελ^-1/2(k^-β + δ^2 β^-1 k^1-β)


Suppose first that we choose δ≲ k^-1/2. Then, the latter sufficient condition becomes ε≳λ^-1/2 k^-β. Hence, we need k to be at least λ^-1/2 βε^-1/β, which implies m_S ≳δ^-2 k log(k/η). Since δ≲ k^-1/2, we get that m_S must be greater than k^2 log(1/η), which further implies


    m_S ≳λ^-1/βε^-2/βlog(1/η).


Now, suppose that we choose δ≳ k^-1/2. Write δ^2 = k^-1+γ, where γ > 0. Since δ < 1, we must have γ∈ (0,1). Further, we need βελ^-1/2δ^2 k^1-β = λ^-1/2 k^γ - β. By assumption, β > 1, hence, γ - β < 0. Hence, the smallest value of k that satisfies the latter inequality is given by 


    k = (ε^-1/βλ^-1/2β)^1/1-γ/β.


On the other hand, the smallest sketching size m_S to achieve an (ε, η) satisfies 


    m_S ≳ k^2-γlog(1/η).


Plugging-in the value of k, we must have 


    m_S ≳(ε^-1/βλ^-1/2β)^2-γ/1-γ/βlog(1/η).


Optimizing over γ∈ (0,1), we finally obtain that the best sufficient sketching size must satisfy


    m_S ≳(ε^-1/βλ^-1/2β)^min(2,β/β-1)log(1/η).


























§ EXTENSION TO THE NON-SMOOTH CASE











Here, we present some results to the case where the function f : ℝ^n →ℝ is proper, convex, but not necessarily smooth. We make the assumption that the function is L-Lipschitz, that is, for any x, y ∈ℝ^n, 


    f(x) - f(y)_2  L x-y_2.


In particular, this implies that the domain of the function f^* is bounded, i.e., for any z ∈dom f^*, it holds that z_2  L.

Let x^* be the solution of the primal program (<ref>), which exists and is unique by strong convexity of the primal objective. 

For a sketching matrix S ∈ℝ^d × m, the sketched primal program (<ref>) admits a solution α^*. Indeed, using arguments similar to the proof of Proposition <ref>, the sketched program is strongly convex over Ker(S)^⊥, and admits a unique solution α^* over that subspace. Further, for any α∈ℝ^m, we can decompose α = α_⊥ + α_∥, with α_⊥∈Ker(S)^⊥ and α_∥∈Ker(S). Then, 


    f(AS α) + λ/2Sα_2^2    = f(AS α_⊥) + λ/2Sα_⊥_2^2
        f(AS α^*) + λ/2Sα^*_2^2.


As for the smooth case, using convex analysis arguments, we obtain that the dual program (<ref>) has a solution z^* which satisfies z^* = ∇ f(Ax^*). Similarly, the sketched dual program (<ref>) has a solution y^* which satisfies y^*= ∇ f(ASα^*). Further, by first-order optimality conditions of x^*, we have x^* = -λ^-1 A^⊤∇ f(Ax^*), i.e., x^* = -λ^-1 A^⊤ z^*.

As for the smooth case, we introduce the candidate approximate solution x, defined as 


    x = -λ^-1 A^⊤∇ f(ASα^*),


where α^* is any minimizer of (<ref>).



For any λ > 0, it holds that


    x - x^*_2 6 L/λ√(σ_1 Z_f).


where σ_1 is the top singular value of A.




Let α^* be any minimizer of the sketched primal program. Following similar lines as in the proof of Theorem <ref> (see Appendix <ref>), it holds that


    A^⊤Δ_2^2 ≤ Z_f^2 Δ_2^2 + z^*^⊤ A P_S^⊥ A^⊤Δ,


where Δ = y^* - z^*. After applying Cauchy-Schwarz and using the definition of Z_f, inequality (<ref>) becomes


    A^⊤Δ_2^2 ≤ Z_f^2 Δ_2^2 + Z_f z^*_2 A^⊤Δ_2.


Using the fact that √(w + w^')≤√(w) + √(w^') with w = Z_f^2 Δ_2^2 and w^' = Z_f z^*_2 A^⊤Δ_2, along with the inequality A^⊤Δ_2 ≤σ_1 Δ_2, we obtain


    A^⊤Δ_2 ≤√(Z_f)(√(Z_f)Δ_2 + √(z^*_2 Δ_2 σ_1)).


Using the inequality 2ww^'≤ w^2 + w^'^2 and the fact that Z_f ≤σ_1, it follows that 


    A^⊤Δ_2  2 √(Z_f σ_1)( Δ_2 + z^*_2)


Dividing by λ and using the fact that x - x^* = -λ^-1 A^⊤Δ, we get


    x - x^*_2 ≤2/λ√(Z_f σ_1)( y^*-z^*_2 + z^*_2 ).


Using the fact that y^*_2, z^*_2  L, we obtain the desired inequality (<ref>).



As for the smooth-case, high-probability bounds follow from the previous deterministic bound on the relative error.






