This paper studies an intriguing phenomenon related to the good generalization performance of estimators obtained by using large learning rates within gradient descent algorithms.
First observed in the deep learning literature, we show that
such a phenomenon can be precisely characterized in the context of kernel methods, even though the resulting optimization problem is convex. Specifically, we consider the minimization of a quadratic objective in a separable Hilbert space, and show that with early stopping, the choice of learning rate influences the spectral decomposition of the obtained solution on the Hessian's eigenvectors. This extends an intuition described by Nakkiran, 2020 on a two-dimensional toy problem to realistic learning scenarios such as kernel ridge regression. While large learning rates may be proven beneficial as soon as there is a mismatch between the train and test objectives, we further explain why it already occurs in classification tasks without assuming any particular mismatch between train and test data distributions.

Optimization; Statistical Learning; Kernel methods

§ INTRODUCTION

Gradient descent methods are omnipresent in machine learning, and a lot of effort has been devoted to better understand their theoretical properties. Optimal rates of convergence have been well characterized for minimizing convex functions in various contexts, including, for instance, stochastic optimization [Nemirovski et al., 2009].
For supervised learning, one is however more interested in the statistical optimality of the resulting estimator rather than in the ability to quickly optimize a training objective [Bottou and Bousquet, 2007]. When considering both optimization and statistical questions, gradient descent methods were proven to be optimal under many assumptions [Yao et al., 2007, Pillaud-Vivien et al., 2018].

An important observation for this paper is that gradient descent algorithms typically require to tune some learning rate, or step size, to achieve the best performance. This has been thoroughly investigated in the optimization literature. For convex smooth problems in particular, the influence of step size on convergence rates is well understood [Nesterov, 2018]. However, recent empirical studies have highlighted a surprising aspect of this parameter: when using gradient descent methods on neural networks, large learning rates were found to be useful for obtaining good generalization properties, or in other words, good statistical performance [Jastrzebski et al., 2021], even though they may be sub-optimal from an optimization point of view.

This paper aims at understanding this phenomenon from a broad but simple perspective, where both the function $\LEmp$ we optimize and the function $\LPop$ used to evaluate the statistical performance are quadratic forms of some separable Hilbert space $\HH$. Specifically, we assume that

\begin{equation}\label{eq:intro_model}
\forall \theta \in \HH, ~~ \LEmp(\theta) = \frac{1}{2} \norm{\theta - \optEmp}_\TEmp^2 + \minLEmp, ~~~\text{and}~~~ \LPop(\theta) = \frac{1}{2} \norm{\theta - \optPop}_\TPop^2,
\end{equation}

where $\norm{\cdot}_A$ denotes the norm on $\HH$ induced by a positive definite operator $A$, i.e., $\norm{\theta}_A^2 = \dotprod{\theta}{A \theta}$ for any $\theta$ in $\HH$. With <ref>, $\LEmp$ and $\LPop$ are characterized by positive definite operators $\TEmp$ and $\TPop$, along with their minimizers denoted by $\optEmp$ and $\optPop$, respectively. The constant value $\minLEmp$ does not affect the optimization problem and can be safely ignored in the rest of this presentation.

The model from <ref> captures a large class of problems such as learning with kernels, detailed in <ref>, but we give here a simple example with ridge regression.

[$\TEmp$ and $\TPop$ with Ridge Regression.]
Let $x_1, \dots, x_n$ be data points in $\RR^d$, and $y_1, \dots, y_n$ prediction variables, with $n \geq d$. Define $X \in \RR^{n \times d}$ the data matrix. We consider the ridge regression estimator with regularization $\lambda>0$, which is defined as the minimum of

\begin{equation*}
\forall \theta \in \RR^d, ~~ \LEmp(\theta) = \frac{1}{n} \sum_{i=1}^n \frac{1}{2} (\theta^\top x_i - y_i)^2 + \frac{\lambda}{2} \norm{\theta}^2 = \frac{1}{2n}\norm{X \theta - y}^2 + \frac{\lambda}{2} \norm{\theta}^2.
\end{equation*}

$\LEmp$ is a quadratic function of $\theta$, which can be rewritten as in <ref> with

\begin{equation*}
\forall \theta \in \RR^d, ~~ \LEmp(\theta) = \frac{1}{2} \norm{\theta - \optPop}_{\TEmp}^2 + \minLEmp, ~~ \text{with} ~~ \begin{cases}
\optPop & = \frac{1}{n} \p{\frac{1}{n} X^\top X + \lambda \II_d}^{-1} X^\top y \\
\TEmp   & = \frac{1}{n} (X^\top X + \lambda \II_d).                            \\
% \minLEmp &= \frac{1}{2n} y^\top \csb{\II_n - \frac{1}{n}X \p{\frac{1}{n} X^\top X + \lambda}^{-1} X^\top}  y.
\end{cases}
\end{equation*}

Assuming that the output can be written $y_i = x_i^\top \nu^\star + \epsilon_i$ with $\nu^\star \in \RR^d$ and $\epsilon_i$ some independent, zero-mean noise, then the population loss is $\calP(\theta) = \EE \nicefrac12 \p{\theta^\top x_i - y_i}^2$, and the excess risk defined by $\LPop(\theta) = \calP(\theta) - \inf_{\nu} \calP(\nu)$ is given with

\begin{equation*}
\LPop(\theta) = \EE \csb{\frac{1}{2} \p{(\theta-\nu^\star)^\top x}^2} = \frac{1}{2} \norm{\theta - \optPop}_\TPop^2, ~~ \text{with} ~~ \TPop = \EE \left[x x^\top\right].
\end{equation*}

In this example, a discrepancy between train and test losses (between $\TEmp$ and $\TPop$) may occur in particular situations (e.g., presence of data augmentation during training, or simply mismatch between train and test distributions). The next example shows that such a mismatch may be in fact frequent for classification problems, even when train and test distributions do match.

[Discrepancy between train and test losses in classification with separable classes.]
The scenario described in <ref> is particularly evident in the context of binary classification, when the classes are separable by a non-zero margin.
This is considered a typical situation in many learning scenarios of interest, as classification over natural images—motivating the wide use of large-margin based classifiers in the field [Rawat and Wang, 2017]. We highlight here, that in this context, the loss we are using for training is not the best loss to consider for the test error, as discussed next.

More precisely, consider a classification problem with two classes with non-zero margin. Let $\XX \subseteq \mathbb{R}^d$ and $\YY =\{-1,1\}$ be the input and the output space. Denote by $\rho(x,y) = \rho_\XX(x) \rho(y|x)$ the probability distribution describing the classification problem, where $\rho_\XX$ is the marginal probability over $\XX$, while $\rho(y|x)$ is the conditional probability of $y$ given $x$. The error that we would like to minimize is the binary error on the population, i.e. $B(\theta) = \mathbb{P}[ \sign[\theta(x)] \neq y]$ for a model $\theta$. Let $\nu^\star$ be a function minimizing the binary error and ${\cal H}$ be the class of models under consideration. Assume, for simplicity, that ${\cal H}$ is a RKHS with norm  $\|\cdot\|$ [Aronszajn, 1950], i.e., there exists a map $\phi:\XX \to {\cal H}$ such that the function in ${\cal H}$ are characterized as $\theta(x) = \langle \phi(x), \theta\rangle$, for any $\theta$ in ${\cal H}$.

It has been shown by Pillaud-Vivien et al., 2018 (in particular, Lemma 1 and Appendix A, Theorem 13), that in the context of two classes separated by a non-zero margin and whose conditional probability is regular enough, then $\nu^\star$ is in $\HH$ and moreover $B(\theta) - B(\nu^\star) \leq e^{-c/\|\theta - \nu^\star\|}$ for some constant $c$.
Therefore, the binary error decreases exponentially in terms of the Hilbert norm $\|\cdot\|$. On the other hand, the norm minimized at training time is some smooth convex surrogate of the binary loss, as, for example, the quadratic loss.

In this case, the trained vector may be obtained by minimizing the population loss (in fact, a regularized empirical version, but we omit this fact here for simplicity) such that $F(\theta) = \mathbb{E} (\theta(x) - y)^2$.
Noting that $\int (\theta(x) - \optEmp(x))^2 \dd\rho(x) = \|\TEmp^{1/2}(\theta - \optEmp)\|^2 = \|\theta - \optEmp\|_\TEmp$ for $\TEmp = \int \phi(x) \phi(x)^\top \dd \rho_\XX(x)$, and $\optEmp = \nu^\star$ under the considered conditions (see the same paper), we have

\begin{equation*}
F(\theta) - F(\optEmp) = \int (\theta(x) - \optEmp(x))^2 \dd\rho_\XX(x) = \|\theta - \optEmp\|^2_{\TEmp}.
\end{equation*}

This is a typical case, where there is a discrepancy between the error of interest

\begin{equation*}
R(\theta) - R(\nu^\star) = \|\theta - \nu^\star\|^2,
\end{equation*}

which, if optimized, would lead to an exponential decrease of the classification error, and the loss $F$ that is instead optimized by the algorithm at training time, for which we have only the slower rate $B(\theta) - B(\nu^\star) \leq (F(\theta) - F(\optPop))^\alpha$, with $\alpha \in (1/2, 1)$ (see, e.g. Audibert, 2004 or Audibert and Tsybakov, 2007 for what concerns the CAR assumption).

In this paper, we are interested in understanding in which regime large learning rates with early stopping could be useful for kernel methods, even if they are suboptimal from an optimization point of view.
We consider indeed the optimization of $\LEmp$ in <ref> with plain gradient descent, starting from a vector $\theta_0$ in $\HH$ with step-size $\eta$, and

we distinguish between two cases: having a small learning rate $\eta_s$ or a large learning rate $\eta_b$, the range of both is to be detailed later.

A simple intuition was suggested by Nakkiran, 2020 on a two-dimensional toy problem, showing that large learning rates may be beneficial as soon as there is a mismatch between $F$ and $R$ (meaning, what we train on does not correspond to what we test on). We show that such an insight can be extended beyond toy problems to realistic scenarios with traditional kernel methods, and that, perhaps surprisingly, this phenomenon occurs already in simple classification tasks.

Under a few assumptions described later in this paper, consider the target accuracy $\alpha$ and large and small step sizes $\theta_b$ and $\theta_s$ (these quantities being defined in the aforementioned assumptions).
Consider the gradient descent iterations $ \theta_{t+1} = \theta_{t} - \eta \TEmp (\theta_{t} - \optEmp)$ either with step size $\eta=\eta_b$ or $\eta=\eta_s$, and stop the procedure as soon as $F(\theta_{t+1}) \leq \alpha$, resulting in two estimators $\theta_b$ or $\theta_s$. Then,

\begin{equation}\label{eq:informal_thm_bound}
\LPop(\theta_b) - \LPop(\optPop) ~\leq~ 34 \frac{\kappa_{\TPop}}{\kappa_\TEmp} (\LPop(\theta_s)- \LPop(\optPop) ),
\end{equation}

where $\kappa_{\TPop}$ and $\kappa_{\TEmp}$ are the condition numbers of the operators $\TPop$ and $\TEmp$, respectively, restricted to $\HH_n$.

Note that $\LPop(\optPop) = 0$ by definition of $R$; we made this quantity explicit in the bound for clarity purposes. The main conclusion from the theorem is that with early stopping (and with a target accuracy that is reasonable according to statistical learning theory, as discussed later), large learning rates can provide better estimators than small ones, even though the quantity $\eta_s$ may yield much faster convergence for minimizing the objective function $F$ than $\eta_b$ (a fact also discussed later).
This phenomenon occurs when the condition number $\kappa_\TEmp$ is much larger than $\kappa_\TPop$, which may already arise in classification tasks, as mentioned earlier in Example <ref> where $\kappa_\TPop = 1$.

Note that <ref> raises several questions and could be easily misinterpreted, since such a relation may suggest that an arbitrarily small risk $R(\theta_b)$ could be obtained by considering minimization problems that are arbitrarily badly conditioned. Unfortunately, but not surprisingly, there is however no free lunch here, as discussed in the next remark.

A naive observation is that $R(\theta_b)$ could be arbitrarily small by making the problem more ill-conditioned. However, the bound on $\LPop(\theta_b)$ in <ref> is relative to $\LPop(\theta_s)$. Notably, a careful reading of the proof shows that $\LPop(\theta_s)$ is an increasing function of the conditioning number, for the chosen level sets $\alpha$.

