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 