---
abstract: 'We develop a novel preconditioning method for ridge regression, based on recent linear sketching methods. By equipping Stochastic Variance Reduced Gradient (SVRG) with this preconditioning process, we obtain a significant speed-up relative to fast stochastic methods such as SVRG, SDCA and SAG.'
author:
- 'Alon Gonen[^1]'
- 'Francesco Orabona[^2]'
- 'Shai Shalev-Shwartz[^3]'
bibliography:
- 'bib.bib'
title: |
    Solving Ridge Regression using\
    Sketched Preconditioned SVRG
---

Introduction {#sec:intro}
============

Consider the *ridge regression* problem: $$\label{eq:ridge}
\min_{w \in \reals^d}\  \left\{L(w) = \frac{1}{n}\sum_{i=1}^n \frac{1}{2}(w^\top x_i-y_i)^2 + \frac{\lambda}{2} \|w\|^2\right\},$$ where $\lambda>0$ is a regularization parameter, $x_i\in{\mathbb{R}}^d$ and $y_i\in {\mathbb{R}}$ for $i=1,\cdots,n$ the training data. We focus on the large scale regime, where both $n$ and $d$ are large. In this setting, stochastic iterative methods such as SDCA [@shalev2013stochastic], SVRG [@johnson2013accelerating], and SAG [@roux2012stochastic] have become a standard choice for minimizing the objective $L$. Specifically, the overall complexity of a recent improved variant of SVRG due to [@xiao2014proximal] depends on the average condition number, which is defined as follows. Denote the empirical correlation matrix and its eigenvalue decomposition by $$\label{eq:covSpec}
C := \frac{1}{n} \sum_{i=1} ^n x_i x_i^\top = \sum_{i=1} ^d \lambda_i u_i u_i^\top~.$$ The average condition number of $C+\lambda I$ is defined as the ratio between the trace of the Hessian of $L$ and its minimal eigenvalue: $$\label{eq:initialCondition}
{\hat{\kappa}}:={\hat{\kappa}}(C+\lambda I) = \frac{\tr(C+\lambda I)}{\lambda_d(C+\lambda I)} = \sum_{i=1}^d \frac{\lambda_i+\lambda}{\lambda_d+\lambda}  ~.$$ The mentioned variant of SVRG finds an $\epsilon$-approximate minimizer of $L$ in time $\tilde{O}(({\hat{\kappa}}+n)d \log(1/\epsilon))$. Namely, the output of the algorithm, denoted $\hat{w}$, satisfies $\bE[L(\hat{w})]-L(w^\star) \le \epsilon$, where the expectation is over the randomness of the algorithm. For an accelerated version of the algorithm, we can replace ${\hat{\kappa}}$ by $\sqrt{n {\hat{\kappa}}}$ [@shalev2014accelerated; @lin2015universal].

The regularization parameter, $\lambda$, increases the smallest eigenvalue of $C+\lambda I$ to be at least $\lambda$, thus improves the condition number and makes the optimization problem easier. However, to control the under/over fitting tradeoff, $\lambda$ has to decrease as $n$ increases [@shalev2014understanding]. Moreover, in many machine learning applications $\lambda_d$ approaches zero and it is usually smaller than the value of $\lambda$. Overall, this yields a large condition number in most of the interesting cases.

A well-known approach for reducing the average condition number is *preconditioning*. Concretely, for a (symmetric) positive definite (pd) matrix $P \in \reals^{d \times d}$, we define the preconditioned optimization problem as $$\label{eq:condRidge}
\min_{{\tilde{w}}\in \reals^d} {\tilde{L}}({\tilde{w}}):= L(P^{-1/2} {\tilde{w}}) ~.$$ Note that ${\tilde{w}}$ is an $\epsilon$-approximate minimizer of ${\tilde{L}}$ if and only if $w = P^{-1/2} {\tilde{w}}$ forms an $\epsilon$-approximate minimizer of $L$. Hence, we can minimize rather than . As we shall see, the structure of the objective allows us to apply the preconditioning directly to the data (as a preprocessing step) and consequently rewrite the preconditioned objective as a ridge regression problem with respect to the preconditioned data (see ). For a suitable choice of a matrix $P$, the average condition number is significantly reduced. Precisely, as will be apparent from the analysis, the pd matrix that minimizes the average condition number is $P=C+\lambda I$, and the corresponding average condition number is $d$. However, we note that such preconditioning process would require both the computation of $P^{-1/2}$ and the computation of $P^{-1/2} x_i$ for each $i \in [n]$. By first order conditions, computing $(C+\lambda I)^{-1/2}$ is equivalent to solving the original problem in , rendering this “optimal” preconditioner useless. Yet, the optimal preconditioner might not needed in many cases. In fact, a common empirical observation (see ) is that (high-dimensional) machine learning problems tend to have few dominant features, while the other coordinates are strongly correlated with the stronger features. As a result, the spectrum of the correlation matrix decays very fast. Hence, it is natural to expect to gain a lot from devising preconditioning methods that focus on the stronger directions of the data.

Our contributions are as follows. We develop a relatively cheap preconditioning method that, coupled with SVRG, assures to speed-up the convergence in practical applications while having a computational cost comparable to SVRG alone. In order to approximately extract the stronger directions while incurring a low computational cost, we rely on a variant of the Block Lanczos method due to [@musco2015stronger] in order to compute an approximated truncated SVD (Singular Value Decomposition) of the correlation matrix $C$. Finally, by equipping SVRG with this preconditioner, we obtain our main result.

Main Result {#sec:mainResult}
===========

\[thm:main\] Let $k \in [d]$ be a given parameter and assume that the regularization parameter, $\lambda$, is larger than $\lambda_d$. Our preconditioning process runs in time $O(ndk \log(n))$. By equipping the SVRG of [@xiao2014proximal] with this preconditioner, we find an $\epsilon$-approximate minimizer for (with probability at least $9/10$) in additional runtime of $O ((\tilde{\kappa}  + n + d )d \log(1/\epsilon) )$, where $\tilde{\kappa}=\frac{k \lambda_k+ \sum_{i>k} \lambda_i}{\lambda}$ or $\tilde{\kappa} = \left(\frac{n(k \lambda_k+ \sum_{i>k} \lambda_i)}{\lambda}\right)^{1/2}$ if we use accelerated SVRG.

When the runtimes of both the (accelerated) SVRG and our preconditioned (accelerated) SVRG are controlled by the average condition number (and both runtimes dominate $ndk$), then ignoring logarithmic dependencies, we obtain a speed-up of order $$\label{eq:ratioRidge}
\textrm{ratio} = \frac{\sum_{i=1}^d \lambda_i}{k\,\lambda_k+\sum_{i>k}\lambda_i}
= \frac{\sum_{i=1}^k \lambda_i +\sum_{i>k}\lambda_i }{k\,\lambda_k ~~~~+\sum_{i>k}\lambda_i}~.$$ (or $\sqrt{\sum_{i=1}^d \lambda_i/(\lambda_k k+\sum_{i>k}\lambda_i)}$ if acceleration is used) over SVRG. If the spectrum decays fast then $k \, \lambda_k \ll \sum_{i=1}^k \lambda_i$ and $\sum_{i>k}\lambda_i \ll k \, \lambda_k$. In this case, the ratio will be large. Indeed, as we show in the experimental section, this ratio is often *huge* for relatively small $k$.

Main challenges and perspective
-------------------------------

While the idea of developing a preconditioner that focuses on the stronger directions of the data matrix sounds plausible, there are several difficulties that have to be solved.

-   First, since a preconditioner must correspond to an invertible transformation, it is not clear how to form a preconditioner based on a low rank approximation and, in particular, how should we treat the non-leading components.

-   One of the main technical challenges in our work is to translate the approximation guarantees of the Lanczos method into a guarantee on the resulted average condition number. The standard measures of success for low-rank approximation are based on either Frobenius norm or spectral norm errors. As will be apparent from the analysis (see ), such bounds do not suffice for our needs. Our analysis relies on stronger per vector error guarantees due to [@musco2015stronger].

It should be emphasized that while we use a variant of SVRG due to [@xiao2014proximal], we could equally use a variant of SDCA [@shalev2016sdca] or develop such a variant for SAG or SAGA. Furthermore, while we focus on the quadratic case, we believe that our ideas can be lifted to more general setting. For example, when applied to self-concordant functions, each step of Newton’s method requires the minimization of a quadratic objective. Therefore, it is natural to ask if we can benefit from applying our method for approximating the Newton step.

Bias-complexity tradeoff {#sec:reduceReg}
------------------------

As we mentioned above, $\lambda$ controls a tradeoff between underfitting and overfitting. In this view, we can interpret our result as follows. Assuming for simplicity that $n \ge d$ and ignoring logarithmic dependencies, we note that if $$\label{eq:decay}
\lambda = \frac{k \lambda_k+\sum_{i>k}\lambda_i}{nk}~,$$ then the runtime of our preconditioned SVRG is $\tilde{O}(ndk)$. For comparison, the runtime of (unconditioned) SVRG is $\tilde{O}(ndk)$ if $$\label{eq:decaySVRG}
\lambda = \frac{\sum_{i=1}^d \lambda_i}{nk}~.$$ The ratio between the RHS of and is the ratio given in . Hence, for a given “runtime budget” of order $\tilde{O}(ndk)$, we can set the regularization parameter of the preconditioned SVRG to be smaller by this ratio. Similar interpretation holds for the accelerated versions.

Related Work {#sec:related}
============

#### Existing algorithms and their complexities: 

Since minimizing is equivalent to solving the system $(C+\lambda I)w=\frac{1}{n} \sum_{i=1}^n y_i x_i$, standard numerical linear algebra solvers such as Gaussian elimination can be used to solve the problem in time $O(nd^2)$.

Iterative deterministic methods, such as Gradient Descent (GD), finds an $\epsilon$-approximate minimizer in time $nd \kappa \log(1/\epsilon)$, where $\kappa=\frac{\lambda_1(C+\lambda I)}{\lambda_d(C+\lambda I)}$ is the condition number of $C+\lambda I$ (see Theorem 2.1.15 in [@nesterov2004introductory]). The Kaczmarz algorithm [@kaczmarz1937angenaherte] has an identical complexity. Both the Conjugate Gradient (CG) method [@hestenes1952methods] and the Accelerated Gradient Descent (AGD) algorithm of [@nesterov1983method] enjoy a better runtime of $nd \sqrt{\kappa}\log(1/\epsilon)$. In fact, CG has a more delicate analysis (see Corollary 16.7 in [@vishnoi2012laplacian]): If all but $c \in [d]$ eigenvalues of $C+\lambda I$ are contained in a range $[a,b]$, then the runtime of CG is at most $nd(c+\sqrt{b/a} \log(1/\epsilon))$. In particular, CG’s runtime is at most $O(nd^2)$. Furthermore, following the interpretation of our main result in , we note that for a “runtime budget” of $\tilde{O}(ndk)$, we can set the regularization parameter of CG to be of order $\lambda_k/k^2$ (which is usually much greater than the RHS of ).

#### Linear Sketching: 

Several recently developed methods in numerical linear algebra are based on the so-called *sketch-and-solve* approach, which essentially suggests that given a matrix $A$, we first replace it with a smaller random matrix $AS$, and then perform the computation on $AS$ [@woodruff2014sketching; @clarkson2013low; @sarlos2006improved]. For example, it is known that if the entries of $S$ are i.i.d. standard normal variables and $S$ has $p=\Omega(k/\epsilon)$ columns, then with high probability, the column space of $AS$ contains a $(1+\epsilon)$ rank-$k$ approximation to $A$ with respect to the Frobenius norm. This immediately yields a fast PCA algorithm (see Section 4.1 in [@woodruff2014sketching]).

While the above sketch-and-solve approach sounds promising for this purpose, our analysis reveals that controlling the Frobenius norm error does not suffice for our needs. We need spectral norm bounds, which are known to be more challenging [@witten2013randomized]. Furthermore, as mentioned above, the success of our conditioning method heavily depends on the stronger per vector error guarantees obtained by [@musco2015stronger] which are not obtained by simpler linear sketching methods.

#### Sketched preconditioning: 

Recently, subspace embedding methods were used to develop cheap preconditioners for linear regression with respect to the squared loss [@woodruff2014sketching]. Precisely, [@clarkson2013low] considered the case $\lambda=0$ (i.e, standard least-squares) and developed a preconditioning method that reduces the average condition number to a constant. Thereafter, they suggest applying a basic solver such as CG. The overall running time is dominated by the preconditioning process which runs in time $\tilde{O}(d^3+nd)$. Hence, a significant improvement over standard solvers is obtained if $n \gg d$.

The main shortcoming of this method is that it does not scale well to large dimensions. Indeed, when $d$ is very large, the overhead resulted from the preconditioning process can not be afforded.

#### Efficient preconditioning based on random sampling: 

While we focus on reducing the dependence on the dimensionality of the data, other work investigated the gain from using only a random subset of the data points to form the conditioner [@yang2014data]. The theoretical gain of this approach has been established under coherence assumptions [@yang2014data].

Preliminaries {#sec:pre}
=============

Additional notation and definitions {#sec:notation}
-----------------------------------

Any matrix $B \in \reals^{d \times n}$ of rank $r$ can be written in (thin) SVD form as $B = U \Sigma V^\top = \sum_{i=1}^r \sigma_i(B) u_i v_i^\top$. The singular values are ordered in descending order. The spectral norm of $B$ is defined by $\|B\|=\sigma_1(B)$. The spectral norm is submultiplicative, i.e., $\|AB\| \le \|A\|\|B\|$ for all $A$ and $B$. Furthermore, the spectral norm is unitary invariant, i.e., for all $A$ and $U$ such that the columns of $U$ are orthonormal, $\|UA\| = \|A\|$. For any $k \in [r]$, it is well known that the truncated SVD of $B$, $B_k := U_k \Sigma_k V_k =  \sum_{i=1} ^{k} \sigma_i(B) u_i v_i^\top$, is the best rank-$k$ approximation of $B$ w.r.t. the spectral norm [@trefethen1997numerical]. A twice continuously differentiable function $f: \reals^d \rightarrow \reals$ is said to be $\beta$-smooth if $\|\nabla^2 f(w)\| \leq \beta$ for all $w$, where $\nabla^2 f(w)$ is the Hessian of $f$ at $w$. $f$ is said to be $\alpha$-strongly convex if $\lambda_d(\nabla^2 f(w)) \ge \alpha$ for all $w$. If $g$ is convex and $f$ is $\alpha$-strongly convex, then $f+g$ is $\alpha$-strongly convex.

Stochastic Variance Reduced Gradient (SVRG) {#sec:svrg}
-------------------------------------------

We consider a variant of the Stochastic Variance Reduced Gradient (SVRG) algorithm of [@johnson2013accelerating] due to [@xiao2014proximal]. The algorithm is an epoch-based iterative method for minimizing an average, $F(w) = \frac{1}{N} \sum_{i=1}^N f_i(w)$, of smooth functions. It is assumed that each $f_i:\reals^d \rightarrow \reals$ is convex and $\beta_i$-smooth. The entire function $F$ is assumed to be $\alpha$-strongly convex. The algorithm is detailed in . Its convergence rate depends on the averaged smoothness of the individual functions and the average condition number of $F$, defined as $$\label{eq:conditionSVRG}
\hat{\beta} = \frac{1}{N} \sum_{i=1}^N \beta_i ~~~;~~~ {\hat{\kappa}}_F = \frac{\hat{\beta}}{\alpha} ~.$$

**[@xiao2014proximal]** \[thm:svrg\] Fix $\epsilon>0$. Running SVRG () with any $w_0$, $S \ge \log((F(w_0) - \min_{w \in \reals^d} F(w))/\epsilon)$, $m = \lceil {\hat{\kappa}}_F \rceil$, and $\eta = 0.1/\hat{\beta}$ yields an $\epsilon$-approximate minimizer of $F$. Furthermore, assuming that each single gradient $\nabla f_i(w)$ can be computed in time $O(d)$, the overall runtime is $O(({\hat{\kappa}}_F+N)d \log(\epsilon_0/\epsilon))$.

**Input:** Functions $f_1,\ldots,f_n, \beta_1,\ldots,\beta_n$ **Parameters:** $\bar{w}_0 \in \reals^d$, $m$, $\eta$, $S \in \mathbb{N}$ $\bar{w} = \bar{w}_{s-1}$ $\bar{v} = \nabla F(\bar{w})$ $w_0 = \bar{w}$ Pick $i_t \in [N]$ with probability $q_{i_t}=\beta_{i_t}/\sum \beta_j$ $v_t = (\nabla f_{i_t}(w_{t-1}) -\nabla f_{i_t}(\bar{w}))/q_{i_t}+\bar{v}$ $w_t = w_{t-1} - \eta v_t$ $\bar{w}_s = \frac{1}{m} \sum_{t=1}^m w_t$ **Output:** the vector $\bar{w}_S$

In the original definition of SVRG [@johnson2013accelerating], the indices $i_t$ are chosen uniformly at random from $[n]$, rather than proportional to $\beta_i$. As a result, the convergence rate depends on the maximal smoothness, $\max \{\beta_i\}$, rather than the average, $\hat{\beta}$. It will be apparent from our analysis (see ) that in our case, $\max \{\beta_i\}$ is proportional to the maximum norm of any preconditioned $x_i$. Since we rely on the improved variant of [@xiao2014proximal], our bound depends on the average of the $\beta_i$’s, which scale with the average norm of the preconditioned $x_i$’s. To simplify the presentation, in the sequel we refer to as SVRG.

Randomized Block Lanczos
------------------------

A randomized variant of the Block Lanczos method due to [@musco2015stronger] is detailed[^4] in . Note that the matrix ${\tilde{U}}_k {\tilde{\Sigma}}_k {\tilde{V}}_k^\top$ forms an SVD of the matrix ${\tilde{A}}_k := Q(Q^\top A)_k = {\tilde{U}}_k {\tilde{U}}_k^\top A$.

**Input:** $A \in \reals^{d \times n}, k \le d, \epsilon' \in (0,1)$ $q= \Theta \left(\frac{\log(n)}{\sqrt{\epsilon}} \right)$, $p=qk$, $\Pi \sim \cN(0,1)^{n \times k}$ Compute $K = [A \Pi,(AA^\top) A \Pi,\ldots, (AA^\top)^{q-1} A \Pi]$Orthonormalize $K$’s columns to obtain $Q \in \reals^{d \times qk}$ Compute the truncated SVD $(Q^\top A)_k = {\tilde{W}}_k {\tilde{\Sigma}}_k {\tilde{V}}_k^\top$ Compute ${\tilde{U}}_k = Q {\tilde{W}}_k$ **Output:** the matrices ${\tilde{U}}_k,{\tilde{\Sigma}}_k,{\tilde{V}}_k$

\[thm:musco\] **[@musco2015stronger]** Consider the run of and denote ${\tilde{A}}_k = {\tilde{U}}_k {\tilde{\Sigma}}_k {\tilde{V}}_k = \sum_{i=1}^k {\tilde{\sigma}}_i {\tilde{u}}_i {\tilde{v}}_i^\top$. Denote the SVD of $A$ by $A = \sum_{i=1}^d \sigma_i v_i u_i^\top$. The following bounds hold with probability at least $9/10$: $$\|A - {\tilde{A}}_k\| \le (1+\epsilon') \|A-A_k\| \le (1+\epsilon') \sigma_k$$ $$\begin{aligned}
\forall i \in [k],~~|z_i^\top AA^\top z_i - u_i^\top AA^\top u_i| &=  |{\tilde{\sigma}}_i^2-\sigma_i^2 | \notag \\
&\le \epsilon' \sigma_{k+1}^2 ~. \label{eq:perVec}\end{aligned}$$ The runtime of the algorithm is $O\left(\frac{nd k \log(n)}{ \sqrt{\epsilon'}}+ \frac{k^2(n+d)}{\epsilon'} \right)$.

Sketched Conditioned SVRG {#sec:scsvrg}
=========================

In this section we develop our sketched conditioning method. By analyzing the properties of this conditioner and combining it with SVRG, we will conclude .

Recall that we aim at devising cheaper preconditioners that lead to a significant reduction of the condition number. Specifically, given a parameter $k \in [d]$, we will consider only preconditioners $P^{-1/2}$ for which both the computation of $P^{-1/2}$ itself and the computation of the set $\{P^{-1/2}x_i, \ldots, P^{-1/2} x_n\}$ can be carried out in time $\tilde{O}(ndk)$. We will soon elaborate more on the considerations when choosing the preconditioner, but first we would like to address some important implementation issues.

Preconditioned regularization {#sec:regCond}
-----------------------------

In order to implement the preconditioning scheme suggested above, we should be able to find a simple form for the function $\tilde{L}$. In particular, since we would like to use SVRG, we should write $\tilde{L}$ as an average of $n$ components whose gradients can be easily computed. Denote by ${\tilde{x}}_i = P^{-1/2}x_i$ for all $i \in [n]$. Since for every $i \in [n]$, $((P^{-1/2} w)^\top x_i-y_i)^2=(w^\top {\tilde{x}}_i  - y_i)^2$, it seems natural to write $\tilde{L}(w)=L(P^{-1/2}w)$ as follows: $$\tilde{L}(w) =  \frac{1}{n} \sum_{i=1}^n \underbrace{\frac{1}{2}(w^\top {\tilde{x}}_i-y_i)^2}_{=:\tilde{\ell}_i}  +\frac{\lambda}{2} \|P^{-1/2}w\|^2~.$$ Assume momentarily that $\lambda=0$. Note that the gradient of $\tilde{\ell}_i$ at any point $w$ is given by $\nabla \tilde{\ell}_i(w_t)=(w^\top {\tilde{x}}_i-y_i){\tilde{x}}_i$. Hence, by computing all the ${\tilde{x}}_i$’s in advance, we are able to apply SVRG directly to the preconditioned function and computing the stochastic gradients in time $O(d)$.

When $\lambda>0$, the computation of the gradient at some point $w$ involves the computation of $P^{-1} w$. We would like to avoid this overhead. To this end, we decompose the regularization function as follows. Denote the standard basis of $\reals^d$ by $e_1,\ldots,e_d$. Note that the function $L$ can be rewritten as follows: $$\begin{aligned}
 L(w) = \frac{1}{n+d} \sum_{i=1}^{n+d} \ell_i(w)~,\end{aligned}$$ where $\ell_i(w) = \frac{n+d}{n} \frac{1}{2} (w^\top x_i-y_i)^2$ for $i=1,\ldots,n$ and $\ell_{n+i}(w)=\lambda(n+d) \frac{1}{2}(w^\top e_i)^2$ for $i=1,\ldots,d$. Finally, denoting $b_i = P^{-1/2} e_i$ for all $i$, we can rewrite the preconditioned function $\tilde{L}$ as follows: $$\begin{aligned}
 {\tilde{L}}(w) = \frac{1}{n+d} \sum_{i=1}^{n+d} {\tilde{\ell}}_i(w)~,\end{aligned}$$ where ${\tilde{\ell}}_i(w) = \frac{n+d}{n} \frac{1}{2} (w^\top {\tilde{x}}_i-y_i)^2$ for $i=1,\ldots,n$ and ${\tilde{\ell}}_{n+i}(w)=\lambda(n+d) \frac{1}{2}(w^\top b_i)^2$ for $i=1,\ldots,d$. By computing the ${\tilde{x}}_i$’s and the $b_i$’s in advance, we are able to apply SVRG while computing stochastic gradients in time $O(d)$.

The effect of conditioning
--------------------------

We are now in position to address the following fundamental question: How does the choice of the preconditioner, $P^{-1/2}$, affects the resulted average condition number of the function ${\tilde{L}}$ (\[eq:conditionSVRG\])? The following lemma upper bounds ${\hat{\kappa}}_{{\tilde{L}}}$ by the average condition number of the matrix $P^{-1/2}(C+\lambda I)P^{-1/2}$, which we denote by ${\tilde{\kappa}}$ (when the identity of the matrix $P$ is understood).

\[thm:condEffect\] Let $P^{-1/2}$ be a preconditioner. Then, the average condition number of $\tilde{L}$ is upper bounded by $${\hat{\kappa}}_{\tilde{L}} \le {\tilde{\kappa}}= \frac{\tr (P^{-1/2} (C + \lambda I ) P^{-1/2} ) } {\lambda_d(P^{-1/2} (C + \lambda I ) P^{-1/2})} ~.$$

The proof is in the appendix. Note that an optimal bound of $O(d)$ is attained by the whitening matrix $P^{-1/2}  = (C+\lambda I)^{-1/2}$.

Exact sketched conditioning
---------------------------

Our sketched preconditioner is based on a random approximation of the best rank-$k$ approximation of the data matrix. It will be instructive to consider first a preconditioner that is based on an exact rank-$k$ approximation of the data matrix. Let $X \in \reals^{d \times n}$ be the matrix whose $i$-th columns is $x_i$ and let $\bar{X}=n^{-1/2} X$. Denote by $\bar{X}= \sum_{i=1}^{\textrm{rank}(\bar{X})} \sigma_i u_i v_i^\top = U \Sigma V^\top$ the SVD of $\bar{X}$ and recall that $\bar{X}_k= \sum_{i=1}^k \sigma_i u_i v_i^\top$ is the best $k$-rank approximation of $\bar{X}$. Note that $\bar{X}\bar{X}^\top=C$ and therefore $\sigma_i^2 = \lambda_i(C)=\lambda_i$. Furthermore, the left singular vectors of $\bar{X}$, $u_1,\ldots,u_k$, coincide with the $k$ leading eigenvectors of the matrix $C$. Consider the preconditioner, $$P^{-1/2}  = \sum_{i=1}^k \frac{u_i u_i^\top}{\sqrt{\lambda_i+\lambda}} + \frac{I-\sum_{i=1}^k  u_i u_i^\top}{\sqrt{\lambda_k+\lambda}}~,$$ where $u_{k+1},\ldots, u_d$ are obtained from a completion of $u_1,\ldots,u_k$ to an orthonormal basis.

\[lem:exactCondSketch\] Let $k \in [d]$ be a parameter and assume that the regularization parameter, $\lambda$, is larger than $\lambda_d$. Using the exact sketched preconditioner, we obtain $$\label{eq:exactCondSketchCond}
{\hat{\kappa}}_{\tilde{L}} \le \frac{k \lambda_k+ \sum_{i>k} \lambda_i}{\lambda} + d ~.$$

A simple calculation shows that for $i=1,\ldots,k$, $$\lambda_{i} (P^{-1/2} (C + \lambda I ) P^{-1/2} )  = \frac{\lambda_i+\lambda}{\lambda_i+\lambda} = 1 ~.$$ Similarly, for $i=k+1,\ldots, d$, $$\lambda_{i} (P^{-1/2} (C + \lambda I ) P^{-1/2} ) = \frac{\lambda_i+\lambda}{\lambda_k+\lambda} ~.$$ Finally, $$\lambda_{d} (P^{-1/2} (C + \lambda I ) P^{-1/2}) \ge \frac{\lambda}{\lambda_k+\lambda}~.$$ Combining the above with , we obtain that $$\begin{aligned}
 {\hat{\kappa}}_{\tilde{L}} &\le \frac{\tr(P^{-1/2} (C + \lambda I ) P^{-1/2})}{\lambda_d(P^{-1/2} (C + \lambda I ) P^{-1/2})} \\
&\le  k\frac{\lambda_k+\lambda}{\lambda} + \sum_{i=k+1}^d \frac{\lambda_i+\lambda}{\lambda}  \\
&= \frac{k \lambda_k + \sum_{i>k} \lambda_i}{\lambda} + d ~. \qedhere\end{aligned}$$

Sketched conditioning {#sec:lowCond}
---------------------

An exact computation of the SVD of the matrix $\bar{X}$ takes $O(nd^2)$. Instead, we will use the Block Lanczos method in order to approximate the truncated SVD of $\bar{X}$. Specifically, given a parameter $k \in [d]$, we invoke the Block Lanczos method with the parameters $\bar{X},k$ and $\epsilon'=1/2$. Recall that the output has the form $\tilde{X}_k = {\tilde{U}}_k {\tilde{\Sigma}}_k {\tilde{V}}_k^\top= \sum_{i=1}^k {\tilde{\sigma}}_i{\tilde{u}}_i {\tilde{v}}_i^\top$. Analogously to the exact sketched preconditioner, we define our sketched preconditioner by $$\label{eq:condSketch}
P^{-1/2}  = \sum_{i=1}^k \frac{{\tilde{u}}_i {\tilde{u}}_i^\top}{\sqrt{{\tilde{\sigma}}_i^2+\lambda}}  + \frac{I-\sum_{i=1}^k  {\tilde{u}}_i {\tilde{u}}_i^\top}{\sqrt{{\tilde{\sigma}}_k^2+\lambda}} ~.$$

\[thm:sketchCond\] Let $k \in [d]$ be a parameter and assume that the regularization parameter, $\lambda$, is larger that $\lambda_d$. Using the sketched preconditioner defined in , up to a multiplicative constant, we obtain the bound on the average condition number with probability at least $9/10$.

The rest of this section is devoted to the proof of . We follow along the lines of the proof of . Up to a multiplicative constant, we derive the same upper and lower bounds on the eigenvalues of $P^{-1/2} (C+\lambda I)P^{-1/2}$.

From now on, we assume that the bounds in (where $\epsilon'=1/2$) hold. This assumption will be valid with probability of at least $9/10$. We next introduce some notation. We can rewrite $P^{-1/2}  = {\tilde{U}}({\tilde{\Sigma}}^2+\lambda I)^{-1/2} {\tilde{U}}^\top$ where ${\tilde{\Sigma}}$ is a diagonal $d \times d$ with ${\tilde{\Sigma}}_{i,i} = {\tilde{\sigma}}_i$ if $i \le k$ and ${\tilde{\Sigma}}_i = {\tilde{\sigma}}_k$ if $i>k$. and the columns of ${\tilde{U}}$ are a completion of ${\tilde{u}}_1,\ldots,{\tilde{u}}_k$ to an orthonormal basis. Recall that the SVD of $\bar{X}$ is denoted by $\bar{X} = \sum_{i=1}^d \sigma_i u_i v_i^\top = U \Sigma V^\top$.

**(Upper bound on the leading eigenvalue)** We have $$\lambda_1(P^{-1/2} (C+\lambda I)P^{-1/2})  \le 17~.$$

Since $\lambda_1(P^{-1/2} (C+\lambda I)P^{-1/2}) = \|P^{-1/2} (C+\lambda I)P^{-1/2}\| = \|P^{-1/2} C P^{-1/2} + \lambda P^{-1}\|$, using the triangle inequality we have that $$\lambda_1(P^{-1/2} (C+\lambda I)P^{-1/2}) ~\le~ \|P^{-1/2} C P^{-1/2}\| + \lambda \|P^{-1}\| ~.$$ By the definition of $P$ we have that $\|P^{-1}\| = \frac{1}{{\tilde{\sigma}}_k^2 + \lambda}$ and therefore the second summand on the right hand side of the above is at most $\frac{\lambda}{{\tilde{\sigma}}_k^2+\lambda} \le 1$. As to the first summand, recall that $C = \bar{X} \bar{X}^\top$ and therefore $\|P^{-1/2} C P^{-1/2}\| = \|\bar{X}^\top P^{-1/2} \|^2$. We will show that $\|\bar{X}^\top P^{-1/2} \| \le 4$ which will imply that $\|P^{-1/2} C P^{-1/2}\| \le 16$. To do so, we first apply the triangle inequality, $$\begin{aligned}
&\|\bar{X}^\top P^{-1/2}\|  = \| ({\tilde{X}}_k + (\bar{X} - {\tilde{X}}_k))^\top P^{-1/2} \| \\
&\quad \le \| {\tilde{X}}_k^\top P^{-1/2} \| + \|(\bar{X} - {\tilde{X}}_k)^\top P^{-1/2} \| ~.\end{aligned}$$ Let us consider one term at the time. Recall that ${\tilde{X}}_k = {\tilde{U}}_k {\tilde{\Sigma}}_k {\tilde{V}}_k^\top$. Since ${\tilde{U}}_k^\top {\tilde{U}}\in \reals^{k,d}$ is a diagonal matrix with ones on the diagonal, and since the spectral norm is invariant to multiplication by unitary matrices, we obtain that $$\begin{aligned}
&\|{\tilde{X}}_k^\top P^{-1/2}\| = \|{\tilde{V}}_k {\tilde{\Sigma}}_k {\tilde{U}}_k^\top {\tilde{U}}({\tilde{\Sigma}}^2+\lambda I)^{-1/2} {\tilde{U}}^\top\| \\
&\quad = \|{\tilde{\Sigma}}_k {\tilde{U}}_k^\top {\tilde{U}}({\tilde{\Sigma}}^2+\lambda I)^{-1/2} \| \\
&\quad = \max_{i \in [k]} \frac{{\tilde{\sigma}}_i}{ \sqrt{{\tilde{\sigma}}_i^2+\lambda}} 
\le \max_{i \in [k]} \frac{{\tilde{\sigma}}_i}{{\tilde{\sigma}}_i+\sqrt{\lambda}} 
\le 1 ~.\end{aligned}$$ Next, by the submutiplicativity of the spectral norm, $$\|(\bar{X} - {\tilde{X}}_k)^\top P^{-1/2} \| \le \|\bar{X}-{\tilde{X}}_k\| \cdot \|P^{-1/2}\| ~.$$ implies that $\|\bar{X}-{\tilde{X}}_k\| \le \tfrac{3}{2} \sigma_k$ and $$\begin{aligned}
&\|P^{-1/2}\| = \frac{1}{ \sqrt{{\tilde{\sigma}}_k^2+\lambda}} 
\le \frac{1}{\sqrt{{\tilde{\sigma}}_k^2}} \le \frac{1}{ \sqrt{\sigma_k^2 - (1/2) \sigma_{k+1}^2}} \\
&\quad \le \frac{1}{ \sigma_{k} \sqrt{\tfrac{1}{2}} } = \frac{ \sqrt{2}}{ \sigma_k} 
<\frac{ 2}{\sigma_k}~.\end{aligned}$$ Hence, $\|\bar{X}-{\tilde{X}}_k\| \cdot \|P^{-1/2}\| \le 3$. Combining all of the above bounds concludes our proof.

**(Refined upper bound on the last $d-k$ eigenvalues)** For any $i \in \{k+1,\ldots,d\}$, $$\lambda_i \left(P^{-1/2} (C+ \lambda I)P^{-1/2} \right) \le \frac{2 (\lambda_i+\lambda)}{\lambda_k+\lambda}~.$$

Using the Courant minimax principle [@bhatia2013matrix], we obtain the following bound for all $i \in \{k+1,\ldots,d\}$: $$\begin{aligned}
&\lambda_i \left(P^{-1/2} (C+ \lambda I ) P^{-1/2} \right) \\
&= \max_{\substack{\mathcal{M} \subseteq \reals^d:\\ \dim(\mathcal{M}) =i}}\min_{\substack{x \in \cM:\\x \neq 0}} \frac{  x^\top P^{-1/2} (C+\lambda I) P^{-1/2} x}{\|x\|^2} \\
& = \max_{\substack{\mathcal{M} \subseteq \reals^d:\\ \dim(\mathcal{M}) =i}}\min_{\substack{x \in \cM:\\ x \neq 0}} \frac{  x^\top P^{-1/2} (C+\lambda I) P^{-1/2} x}{\|P^{-1/2}x\|^2} \cdot \frac{\|P^{-1/2}x\|^2} {\|x\|^2} \\
& \le \left(\max_{\substack{\mathcal{M} \subseteq \reals^d:\\ \dim(\mathcal{M}) =i}}\min_{\substack{x \in \cM: \\x \neq 0}} \frac{  x^\top P^{-1/2} (C+\lambda I) P^{-1/2} x}{\|P^{-1/2}x\|^2}\right) \times\\
& \qquad \left( \max_{\substack{x \in \reals^d:\\ x \neq 0}} \frac{x^\top P^{-1} x}{\|x\|^2} \right)\\
&= \lambda_i \left(C+\lambda  I \right) \cdot \lambda_1(P^{-1})  = (\lambda_i+\lambda) \cdot (\tilde{\sigma}_k^2+\lambda)^{-1} ~.
$$ Finally, using we have that ${\tilde{\sigma}}_k^2 \ge \sigma_k^2 - \frac{1}{2} \sigma_{k+1}^2 \ge \tfrac{1}{2} \sigma_k^2 = \tfrac{1}{2} \lambda_k$ and therefore, $$(\tilde{\sigma}_k^2+\lambda)^{-1} ~\le~ (\tfrac{1}{2} \lambda_k +\lambda)^{-1} ~\le~ 2\, (\lambda_k+\lambda)^{-1} ~. \qedhere$$

\[lem:lowerSmall\] **(Lower bound on the smallest eigenvalue)** $$\lambda_d( P^{-1/2} C P^{-1/2}) \ge \frac{\lambda}{19 (\lambda_k+\lambda)}~.$$

Note that $$\label{eq:digging0}
\lambda_d (P^{-1/2} (C+\lambda I) P^{-1/2}) = \frac{1}{\|P^{1/2} (C+\lambda I)^{-1} P^{1/2}\|}~,$$ so we can derive an upper bound on $\|P^{1/2} (C+\lambda I)^{-1} P^{1/2}\|$. Consider an arbitrary completion of $\tilde{v}_1,\ldots,\tilde{v}_k$ to an orthonormal set, $\tilde{v}_1,\ldots,\tilde{v}_d \in \reals^n$. Let $\tilde{V} \in \reals^{n \times d}$ be the matrix whose $i$-th column is ${\tilde{v}}_i$. Since the spectral norm is unitary invariant and both ${\tilde{U}}$ and ${\tilde{V}}$ have orthonormal columns, $$\begin{aligned}
 &\|P^{1/2} (C+\lambda I)^{-1} P^{1/2}\| \notag\\
&= \|{\tilde{U}}({\tilde{\Sigma}}^2+\lambda I)^{1/2} {\tilde{U}}^\top (C+\lambda I)^{-1} {\tilde{U}}({\tilde{\Sigma}}^2+\lambda I)^{1/2} {\tilde{U}}^\top\| \notag \\
&= \|{\tilde{V}}({\tilde{\Sigma}}^2+\lambda I)^{1/2} {\tilde{U}}^\top (C+\lambda I)^{-1} {\tilde{U}}({\tilde{\Sigma}}^2+\lambda I)^{1/2} {\tilde{V}}^\top\| ~.\end{aligned}$$ Denote by ${\tilde{Z}}= {\tilde{U}}({\tilde{\Sigma}}^2+\lambda I)^{1/2} {\tilde{V}}^\top$. By the triangle inequality and the submutiplicativity of the spectral norm, $$\begin{aligned}
&\|{\tilde{Z}}^\top(C+\lambda I)^{-1} {\tilde{Z}}\| \le \|\bar{X}^\top (C+\lambda I)^{-1} \bar{X}\| \notag\\
&\qquad + \|({\tilde{Z}}-\bar{X})^\top (C+\lambda I)^{-1} ({\tilde{Z}}-\bar{X})\| \notag \\
& \le \|\bar{X}^\top (C+\lambda I)^{-1} \bar{X}\|+  \|{\tilde{Z}}-\bar{X}\|^2 \|(C+\lambda I)^{-1}\|~. \label{eq:digging15}\end{aligned}$$ To bound the first summand of , we use the unitary invariance to obtain $$\begin{aligned}
 &\|\bar{X}^\top (C+\lambda I)^{-1} \bar{X} \| = \|V \Sigma U^\top U (\Sigma^2 +\lambda I)^{-1} U^\top U \Sigma V^\top\| \notag\\
&\quad = \|\Sigma  (\Sigma^2 +\lambda I)^{-1} \Sigma \| = \max_i \frac{\lambda_i^2}{\lambda_i^2+\lambda} \le 1~.\end{aligned}$$ For the second summand of , note that $\|(C+ \lambda I)^{-1}\|= \frac{1}{\lambda_d+\lambda}$ and that, using the triangle inequality, $$\begin{aligned}
&\|{\tilde{Z}}- \bar{X}\| = \|({\tilde{U}}{\tilde{\Sigma}}{\tilde{V}}^\top - \bar{X}) + ({\tilde{Z}}-  {\tilde{U}}{\tilde{\Sigma}}{\tilde{V}}^\top) \|\\
&\quad \le \|{\tilde{U}}{\tilde{\Sigma}}{\tilde{V}}^\top - \bar{X}\| + \|{\tilde{U}}(({\tilde{\Sigma}}^2+\lambda I)^{1/2}- {\tilde{\Sigma}}) {\tilde{V}}^\top\|~.\end{aligned}$$ By using unitary invariance together with the inequality $\sqrt{{\tilde{\sigma}}_i^2+\lambda} - {\tilde{\sigma}}_i \le \sqrt{\lambda}$ (which holds for every $i$), we get $$\|{\tilde{U}}(({\tilde{\Sigma}}^2+\lambda I)^{1/2}- {\tilde{\Sigma}}) {\tilde{V}}^\top\| = \|({\tilde{\Sigma}}^2+\lambda I)^{1/2}- {\tilde{\Sigma}}\| \le \sqrt{\lambda}~.$$ Hence, using the inequality $(x+y)^2 \le 2x^2 + 2y^2$, we obtain $$\begin{aligned}
  \|{\tilde{Z}}- \bar{X}\|^2 &\le  2 \|{\tilde{U}}{\tilde{\Sigma}}{\tilde{V}}^\top - \bar{X}\|^2 + 2\lambda ~.\end{aligned}$$ We next derive an upper bound on $\|{\tilde{U}}{\tilde{\Sigma}}{\tilde{V}}^\top - \bar{X}\|$. Since ${\tilde{U}}{\tilde{\Sigma}}{\tilde{V}}^\top = \tilde{X}_k + {\tilde{\sigma}}_k \sum_{i=k+1}^d {\tilde{u}}_i {\tilde{v}}_i^\top$, $$\|{\tilde{U}}{\tilde{\Sigma}}{\tilde{V}}^\top-\bar{X}\| \le \|\tilde{X}_k - \bar{X}\|+{\tilde{\sigma}}_k \left\|\sum_{i=k+1}^d {\tilde{u}}_i {\tilde{v}}_i^\top\right\| ~.$$ Using we know that $\|\tilde{X}_k - \bar{X}\| \le 1.5 \, \sigma_k$ and that ${\tilde{\sigma}}_k \le \sqrt{\sigma_k^2+0.5 \, \sigma_{k+1}^2} \le  1.5 \, \sigma_k$. Combining this with the fact that $\|\sum_{i=k+1}^d {\tilde{u}}_i {\tilde{v}}_i^\top\| = 1$, we obtain $$\begin{aligned}
 \|{\tilde{U}}{\tilde{\Sigma}}{\tilde{V}}^\top-\bar{X}\| \le 3\,\sigma_k ~.\end{aligned}$$ Combining the above inequalities, we obtain $$\begin{aligned}
\|P^{1/2} (C+\lambda I)^{-1} P^{1/2}\| 
&\le  1 + \frac{2 \cdot(3 \sigma_k)^2+2\lambda}{\lambda_d+\lambda} \\
&\le \frac{19(\lambda_k+\lambda)}{\lambda} ~,\end{aligned}$$ and using we conclude our proof.

**(of )** The three last lemmas imply that the inequalities derived during the proof of remain intact up to a multiplicative constant. Therefore, the bound on the condition number also holds up to a multiplicative constant. This completes the proof.

Sketched Preconditioned SVRG
----------------------------

By equipping SVRG with the sketched preconditioner , we obtain the Sketched Preconditioned SVRG (see ).

**(of )** The theorem follows from and .

**Input:** $x_1,\ldots, x_n \in \reals^d, y_1,\ldots,y_n \in \reals, \epsilon>0$ **Parameters:** $\lambda>0, k \in [d]$ Let $\bar{X} \in \reals^{d,n}$ be the matrix whose $i$’th column is $(1/n) x_i$ Run the Block Lanczos method () with the input $\bar{X}, k, \epsilon'=1/2$ to obtain ${\tilde{X}}_k = {\tilde{U}}_k {\tilde{\Sigma}}_k {\tilde{V}}_k$ Let ${\tilde{u}}_i$ be the columns of ${\tilde{U}}_k$ and $\tilde{\sigma}_i$ be the diagonal elements of ${\tilde{\Sigma}}_k$ Form the preconditioner $P^{-1/2}$ according to Compute ${\tilde{x}}_i = P^{-1/2} x_i, b_i = P^{-1/2} e_i$ Let $\ell_i(w) = \frac{n+d}{n} \frac{1}{2} (w^\top {\tilde{x}}_i-y_i)^2$ for $i=1,\ldots,n$ and $\ell_i(w) = \lambda(n+d) (w^\top b_i)^2$ for $i=n+1,\ldots,n+d$ Let $\beta_i = \frac{n+d}{n} \|{\tilde{x}}_i\|^2$ for $i=1,\ldots,n$ and $\beta_i  = \lambda(n+d) \|b_i\|$ for $i=n+1,\ldots,n+d$. Let $\hat{\beta} = \frac{1}{n} \sum_{i=1}^{n+d} \beta_i$ Run SVRG () Return $\hat{w} = P^{1/2} \tilde{w}$

The Empirical Gain of Sketched Preconditioning {#sec:empirical}
==============================================

\[b\] [0.23 ]{}

\[b\] [0.23]{}

\[b\] [0.23]{}

\[b\] [0.23]{}

[0.23 ]{}

[0.23 ]{}

[0.23 ]{}

[0.23 ]{}

[0.23 ]{}

[0.23 ]{}

In this section we empirically demonstrate the gain of our method. We consider both regression problems and binary classifications tasks, where the square loss serves as a surrogate for the zero-one loss. We use the following datasets:

-   *Synthetic*: We draw two random $5000 \times 20000$ matrices, $X^{(1)}$ and $X^{(2)}$, whose singular vectors are drawn uniformly at random and the $q$-th singular value is $1/q$ and $1/q^2$, respectively. We then normalize the columns. For each $X=X^{(j)}$, we consider a regression problem, where the labels are generated as follows: we first draw a vector $w^\star \in \cN(0,1)^{5000}$ and then set $y_i = {w^\star}^\top X_{\cdot,i}+z_i$, where $z_i \sim \cN(0,0.1)$.

-   *MNIST*:[^5] A subset of MNIST, corresponding to the digits $4$ and $7$, where the task is to distinguish between the two digits. Here, $n=12107, d=784$.

-   *RCV1*:[^6] The Reuters RCV1 collection. Here, $n=20242, d=47236$ and we consider a standard binary document classification task.

-   *CIFAR-10*:[^7] Here, $n=50000, d=3072$. Following [@frostig2015regularizing], the classification task is to distinguish between the animal categories to the automotive ones.

-   *real-sim*:[^8] Here, $n=72309, d=20958$, and we consider a standard binary document classification task.

Inspecting our theoretical speed-up {#sec:ratioEmp}
-----------------------------------

Recall that the ratio quantifies our theoretical speedup. Hence, we first empirically inspect the prefixes of the corresponding quantities (as a function of $k$) for each of the datasets (see ). We can see that while in MNIST and CIFAR-10 the ratio is large for small values of $k$, in RCV1 and real-sim the ratio increases very slowly (note that for the former two datasets we use logarithmic scale).

Empirical advantage of Sketched Preconditioned SVRG {#sec:emp}
---------------------------------------------------

We now evaluate and compare it to the SVRG algorithm of [@xiao2014proximal]. To minimally affect the inherent condition number, we added only a slight amount of regularization, namely, $\lambda=10^{-8}$. The loss used is the square loss. The step size, $\eta$, is optimally tuned for each method. Similarly to previous work on SVRG [@xiao2014proximal; @johnson2013accelerating], the size of each epoch, $m$, is proportional to the number of points, $n$. We minimally preprocessed the data by average normalization: each instance vector is divided by the average $\ell_2$-norm of the instances. The number of epochs is up to $60$. Note that in all cases we choose a small preconditioning parameter, namely $k=30$, so that the preprocessing time of is negligible. There is a clear correspondence between the ratios depicted in and the actual speedup. In other words, the empirical results strongly affirm our theoretical results.

Acknowledgments {#acknowledgments .unnumbered}
===============

We thank Edo Liberty for helpful discussions. The work is supported by ICRI-CI and by the European Research Council (TheoryDL project).

Omitted Proofs
==============

**(of )** We first show that the average smoothness of $\tilde{L}$ is bounded by $$\label{eq:condEffectSmooth}
\frac{1}{n+d} \sum_{i=1}^{n+d} \tilde{\beta}_i \le  \tr \left(P^{-1/2} \left(C + \lambda I \right) P^{-1/2} \right)~.$$ Note that for any $w$, $$\nabla^2 \tilde{\ell}_i(w) = \begin{cases} \frac{n+d}{n}{\tilde{x}}_i {\tilde{x}}_i^\top  & 1 \le i \le n, \\ \lambda (n+d) b_{i-n} b_{i-n}^\top  & n < i \le n+d~.  \end{cases}$$ Therefore, using the fact that the spectral norm of a rank-$1$ psd matrix is equal to its trace, we obtain $$\begin{aligned}
\frac{1}{n+d} \sum_{i=1}^n \tilde{\beta}_i &= \frac{1}{n+d} \frac{n+d}{n} \sum_{i=1}^n \|{\tilde{x}}_i {\tilde{x}}_i^\top\| + \frac{1}{n+d}\lambda (n+d) \sum_{j=1}^d \|b_j b_j^\top\| \\
& = \frac{1}{n} \sum_{i=1}^n \tr({\tilde{x}}_i {\tilde{x}}_i^\top) + \lambda \sum_{j=1}^d \tr(b_i b_i^\top) \\
&= \frac{1}{n}\tr( \sum_{i=1}^n P^{-1/2} x_i x_i^\top P^{-1/2}) + \lambda \, \tr \sum_{j=1}^d (P^{-1/2} e_i e_i^\top P^{-1/2}) \\
&=\tr(P^{-1/2} ( C + \lambda I)P^{-1/2})~.\end{aligned}$$ Hence, we deduce (\[eq:condEffectSmooth\]).

We will conclude the theorem by showing that $\tilde{L}$ is $\lambda_d (P^{-1/2} (C + \lambda I)  P^{-1/2})$-strongly convex. Indeed, a similar calculation shows that the Hessian of $L$ at any point $w$ is given by $$\begin{aligned}
\nabla^2 \tilde{L}(w) &=  P^{-1/2} (C+ \lambda I) P^{-1/2} ~.\end{aligned}$$ Hence, we conclude the claimed bound.

[^1]: School of Computer Science, The Hebrew University, Jerusalem, Israel

[^2]: Yahoo Labs, New York, NY, USA

[^3]: School of Computer Science, The Hebrew University, Jerusalem, Israel

[^4]: More precisely, Algorithm 2 in [@musco2015stronger] returns the projection matrix ${\tilde{U}}_k {\tilde{U}}_k^\top$, while we also compute the SVD of ${\tilde{U}}_k {\tilde{U}}_k^\top A$. The additional runtime is negligible.

[^5]: http://yann.lecun.com/exdb/mnist/

[^6]: https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/

[^7]: http://www.cs.toronto.edu/ kriz/cifar.html

[^8]: https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/