Summary of contributions.
Our first contribution is the relation described in <ref>, highlighting potential benefits of large learning rate strategies when the training objective has a worse condition number than the one used to evaluate the quality of the estimator. This is illustrated in <ref>, a figure inspired by Nakkiran, 2020.
Our second contribution is to show that such a mismatch systematically occurs in simple classification scenarios with low noise, where the quantity of interest to minimize may not be the population risk, as discussed earlier.
Overall, this allows us to show that the previous phenomenon occurs in realistic learning scenarios with kernels, which we also check in practice through numerical experiments.

[->] (-5,0) – (5,0) node[below right] $e_2$;
[->] (0,-2.5) – (0,2.5) node[left] $e_1$;
[black] (0,0) circle (1pt) node[below left] $\optEmp$;

[dashed,rotate around=60:(,)] (,)
ellipse (sqrt(/) and sqrt(/));
[dashed,rotate around=60:(,)] (,)
ellipse (sqrt(/) and sqrt(/));
[black] (,) circle (1pt) node[above right] $\optPop$;

[red] (0, 0) ellipse (sqrt(/) and sqrt(/));
(0, 0) ellipse (sqrt(/) and sqrt(/));
(0, 0) ellipse (sqrt(/) and sqrt(/));

[black] (, ) circle (1pt) node[above left] $\theta_0$;

= (1 - * )^* ;
= (1 - * )^* ;
= (1 - * )^(- 1) * ;
= (1 - * )^(- 1) * ;
] in 1,...,30
[blue] (, ) circle (1pt);
[blue, dotted] (, ) – (, );

= (1 - * )^* ;
= (1 - * )^* ;
= (1 - * )^(- 1) * ;
= (1 - * )^(- 1) * ;
] in 1,...,19
[brown] (, ) circle (1pt);
[brown, dotted] (, ) – (, );

We optimize the quadratic $\LEmp$ (level sets are filled lines, centered in $\optEmp$) with gradient descent, starting from $\theta_0$ until we reach the level sets $\alpha$ (filled line, red). However, we evaluate the quality of the estimator through $\LPop$ (level sets are dashed lines, centered in $\optPop$). Doing small step size (blue dots) optimizes the direction $e_1$ first, and yields an estimate which is far from $\optPop$ in $\TPop$ norm; doing big step size (brown dots) oscillates in the direction $e_1$, but ultimately yields an estimator which is close to $\optPop$ in $\TPop$ norm.

§ RELATED WORK

Our main motivation is to better understand the role of learning rate in obtaining good generalization for supervised learning. Even though empirical benefits of large learning rates were first described for neural networks, a few recent works have studied this phenomenon for convex problems. We review here some relevant work.

Setting the learning rate in neural networks.

Stochastic gradient descent has become the standard tool for optimizing neural networks. When the learning rate is very small, the network evolves in a so-called “lazy-regime” where its dynamics are well understood [Chizat et al., 2019, Jacot et al., 2018] but which fails to capture the good generalization performance observed with large learning rate.

Specifically, this phenomenon has been empirically observed numerous times <cit.>; common strategies consist of using first a large learning rate, before annealing it to a smaller value. As a first step towards proving theoretically the effect of choosing large learning rates for training neural networks, Li et al., 2019 devise a two-layer neural network model with different set of features where the order in which they are learnt matters, where the previous annealing strategy could be shown to be useful in theory.

A convex perspective.
Recently, different papers tried to reproduce this phenomenon in convex settings. This is probably thanks to the observation made by [Nakkiran, 2020], where a toy dataset is exhibited, which was the main motivation for this work. However, it fails to capture realistic scenarios where the data distribution is not isotropic, or with non linear data embeddings. Wu et al., 2021 aimed at filling these gaps but again, relies on the data distribution to be linear, isotropic, with the number of dimension going to infinity in order to have all data points approximately orthogonal; we do not make any of those assumptions. The bound with the condition number we obtain in <ref> is totally new. Finally, we highlight that we use plain gradient descent, and do not need stochasticity to exhibit the big learning rate phenomenon. This is consistent with recent work [Geiping et al., 2021] which shows that SGD is not necessary to obtain state of the art performances, and that GD simply needs a better fine tuning of hyper parameters.

§ MAIN RESULT

In this section, we show that by performing standard gradient descent on the empirical loss $\LEmp$, choosing a big learning rate will first optimize the smallest eigencomponent of $\TEmp$. That is, the resulting estimator is mostly located on the biggest eigenvector of $\TEmp$. On the other hand, the smaller the learning rate, the more will the solution be located on the small eigencomponents, with biggest eigenvectors of $\TEmp$ being learnt first.

§.§ Settings and notations

Gradient descent updates.
We perform standard gradient descent on the empirical loss $\LEmp$, starting from some $\theta_{0} \in \HH$, with step size $\eta$. We obtain

\begin{equation}\label{eq:gradient_descent_update}
\forall t \geq 0, ~~ \theta_{t+1} = \theta_{t} - \eta \TEmp (\theta_{t} - \optEmp), ~~ \text{thus} ~~ \theta_{t} - \optEmp = (\II - \eta \TEmp)^t (\theta_{0} - \optEmp).
\end{equation}

This enables a very simple analysis of the training in the eigenbasis of $\TEmp$. We now give a more precise definition of the model in <ref>.

[Representer theorem assumption]
There is a $n$-dimensional subspace $\HH_n \subseteq \HH$ that is invariant by $T$—that is, $T \theta$ is in $\HH_n$ for all $\theta$ in $\HH_n$ and such that $\optPop$ is in $\HH_n$.

We denote by $(\sigma_i, e_i)$ the eigenbasis of the p.d. operator $\TEmp$ restricted to $\HH_n$, with $\sigma_1 > \dots > \sigma_n > 0$, assuming eigenvalues are distinct from each other, and we call $\kappa_\TEmp = \sigma_1/\sigma_n$ the condition number.
Similarly, the restriction of $\TPop$ to $\HH_n$ is a positive definite operator whose spectrum is $\varsigma_1 > \dots > \varsigma_n$, with condition number $\kappa_\TPop = \varsigma_1 / \varsigma_n$. Since, rescaling the objectives $F$ and $R$ by constant factors does not change their minimizers, we also safely assume that $\sigma_1 = \varsigma_1 = 1$.

The model described by <ref> is quite natural, and ensures that a representer theorem holds when learning on a finite training set of $n$ points. It is notably satisfied in classical learning formulations with kernels.

With the notations of <ref>, we can now rewrite the update of <ref> along a specific direction $e_i$:

\begin{equation}\label{eq:gradient_descent_update_eigen}
\forall t \geq 0, ~~ \dotprod{\theta_{t} - \optEmp}{e_i}_\HH = (1 - \eta \sigma_i)^t \dotprod{\theta_{0} - \optEmp}{e_i}_\HH.
\end{equation}

Consider the quantity $\absv{1 - \eta \sigma_i}$. The closer to $0$, the smaller will the $i$-th component of $\theta_t - \optEmp$ on the eigenbasis be when the number of steps $t$ increases. We plot $\absv{1 - \eta \sigma_i}$ in <ref>.

xmin = 0, xmax = 2.1,
ymin = 0, ymax = 1.1,
ytick distance = 0.5,

minor tick num = 1,
width = .7,
height = 0.3,
xlabel = $\eta$,
ylabel = $\absv{1- \eta \sigma_i}$,
legend style=
anchor=south west,

nodes=scale=0.5, transform shape,
ticklabel style=font=

domain = 0:2.1,
samples = 200,
] sqrt((1 - x)^2);

domain = 0:2.1,
samples = 200,
] sqrt((1 - .9 * x)^2);

domain = 0:2.1,
samples = 200,
] sqrt((1 - .3 * x)^2);

domain = 0:2.1,
samples = 200,
] sqrt((1 - .2 * x)^2);

coordinates (2, 0) (2, 1.1);

coordinates (2/(1+.2), 0) (2/(1+.2), 1.1);

$\sigma_1 = 1.$,
$\sigma_2 = 0.9$,
$\sigma_{n-1} = 0.3$,
$\sigma_n = 0.2$,

Attenuation coefficient $\absv{1 - \eta \sigma_i}$ function of the step size $\eta$, for 4 eigenvalues. For $1 \leq i \leq n$, the attenuation is the quantity by which decays the projection of $\theta - \optEmp$ on $e_i$ at each step. The closer to $0$, the faster will the direction $e_i$ of $\TEmp$ be learnt. In our analysis, the learning rate must satisfy $ \eta_s < 2/(\sigma_1 + \sigma_n) < \eta_b < 2/\sigma_1$.

Specifically, two ranges of learning rate naturally appear. With a small learning rate satisfying $\eta_s < 2/(\sigma_1 + \sigma_n)$, we see on <ref> that the attenuation $\absv{1 - \eta \sigma_i}$ is biggest for the smallest eigenvalue. On the other hand, with a big learning rate satisfying $\eta_b > 2/(\sigma_1 + \sigma_n)$, we see that the attenuation is biggest for the biggest eigenvalue. This motivates the next assumption.

[Learning rate]
The learning rates satisfy

\begin{equation}
0 < \eta_s < \frac{2}{\sigma_1 + \sigma_n} < \eta_b < \frac{2}{\sigma_1}.
\end{equation}

Note that the quantity $\eta = \frac{2}{\sigma_1 + \sigma_n}$ which naturally appears in <ref> is a classical upper bound for proving the convergence of $\sigma_1$-smooth and $\sigma_n$-strongly convex function, of which $\LEmp$ belongs to, see e.g. Thm 2.1.15 in [Nesterov, 2018]. The rate $1/\sigma_1$ is the classical one when we do not have a strong convexity assumption.
This means that in our model, the concept of “small” learning rate simply means being of the order of the best possible learning rates available from an optimization point of view, while the concept of “large” means being close to values leading to diverging algorithms.

In a recent work, Cohen et al., 2021 observed that neural networks trained with gradient descent and “good” constant step size $\eta$ were often in a regime where $\sigma_1$ – the maximum value of the Hessian of the loss – hovered just around $2/\eta$. It is surprisingly analogous to our model: the range of learning rate we consider for big step sizes in Assumption <ref> enforces $2/\eta_b$ to be close to $\sigma_1$.

Geometrically, having $\eta > 1/ \sigma_1$ means that the estimator will oscillate along the direction $e_1$, i.e. $\dotprod{\theta_{t} - \optEmp}{e_1}$ will change sign at each iteration. Such behaviour was observed for neural networks trained with classical learning rate strategies, where the weights' sign change in the early phase of training [Xing et al., 2018].

Then, the following technical assumption is needed to ensure that there is a signal on the lowest and biggest eigendirection. It is satisfied e.g. as soon as the initialization is chosen at random.

We assume

\begin{align*}
\dotprod{\theta_0 - \optEmp}{e_1}_\HH \neq 0, ~~ \dotprod{\theta_0 - \optEmp}{e_n}_\HH \neq 0.
\end{align*}

Finally, we assume that the target accuracy in terms of optimization is not too small compared to the model error $R(\optPop)$:
[Target accuracy and model error]
Consider some learning rates $\eta_b$ and $\eta_s$ chosen in the range of Assumption <ref>. We assume that the target accuracy $\alpha$ satisfies $\alpha \leq \alpha_1$, where $\alpha_1$ is given in Definition <ref> and only depends on the spectrum of $\TEmp$ and the learning rates. Furthermore, we assume that

\begin{equation}\label{eq:asmpt_target_bound}
\frac{\LPop(\optEmp)}{\alpha} \leq \min \cb{\frac{1}{4}, \frac{\kappa_\TEmp}{72\kappa_\TPop}}.
\end{equation}

This assumption is twofold, providing on the target accuracy $\alpha$ both an upper bound (with $\alpha \leq \alpha_1$) and a lower bound (with <ref>). The lower bound not being satisfied amounts to $F$ being a poor approximation of $R$: $R(\theta^\star)$ is too big, or the two ellipsoids' centers are far apart in <ref>. Usual machine learning settings often involve large amounts of data, where the limiting factor for good generalization is poor optimization rather than lacking information, so we can expect this assumption to often hold in practice.
The upper bound stems from the proof technique of our main result in Theorem <ref>, a brief sketch of which is available in <ref>. The proof relies on the fact that sufficiently many steps $t_s$ (resp. $t_b$) are made before the gradient descent is stopped, so that the biggest (resp. the smallest) eigen component is attenuated enough. To ensure this, we can (i) either make the learning rate smaller[Informally, by making $\eta_b \to 2/\sigma_1$ (resp. $\eta_s \to 0$) we need bigger $t_b$ (resp. $t_s$) to achieve optimization error $\alpha$. Thus, $\alpha < \alpha_1$ can be replaced with $\eta_b > 2/\sigma_1 - \epsilon$ (resp. $\eta_s < \epsilon$).], or (ii) take the target error $\alpha$ sufficiently small. We choose the latter, hence the assumption $\alpha \leq \alpha_1$.

§.§ The Main Theorem

We now give our main theorem, whose proof is given in <ref>.

Consider the different quantities defined in Assumptions <ref>, <ref>, <ref> and <ref>.
Then, perform the gradient descent updates of <ref>, with either small step size $\eta_s$ or big step size $\eta_b$, and stop as soon as $F(\theta_{t}) \leq \alpha$, assuming that the resulting estimators satisfy $F(\theta_s) \geq \alpha/2$ and $F(\theta_b) \geq \alpha/2$[This is a mild assumption that could be removed at the price of cumbersome technical details.].

\begin{equation}\label{eq:thm_body_main_bound}
\LPop(\theta_b) - \LPop(\nu^\star) ~\leq~ 34 \frac{\kappa_{\TPop}}{\kappa_\TEmp} (\LPop(\theta_s)-\LPop(\nu^\star)). %~~ \text{with} ~~ c_\alpha = \frac{1 + 2\frac{\sigma_1}{\varsigma_1} \frac{\LPop(\optEmp)}{\alpha}}{\p{1 - \sqrt{18 \frac{\sigma_n}{\varsigma_n} \frac{\LPop(\optEmp)}{\alpha}}}_+}.
\end{equation}

Recall that $\LPop$ has minimum $0$, so <ref> essentially guarantees better performance of $\theta_b$.

The estimator obtained with big step size is better that the one obtained with small step size as soon as the operator $\TEmp$ is ill-conditioned. This is notably the case when doing classification with kernel methods with the ridge estimator, as we discuss in <ref>.

§.§ Discussion

Implications in classification with separable classes. In the context of the example discussed in the introduction, we see clearly that, under the simplifying hypothesis that the population error behaves similarly to the empirical error, the choice of the learning step has the unexpected impact of reducing the Hilbert norm by a multiplicative constant that can be significantly smaller than $1$, leading to an exponential improvement in the classification error.

Comparison with analysis techniques based on learning rate annealing.
Most recent approaches to explain the role of the learning rate in the generalization [Li et al., 2019, Wu et al., 2021] rely on annealing the learning rate: the first phase of the training is carried out with a large step size, before it is discounted to a lower value. We do not need such mechanism in our theoretical analysis, which turns to be simpler with a unique value for the step size. With annealing, our analysis could sum up the following way: do $t$ steps with learning rate greater or equal to $2/\sigma_1$, so that all attenuation coefficient $\absv{1 - \eta_b \sigma_i}$ in <ref> are smaller than $1$ except for the first (few) eigencomponent. Doing so, all eigendirections would be optimized except for the first (few) ones. Then, anneal the learning rate until the $\alpha$ level set of the loss are reached.

Discussion on the complexity.
The fact that the result of Theorem <ref> relies on the condition number can be somewhat surprising. Indeed, we may wonder why the other eigenvalues of the spectrum do not play a role in the result. This is in fact due to the proof technique, which relies on comparing the estimator mostly located on $e_1$ (for big learning rates) and the one mostly on $e_n$ (for the small learning rates). Thus, the distance $\sigma_{n-1} - \sigma_n$ and $\sigma_1 - \sigma_2$ play a role in the complexity of the gradient descent, which is highlighted by the next lemma, which is a consequence of Lemmas <ref> and <ref> in <ref>.

Under the settings of Theorem <ref>, denote $t_s$ (resp. $t_b)$ the number of steps necessary to obtain $\theta_s$ (resp. $\theta_b$). Then,

\begin{equation}
% t_s \geq \frac{1}{2} \frac{\log \csb {\frac{\norm{\theta_0 - \optEmp}_\HH^2}{\iota_n^2} \max \cb{16n\kappa_{\TPopProj}, 4\kappa_\TEmp}}}{\log \frac{1 - \eta_s\sigma_{n}}{1 - \eta_s\sigma_{n-1}}},
% ~~~
% t_b \geq \frac{1}{2}\frac{\log \frac{\norm{\theta_0 - \optEmp}_\HH^2}{\iota_1^2} 4n\kappa_{\TPopProj}}{\log \frac{\absv{\eta_b \sigma_1 - 1}}{\max \cb{\eta_b \sigma_2 - 1, 1 - \eta_b \sigma_n}}}.
t_s \geq O\p{\log \frac{1 - \eta_s\sigma_{n}}{1 - \eta_s\sigma_{n-1}}},
t_b \geq O\p{\log \frac{\absv{\eta_b \sigma_1 - 1}}{\max \cb{\eta_b \sigma_2 - 1, 1 - \eta_b \sigma_n}}}.
\end{equation}

In particular, if $s= \sigma_{n-1} - \sigma_n$, then $t_s = O(1/s)$, and the same holds for $t_b$ with $s =\sigma_1 - \sigma_2$.

We emphasize that the complexity incurred by using learning rates satisfying Assumption <ref> may be very large. However, the motivation of our work is to study an existing practice in deep learning (using large step sizes that are not optimal in terms of optimization of the training loss, but which are better in terms of test loss). We study this phenomenon from a kernel perspective to be able to leverage theoretical tools. Analysing the computational complexity to obtain a given statistical accuracy in convex problems is a separate and well-documented problem, and we do not necessarily advocate the use of very large step sizes for kernel regression.

Nyström projections.
In practice projections techniques (known as Nyström projections) are used to reduce the dimension and avoid a cost quadratic in $n$ for gradient descent. Given some $m \ll n$, this amounts to choosing $m$ anchor points among the data and approximating the kernel matrix $K$ with a rank-$m$ matrix $\tilde{K} = K_{nm}K_m^{-1}K_{nm}^\top$. Fortunately, this is still encompassed in the framework of the model in <ref>. Indeed, the resulting problem is still a quadratic problem, for which a representer theorem holds – the solution lies on the span of the anchor points.

Beyond the square loss.
The result of <ref> relies on using the square loss for obtaining a closed-form expression of the gradient descent update. This is fairly common in learning theory and already provides interesting hindsight. A natural extension to other loss functions would be to consider local quadratic approximations using the Hessian. This is for instance what has been done for kernel ridge regression and self-concordant loss functions by [Bach, 2010].

§.§ Sketch of proof for the main result

The detailed proof is delayed to <ref>. The idea is the following: by tuning the number of steps, we can have the estimator trained with small (resp. big) step size mostly aligned with the smallest (resp. biggest) eigenvector.

The directional bias induced by the step size.
As shown in <ref>, having the learning rates satisfying Assumption <ref> ensures that the quantities

\begin{equation}
\epsilon_b^2 =
\frac{\sum_{2 \leq i \leq n} \dotprod{\theta_b - \optEmp}{e_i}^2}{\dotprod{\theta_b - \optEmp}{e_1}^2},
\epsilon_s^2 =
\frac{\sum_{1 \leq i \leq n-1} \dotprod{\theta_s - \optEmp}{e_i}^2}{\dotprod{\theta_s - \optEmp}{e_n}^2},
\end{equation}

can be made arbitrarily small, while Assumption <ref> ensures they are well defined. $\epsilon_b$ (resp. $\epsilon_s$) quantifies to what extent is $\theta_b - \optEmp$ (resp. $\theta_s - \optEmp$) mostly on $e_1$ (resp. $e_n$). For instance, in the extreme case where $\epsilon_b = 0$ (resp. $\epsilon_s = 0$), then $\theta_b - \optEmp = x e_1, x\in \RR$ (resp. $\theta_s - \optEmp = y e_n, y \in \RR$). To see why it can be made small, refer to the closed-form expression of $\theta_{(\eta, t)}$ in <ref>, and assume for simplification that $\dotprod{\theta_0 - \optEmp}{e_i} = c_0$ for all $i$, with $c_0$ some constant factor. This gives

\begin{equation}
\epsilon_b^2 =
\frac{\sum_{2 \leq i \leq n} (1 - \eta_b \sigma_i)^{2t}}{(1-\eta_b \sigma_1)^{2t}},
\epsilon_s^2 =
\frac{\sum_{1 \leq i \leq n-1} (1 - \eta_s \sigma_i)^{2t}}{(1-\eta_s \sigma_n)^{2t}}.
\end{equation}

Following the discussion of Assumption <ref>, we have that $|1-\eta_b \sigma_1| > |1 - \eta_b \sigma_i|$ for all $i > 1$, and $|1 - \eta_s \sigma_n| > |1 - \eta_s \sigma_n|$ for all $i < n$. Thus, we have that for any $\delta > 0$,

\begin{equation}
t_b \geq \frac{1}{2} \frac{\log 1/\delta^2}{\log \frac{\eta_b \sigma_1 - 1}{\max_{2 \leq i \leq n} |1 - \eta_b \sigma_i|}} \implies \epsilon_b^2 \leq \delta^2,
t_s \geq \frac{1}{2} \frac{\log 1/\delta^2}{\log \frac{1 - \eta_s \sigma_n}{\max_{1 \leq i \leq n-1} |1 - \eta_s \sigma_i|}} \implies \epsilon_s^2 \leq \delta^2.
\end{equation}

Different risk on the $\alpha$-level sets.
In this paragraph, assume (A) $\theta_b = x e_1$, and $\theta_s = y e_n$, that is $\epsilon_b = \epsilon_s = 0$, (B) that $\LPop(\optEmp) = 0$ and (C) that

\begin{equation}\label{eq:sketch_proof_level_sets}
\alpha/2 \leq \LEmp(\theta_b) \leq \alpha, ~~ \alpha/2 \leq \LEmp(\theta_s) \leq \alpha.
\end{equation}

Then, we have for $\theta_b$ that

\begin{equation}\label{eq:sketch_proof_big_lr}
\LPop(\theta_b) = \frac{1}{2} \norm{x e_1}_\TPop^2 \leq \frac{1}{2} \varsigma_1 x^2 \leq \frac{\alpha}{2} \frac{\varsigma_1}{\sigma_1},
\end{equation}

where we used <ref> to bound $\alpha \geq \LEmp(\theta_b) = 1/2 \cdot \sigma_1 x^2$. We do the same with $\theta_s$, this time using $1/2 \cdot \sigma_n y^2 = \LEmp(\theta_s) \leq \alpha$ to obtain

\begin{equation}\label{eq:sketch_proof_small_lr}
\LPop(\theta_s) = \frac{1}{2} \norm{y e_n}_\TPop^2 \geq \frac{1}{2} \varsigma_n y^2 \geq \alpha \frac{\sigma_n}{\varsigma_n},
\end{equation}

Finally, combining <ref> with <ref>, we obtain

\begin{equation*}
\LPop(\theta_b) \leq 2 \frac{\kappa_\TPop}{\kappa_\TEmp} \LPop(\theta_s).
\end{equation*}

Ensuring both conditions can be met together.
We now point out the main differences between this simplified sketch of proof and the rigorous proof in <ref>. First of all, we do not have (A) but rather an approximation of it, with $\epsilon_s \leq \delta$ and $\epsilon_b \leq \delta$. Second, we do not have (B) and rather take into account the error $\LPop(\optEmp)$ to derive Theorem <ref>. Finally, and most importantly, we check that we can have low $\epsilon$ and $F(\theta) \geq \alpha/2$ at the same time. Indeed, we need a big number of iterations to achieve low $\epsilon_b$ or $\epsilon_s$. This implies better optimization of the objective function $\LEmp$. To prevent this, we can either tune the learning rate (having $\eta_s$ close to $0$ and $\eta_b$ close to $2 / \sigma_1$) or provide an upper bound on $\alpha$. We choose the later, hence the hypothesis $\alpha \leq \alpha_1$ in Assumption <ref>.

§ COMPARISON WITH RESULTS IN KERNEL REGRESSION

To provide some intuition over the result of Theorem <ref>, we consider its implications in a supervised learning setting, specifically classification on a low-noise dataset.

§.§ Background on the kernel ridge regression estimator

We consider standard settings. We denote by $\XX \subseteq \RR^d$ the input space, and $\YY = \cb{-1, 1}$ the output space. We draw $n$ i.i.d samples $(x_i, y_i)_{1 \leq i \leq n}$ from an unknown distribution $\rho$ on $\XY$, and we search a prediction function $\theta \in \HH$, where $\HH$ is a RKHS with feature map $\phi: \XX \to \HH$. We assume the kernel to be bounded by a constant $C_K$. See [Aronszajn, 1950] for a precise account on RKHS. We use the square loss as loss function. In order to find a function $\theta$ which maps elements of $\XX$ to $\YY$, we optimize the (regularized) empirical risk $\LEmp$, defined for all $\theta \in \HH$ and $\lambda \geq 0$ a regularization parameter with

\begin{equation}\label{eq:def_empirical_risk}
\LEmp\p{\theta} = \frac{1}{n}\sum_{i=1}^n \frac{1}{2} \p{\theta(x_i) - y_i}^2 + \frac{\lambda}{2} \norm{\theta}^2.
\end{equation}

The minimizer $\optEmp$ of $\LEmp$ is always well defined as soon the training samples $x_i$ are distinct (in the case $\lambda = 0$), which we assume now. We will be minimizing $\LEmp$ with gradient descent when we are in fact interested in minimizing the test error

\begin{equation}\label{eq:def_test_error}
B(\theta) = {\mathbb P}[\sign[\theta(x) \neq y]].
\end{equation}

We will relate the test error and the empirical risk to quadratics forms in $\HH$ by means of other quantities. To do that, we first define the population loss along with its regularized version with

\begin{equation}\label{eq:def_pop_loss}
\forall \theta, ~~ \mathcal{P} (\theta) = \int_{\XY} \frac{1}{2} \p{\theta(x) - y}^2 \dd\rho(x,y), ~~ \mathcal{P}_\lambda(\theta) = \mathcal{P}(\theta) + \frac{\lambda}{2} \norm{\theta}^2.
\end{equation}

The minimizer of $\mathcal{P}$ on $\LL_2(\rho_x)$ is the regression function $\regFunc(x) = \EE[y|x]$. It is an element of $\LL_2(\rho_x)$ but not necessary of $\HH$. We denote by $\optPop$ the minimizer of $\mathcal{P}_\lambda$ on $\HH$. If $\lambda>0$, it is always well defined; otherwise, with $\inclOp : \HH \to \LL_2(\rho_x)$ the inclusion operator, $\optPop$ exists as soon as the projection of the regression function on the closure of the range of $\inclOp$ belongs to the range of $\inclOp$. See [Vito et al., 2005] for a precise account.

In the following, we assume $\lambda \geq 0$ and that $\optPop$ is well defined, and we consider specific assumptions on $\rho$, via assumptions on $\optPop$ and $\regFunc$.

§.§ Relating supervised learning with quadratic forms in H

To relate the problems of <ref> to quadratic forms in the RKHS, we simply need to introduce the empirical covariance operator, with

\begin{equation}\label{eq:def_empirical_covariance}
\TEmp = \frac{1}{n} \sum_{i=1}^n \phi(x_i) \otimes \phi(x_i).
\end{equation}

Then, as optimizing the Hilbert norm is a good proxy for optimizing the test error (following Example <ref> and Lemma <ref> in <ref>), we define

\begin{equation}\label{eq:def_empirical_risk_quadratic}
\LEmp\p{\theta} = \frac{1}{2} \norm{\theta - \optEmp}_{\TEmp}^2 + \minLEmp, ~~ \LPop(\theta) = \frac{1}{2} \norm{\theta - \optPop}_\HH^2, ~~ \text{with} ~~ \minLEmp =  \frac{1}{2n} y^\top \csb{\II_n - \frac{K}{n}\p{\frac{K}{n} + \lambda}^{-1}}  y.
\end{equation}

$K$ is the kernel matrix ($K/n$ shares the same spectrum than $\TEmp$), and the minimum $\LEmp(\optEmp) = \minLEmp$ is $0$ when $\lambda = 0$. We are in the settings of the model of <ref>, and we can readily apply Theorem <ref> with $\TPop$ being the identity operator $\II_\HH$.

Under the settings of Theorem <ref> with the additional assumption that $\TPop = \II_{\HH_n}$, we have

\begin{equation*}
\LPop(\theta_b)- \LPop(\nu^\star) ~\leq~ \frac{34}{\kappa_\TEmp} (\LPop(\theta_s) - \LPop(\nu^\star)).
\end{equation*}

§ EXPERIMENTS

We evaluate the claims of <ref> on CKN-MNIST, a dataset consisting of the MNIST dataset embedded by a convolutional kernel network [Mairal, 2016]. It allows for a realistic use-case, with classification accuracy close to $99\%$, by necessitating a reasonable number of samples $n = 1000$.

On CKN-MNIST, we achieve $98.5\%$ test accuracy with the Gaussian kernel with scale parameter $30$ and no regularization. Adding regularization only improves the test accuracy by $0.04\%$.

(Left) Test accuracy (blue) function of step size $\eta$ for CKN-MNIST. As the learning rate increases, the projection on the first component (orange) increases, which makes the Hilbert norm (green) decreases. This results in predictions closest in $\LL_\infty$ norm to the optimum.

(Right) Test accuracy function of level set $\alpha$ for CKN-MNIST. As we optimize more, the better performances of big step size (orange) compare to optimal step size (blue) vanish to reach the prediction of the optimum of $\LEmp$ (green, dashed).

Shaded areas show standard deviation (train set and initialization) over 10 runs.

[$s=30$, $\kappa_\TEmp=1.3 \cdot 10^4$][t]

[$s=60$, $\kappa_\TEmp=2.8 \cdot 10^5$][t]

[$s=90$, $\kappa_\TEmp=1.4 \cdot 10^6$][t]

[$s=120$, $\kappa_\TEmp=4.0 \cdot 10^7$][t]

Theorem <ref> predicts that the improvement in taking big rather than small step size increases with the condition number of $\TEmp$. We test this claim on our dataset: larger kernel scale $s$ makes the condition number of the kernel matrix $\kappa_\TEmp$ increases, which results in a larger margin between the test accuracy of $\theta_b$ and the one of $\theta_s$.

Test accuracy function of the step size.
We define some final level set $\alpha$. We plot three quantities function of the learning rate $\eta$. The projection on the first component $\dotprod{\theta - \optEmp}{e_1}_\HH$ must be small for moderate learning rate and big for learning rate close to $2/\sigma_1$, as the attenuation satisfies $\absv{1 - \eta \sigma_1} \to 1$ when $\eta \to 2/\sigma_1$, see <ref>. This makes the Hilbert norm $\LPop(\theta(\eta))$ decreases as $\eta$ increases, as predicted by Corollary <ref>. This results in better test accuracy $1-B$ for big learning rate, following Lemma <ref>. This is summed up in <ref>, left. Note that the range for big step size is narrow. This is consistent with Assumption <ref>, which requires a range equal to $2 \sigma_n / (\sigma_1 + \sigma_n) \approx 2\kappa_T^{-1}$. This is also in line with practices in deep learning where the best performance in generalization is often obtained by choosing the largest possible learning rate such that the model does not diverge.

Test accuracy function of optimization.
Our main bound in Theorem <ref> relies on Assumption <ref>: in learning settings, it means that the optimization error $\alpha$ must be greater than a constant times the statistical error $\LPop(\optEmp)$ in order to observe improvements with big step sizes. This is shown in <ref>, right. Additional details on the plot is available in <ref>.

Scale of the kernel.
In <ref>, the scale of the Gaussian kernel is set to $s=30$, with which we obtained the best results on the test set. It is worth noting though that the scale of the kernel directly impacts the conditioning of the matrix. Notably, when $s \to 0$, the kernel matrix $K$ tends to the identity (hence $\kappa_\TEmp \to 1$) while when $s \to \infty$, $K$ tends to a rank-1 operator (hence $\kappa_\TEmp \to \infty$). A core result of Theorem <ref> is that we have bigger improvement for bigger $\kappa_\TEmp$. We reproduce the experiment on the test accuracy with different scale in <ref>, and notice indeed bigger improvements for larger scale $s$, hence worse conditioning.

§ CONCLUSION

A large class of learning problems can be formulated as optimizing a function $\LEmp$ with gradient descent while we are interested in optimizing another function $\LPop$. Using simple quadratic forms to model this mismatch is a natural thing to do, while already providing lots of insight. We indeed show that the choice of large step sizes that may be suboptimal from an optimization point of view may provide better estimators than small/medium step sizes. In particular, we show that this phenomenon occurs in realistic classification tasks with low noise when learning with kernel methods.
In future work, we are planning to study other variants of gradient-based algorithms, which may be stochastic, or accelerated, and perhaps exploit the insight developed in our work to design new algorithms, which would focus on statistical efficiency, exploiting prior knowledge on the test loss, rather than on optimization of the training objective.

A.R. acknowleges support of the French government under management of Agence Nationale de la Recherche as part of the “Investissements d'avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). A.R. acknowledges support of the European Research Council (grant REAL 947908).  J. Mairal was supported by the ERC grant number 714381 (SOLARIS project) and by ANR 3IA MIAI@Grenoble Alpes, (ANR19-P3IA-0003).

[Aronszajn, 1950]
Nachman Aronszajn.
Theory of reproducing kernels.
Transactions of the American Mathematical Society, 680
(3):0 337–404, 1950.

[Audibert, 2004]
J-Y Audibert.
Classification under polynomial entropy and margin assumptions and
randomized estimators.

[Audibert and Tsybakov, 2007]
Jean-Yves Audibert and Alexandre B Tsybakov.
Fast learning rates for plug-in classifiers.
The Annals of statistics, 350 (2):0 608–633,

[Bach, 2010]
F. Bach.
Self-concordant analysis for logistic regression.
Electronic Journal of Statistics, 4:0 384 – 414,
URL <https://doi.org/10.1214/09-EJS521>.

[Bauer et al., 2007]
Frank Bauer, Sergei Pereverzev, and Lorenzo Rosasco.
On regularization algorithms in learning theory.
Journal of Complexity, 230 (1):0 52–72, 2007.

[Blanchard and Mücke, 2016]
Gilles Blanchard and Nicole Mücke.
Optimal rates for regularization of statistical inverse learning
Foundations of Computational Mathematics, 18:0
971–1013, 2016.

[Bottou and Bousquet, 2007]
Léon Bottou and Olivier Bousquet.
The tradeoffs of large scale learning.
Advances in Neural Information Processing Systems (NIPS),

[Chizat et al., 2019]
Lénaïc Chizat, Edouard Oyallon, and Francis Bach.
On lazy training in differentiable programming.
In Advances in Neural Information Processing Systems
(NeurIPS), 2019.

[Cohen et al., 2021]
Jeremy Cohen, Simran Kaur, Yuanzhi Li, J Zico Kolter, and Ameet Talwalkar.
Gradient descent on neural networks typically occurs at the edge of
In International Conference on Learning Representations
(ICLR), 2021.

[Geiping et al., 2021]
Jonas Geiping, Micah Goldblum, Phillip E. Pope, Michael Moeller, and Tom
Stochastic training is not necessary for generalization.
arXiv preprint arXiv:2109.14119, 2021.

[Gerfo et al., 2008]
L. Lo Gerfo, L. Rosasco, F. Odone, E. De Vito, and A. Verri.
Spectral algorithms for supervised learning.
Neural Computation, 200 (7):0 1873–1897,

[Jacot et al., 2018]
Arthur Jacot, Franck Gabriel, and Clement Hongler.
Neural tangent kernel: Convergence and generalization in neural
In Advances in Neural Information Processing Systems
(NeurIPS), 2018.

[Jastrzebski et al., 2020]
Stanislaw Jastrzebski, Maciej Szymczak, Stanislav Fort, Devansh Arpit, Jacek
Tabor, Kyunghyun Cho, and Krzysztof Geras.
The break-even point on optimization trajectories of deep neural
In International Conference on Learning Representations
(ICLR), 2020.

[Jastrzebski et al., 2021]
Stanislaw Jastrzebski, Devansh Arpit, Oliver Astrand, Giancarlo B Kerg, Huan
Wang, Caiming Xiong, Richard Socher, Kyunghyun Cho, and Krzysztof J Geras.
Catastrophic Fisher explosion: Early phase Fisher matrix impacts
In International Conference on Machine Learning (ICML), 2021.

[Li et al., 2019]
Yuanzhi Li, Colin Wei, and Tengyu Ma.
Towards explaining the regularization effect of initial large
learning rate in training neural networks.
In Advances in Neural Information Processing Systems
(NeurIPS), 2019.

[Mairal, 2016]
Julien Mairal.
End-to-End Kernel Learning with Supervised Convolutional Kernel
In Advances in Neural Information Processing Systems (NIPS),

[Nakkiran, 2020]
Preetum Nakkiran.
Learning rate annealing can provably help generalization, even for
convex problems.
arXiv preprint arXiv:2005.07360, 2020.

[Nemirovski et al., 2009]
Arkadi Nemirovski, Anatoli Juditsky, Guanghui Lan, and Alexander Shapiro.
Robust stochastic approximation approach to stochastic programming.
SIAM Journal on optimization, 190 (4):0
1574–1609, 2009.

[Nesterov, 2018]
Yurii Nesterov.
Lectures on Convex Optimization.
Springer, 2018.

[Pillaud-Vivien et al., 2018]
Loucas Pillaud-Vivien, Alessandro Rudi, and Francis Bach.
Statistical Optimality of Stochastic Gradient Descent on Hard
Learning Problems through Multiple Passes.
In Neural Information Processing Systems (NeurIPS),

[Pillaud-Vivien et al., 2018]
Loucas Pillaud-Vivien, Alessandro Rudi, and Francis Bach.
Exponential convergence of testing error for stochastic gradient
In Conference On Learning Theory (COLT), 2018b.

[Rawat and Wang, 2017]
Waseem Rawat and Zenghui Wang.
Deep convolutional neural networks for image classification: A
comprehensive review.
Neural computation, 290 (9):0 2352–2449,

[Smith et al., 2021]
Samuel L Smith, Benoit Dherin, David Barrett, and Soham De.
On the origin of implicit regularization in stochastic gradient
In International Conference on Learning Representations
(ICLR), 2021.

[Vito et al., 2005]
Ernesto De Vito, Lorenzo Rosasco, Andrea Caponnetto, Umberto De Giovannini, and
Francesca Odone.
Learning from examples as an inverse problem.
Journal of Machine Learning Research, 60
(30):0 883–904, 2005.

[Wu et al., 2021]
Jingfeng Wu, Difan Zou, Vladimir Braverman, and Quanquan Gu.
Direction matters: On the implicit bias of stochastic gradient
descent with moderate learning rate.
In International Conference on Learning Representations
(ICLR), 2021.

[Xing et al., 2018]
Chen Xing, Devansh Arpit, Christos Tsirigotis, and Yoshua Bengio.
A walk with SGD.
arXiv preprint arXiv:1802.08770, 2018.

[Yao et al., 2007]
Yuan Yao, Lorenzo Rosasco, and Andrea Caponnetto.
On early stopping in gradient descent learning.
Constructive Approximation, 260 (2):0
289–315, 2007.


#  Spin Hall response at finite wave vector in ferromagnets

Gen Tatara RIKEN Center for Emergent Matter Science (CEMS) and RIKEN Cluster
for Pioneering Research (CPR)RIKEN Center for Emergent Matter Science (CEMS)
and RIKEN Cluster for Pioneering Research (CPR) 2-1 Hirosawa 2-1 Hirosawa Wako
Wako Saitama Saitama 351-0198 Japan 351-0198 Japan

###### Abstract

Spin Hall effect at finite wave vector in a ferromagnetic conductor is
theoretically studied by calculating the spin density as the linear response
to an applied electric field. The cases of a spin-orbit interaction due to
random impurities and a localized Rashba interaction are considered. It is
shown that the spin Hall effect has a broad response for the wave vector
$q\lesssim 2{k_{F}}$ where ${k_{F}}$ is the Fermi wave vector. This fact
confirms the local nature of the spin-charge conversion effects.

##  1 Introduction

Spin Hall effect is widely used as a method to generate spin current by
applying an electric field to nonmagnetic conductors [1, 2]. It induces a spin
accumulation at edges and interfaces and a torque in the case of a junction
with a ferromagnet [3]. Conventionally, spin Hall effect has been
theoretically argued in terms of spin Hall conductivity which is the ratio of
the spin current and the applied electric field [2]. Although this description
is convenient as an analogy with electric conductivity is applicable, it has a
fundamental ambiguity of definition of spin current because of the
nonconservation of spin current.

The ambiguity can be avoided by discussing spin density, which is a physical
observable, as in the original proposal of the spin Hall effect by Dyakonov
and Perel [1]. The spin Hall effect in the spin representation in the linear
response regime is described by a response function of spin to an applied
electric field, namely, a correlation function of spin and electric current
[4]. A mechanism to couple the spin and the electric current, such as the
spin-orbit interaction and a magnetic structure, is essential to have a finite
spin-electric current correlation function. Considering a spin-orbit
interaction keeping the inversion symmetry, the spin-electric current
correlation has vanishing uniform component, as it is the average of
$s_{i}k_{j}$, a product of spin ${{\bm{s}}}$ and wave vector ${\bm{k}}$ ($i$
and $j$ represent directions), which is odd in inversion. This is in contrast
to the spin Hall conductivity, which is inversion-even and has a uniform
component. The result is natural physically as the spin accumulation generated
by the spin Hall effect is not uniform, but is localized at the edges within
the penetration length, determined by the rates of the spin relaxation and
diffusion [1]. It was recently pointed out theoretically that the leading
contribution of the spin-electric current correlation is the term linear in
the external wave vector ${\bm{q}}$ (or spatial derivative), which corresponds
to the inhomogeneous spin accumulation at edges [4]. When electron diffusion
is neglected, the induced spin density ${{\bm{s}}}$ was shown to be
proportional to the vorticity of the electric current or electric field
${\bm{E}}$ as

$\displaystyle{{\bm{s}}}=\lambda_{\rm sh}(\nabla\times{\bm{E}}),$ (1)

where $\lambda_{\rm sh}$ is a coefficient proportional to the strength of the
spin orbit interaction and to the spin Hall angle of the conventional
description. As was argued in Ref. [5], the spin Hall effect is a direct
consequence of the spin-vorticity coupling [6].

The expression (1) clearly indicates that the spin Hall effect in terms of
spin accumulation has vanishing bulk impact in the long wavelength limit,
${\bm{q}}=0$. At large ${\bm{q}}$ comparable to the Fermi wave vector $k_{\rm
F}$, the correlation function is also small in the nonmagnetic case because
the electron excitation is away from the Fermi energy in the spin-unpolarized
case. The situation is different if electrons have different Fermi surfaces,
like in ferromagnetic metal. There, the Fermi surfaces for the electrons with
spin up and down do not overlap and low energy excitation becomes possible at
high ${\bm{q}}$, in the broad regime of $|k_{{\rm F}+}-k_{{\rm F}-}|<q<k_{{\rm
F}+}+k_{{\rm F}-}$ for the case of spherical Fermi surface with Fermi wave
vectors $k_{{\rm F}\pm}$ for spin $\pm$. Such a low energy excitation with
spin flip, known as the Stoner excitation, leads to a large weight for the
magnetic susceptibility at large $q$. The low energy excitation is expected to
lead to an enhancement of the spin Hall response function in the Stoner regime
in ferromagnets.

Spin response function for the applied electric current was discussed in the
context of spin-orbit torque [3] phenomenologically in Ref. [7] and density-
functional calculation has been carried out for bilayers of a ferromagnet and
a heavy metal [8, 9, 10]. The torque is induced by the spin accumulation at
the interface generated by the spin Hall effect in the heavy metal, and thus
the theories calculate a spin response function for the applied electric
current. In the present paper, we explore analytically the finite-$q$ regime
of the response function in a ferromagnetic metal with a quadratic free
electron dispersion.

##  2 Linear response theory : Random impurity spin-orbit model

Figure 1:  Feynman diagrams representing the response function of spin (left
vertex) and electric current (right vertex). Solid lines are electron Green’s
functions ($\sigma=\pm$ represents spin) and dotted lines denote impurity
effects ($v_{\rm i}$ is the impurity potential) and the spin-orbit interaction
is denoted by $\lambda$. Diagram (b) is the dominant contributions at $q\sim
0$ as has been known in the context of the anomalous Hall [11] and spin Hall
effects [4], while diagram (a) is dominant for large $q$.

The spin Hall response function is a correlation function of spin and electric
current, consisting of is a spin vertex and a current vertex connected by the
spin-orbit interaction, depicted in Fig. 1 [4]. We consider ferromagnetic
conductor with spin splitting $M$ with spin polarization along a unit vector
$\hat{{\bm{M}}}$. The electron energy for wave vector ${\bm{k}}$ and spin
$\sigma=\pm$ measured from the Fermi energy ${\epsilon_{F}}$ is
$\epsilon_{{\bm{k}}\sigma}=\frac{k^{2}}{2m}-\sigma M-{\epsilon_{F}}$, $m$
being the mass. The spin-orbit interaction we consider is the one induced by
impurities, represented by

$\displaystyle H_{\rm so}$ $\displaystyle=\lambda\int
d^{3}rc^{\dagger}({\bm{r}})[(\nabla
v({\bm{r}})\times\hat{{\bm{p}}})\cdot{\bm{\sigma}}]c({\bm{r}})$
$\displaystyle=i\lambda\sum_{{\bm{k}}{\bm{k}}^{\prime}}v_{{\bm{k}}^{\prime}-{\bm{k}}}({\bm{k}}^{\prime}\times{\bm{k}})\cdot
c^{\dagger}_{{\bm{k}}^{\prime}}{\bm{\sigma}}c_{{\bm{k}}}$ (2)

where $\lambda$ is a coupling constant, $v({\bm{r}})=v_{\rm
i}\sum_{{\bm{R}}_{n}}\delta({\bm{r}}-{\bm{R}}_{n})$ is an impurity potential,
where $v_{\rm i}$ and ${\bm{R}}_{n}$ are the strength of the impurity
potential and the position of $n$-th impurity, respectively, and the impurity
scattering potential is $H_{\rm i}=\int
d^{3}rv({\bm{r}})c^{\dagger}({\bm{r}})c({\bm{r}})$. The impurity scattering
induces an electron lifetime of elastic scattering, $\tau$, given by
$\tau^{-1}=2\pi\nu n_{\rm i}{v_{\rm i}}^{2}$, where $n_{\rm i}$ is the
impurity concentration and $\nu$ is the density of states of electron at the
Fermi energy.

We consider a spin density ${{\bm{s}}}$ induced by an applied electric field
${\bm{E}}$ in the linear response regime. The field is temporary static and
spatially nonuniform as represented by a finite external wave vector
${\bm{q}}$. We first revisit the nonmagnetic case argued in Ref. [4], where
the spin density response, $s^{\alpha}\equiv s_{i}^{{\rm sh},\alpha}E_{i}$
($\alpha$ is spin direction and $i$ is spatial direction), was studied at
$q\sim 0$ for spin-unpolarized case. Near $q=0$, the spin response function
was shown to be described by the process Fig. 1(b), which reads

$\displaystyle s_{i}^{{\rm sh},\alpha}({\bm{q}})$
$\displaystyle=\frac{2i}{\pi}\frac{\lambda n_{\rm i}{v_{\rm
i}}^{3}}{m}\nu\epsilon_{\alpha jk}{\rm
Im}[J_{k}({\bm{q}})\Sigma_{ij}({\bm{q}})]$ (3)

where

$\displaystyle J_{i}({\bm{q}})$
$\displaystyle=\frac{1}{V}\sum_{{\bm{k}}}k_{i}g_{{\bm{k}}+\frac{{\bm{q}}}{2}}^{\rm
r}g_{{\bm{k}}-\frac{{\bm{q}}}{2}}^{\rm a}$ (4)

is the spin part (the factor $k_{i}$ arises from the spin-orbit interaction)
and

$\displaystyle\Sigma_{ij}({\bm{q}})\equiv\frac{1}{V}\sum_{{\bm{k}}}k_{i}k_{j}g_{{\bm{k}}+\frac{{\bm{q}}}{2}}^{\rm
r}g_{{\bm{k}}-\frac{{\bm{q}}}{2}}^{\rm a}$ (5)

is a current-current response function ($V$ is a dimensionless volume). The
retarded and advanced components of the electron free Green’s functions at
zero angular frequency are denoted by $g_{{\bm{k}}}^{\rm
r}\equiv[-\epsilon_{{\bm{k}}}+\frac{i}{2\tau}]^{-1}$ and $g_{{\bm{k}}}^{\rm
a}$, respectively. Near $q\sim 0$ [4], we have
$J_{i}({\bm{q}})=-i\frac{\pi}{3m}\nu{k_{F}}^{2}\tau^{2}q_{i}+O(q^{3})$ and
$\Sigma_{ij}=\delta_{ij}\frac{a^{3}m^{2}}{e^{2}}\sigma_{{\rm e}}+O(q^{2})$,
where ${k_{F}}$ is the Fermi wave vector, $\sigma_{{\rm e}}=e^{2}n\tau/m$ is
the Boltzmann conductivity, $n$ is the electron density, $a\sim{k_{F}}^{-1}$
is the lattice constant. Thus we obtain ${{\bm{s}}}=\lambda_{\rm
sh}(\nabla\times{\bm{E}})$ with $\lambda_{\rm
sh}=\frac{1}{9}\frac{\epsilon_{\rm
so}}{{\epsilon_{F}}}({\epsilon_{F}}\tau)^{2}\frac{1}{{\epsilon_{F}}{k_{F}}^{2}}$
in the nonmagnetic case (assuming $\nu\sim{\epsilon_{F}}^{-1}$), where
$\epsilon_{\rm so}\equiv\lambda{v_{\rm i}}{k_{F}}^{2}$ and ${\epsilon_{F}}$
are the spin-orbit energy and the Fermi energy, respectively.

Now we study the ferromagnetic case. We evaluate the electron spin density,

$\displaystyle s^{\alpha}({\bm{r}},t)$
$\displaystyle=-i\tr[\sigma_{\alpha}{\cal G}^{<}({\bm{r}},t,{\bm{r}},t)],$ (6)

where $\tr$ is the trace over spin, $\sigma_{\alpha}$ is the Pauli matrix with
direction $\alpha$, ${\cal G}$ is the full electron Green’s function and <
denotes the lesser component. To explore spin-split case and general
${\bm{q}}$, we go beyond the Fermi surface approximation employed in Ref. [4],
where the excitation represented by the product of the retarded and advanced
Green’s functions were considered. In fact, the expression Eq. (5), which is
commonly used for conductivity, diverges at large $k$ in the continuum model,
but it gives a good approximation concerning Fermi surface effects near $q=0$.
The correct expression applicable to general cases is obtained
straightforwardly by evaluating the lesser component in Eq. (6). The
contribution linear in a static applied electric field is written in the
momentum representation as

$\displaystyle s^{\alpha}({\bm{q}})$
$\displaystyle=i\frac{e}{V}\sum_{{\bm{k}}}\int\\!\frac{d\omega}{2\pi}\sum_{k}k_{k}\lim_{\Omega\rightarrow
0}\tr[\sigma_{\alpha}G_{{\bm{k}}+\frac{{\bm{q}}}{2},\omega+\Omega}G_{{\bm{k}}-\frac{{\bm{q}}}{2},\omega}]^{<}A_{k}({\bm{q}},\Omega)$
(7)

where ${\bm{A}}$ is a vector potential related to the electric field as
${\bm{E}}=-\dot{{\bm{A}}}$ ($i\Omega A_{k}({\bm{q}},\Omega)=E_{k}({\bm{q}})$)
and the Green’s function $G$ is without the external field. The linear
response coefficient defined as $s^{\alpha}=s_{k}^{\alpha}E_{k}$ is therefore

$\displaystyle s_{k}^{\alpha}({\bm{q}})$
$\displaystyle=i\frac{e}{V}\sum_{{\bm{k}}}\int\\!\frac{d\omega}{2\pi}\sum_{k}k_{k}\lim_{\Omega\rightarrow
0}\frac{1}{\Omega}\tr[\sigma_{\alpha}G_{{\bm{k}}+\frac{{\bm{q}}}{2},\omega+\Omega}G_{{\bm{k}}-\frac{{\bm{q}}}{2},\omega}]^{<}$
(8)

The static limit in the case of static interactions is given by the
contribution of

$\displaystyle[G_{{\bm{k}}+\frac{{\bm{q}}}{2},\omega+\Omega}G_{{\bm{k}}-\frac{{\bm{q}}}{2},\omega}]^{<}=\Omega
f^{\prime}(\omega)\left[G_{{\bm{k}}+\frac{{\bm{q}}}{2},\omega}^{\rm
r}G_{{\bm{k}}-\frac{{\bm{q}}}{2},\omega}^{\rm
a}-\frac{1}{2}(G_{{\bm{k}}+\frac{{\bm{q}}}{2},\omega}^{\rm
r}G_{{\bm{k}}-\frac{{\bm{q}}}{2},\omega}^{\rm
r}+G_{{\bm{k}}+\frac{{\bm{q}}}{2},\omega}^{\rm
a}G_{{\bm{k}}-\frac{{\bm{q}}}{2},\omega}^{\rm a})\right]+O(\Omega^{2})$ (9)

where $f(\omega)\equiv(e^{\beta\omega}+1)^{-1}$ is the Fermi distribution
function ($\beta\equiv({k_{B}}T)^{-1}$ is the inverse temperature, ${k_{B}}$
being the Boltzmann constant). The Green’s function $G$ is calculated
perturbatively taking account of the spin-orbit interaction. The impurity
scattering is included to compensate the momentum change due to the spin-orbit
interaction. In the case of uniform ($q=0$) response function, the diagram of
Fig. 1(b) including the impurity to the third order is dominant [11, 4], while
the second-order diagram is also important at finite $q$. The second-order
contribution is

$\displaystyle s_{k}^{(2)\alpha}({\bm{q}})$
$\displaystyle=-i\frac{e}{\pi}\frac{\lambda n_{\rm i}{v_{\rm
i}}^{2}}{m}\frac{1}{V^{2}}\sum_{{\bm{k}}{\bm{k}}^{\prime}}k^{\prime}_{k}\left[\left({\bm{k}}^{\prime}+\frac{{\bm{q}}}{2}\right)\times\left({\bm{k}}+\frac{{\bm{q}}}{2}\right)\right]_{\beta}$
$\displaystyle\times{\rm
Re}\tr\biggl{[}\sigma_{\alpha}g_{{\bm{k}}+\frac{{\bm{q}}}{2}}^{\rm
r}\sigma_{\beta}g_{{\bm{k}}^{\prime}+\frac{{\bm{q}}}{2}}^{\rm
r}g_{{\bm{k}}^{\prime}-\frac{{\bm{q}}}{2}}^{\rm
a}g_{{\bm{k}}-\frac{{\bm{q}}}{2}}^{\rm
a}-\sigma_{\alpha}g_{{\bm{k}}+\frac{{\bm{q}}}{2}}^{\rm
r}\sigma_{\beta}g_{{\bm{k}}^{\prime}+\frac{{\bm{q}}}{2}}^{\rm
r}g_{{\bm{k}}^{\prime}-\frac{{\bm{q}}}{2}}^{\rm
r}g_{{\bm{k}}-\frac{{\bm{q}}}{2}}^{\rm r}\biggr{]}$ (10)

where $g_{{\bm{k}}\sigma}^{\rm
r}\equiv[-\epsilon_{{\bm{k}}\sigma}+\frac{i}{2\tau}]^{-1}$ and
$g_{{\bm{k}}\sigma}^{\rm a}$ are the Free Green’s functions at zero-frequency.

The trace for spin is calculated using

$\displaystyle\tr[\sigma^{\alpha}A\sigma^{\beta}B]$
$\displaystyle=(\delta_{\alpha\beta}-\delta_{\alpha z}\delta_{\beta
z})\sum_{\sigma}A_{\sigma}B_{-\sigma}+\delta_{\alpha z}\delta_{\beta
z}\sum_{\sigma}A_{\sigma}B_{\sigma}-i\epsilon_{\alpha\beta
z}\sum_{\sigma}\sigma A_{\sigma}B_{-\sigma}$ (11)

for $2\times 2$ diagonal matrices, $A\equiv\left(\begin{array}[]{cc}A_{+}&0\\\
0&A_{-}\end{array}\right)$ and $B\equiv\left(\begin{array}[]{cc}B_{+}&0\\\
0&B_{-}\end{array}\right)$. As was argued in Ref. [5], we have three
components of spin accumulation in ferromagnets, $s_{k}^{\alpha}=s_{k}^{{\rm
sh},\alpha}+s_{k}^{{\perp},\alpha}+s_{k}^{{\rm ad},\alpha}$, where
$s_{k}^{{\rm sh},\alpha}$ is the spin Hall contribution proportional to
${\bm{q}}\times{\bm{E}}$ (vorticity in the ${\bm{q}}$-representation),
$s_{k}^{{\perp},\alpha}$ is the perpendicular component proportional to
${\bm{M}}\times({\bm{q}}\times{\bm{E}})$ and $s_{k}^{{\rm ad},\alpha}$ is the
adiabatic component parallel to ${\bm{M}}$. (Note that
${\bm{q}}\times{\bm{E}}$ at finite $q$ regime with extra dependence on
${\bm{q}}$ does not exactly corresponds to the vorticity in the real space,
$\nabla\times{\bm{E}}$.) We focus on the nonadiabatic components $s_{k}^{{\rm
sh},\alpha}$ and $s_{k}^{{\perp},\alpha}$. The contribution second order of
impurity, Eq. (10), is

$\displaystyle s_{k}^{{\rm sh}(2)\alpha}({\bm{q}})$
$\displaystyle=-i\frac{e}{\pi}\frac{\lambda n_{\rm i}{v_{\rm
i}}^{2}}{m}\epsilon_{\alpha ij}\sum_{\sigma}{\rm
Re}\left[J_{j}^{{\rm(ra)},\sigma,-\sigma}({\bm{q}})\Sigma_{ik}^{{\rm(ra)},-\sigma,-\sigma}-J_{j}^{{\rm(rr)},\sigma,-\sigma}({\bm{q}})\Sigma_{ik}^{{\rm(rr)},-\sigma,-\sigma}\right]$
$\displaystyle s_{k}^{{\perp}(2)\alpha}({\bm{q}})$
$\displaystyle=-i\frac{e}{\pi}\frac{\lambda n_{\rm i}{v_{\rm
i}}^{2}}{m}\epsilon_{\alpha\beta z}\epsilon_{\beta ij}\sum_{\sigma}\sigma{\rm
Im}\left[J_{j}^{{\rm(ra)},\sigma,-\sigma}({\bm{q}})\Sigma_{ik}^{{\rm(ra)},-\sigma,-\sigma}-J_{j}^{{\rm(rr)},\sigma,-\sigma}({\bm{q}})\Sigma_{ik}^{{\rm(rr)},-\sigma,-\sigma}\right]$
(12)

Here (A, B=r, a)

$\displaystyle J_{j}^{{\rm(AB)},\sigma,-\sigma}({\bm{q}})$
$\displaystyle\equiv\frac{1}{V}\sum_{{\bm{k}}}\left({\bm{k}}+\frac{{\bm{q}}}{2}\right)_{j}g_{{\bm{k}}+\frac{{\bm{q}}}{2},\sigma}^{\rm
A}g_{{\bm{k}}-\frac{{\bm{q}}}{2},-\sigma}^{\rm B}$
$\displaystyle\Sigma_{ik}^{{\rm(AB)},\sigma,\sigma}$
$\displaystyle\equiv\frac{1}{V}\sum_{{\bm{k}}^{\prime}}k^{\prime}_{k}\left({\bm{k}}^{\prime}+\frac{{\bm{q}}}{2}\right)_{i}g_{{\bm{k}}^{\prime}+\frac{{\bm{q}}}{2},\sigma}^{\rm
A}g_{{\bm{k}}^{\prime}-\frac{{\bm{q}}}{2},\sigma}^{\rm B},$ (13)

are extensions of Eqs. (4) (5) to spin-polarized case with finite $q$. As is
easily seen from symmetry,

$\displaystyle J_{j}^{{\rm(AB)},\sigma,-\sigma}$
$\displaystyle=\hat{q}_{j}\overline{J}^{{\rm(AB)},\sigma,-\sigma}$
$\displaystyle\Sigma_{ik}^{{\rm(AB)},\sigma,\sigma}$
$\displaystyle=\delta_{ik}\overline{\Sigma}^{{\rm(AB)},\sigma,\sigma}+\hat{q}_{i}\hat{q}_{k}\overline{\Sigma}^{{\rm(AB)}(2),\sigma,\sigma}$
(14)

where

$\displaystyle\overline{J}^{{\rm ra},\sigma,-\sigma}$
$\displaystyle\equiv\frac{1}{V}\sum_{{\bm{k}}}\left(k_{z}+\frac{q}{2}\right)g_{{\bm{k}}+\frac{{\bm{q}}}{2},\sigma}^{\rm
r}g_{{\bm{k}}-\frac{{\bm{q}}}{2},-\sigma}^{\rm a}$
$\displaystyle\overline{\Sigma}^{{\rm ra},\sigma,\sigma}$
$\displaystyle\equiv\frac{1}{V}\sum_{{\bm{k}}^{\prime}}\frac{(k^{\prime})^{2}}{3}g_{{\bm{k}}^{\prime}+\frac{{\bm{q}}}{2},\sigma}^{\rm
r}g_{{\bm{k}}^{\prime}-\frac{{\bm{q}}}{2},\sigma}^{\rm a},$
$\displaystyle\overline{\Sigma}^{{\rm ra}(2),\sigma,\sigma}$
$\displaystyle\equiv\frac{1}{V}\sum_{{\bm{k}}^{\prime}}\left[k^{\prime}_{z}\left(k^{\prime}_{z}+\frac{q}{2}\right)-\frac{(k^{\prime})^{2}}{3}\right]g_{{\bm{k}}^{\prime}+\frac{{\bm{q}}}{2},\sigma}^{\rm
r}g_{{\bm{k}}^{\prime}-\frac{{\bm{q}}}{2},\sigma}^{\rm a},$ (15)

and similarly for r and a components, where the $z$ axis of ${\bm{k}}$ here
are chosen along $\hat{{\bm{q}}}$. Eq. (12) therefore read

$\displaystyle s_{k}^{{\rm sh}(2)\alpha}({\bm{q}})$
$\displaystyle=-i\overline{\lambda}_{\rm sh}\epsilon_{\alpha
kj}\hat{q}_{j}\overline{s}^{{\rm sh}(2)}(q)$ $\displaystyle
s_{k}^{{\perp}(2)\alpha}({\bm{q}})$ $\displaystyle=-i\overline{\lambda}_{\rm
sh}\epsilon_{\alpha\beta z}\epsilon_{\beta
kj}\hat{q}_{j}\overline{s}^{{\perp}(2)}(q)$ (16)

where

$\displaystyle\overline{\lambda}_{\rm sh}\equiv$
$\displaystyle\frac{e}{\pi}\frac{\lambda n_{\rm i}{v_{\rm
i}}^{2}}{m}\frac{{k_{F}}^{3}}{{\epsilon_{F}}^{4}}=\frac{e}{\pi}\frac{\epsilon_{\rm
so}}{{\epsilon_{F}}}\frac{n_{\rm i}v_{\rm
i}}{{\epsilon_{F}}}\frac{1}{{k_{F}}{\epsilon_{F}}}$ (17)

and

$\displaystyle\overline{s}^{{\rm sh}(2)}(q)\equiv$
$\displaystyle\frac{{\epsilon_{F}}^{4}}{{k_{F}}^{3}}\sum_{\sigma}{\rm
Re}\left[\overline{J}^{{\rm(ra)}\sigma,-\sigma}\overline{\Sigma}^{{\rm(ra)}-\sigma,-\sigma}-\overline{J}^{{\rm(rr)}\sigma,-\sigma}\overline{\Sigma}^{{\rm(rr)}-\sigma,-\sigma}\right]$
$\displaystyle\overline{s}^{{\perp}(2)}(q)\equiv$
$\displaystyle\frac{{\epsilon_{F}}^{4}}{{k_{F}}^{3}}\sum_{\sigma}\sigma{\rm
Im}\left[\overline{J}^{{\rm(ra)}\sigma,-\sigma}\overline{\Sigma}^{{\rm(ra)}-\sigma,-\sigma}-\overline{J}^{{\rm(rr)}\sigma,-\sigma}\overline{\Sigma}^{{\rm(rr)}-\sigma,-\sigma}\right]$
(18)

are dimensionless coefficients. Equation (16) indicates that directions of
spin polarizations are ${{\bm{s}}}^{{\rm
sh}(2)}({\bm{q}})\propto({\bm{q}}\times{\bm{E}})$ and
${{\bm{s}}}^{{\perp}(2)}({\bm{q}})\propto\hat{{\bm{M}}}\times({\bm{q}}\times{\bm{E}})$.

A key factor determining the spin Hall response is the spin part function
$\overline{J}^{{\rm AB},\sigma,-\sigma}$, that involves the two opposite
spins. This function has a large contribution from ${\bm{k}}$ satisfying
$\epsilon_{{\bm{k}},\sigma}=0$ (redefining ${\bm{k}}+\frac{{\bm{q}}}{2}$ as
${\bm{k}}$ in Eq. (15)) and $\epsilon_{{\bm{k}}+{\bm{q}},-\sigma}=0$, i.e., at
$k=k_{{\rm F}\sigma}$ and $|{\bm{k}}+{\bm{q}}|=k_{{\rm F},-\sigma}$. These
conditions are satisfied when $|k_{{\rm F}+}-k_{{\rm F}-}|\leq q\leq k_{{\rm
F}+}+k_{{\rm F}-}$, namely, in the Stoner regime for the external wave vector
$q$. In the regime, the amplitude of the function is the order of
${k_{F}}/{\epsilon_{F}}^{2}$ regardless of the value of
$\eta\equiv(2\tau)^{-1}$. Considering $q=O({k_{F}})$ in the Stoner regime, the
spin density has therefore an enhancement factor of the order of
${k_{F}}\ell_{s}$ compared to the conventional spin Hall effect near $q=0$
(Eq. (1)), where $\ell_{s}$ is the length of spin propagation (spin diffusion
length), related to $q$ as $q(=-i\nabla)\sim 1/\ell_{s}$ (Fig. 2). In
contrast, function $\overline{\Sigma}^{{\rm AB},\sigma,\sigma}$ is not
sensitive to spin-splitting. The real part of $\overline{\Sigma}^{{\rm
ra},\sigma,\sigma}$ has a sharp peak at $q=0$ of magnitude and width
proportional to $\tau$ and $\tau^{-1}$, respectively, while the imaginary part
is small.

The contribution third order of the impurity, Fig. 1(b), is similarly obtained
as

$\displaystyle s_{k}^{{\rm sh}(3)\alpha}({\bm{q}})$
$\displaystyle=-i\overline{\lambda}_{\rm sh}\pi\frac{{v_{\rm
i}}}{{\epsilon_{F}}}\epsilon_{\alpha ij}\hat{q}_{j}\overline{s}^{{\rm
sh}(3)}(q)$ $\displaystyle s_{k}^{{\perp}(3)\alpha}({\bm{q}})$
$\displaystyle=-i\overline{\lambda}_{\rm sh}\pi\frac{{v_{\rm
i}}}{{\epsilon_{F}}}\epsilon_{\alpha\beta z}\epsilon_{\beta
ij}\hat{q}_{j}\overline{s}^{{\perp}(3)}(q)$ (19)

where

$\displaystyle\overline{s}^{{\rm sh}(3)}(q)\equiv$
$\displaystyle-\frac{{\epsilon_{F}}^{5}}{{k_{F}}^{3}}\sum_{\sigma}\nu_{-\sigma}{\rm
Im}\left[\overline{J}^{{\rm(ra)}\sigma,-\sigma}\overline{\Sigma}^{{\rm(ra)}-\sigma,-\sigma}+\overline{J}^{{\rm(rr)}\sigma,-\sigma}\overline{\Sigma}^{{\rm(rr)}-\sigma,-\sigma}\right]$
$\displaystyle\overline{s}^{{\perp}(3)}(q)\equiv$
$\displaystyle\frac{{\epsilon_{F}}^{5}}{{k_{F}}^{3}}\sum_{\sigma}\sigma\nu_{-\sigma}{\rm
Re}\left[\overline{J}^{{\rm(ra)}\sigma,-\sigma}\overline{\Sigma}^{{\rm(ra)}-\sigma,-\sigma}+\overline{J}^{{\rm(rr)}\sigma,-\sigma}\overline{\Sigma}^{{\rm(rr)}-\sigma,-\sigma}\right]$
(20)

Compared to the second-order process, the third-order one includes one more
self-energy, which is dominantly pure imaginary ($\pm i\pi\nu$), and thus the
roles of the real and imaginary parts of the response function are
interchanged [11]. At $q=0$, the current-current correlation function
$\Sigma_{ik}^{{\rm(ra)},\sigma,\sigma}$ has large real part with vanishing
imaginary part, while the current function $J_{j}^{{\rm(ra)},\sigma,-\sigma}$
is pure imaginary proportional to $q_{j}$. We have therefore large imaginary
part of the product of the two functions, resulting in a large third-order
contribution, $\overline{s}^{{\rm sh}(3)}$, at $q\sim 0$ (Fig. 2). (The same
mechanism works for the anomalous Hall effect [11].) The feature is different
at finite $q$ in the ferromagnetic case, and the second-order diagram becomes
dominant.

The plot of the amplitudes $\overline{s}^{{\rm sh}(n)}({\bm{q}})$ and
$\overline{s}^{{\perp}(n)}({\bm{q}})$ ($n=2,3$) are shown in Fig. 2. The
perpendicular components $\overline{s}^{{\perp}(n)}({\bm{q}})$ are small
compared to $\overline{s}^{{\rm sh}(n)}({\bm{q}})$. It is seen that the
second-order spin Hall contribution $\overline{s}^{{\rm sh}(2)}$ is dominant
for a broad regime of $q\lesssim 2{k_{F}}$ with a height insensitive to
$\eta$. The spin Hall effect has thus a broad response up to a high-wave
vector regime in ferromagnets. This feature indicates that the induced spin
accumulation has a localized nature, namely, a local nature of the spin-charge
conversion. Conventional spin Hall effect in nonmagnets corresponds to
$\overline{s}^{{\rm sh}(3)}$ near $q=0$ with a high peak $\propto\tau$.

Figure 2: The amplitude of spin Hall response functions $\overline{s}^{{\rm
sh}(n)}$ and $\overline{s}^{{\perp}(n)}$ for the second and third-order of
impurity ($n=2,3$ respectively) plotted as functions of wave vector $q$
normalized by ${k_{F}}$. Spin splitting is $M=0,0.4,0.8$ in unit of
${\epsilon_{F}}$ and the damping constant is $\eta/{\epsilon_{F}}=0.1$. Shaded
regime is the contribution for the uniform spin Hall effect for the case of
spin diffusion length of $\ell_{s}=10/{k_{F}}$.

###  2.1 Discussion

In Ref. [10], the spin magnetoelectric susceptibility $\chi_{s}$ in a junction
was defined as $s=\chi_{s}E$ and its $q=0$ contribution was theoretically
discussed. Their spin density proportional to ${\bm{E}}\times\bm{u}$ (called
${\bm{E}}_{\perp}$ component), where $\bm{u}$ is the normal vector to the
interface, corresponds to the $q\rightarrow 0$ limit of our ${{\bm{s}}}^{{\rm
sh}}$, and density proportional to ${\bm{M}}\times({\bm{E}}\times\bm{u})$
(called ${\bm{M}}_{\perp}$) corresponds to our ${{\bm{s}}}^{{\perp}}$. Thus
the susceptibilities in our notation read

$\displaystyle\chi_{s}^{{\bm{E}}_{\perp}}$
$\displaystyle=\overline{\lambda}_{\rm sh}\lim_{q\rightarrow
0}\frac{1}{q}\left[\overline{s}^{{\rm sh}(2)}(q)+\frac{\pi{v_{\rm
i}}}{{\epsilon_{F}}}\overline{s}^{{\rm sh}(3)}(q)\right]$
$\displaystyle\chi_{s}^{{\bm{M}}_{\perp}}$
$\displaystyle=\overline{\lambda}_{\rm sh}\lim_{q\rightarrow
0}\frac{1}{q}\left[\overline{s}^{{\perp}(2)}(q)+\frac{\pi{v_{\rm
i}}}{{\epsilon_{F}}}\overline{s}^{{\perp}(3)}(q)\right]$

The magnitude of $\frac{e}{\pi}\frac{1}{{k_{F}}{\epsilon_{F}}}$ in
$\overline{\lambda}_{\rm sh}$ is evaluated for ${k_{F}}^{-1}=1\AA$ and
${\epsilon_{F}}=1$ eV as $\frac{e}{\pi}\frac{1}{{k_{F}}{\epsilon_{F}}}=3\times
10^{-11}$ m/V. Considering the numerical factor of $\overline{s}^{{\rm
sh}(2)}\sim 0.02$ in the Stoner regime (Fig. 2), our estimate of the
susceptibility in the Stoner regime is $\chi_{s}\sim 6\times
10^{-13}\times\frac{\epsilon_{\rm so}}{{\epsilon_{F}}}\frac{n_{\rm i}v_{\rm
i}}{{\epsilon_{F}}}$ m/V. At $q=0$, the susceptibility is enhanced by a factor
of $\eta^{-1}$ due to a peak of ${\rm Re}\overline{\Sigma}^{{\rm ra}}$. In the
case of Pt heavy metal, $\chi_{s}\sim 10^{-11}$ m/V at $q=0$ is predicted in
Ref. [10]. Our estimate based on a simplified model therefore agrees by the
order of magnitudes with Ref. [10] if $\frac{\epsilon_{\rm
so}}{{\epsilon_{F}}}\frac{n_{\rm i}v_{\rm i}}{{\epsilon_{F}}}\sim 1$ is
assumed for heavy metals. Considering the fact that spin-orbit torque is
highly efficient for magnetization switching [3], spin Hall effect in the
Stoner regime is expected to induce various finite-$q$ magnetization dynamics.

##  3 Rashba spin-orbit interaction

Spin-orbit interaction due to random impurities conserves momentum, and thus
the wave vector ${\bm{q}}$ driving the finite-${\bm{q}}$ spin Hall effect
needs to be given by the applied electric field. Due to the momentum
conserving response function, the spin accumulation and the electric field
related by the direct and inverse spin Hall effects have the same wave length.
It is not, however, easy to access by electrical methods the Stoner regime,
which corresponds to a high wave vector regime of the order of lattice
constant in ferromagnetic metals. Doped ferromagnetic semiconductors [12]
having an orders of magnitude longer wave length is one promising system for
finite-$q$ spin Hall effect. From the viewpoint of finite $q$, use of surface
and interface is of special interest, as finite wave vector transfer is
allowed due to the inversion symmetry breaking.

We consider the Rashba spin-orbit interaction, whose Hamiltonian is

$\displaystyle H_{\rm R}$ $\displaystyle=\int d{\bm{r}}\alpha_{{\rm
R},i}\epsilon_{ijk}c^{\dagger}\hat{p}_{j}\sigma_{k}c$ (22)

where ${\bm{\alpha}}_{\rm R}$ is the Rashba vector. We consider the case where
${\bm{\alpha}}_{\rm R}$ is along the $z$ axis and localized at the surface
chosen as $z=0$; Namely, ${\bm{\alpha}}_{\rm R}({\bm{r}})=\alpha_{\rm
R}(z)\hat{\bm{z}}$ with a coefficient $\alpha_{\rm
R}(z)=\overline{\alpha}_{\rm R}a\delta(z)$ ($a$ is the lattice constant
introduced for using conventional unit). In the momentum representation, we
have

$\displaystyle H_{\rm R}$ $\displaystyle=-\sum_{{\bm{p}}{\bm{k}}}\alpha_{\rm
R}({\bm{p}})\epsilon_{ijk}k_{j}c^{\dagger}_{{\bm{k}}+\frac{{\bm{p}}}{2}}\sigma_{k}c_{{\bm{k}}-\frac{{\bm{q}}}{2}}$
(23)

The interaction modifies the electron current density as

$\displaystyle j_{\alpha}({\bm{r}})$
$\displaystyle=e^{-i({\bm{q}}-{\bm{p}})\cdot{\bm{r}}}c^{\dagger}_{{\bm{k}}+{\bm{q}}}\left[\frac{k_{\alpha}}{m}\delta_{{\bm{p}},0}-\epsilon_{\alpha
lm}\alpha_{{\rm R},l}({\bm{p}})\sigma_{m}\right]c_{{\bm{k}}}$ (24)

We consider the Rashba interaction to the linear order to evaluate the
response function of electron spin and current for a uniform applied electric
field. The Feynman diagrams representing the response function are in Fig. 3.

Figure 3:  Diagrams contributing to the response function of electron spin and
current at the linear order of the Rashba spin-orbit interaction denoted by
circles. Diagram (b) is the one with a current vertex modified by the Rashba
interaction.

The response functions defined by the ratio of induced spin in the direction
$\alpha$ to the applied electric field in the $i$ direction is ((a) and (b)
correspond to Fig. 3.)

$\displaystyle s_{{\rm R},i}^{{\rm(a)},\alpha}(p)$
$\displaystyle=-2\epsilon_{jl\beta}\sum_{p}\alpha_{{\rm
R},j}(p)\frac{1}{m}\frac{1}{V}\sum_{{\bm{k}}}\left(k-\frac{p}{2}\right)_{l}k_{i}{\rm
Re}\tr\left[\sigma_{\alpha}g^{\rm r}_{{\bm{k}}-{\bm{p}}}\sigma_{\beta}g^{\rm
r}_{{\bm{k}}}g^{\rm a}_{{\bm{k}}}\right]$ $\displaystyle s_{{\rm
R},i}^{{\rm(b)},\alpha}(p)$
$\displaystyle=\epsilon_{ij\beta}\sum_{p}\alpha_{{\rm
R},j}(p)\frac{1}{V}\sum_{{\bm{k}}}{\rm Re}\tr\left[\sigma_{\alpha}g^{\rm
r}_{{\bm{k}}-\frac{{\bm{p}}}{2}}\sigma_{\beta}g^{\rm
a}_{{\bm{k}}+\frac{{\bm{p}}}{2}}\right]$ (25)

Evaluating the trace for the spin, we obtain the spin Hall and anomalous Hall
components, $s_{\rm R}^{{\rm(n)},\alpha}=s_{\rm R}^{{\rm
sh}{\rm(n)},\alpha}+s_{\rm R}^{{\perp}{\rm(n)},\alpha}$ (n=a,b). In the
geometry we consider, $\alpha_{{\rm R},j}\propto\delta_{jz}$, and we obtain
the sum of contributions (a) and (b) as

$\displaystyle s_{{\rm R},i}^{{\rm sh},\alpha}(p)$
$\displaystyle=2\epsilon_{ij\alpha}\sum_{p}\alpha_{{\rm
R},j}(p)\sum_{\sigma}{\rm Re}\biggl{[}J_{\rm R}^{\sigma}(p)+\Sigma_{\rm
R}^{\sigma}(p)\biggr{]}$ $\displaystyle s_{{\rm R},i}^{\perp,\alpha}(p)$
$\displaystyle=2\epsilon_{ij\beta}\epsilon_{\alpha\beta\gamma}\hat{M}_{\gamma}\sum_{p}\alpha_{{\rm
R},j}(p)\sum_{\sigma}\sigma{\rm Im}\biggl{[}J_{\rm R}^{\sigma}(p)+\Sigma_{\rm
R}^{\sigma}(p)\biggr{]}$ (26)

where

$\displaystyle J_{\rm R}^{\sigma}(p)$
$\displaystyle\equiv\frac{1}{V}\sum_{{\bm{k}}}\frac{k^{2}}{3m}g^{\rm
r}_{{\bm{k}}-{\bm{p}},\sigma}g^{\rm r}_{{\bm{k}},-\sigma}g^{\rm
a}_{{\bm{k}},-\sigma}$ $\displaystyle\Sigma_{\rm R}^{\sigma}(p)$
$\displaystyle\equiv\frac{1}{V}\sum_{{\bm{k}}}g^{\rm
r}_{{\bm{k}}-\frac{{\bm{p}}}{2},\sigma}g^{\rm
a}_{{\bm{k}}+\frac{{\bm{p}}}{2},-\sigma}$ (27)

As was shown in Ref. [4], the spin Hall contribution $s_{{\rm R},i}^{{\rm
sh}}(p)$ vanishes linearly at $p=0$ if $M=0$. Numerically calculated response
function is plotted in Fig. 4. The Rashba system has a broad response for
$q\lesssim 2{k_{F}}$, meaning that the real-space distribution of the induced
spin has a localized nature. This fact is consistent with the local spin-
charge conversion picture of the Rashba-Edelstein effect.

The broad spin response of the localized Rashba interaction would be useful
for the emission of finite-$q$ magnons.

Figure 4:  The response functions $s_{{\rm R},i}^{{\rm sh},\alpha}$ and
$s_{{\rm R},i}^{\perp,\alpha}$ of electron spin and current at the linear
order of the Rashba spin-orbit interaction as function of momentum transfer
$p$ normalized by ${k_{F}}$. The response for $M=0$ near $p=0$ is
qualitatively different from those for $M\neq 0$. Damping constant is
$\eta/{\epsilon_{F}}=0.1$.

##  4 Discussion

The electron spin density induced by the spin Hall effect couples to the
magnetization via the exchange coupling (with coupling constant $J_{sM}$). The
coupling leads to a so called the spin-orbit torque [3] and drives
magnetization and induces magnon emission. Writing the transverse induced spin
as ${{\bm{s}}}_{\rm ind}$, i.e., ${{\bm{s}}}_{\rm ind}={{\bm{s}}}^{\rm
sh}+{{\bm{s}}}^{\perp}$ and ${{\bm{s}}}_{\rm ind}={{\bm{s}}}^{\rm sh}_{\rm
R}+{{\bm{s}}}^{\perp}_{\rm R}$ for spin-orbit interaction due to random
impurities and the Rashba interaction, respectively, the exchange coupling
Hamiltonian is $H_{sM}\equiv-J_{sM}{{\bm{s}}}_{\rm ind}\cdot{\bm{M}}$. In
terms of the Holstein-Primakoff boson operators, $a$ and $a^{\dagger}$, the
magnetization (with quantum axis along $\hat{\bm{z}}$) is
${\bm{M}}=M\hat{\bm{z}}+\sqrt{\frac{M}{2}}(a^{\dagger}+a,i(a^{\dagger}-a),0)+O((a,a^{\dagger})^{2})$
to the linear order. The induced transverse spin thus couples linearly to the
magnon operators as

$\displaystyle H_{sM}=-J_{sM}\sqrt{\frac{M}{2}}(s^{+}_{\rm
ind}a^{\dagger}+s^{-}_{\rm ind}a)$ (28)

where $s^{\pm}_{\rm ind}\equiv s^{x}_{\rm ind}\pm is^{y}_{\rm ind}$. The wave
vector of the induced spin is therefore transferred directly to magnons, and
this feature would be useful for magnon manipulations [13]. Experimentally,
magnon energy can be chosen by use of an alternating current, and the magnon
emission amplitude in this case is expected to be proportional to
$s^{\pm}_{\rm ind}(q_{\rm m})$, where $q_{\rm m}$ is the wave vector for the
excited magnon.

##  5 Summary

We have studied the spin Hall effect in ferromagnets at finite wave vector by
calculating the linear response function of spin and electric current for the
spin-orbit interactions due to random impurity and the Rashba interaction. It
was shown that the response function shows a broad response in the regime
$q\lesssim 2{k_{F}}$ due to the gapless Stoner excitation, suggesting a local
nature of the spin-charge conversion effects. For experimental study of
finite-$q$ regime, dilute magnetic semiconductors or use of STM would be
suitable.

###### Acknowledgements.

The author thank E. Saitoh for valuable discussion. This study was supported
by a Grant-in-Aid for Scientific Research (B) (No. 17H02929 and No. 21H01034)
from the Japan Society for the Promotion of Science.

## References

  * [1] M. Dyakonov and V. I. Perel: Sov. Phys. JETP Lett. 13 (1971) 467.
  * [2] J. E. Hirsch: Phys. Rev. Lett. 83 (1999) 1834.
  * [3] A. Manchon, J. Železný, I. M. Miron, T. Jungwirth, J. Sinova, A. Thiaville, K. Garello, and P. Gambardella: Rev. Mod. Phys. 91 (2019) 035004.
  * [4] G. Tatara: Phys. Rev. B 98 (2018) 174422.
  * [5] G. Tatara: Phys. Rev. B 104 (2021) 184414.
  * [6] M. Matsuo, J. Ieda, E. Saitoh, and S. Maekawa: Phys. Rev. B 84 (2011) 104410.
  * [7] K. M. D. Hals and A. Brataas: Phys. Rev. B 88 (2013) 085423.
  * [8] F. Freimuth, S. Blügel, and Y. Mokrousov: Journal of Physics: Condensed Matter 26 (2014) 104202.
  * [9] F. Freimuth, S. Blügel, and Y. Mokrousov: Phys. Rev. B 92 (2015) 064415\.
  * [10] L. Salemi, M. Berritta, and P. M. Oppeneer: Phys. Rev. Materials 5 (2021) 074407.
  * [11] A. Crépieux and P. Bruno: Phys. Rev. B 64 (2001) 014416.
  * [12] J. König, H.-H. Lin, and A. H. MacDonald: Phys. Rev. Lett. 84 (2000) 5628.
  * [13] A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands: Nat Phys 11 (2015) 453.  Review.