Year: 2021

Abstract:
Most convergence guarantees for stochastic gradient descent with momentum (SGDm) rely on iid sampling. Yet, SGDm is often used outside this regime, in settings with temporally correlated input samples such as continual learning and reinforcement learning. Existing work has shown that SGDm with a decaying stepsize can converge under Markovian temporal correlation. In this work, we show that SGDm under covariate shift with a fixed step-size can be unstable and diverge. In particular, we show SGDm under covariate shift is a parametric oscillator, and so can suffer from a phenomenon known as resonance. We approximate the learning system as a time varying system of ordinary differential equations, and leverage existing theory to characterize the system’s divergence/convergence as resonant/nonresonant modes. The theoretical result is limited to the linear setting with periodic covariate shift, so we empirically supplement this result to show that resonance phenomena persist even under non-periodic covariate shift, nonlinear dynamics with neural networks, and optimizers other than SGDm.

heading: 1 INTRODUCTION
text: Stochastic gradient descent (SGD) (Robbins & Monro, 1951) – and its variants such as Adagrad (Duchi et al., 2011), ADAM (Kingma & Ba, 2014) and RMSprop (Hinton et al., 2012) – are very widely used optimization algorithms across machine learning. SGD is conceptually straightforward, easy to implement, and often performs well in practice. Among the variants of SGD, accelerated versions based on Polyak’s or Nesterov’s acceleration (Polyak, 1964; Nesterov, 1983), known generally as Stochastic Gradient Descent with Momentum (SGDm), are used widely due to the improvements in convergence rate they offer. SGDm can give up to a quadratic speedup to SGD on many functions, and is in fact optimal among all methods having only information about the gradient at consecutive iterates for convex and Lipschitz continuous optimization problems (Nesterov, 2004; Goh, 2017). SGDm has the same computational complexity as SGD, but exhibits superior convergence rates under reasonable assumptions (Su et al., 2014).
These convergence results for SGDm, however, rely on independent and identically distributed (iid) sampling. Little is known about the convergence properties of SGDm under non-iid sampling, yet the non-iid setting is critical. In many machine learning problems it is expensive or impossible to obtain iid samples. In online learning (Rakhlin et al., 2010) and reinforcement learning (RL) (Sutton & Barto, 2018), the data becomes available in a sequential order and there is a temporal dependence among the samples. There is a particularly strong temporal dependence in RL, where observed states are sampled according to the transition dynamics of the underlying Markov decision process (MDP). Federated learning (Hsieh et al., 2020) and time-series learning (Kuznetsov & Mohri, 2014) provide further examples of when non-iid sampling is essential to the learning problem.
Without momentum, SGD’s convergence rate has been examined under non-iid sampling with the stochastic approximation framework in (Benveniste et al., 2012; Kushner & Yin, 2003), and more recently under specific assumptions of ergodicity (Duchi et al., 2012) or Markovian sampling (Nagaraj et al., 2020; Doan et al., 2020b; Sun et al., 2018). Convergence rates under Markovian sampling are also known for ADAM-type algorithms when applied to policy gradient and temporal difference learning (Xiong et al., 2020). To our knowledge, however, there has been little work on providing convergence rates or guarantees for SGDm under non-iid sampling. In (Doan et al., 2020a), a progress bound is provided for SGDm under Markovian sampling based on mixing time—the time required
for a distribution’s convergence toward its stationary distribution—along with a convergence rate guarantee under decaying step-sizes. In this work, we assume a fixed step-size, since it is a common choice in the online setting, especially when the practitioner is unsure of mixing rate or stationarity.
There is a broad literature using ordinary differential equations (ODEs) to analyze gradient descent methods by approximating descent updates as continuous-time flows. Early work is comprehensively discussed in (Kushner & Yin, 2003; Benveniste et al., 2012). Despite the age and establishment of the linear setting, the gradient flow lens continues to reveal new insights (e.g. implicit rank reduction (Arora et al., 2019)) and new perspectives on old insights (e.g. regularization of early stopping (Ali et al., 2019)). Recent work has paid particular attention to the flow induced by momentum accelerated methods (Su et al., 2014; Wibisono et al., 2016; Wilson et al., 2016; Scieur et al., 2017; Muehlebach & Jordan, 2021; Diakonikolas & Orecchia, 2019; Muehlebach & Jordan, 2019; Li et al., 2017; Simsekli et al., 2020; Betancourt et al., 2018; Kovachki & Stuart, 2019; Attouch et al., 2018; Shi et al., 2019; Siegel, 2019; Zhang et al., 2018; Berthier et al., 2021; Kovachki & Stuart, 2021). In these works, it is demonstrated that linear regression under iid sampling with SGDm can be represented as an ODE resembling a harmonic oscillator, sometimes with a time-decaying damping coefficient.
Contributions: In this work, we show that non-iid sampling, due to covariate shift, induces a related system: the parametric oscillator, a harmonic oscillator having coefficients which vary over time in a manner capable of exponentially exciting the system (Mumford, 1960; Csörgő & Hatvani, 2010; Halanay, 1966). We characterize the underlying time-varying ODE, which requires a novel approach to incorporate covariate shift, and develop specific conditions on covariate shift that lead to divergence in SGDm. We provide an empirical design, with synthetic data, to systematically test the response of the learning system to different covariate shift frequencies. We empirically validate that resonancedriven divergence occurs in a learning system which aligns well with theoretical assumptions. We follow with similar demonstrations on learning systems which progressively relax further and further away from our theoretical assumptions, including non-periodic covariate shift, learning with neural networks and optimizers other than SDGm. The implications of this work are that resonance can occur in learning systems with momentum, causing instability or poor accuracy, simply due to certain temporal correlations in the input sequence. Our results provide a novel direction to better understand the properties of our learning systems for time varying inputs—namely beyond iid sampling—and point to a gap in the robustness of our algorithms that merits further investigation.

heading: 2 PROBLEM SETTING
text: We investigate the effect of non-iid sampling on a model optimized using SGDm. We assume labelled training data sampled from discreteindexed, real-valued stochastic processes {Xk}k∈N and {Yk}k∈N such that {Xk}k∈N has a unique stationary distribution Π. Yk is a function of Xk with zero-mean observation noise, Yk = f(Xk) + k. The learning algorithm does not have access to iid samples from Π. Instead, the learning algorithm may only update parame-
ters θk at time k using samples zk = (xk, yk), where Xk’s marginal over time converges to Π. The goal is to train a model with parameters θ to minimize an objective function L(z; θ) with respect to the stationary distribution Π: θ∗ = arg minθ∈Rd EΠ[L(z; θ)] = ∫ R L(z; θ)dΠ(x). An example of this problem setting is illustrated in Figure 1.
This setting is identical to that explored in prior work on Markovian sampling (Sun et al., 2018; Xiong et al., 2020; Nagaraj et al., 2020), but we do not require that our stochastic processes adhere to the Markov property. Note, also, that the underlying functional relationship between the inputs, Xk, and the targets, Yk, does not change over time. That is, the non-iid sampling is a result of covariate shift rather than a changing relationship between the inputs and targets.
We consider Polyak’s Heavy Ball method (Polyak, 1964) with the formulation from Sutskever et al. (2013), namely SGD with momentum (SGDm)
vk+1 = µvk − η∇θL(zk; θk) θk+1 = θk + vk+1 (1)
where η ≥ 0 is the learning rate and µ ∈ (0, 1) the momentum coefficient. Our goal is to understand the behavior of SGDm when learning on this non-iid sequence zk.

heading: 3 COVARIATE SHIFT AS A DRIVING FORCE
text: In this section we characterize SGDm as a discretization of a particular parametric oscillator in continuous time. A parametric oscillator resembles a harmonic oscillator ODE, but has time-varying system coefficients which are capable of driving the system. It is well understood that parametric oscillators can suffer from global solution instability due to coefficients oscillating at particular frequencies (Halanay, 1966), a condition known as parametric resonance. We show parametric resonance conditions sufficient to induce exponential divergence in SGDm. We provide proof sketches and defer complete proofs to Appendix A.3.
Notation: Capital letters are matrices and random variables. {Xk} is a stochastic process, shorthand for {Xk}k∈N. k and t index discrete and continuous time, respectively. When a discrete sequence and a function approximate each other, they share the same symbol and are differentiated by subscript k and function argument t, e.g. the sequence {θk} and the function θ(t). Pk is the joint distribution of Xk and Yk. For loss function L(z; θ) on training pair z = (x, y) and weights θ, we denote the time-varying expected gradient gk(θ) := EPk [∇θL(z; θ)]

heading: 3.1 COVARIATE SHIFT INDUCES A LINEAR TIME-VARYING EXPECTED GRADIENT
text: We begin by showing that linear least squares regression with covariate shift and a fixed target induces a linear time-varying expected loss gradient. That is, the gradient function gk(θ) = Bk(θ − θ∗) for some matrix Bk that varies with time, due to covariate shift.
Assumption 1. The covariate generating process {Xk}k∈N is not identically distributed: Corr(Xk1 , Xk1) 6= Corr(Xk2 , Xk2) for some k1, k2 ∈ N. Assumption 2. The targets Yk are a fixed linear function of Xk with iid zero-mean observation noise, k such that E[ k] = 0. That is, Yk = 〈θ∗, Xk〉+ k where θ∗ ∈ Rd is fixed for all k.
Proposition 1. Under Assumptions 1, 2, a linear model with weights θk making predictions Ŷk = 〈θk, Xk〉 with a mean squared error (MSE) objective will induce time-varying linear expected loss gradients gk(θk) = Bk(θk − θ∗) for all k ∈ N, where Bk ∈ Rd×d and Bk1 6= Bk2 for some k1, k2.
Proof Sketch: We fix a time step k and take the expectation of the MSE loss over observation noise, and over the distribution of Xk. A time-varying expected loss surface remains, and the linear setting implies the loss gradient is a linear function for all steps k, with Bk depending on Xk:
Bk = 2Corr(Xk, Xk) = 2Cov(Xk, Xk) + 2EPk [Xk]EPk [Xk]T (2)

heading: 3.2 ODE CORRESPONDENCE
text: Next, we show that the iterates {θk} generated by SGDm are a first order numerical integration of a particular ODE. The procedure is similar to (Muehlebach & Jordan, 2021), but due to the time-varying loss gradient, we must pay specific attention to the conditions on the ODE necessary to have integration consistency, . Let B(t) be a matrix-valued function, Lipschitz continuous in t, such that {Bk} are samples spaced √ η apart, i.e. Bk = B( √ ηk). The existence and uniqueness of this B(t) are ensured if we assume the {Xk} are sampled from an underlying continuous-time X(t). Proposition 2. The SGDm iterates {θk} numerically integrate the ODE system in equation 3 with integration step √ η and first order consistency.
θ̈(t) + 1− µ √ η θ̇(t) +B(t)(θ(t)− θ∗) = 0 (3)
Proof Sketch: First, the ODE in equation 3 is converted to first order linear form, and split into a sum of two separate systems. The sum’s terms have implicit and explicit Euler integrations; using operator splitting, these can be composed into a consistent numerical integrator for linear time-varying systems. Then, it is shown that the composed integrator is precisely equivalent to SGDm when the integration time step is √ η.
The first order consistency guarantee in Proposition 2 means that when the step-size η is small and initial conditions agree between continuous and discrete time, i.e. θ0 = θ(0), the continuous and discrete time trajectories approximate each other as θk ≈ θ(k √ η), where the difference between them depends on the integration time step h = √ η (Hairer et al., 2006). Therefore, assuming θ0 = θ(0), the difference accrued in one step of k (local error) is θ1 = θ( √ η) +O(η) and over arbitrarily many steps (global error) is θk = θ( √ ηk) +O(ηk).

heading: 3.3 PARAMETRIC RESONANCE FOR ODE CONVERGENCE AND DIVERGENCE
text: The conditions sufficient for convergence and divergence in equation 3 may be shown using established dynamical systems theory. The conditions sufficient for divergence are precisely the conditions for parametric resonance, when B(t) is periodic. Note that periodicity in the first or second moment of {Xk} is sufficient to induce periodicity in B(t). As per elementary results in linear ODE theory, the system in equation 3 can be transformed into a linear time-varying first order form
ξ̇(t) = A(t)ξ(t) A(t) = [ 0d×d Id×d B(t) 1−µ√η Id×d ] (4)
where solution trajectories θ(t) of equation 3 are embedded in solution trajectories ξ(t) of equation 4. Equation 4 admits a fundamental solution matrix1 ψ(t) such that the spectral radius ρ of ψ(T ) characterizes ODE instability, which implies divergence for SGDm, in Theorem 1. Theorem 1. When B(t) is periodic such that B(t) = B(t+ T ) for some T > 0, the spectral radius ρ of ψ(T ) characterizes the stability of solution trajectories of equation 3 as follows:
• ρ > 1 =⇒ trivial solution θ(t) = θ∗ is unstable. All other solutions diverge as θ(t) → ∞ exponentially with rate ρ. • ρ < 1 =⇒ trivial solution is asymptotically stable, all other solutions converge as θ(t) → θ∗ exponentially with rate ρ.
Proof Sketch: ODE 3’s stability depends on the relationship between the coefficient 1−µ√η and the expected gradient signal because they determine the spectral radius of ψ(T ). The spectral radius ρ is the convergence/divergence rate per period T towards/away from the stationary point θ∗.
The spectral radius conditions characterize when ODE solution trajectories θ(t) will converge or diverge, and Proposition 2 tells us that these solution trajectories are approximations of discrete time SGDm trajectories {θk} under identical initialization. But does convergence or divergence of θ(t) imply the same for SGDm? Experiment 4.1 empirically suggests that the approximation is sufficient, since the boundary at ρ = 1 in Figure 2a agrees with both continuous and discrete time. However, given Theorem 1 and Proposition 2, we have a theoretical guarantee of SGDm’s divergence, but not its convergence. This is because the divergence rate of θ(t) is exponential, and Proposition 2 provides a bound on long tail behaviour as θk = θ( √ ηk) +O(ηk). Since the approximation error is linear in time k, and the divergence rate is exponential, the divergence rate dominates, and we are guaranteed that a diverging θ(t) corresponds to a diverging {θk}. However, the same argument does not imply that a convergent θ(t) corresponds to a convergent {θk}, because the linear error bound technically permits the discrete trajectory to escape θ∗ at a linear rate. Empirically, we see agreement for both convergent and divergent cases, but we defer theoretical proof to future work.
Summary: There is a chain of dependencies starting from the non-i.i.d. sampling in {Xk}, and ending at the monodromy’s spectral radius.
{Xk} → X(t)→ B(t)→ A(t)→ ψ(t)→ ψ(T )→ ρ In order, we have the discrete stochastic process {Xk} from which training inputs are sampled, and its underlying continuous stochastic process X(t). If {Xk} and X(t) were i.i.d., then the matrix B
1A fundamental solution matrix for system equation 4 is any matrix-valued function of time ψ(t) whose columns are linearly independent solutions to equation 4. We choose ψ(t) such that ψ(0) = I2d×2d.
would be constant, but the non-i.d.d. nature means B(t) is a function of time. The matrix A(t) is the matrix describing SGDm’s continuous time dynamics as a linear time-varying ODE ξ̇ = A(t)ξ, and the matrix B(t) is a simply a submatrix of A(t). Since we have a linear ODE, the space of solution trajectories is spanned by the columns of the ODE’s fundamental solution matrix ψ(t). Moreover, since we have a periodic A(t) with period T , the stability of all solutions can be determined simply from the largest eigenvalue of ψ(T ), a.k.a. its spectral radius ρ.
The intuition behind ρ’s importance comes from the following fact: after each elapsed period T , the phase space (including weight space) is subject to the linear transformation ψ(T ). So as time increases, say n periods, the linear transformation is applied iteratively as ψ(T )n, so it is clear that any eigenvalue larger than unity means all solutions eventually diverge exponentially. This is precisely the parametric resonance condition. For a detailed example of all these quantities characterizing the divergence of a simple learning system, see Appendix A.1.
Technical Novelty: There are two key novel technical contributions to obtain the results above. (1) The ODE correspondence for non-iid data in combination with momentum methods is new. The non-iid data, which induces time variation in ODE dynamics, is subtle to address, because elementary ODE correspondence (i.e. consistency guarantees) require time invariant ODE dynamics. Our proof specifically addresses this difficulty using operator splitting theory from recent numerical integration literature. (2) Once we have this ODE, it is straightforward to recognize that it is a parametric oscillator. The utility of this connection, however, is that it opens up entirely new avenues for analysis, since there is deeply established literature on the dynamics of these systems. The primary second technical novelty, therefore, is to leverage this connection, and use Floquet theory to give precise divergence conditions for a given η and µ, under periodic covariate shift.

heading: 4 VALIDATING THEORY & ABLATING TOWARDS CONDITIONS IN THE WILD
text: We now empirically evaluate the effect of resonance in learning problems. First, Experiment 4.1 validates the theoretical predictions by investigating the dynamics of a learning problem which closely matches the theoretical assumptions in Section 3. Specifically, we assess whether or not the spectral radius ρ predicts optimizer convergence or divergence. Then, the remaining Experiments 4.2 - 4.6 test for this phenomena in settings beyond our theory. Each experiment makes a single step away from a theoretical assumption, first relaxing the periodicity assumption, then using stochastic rather than expected gradients and then moving to more widely used nonlinear models (neural networks) and optimizers (ADAM).
Across all experiments, input samples at each training step k are drawn from Gaussian distributions with diagonal covariance matrices, i.e. Xk ∼ N (x̄k, cId×d), and we induce covariate shift by constructing a time-varying mean sequence {x̄k}. Each mean sequence is designed such that (a) its frequency content can be swept with a single parameter f or T and (b) iid sampling is induced by setting f = 0 or T = 0, so that each experiment has an iid baseline for comparison. The specification of x̄k is provided in Table 1 for each experiment. This empirical design allows us to systematically test response to different frequencies. See Appendix A.4 for specification of all experiment details.

heading: 4.1 VALIDATING THEORY
text: We start in a setting as close as possible to the theoretical predictions, with linear regression for a quadratic loss, and covariate shift such that the mean E[Xk] varies as a strict sinusoid. We perform regression in two weights (i.e. inputs Xk and labels Yk are scalar-valued). We draw 20 samples from each Xk, so that the loss gradient for each time step k is close to its expected value.
Since we have only two weights, the learning system’s underlying ODE can be easily specified, such that the fundamental solution matrix can be numerically computed. As per Theorem 1, we can use the spectral radius of the fundamental solution matrix evaluated at time T to make theoretical predictions of where the system should converge or diverge. As depicted in Figure 2a, the theory agrees very well with empirical results. This suggests that parametric resonance is indeed the dominant mechanism behind SGDm divergence under covariate shift. Refer to example A.1 for the procedure used to compute the theoretical predictions (i.e. spectral radii ρ), and Figure 6a for the full surface from which the contour lines in Figure 2a are rendered.

heading: 4.2 ABLATING PERIODICITY
text: We repeat Experiment 4.1, but with the mean of Xk varying stochastically instead of deterministically. The result in Theorem 1 assumes the expected gradient varies periodically over time given fixed weights θ. But even with aperiodic and/or stochastic time variation, our LTV system equation 4 might be similarly susceptible to instability; the theory does not have periodicity as a necessary condition of divergence. To test this, we replace our periodic covariate shift mean with a mean that moves according to an AR(2) process. The AR(2) process is tuned to have a frequency peak exactly at the frequency of the sinusoidal covariate shift of Experiment 4.1, to change only the nature of the shift—aperiodicity—rather than the frequency.
We can see from the heatmap in Figure 2b that under this setting we observe very similar resonance behaviour in the learning system, which aligns well with the identical predicted stability regions. Note that both theoretical and empirical results in Figures 2a and 2b suggest that sufficiently low momentum values µ mitigate the resonance phenomenon, which agrees with the role it plays in the ODE system: decreasing µ increases damping. We also observe that resonance regions shrink as step size η is decreased; see the Appendix for plots demonstrating this trend. Analytically characterizing the bounds of stable µ, η in terms of system properties is an interesting future direction.

heading: 4.3 ABLATING EXPECTED GRADIENT
text: Here we perform linear regression with periodic mean covariate shift, and ablate the number of samples drawn from each Xk to show the effect of increasing noise in the gradient signal. Linear regression is performed mapping R5 → R, with the weights learned via SGDm. The mean sequence x̄k oscillates with strict periodicity between ±x̄, where x̄ is a unit norm 5-vector randomly chosen for each run. i.e. the mean signal is a square wave in R5 with period T .
Linear regression implies a quadratic loss, and periodic covariate shift implies a loss with periodic time variation in expectation, so we are very near to the setting in which Theorem 1 is directly applicable. However, as we ablate from 5 samples to 1 sample drawn from each Xk, we move further away from the expected gradient, towards the fully stochastic gradient setting. As we can see in Figure 3 resonance is dampened by stochasticity in the gradient signal, but nonetheless is still present.

heading: 4.4 ABLATING PERIODICITY FURTHER
text: We further depart from strictly periodic covariate shift by having piece-wise stationary means and nonsmooth changes: we randomly sample a new mean from a normal distribution N (0, v2) after every T update steps. We further investigate the impact of this covariate shift under (a) increasing variance v2 and (b) increasing dimensionality of the inputs.
We observe divergence characteristic of the parametric resonance, with smaller T having a highly divergent response. In Figure 4a we can see that a higher variance results in much more instability, though even for smaller variance the distance to θ∗ is still on the order of 103. This results makes sense, as the variance is connected to the driving signal amplitude, on which parametric resonance has a strong dependence. We also observe a strong dependence on the number of input space dimensions, in Figure 4b. This aligns with the fact that the expected norm of samples drawn from a multivariate Gaussian increases with dimensionality, even with a fixed covariance scale (of 0.1). This result suggests that resonance phenomena may actually be exacerbated for higher-dimensional inputs, which is the setting we expect to see in practice.

heading: 4.5 ABLATING OPTIMIZER LINEARITY
text: So far our experiments have ablated periodicity, input dimensionality, and gradient stochasticity. As per the conditions in Section 3, each of these ablated systems still corresponds to a discretization of an LTV ODE. We now ask whether resonance can be observed in systems which do not correspond to a discretization of our ODE. To investigate this hypothesis, we replace SGDm with ADAM as an update rule for our learning algorithm. We no longer have the ODE representation of this learning system, but can empirically investigate the behaviour of the system under covariate shift.
We use the same stochastic mean switching problem as the previous experiment, where we regress from 5 dimensions to 1 with input samples. Again, we measure the distance ||θk − θ∗||, and we vary the ADAM parameter β1. Similar to the SGDm optimizer, we can see in Figure 4c that there is a band of mean switching intervals T for which convergence is worse. But unlike SGDm, there is no divergence. Proper parametric resonance induces exponential divergence, which is why the previous results were presented with Euclidean distances between weights on the order of 1019. Note that in Figure 4c, the frequency response in weight space distance is instead measured on the order of 100.

heading: 4.6 ABLATING MODEL LINEARITY
text: Our final step away from ODE equation 3 abandons the quadratic loss surface. We fit a nonlinear target function using a fully connected ReLU network, reusing the previous covariate shift. In previous experiments, we measured the distance ||θk − θ∗||. Now a unique target is not guaranteed,
(c) ADAM w/ stochastic x̄k, β1 sensitivity
so frequency response is in terms of loss. In order to amplify frequency response, 10 samples were drawn from each Xk to reduce gradient signal noise.
Figure 5 shows a band of mean switching intervals T which damage convergence, which is characteristic of resonance. Though we do not see divergence, the impact on accuracy is problematic. In this case, models with loss > 0.3 fit extremely poorly, but those with < 0.05 perform well.

heading: 5 LIMITATIONS AND PRACTICALITY
text: This work is a first step towards understanding parametric resonance due to covariate shift. Naturally, there are limitations and questions about broader applicability, which we discuss in this section.
Limitations and Next Steps: In our theoretical analysis, we take expectation over observation noise and over the distribution of inputs Xk, so the gradient coefficient matrix B(t) in the ODE equation 3 is deterministic. This allows us to rigorously charac-
terize when resonance will occur for periodic B(t). Future work might employ stochastic differential to avoid taking expectation. The case when B(t) is aperiodic and stochastic is of particular interest, since Experiments 4.2 and 4.4 demonstrate that resonance still drives divergence. In particular, Experiment 4.2 aligns exactly with the theoretical predictions, suggesting that the resonance is in response to the frequency content of Xt, not its periodicity.
For the sake of consistency across experiments, our non-iid sampling is always induced by Gaussian Xk with fixed diagonal covariance matrix. But Assumption 1 is far more relaxed. Even if one wishes to adhere to the strict periodicity requirements of Theorem 1, it is possible to induce time variation by varying other aspects of Xk. Comparing the effects of these different time variations is an interesting future direction.
When assessing a new phenomenon like resonance, assessment on simple problems is a critical first step, and we have endeavored to be comprehensive in that assessment. But we only scratch the surface of more complex problems, including more advanced optimizers and nonlinear models. Similarly, for the sake of experimental control and interpretability, we assess resonance on synthetic data instead of real world data. While it is obvious that many sources of real data have nontrivial frequency content (e.g. audio samples, machine sensors, etc.) it is not yet clear when realistic frequency content will resonate with the system trying to learn from it. The two most exciting future directions for this work are to address these two gaps: to extend into the frequency response of more complex optimizers and models, and into the frequency content of real data.
Finally, the theory in this work is for linear models. This is the right setting to understand first, and is important in its own right. Our theory applies to any linear models that are augmented with a nonlinear basis, such as a Fourier basis, polynomials or radial basis functions. Such bases are particularly useful in the online and streaming settings. Further, a strategy in the continual learning setting is to slowly adjust the features, potentially produced by a neural network, and learn linear predictions that change more quickly on those slowly changing features. The features could be viewed as stationary relative to the linear prediction weights. Extensions to this setting are a natural direction to pursue next.
Practical Implications: One potential conclusion from our results is that though resonance could potentially arise in our learning systems, it may not be problematic. First, the heatmaps in Figures 2a and 2b indicate that SGDm converges much more often than it diverges, and that we can practically avoid resonance simply by reducing the momentum parameter µ. As per Appendix A.2, we find that sufficiently low µ does mitigate resonance. Second, we may simply be able to avoid covariate shift in the first place, by keeping buffers of data and shuffling samples.
However, the reality is not so simple. Higher levels of momentum are typically more effective in practice. There is a tension between the level of momentum and avoiding resonance, that cannot easily be addressed by simply picking low momentum by default. Instead, smarter approaches which recognize potential resonance and then reduce momentum may be a much more effective way to guard against this issue; the development of such methods can leverage the insights in this work. In addition, relying on well-behaved or shuffled inputs may be impractical in certain settings. Algorithms in signal processing and the streaming literature often process data sequentially. Reinforcement learning approaches use sliding window buffers (in replay), where covariate shift may still persist across buffers. Further, it is even possible that adversarial attacks could be designed to exploit this weakness of SGDm, causing poor performance or even instability from attackers attempting to generate data that results in resonant behavior.
Our work also suggests avenues towards the development of diagnostic tools. Indeed, our experimental approach to investigate this phenomena mimics those in the modal analysis literature (Avitabile, 2017). We treat the input sampling process as an input signal which is noisy and/or stochastic, vary the signal’s frequency content across a wide band, and measure the learning system’s output as a response across input frequencies. This is analogous to a mechanical engineer triggering an oscillatory vibration within a machine’s engine compartment, sweeping across frequencies, and measuring the resulting vibration amplitude in the machine’s housing.

heading: 6 CONCLUSION
text: To our knowledge, this work is the first to investigate SGDm under covariate shift from the ODE perspective. This perspective reveals new insights, in particular the revelation that SGDm under covariate shift numerically integrates a parametric oscillator, rather than the simple undriven harmonic oscillator induced by iid sampling. We leverage dynamical systems theory to analyze the stability of a learning system using SGDm under covariate shift, and are able to provide conditions for exponential divergence due entirely to the frequency content induced by this non-iid sampling. Despite its age and establishment, understanding SGDm is an ongoing process. We contribute to this process in a way that provides physical intuition, rather than a merely algebraic proof. We hope that the connection we have drawn to parametric oscillation in dynamical systems theory will help to further propel understanding of non-iid sampling and SGDm.

heading: 7 ETHICS STATEMENT
text: Our work is theoretical and attempts to provide greater understanding of an optimization algorithm. Due to the remoteness of this work from direct ethical concerns, it is difficult to reason about any ethical implications of the work. We hope that the insights into SGDm could help to make machine learning algorithm more robust, but beyond that we are not aware of any ethical concerns or implications of this work.

heading: 8 REPRODUCIBILITY STATEMENT
text: We have endeavored to be very clear and thorough in our explanation of the theoretical results, and all assumptions and proofs are explicitly stated and contained in either the main body of the paper or in the appendix. The code for all of the experiments is included in the supplementary material and the instructions for reproduction are in the README file. The data used is synthetic and code to generate the data is contained in the supplementary material. All of the experimental details (e.g. hyperparameter selection) are contained in tables in the paper and appendix and the supplementary material. All of the experiments should be easily run on modest commodity hardware in the space of a few days, at most.

heading: A APPENDIX

heading: A.1 A PROCEDURE TO NUMERICALLY DETERMINE SYSTEM STABILITY
text: Here we expand on the example given at the end of Section 3. As a reminder, we consider online linear least squares regression using SGDm with a fixed η, µ under a periodic covariate shift. Data is sampled from the sequence of random variables {Xk, Yk}k∈N where Xk is normally distributed and shifts periodically in expectation. Yk is an unchanging linear function of Xk with zero-mean observation noise. For ease of exposition we consider the one dimensional regression problem with Xk ∈ R and Yk ∈ R, however the theory applies to any learning system representable as an LTV ODE satisfying the assumptions required for Theorem 1.
Under this covaiate shift, x̄k varies sinusoidally in [−a, a] with period T = 1f . Note that although Xk is periodic in expectation it admits a stationary distribution. The full specification of the data generating process is as follows
Xk ∼ N (x̄k, 1) x̄k = a cos(2πfk)
Yk = θ ∗ 1Xk + θ ∗ 2 +
∼ N (0, 1)
where θ∗ = [θ∗1 , θ ∗ 2 ] T are the target weights.
We consider a squared error loss function,
L(z; θ) = (ŷ − y)2
= [〈θ,x〉 − (〈θ∗,x〉+ )]2
= [(θ1x+ θ2)− (θ∗1x+ θ∗2 + )]2
Hence, taking gradients at time k w.r.t. θk we have ∇θkL(zk; θk) = [ 2xk[(θk1xk + θk2)− (θ∗k1xk + θ ∗ k2
+ )] 2[(θk1xk + θk2)− (θ∗k1xk + θ ∗ k2 + )] ] Taking expected gradients,
E[∇θkL(zk; θk)] = [ E[2[(θk1x2k + θk2xk)− (θ∗k1x 2 k + θ ∗ k2 xk + xk)]]
E[2[(θk1xk + θk2)− (θ∗k1xk + θ ∗ k2 + )]] ] = [ 2[(θk1(1 + x̄ 2 k) + θk2 x̄k)− (θ∗k1(1 + x̄ 2 k) + θ ∗ k2 x̄k + E[ ]x̄k)]
2[(θk1 x̄k + θk2)− (θ∗k1 x̄k + θ ∗ k2 + E[ ])] ] = 2 [ (θk1(1 + x̄ 2 k) + θk2 x̄k)− (θ∗k1(1 + x̄ 2 k) + θ ∗ k2 x̄k)
(θk1 x̄k + θk2)− (θ∗k1 x̄k + θ ∗ k2 ) ] = 2 [ θk1(1 + x̄ 2 k)− θ∗k1(1 + x̄ 2 k) + θk2 x̄k − θ∗k2 x̄k
θk1 x̄k − θ∗k1 x̄k + θk2 − θ ∗ k2 ] = 2 [ (1 + x̄2k) x̄k x̄k 1 ] [ θk1 − θ∗k1 θk2 − θ∗k2
] In the notation of Proposition 1, we can write this as
gk(θk) = Bk[θk − θ∗]T with Bk = 2 [ (1 + x̄2k) x̄k x̄k 1 ] If we let x̄(t) = a cos(2πfη− 1 2 t), then we induce an X(t) from which {Xk} is sampled, as per Section A.3.3. This results in Bk = B( √ ηk with the matrix-valued function B(t)
B(t) = 2
[ (1 + x̄(t)2) x̄(t)
x̄(t) 1
]
From Proposition 2, we have that our learning system numerically integrates an LTV of the form
θ̈(t) + 1− µ √ η θ̇(t) +B(t)(θ(t)− θ∗) = 0 with (θ(t)− θ∗) = [ θ1(t)− θ∗1(t) θ2(t)− θ2∗(t) ] (5)
Note that B(t) is piecewise continuous and periodic with period T = 1f .
This second order ODE can be represented as a system of first order ODEs with the transformation
ξ1 = θ1 − θ∗1 ξ2 = θ2 − θ∗2 ξ3 = θ̇1
ξ4 = θ̇2
so that the equation equation 5 can be written in standard first order linear form ξ̇ = A(t)ξ: ξ̇1 ξ̇2 ξ̇3 ξ̇4  =  0 0 1 0 0 0 0 1 −2(1 + x̄2(t)) −2x̄(t) µ−1√η 0 −2x̄(t) −2 0 µ−1√η  ξ1ξ2ξ3 ξ4  (6) The matrix B(t) is piecewise continuous and periodic, which implies the matrix A(t) also shares these properties. The ODE equation 5 satisfies all the assumptions needed to make use of Theorem 1 to determine the stability of solution trajectories. We can obtain a fundamental solution matrix, ψ(t), of equation 6 satisfying initial conditions ψ(0) = I4×4 by numerical computation through the use of an ODE solver.
Once we have this matrix ψ(t), we evaluate it at time T = 1f to obtain the system’s monodromy matrix. Now we only need to determine the spectral radius, ρ, of ψ(T ) to determine the stability of solution trajectories for the system. Details regarding why this is true are given in the proof of Theorem 1.
We repeat this process, sweeping over values of µ, η, and f to determine the stability of solution trajectories for systems to triples of fixed values for these hyperparameters. The results of this process are given in Figure 6 below.

heading: A.2 REDUCING RESONANT RESPONSES
text: In experiments 4.1 and 4.2, we see that as the momentum parameter is reduced, the tendency to resonate is mitigated. This aligns with theoretical predictions, which (in the context of those particular experiments) suggest that a sufficiently low µ makes SGDm convergent across all frequencies in the band we evaluate. Intuitively, this aligns with the role µ plays in the ODE, as it appears in the following coefficient on the first order derivative of system equation 3 (i.e. the system’s ‘damping’ coefficient), repeated below with the damping coefficient α explicitly labelled:
θ̈(t) + αθ̇(t) +B(t)(θ(t)− θ∗) = 0 α := 1− µ√ η
(7)
When α = 0, the system has no damping (i.e. friction is zero, in the physical analogue) so resonant responses are maximized. As α increases, damping increases, and resonant responses are reduced. There are two ways to increase α: reduce µ or reduce η.

heading: A.2.1 REDUCING THE MOMENTUM PARAMETER
text: Figure 6b shows the theoretical heatmap across all possible momentum values, and across a much wider band of frequencies, which suggests a trend: setting µ to a sufficiently low value will completely mitigate resonance, though setting µ too low will worsen convergence rate. For now, we will set aside the observation that µ too low worsens convergence rate, and we will now show that reducing µ will reliably dampen resonance across all other experiments.
In 4.5, ADAM’s parameter β1 was varied, and we see that reducing β1 decreases the tendency to resonate. Since β1 is the nearest parameter to µ in SGDm, the desired trend has already been
demonstrated. For the remainder of this section, we show the resonance-damping behaviour of µ in the remaining experiments: 4.3, 4.4, and 4.6. In particular, from each experiment we choose the configuration which had the highest tendency to resonate, and we modify the experiment by running them with several decreasing values of momentum µ. See Figure 7 for results. In Section 4, these experiments used momentum µ = 0.95, and here we run with µ ∈ [0.85, 0.95], with all other experimental parameters identical.

heading: A.2.2 REDUCING THE STEP-SIZE
text: Another way to reduce tendency to resonate is suggested by the damping coefficient equation 7: reducing step-size η. Here we repeat Experiments 4.1 and 4.2, including two smaller step-sizes η to empirically demonstrate that trend. Specifically, Figure 8 shows the resonance heatmap and spectral radius contour lines for regression in two weights with covariate shift, identically to Experiments 4.1 and 4.2. The left column shows sinusoidal covariate shift, and the right column the AR(2) covariate shift. Each row corresponds to a fixed step-size η ∈ {0.01, 0.005, 0.001}, and it is clear that resonant, diverging regions are significantly reduced in size as step-size is decreased, with a narrower band of frequencies diverging, and the minimum momentum µ required for resonance increasing towards 1. This trend is reflected in both the empirical heatmap results, as well as the theoretically predicted spectral radius contour lines.

heading: A.3 PROOFS

heading: A.3.1 ADDITIONAL NOTATION
text: Wherever we write the training input generating process {Xk}, it is implicit that the last dimension is fixed at 1 if one wishes to describe linear models with a bias term. In this way, the notation 〈θ,Xk〉 accommodates linear models with or without a bias term.
We use ξ to denote the phase space coordinates of weights θ, meaning ξ = [ θ
θ̇
] , where vectors
θ, θ̇ ∈ Rd are stacked so that the resulting ξ ∈ R2d. We denote the d-dimensional zero vector as 0d. For an ODE system ξ̇(t) = f(ξ(t), t) with arbitrary solution trajectory ξ(t) approximated by a sequence of iterates {ξk}, we denote the integration time step as h, and the numerical flow as φh,k, which is the map such that ξk+1 = φh,k(ξk). We denote the k-th discrete timestep as tk, which we always use to refer to the k-th multiple of integration time step h, i.e. tk = hk.
Similar to the main body, when a function of time θ(t) is approximated by a discrete time sequence {θk}, they are denoted with the same symbol, and distinguished by parenthetical argument t and subscript k, respectively. At the risk of abusing notation, we will adhere to this convention and use {θ̇k} to be a discrete sequence approximating the function of time θ̇(t), which refers to the time derivative of θ(t). This means the symbol θ̇k is the k-th iterate of a sequence approximating θ̇(t).
We now restate Assumptions, Propositions, and the Theorem from Section 3, providing complete proofs.

heading: A.3.2 LINEAR TIME-VARYING EXPECTED GRADIENT VIA COVARIATE SHIFT
text: Assumption 1. The covariate generating process {Xk}k∈N is not identically distributed: Corr(Xk1 , Xk1) 6= Corr(Xk2 , Xk2) for some k1, k2 ∈ N.
Assumption 2. The targets Yk are a fixed linear function of Xk with iid zero-mean observation noise, k such that E[ k] = 0. That is, Yk = 〈θ∗, Xk〉+ k where θ∗ ∈ Rd is fixed for all k.
Proposition 1. Under Assumptions 1, 2, a linear model with weights θk making predictions Ŷk = 〈θk, Xk〉 with a mean squared error (MSE) objective will induce time-varying linear expected loss gradients gk(θk) = Bk(θk − θ∗) for all k ∈ N, where Bk ∈ Rd×d and Bk1 6= Bk2 for some k1, k2.
Proof. We will show that the gradient of MSE loss ∇θkL(Z; θk) is a random variable whose expectation will take the linear form Bk(θk − θ∗). We start with gradient for arbitrary time step k
∇θkL(Z; θk) = ∇θk(Ŷk − Yk)2 MSE def’n = ∇θk [(〈θk, Xk〉 − Yk)2] Ŷk def’n from Prop. 1 = ∇θk [(〈θk, Xk〉 − 〈θ∗, Xk〉+ k)2] Yk def’n from Ass. 2 = ∇θk [(〈θk − θ∗, Xk〉+ k)2] 〈·, ·〉 distributivity = 2(〈θk − θ∗, Xk〉+ k)∇θk [〈θk − θ∗, Xk〉+ k] chain rule = 2(〈θk − θ∗, Xk〉+ k)Xk = 2〈θk − θ∗, Xk〉Xk + 2 kXk (*)
Taking expectation with respect to distribution Pk:
gk(θk) = EPk [∇θkL(Z; θk)] g def’n = EPk [2〈θk − θ∗, Xk〉Xk + 2 kXk] by equation * = 2EPk [〈θk − θ∗, Xk〉Xk] + 2EPk [ k]EPk [Xk] E linear, k indep. = 2EPk [〈θk − θ∗, Xk〉Xk] E[ k] = 0 = 2EPk [Xk〈θk − θ∗, Xk〉] scalar and vector commute = 2EPk [Xk〈Xk, θk − θ∗〉] 〈·, ·〉 commutative = 2EPk [Xk(XTk (θk − θ∗))] 〈·, ·〉 matrix form = 2EPk [XkXTk ](θk − θ∗) matrix mult. associative = Bk(θk − θ∗)
In the last step, we have defined the matrix Bk as twice the expected outer product of Xk with itself, which is precisely the autocorrelation matrix, and is a function of Xk’s first two moments as follows.
Bk = 2EPk [XkXTk ] = 2Corr(Xk, Xk)
= 2 ( Cov(Xk, Xk) + EPk [Xk]EPk [Xk]T ) So if there exists k1, k2 ∈ N such that the autocorrelation matrices differ
Corr(Xk1 , Xk1) 6= Corr(Xk2 , Xk2)
then Bk1 6= Bk2 . Assumption 1 provides the necessary k1, k2, so the proposition is proved.

heading: A.3.3 ODE CORRESPONDENCE
text: In Section 3 of the main body, we defined B(t), Lipschitz continuous in t, such that Bk = B( √ ηk). If we assume that {Xk} are sampled at intervals of √ η from an underlying continuous-time stochastic processX(t) for which EP (t)[X(t)] and Cov(X(t), X(t)) are Lipschitz continuous in t, this naturally induces a unique Lipschitz continuous B(t). Note P (t) is the distribution at time t, and we can extend equation 2 as follows.
B(t) = 2Cov(X(t), X(t)) + 2EP (t)[X(t)]EP (t)[X(t)]T (8)
This means that {Bk} is sampled from a continuous B(t), with subsequent Bk spaced √ η apart. In essence, step-size η acts as a scaling constant, such that tne unit of t is ≈ 1√η discrete steps of k. (Exact when 1√η is an integer.)
Proposition 2. The SGDm iterates {θk} numerically integrate the ODE system equation 3 with integration step √ η and first order consistency.
θ̈(t) + 1− µ √ η θ̇(t) +B(t)(θ(t)− θ∗) = 0 (3)
Proof. The proof proceeds two parts. First we derive a first order operator-splitting integrator for the ODE equation 3, then we show that it is equivalent to SGDm equation 1.
(Part 1: Operator-Splitting Integrator)
Let ξ : R 7→ R2d be the vector valued function of time whose first d elements are θ(t), and last d elements are θ̇(t):
ξ(t) = [ θ(t)− θ∗ θ̇(t) ] (9)
i.e. ξ is a phase space transformation allowing us to rewrite ODE equation 3 in the form of equation 4, which can then be split into a sum of separate systems f [1], f [2] as follows
ξ̇(t) = f(ξ(t), t) (10)
ξ̇(t) = [ 0d×d Id×d −B(t) − 1−µ√η Id×d ] ξ(t)
ξ̇(t) = [ 0d×d Id×d 0d×d 0d×d ] ξ(t) + [ 0d×d 0d×d −B(t) − 1−µ√η Id×d ] ξ(t)
ξ̇(t) = f [1](ξ(t), t) + f [2](ξ(t), t)
Let h > 0 be an integration time step, φ[1]h,k(ξk) be the implicit Euler numerical flow of f [1](ξ(t), t):
φ [1] h,k(ξk) = ξk + hf [1](φ [1] h,k(ξk), tk+1)
= ξk + h [ 0d×d Id×d 0d×d 0d×d ] φ [1] h,k(ξk) (11)
Let φ[2]h,k(ξk) be explicit Euler numerical flow of f [2](ξ(t), t):
φ [2] h,k(ξk) = ξk + hf [2](ξk, tk)
= ξk + h [ 0d×d 0d×d −B(tk) − 1−µ√η Id×d ] ξk (12)
The composed flow φh,k := φ [1] h,k ◦ φ [2] h,k (13) is a sequentially split operator, which has splitting error order 1 because equation 10 is time-varying linear system (Faragó et al., 2011). The operators being composed, implicit and explicit Euler, are both order 1 consistent with their respective systems (Iserles, 2008). The overall order of consistency of a split operator is the minimum of splitting error, and the orders of the composed flows (Csomós & Faragó, 2008). So we have that equation 13 approximates equation 10 with order 1 consistency.
(Part 2: SGDm Equivalency) We will now show that equation 13 is equivalent to SGDm when integration timestep h = √ η, where η is the SGDm step-size. We start with the definition of ξk+1 as the numerical flow of ξk: ξk+1 := φh,k(ξk)
= φ [1] h,k(φ [2] h,k(ξk)) by equation 13
= φ [2] h,k(ξk) + h [ 0d×d Id×d 0d×d 0d×d ] φ [1] h,k(φ [2] h,k(ξk)) subst’n from equation 11
= φ [2] h,k(ξk) + h [ 0d×d Id×d 0d×d 0d×d ] ξk+1 def’n of ξk+1
= φ [2] h,k(ξk) + h [ θ̇k+1 0d ] matrix mult.
= ( ξk + h [ 0d×d 0d×d −B(tk) − 1−µ√η Id×d ] ξk ) + h [ θ̇k+1 0d ] subst’n from equation 12
= ( ξk + h [ 0d
−B(tk)(θk − θ∗)− 1−µ√η θ̇k
]) + h [ θ̇k+1 0d ] matrix mult.
ξk+1 = ξk + h
[ θ̇k+1
−B(tk)(θk − θ∗)− 1−µ√η θ̇k
] simplification
Recalling that ξk = [ θk − θ∗ θ̇k ] , the above immediately provides the following two recurrence relations:
θk+1 = θk + hθ̇k+1 θ̇k+1 = θ̇k + h ( −B(tk)(θk − θ∗)−
1− µ √ η θ̇k
)
Now let integration time step h be the square root of SGDm step-size, h = √ η, and define vk = √ ηθ̇k. Substituting h, vk, we proceed via elementary algebra:
θk+1 = θk + √ η ( vk+1√ η ) ( vk+1√ η ) = ( vk√ η ) + √ η ( −B(tk)θk − 1− µ √ η ( vk√ η )) θk+1 = θk + vk+1 vk+1 = vk + η ( −B(tk)(θk − θ∗)−
(1− µ)vk η ) vk+1 = vk − ηB(tk)(θk − θ∗)− (1− µ)vk vk+1 = µvk − ηB(tk)(θk − θ∗)
Note that B(tk) = Bk, and that Bk defines the time-varying gradients induced by covariate shift and linear regression to a linear target (Proposition 1). Hence, under those conditions, we have arrived at SGDm equation 1. The proof is complete.
The reader may observe that, given a solution ξ(t) to the system ξ̇(t) = A(t)ξ(t), the proof above shows that the SGDm iterates {θk} are precisely a first order numerical approximation of the first d dimensions of ξ(t) = [ θ(t)− θ∗ θ̇(t) ] , but the remaining d dimensions are approximated only up to a scale factor √ η by iterates {vk}. However, this is not an issue. Solutions θ(t) to the system equation 3 are embedded within the first d dimensions of ξ(t), so the remaining d dimensions and iterates {vk} do not affect the result.

heading: A.3.4 PARAMETRIC RESONANCE FOR ODE CONVERGENCE AND DIVERGENCE
text: Theorem 1. When B(t) is periodic such that B(t) = B(t + T ), the spectral radius ρ of ψ(T ) characterizes the stability of solution trajectories of equation 3 as follows:
• ρ > 1 =⇒ trivial solution θ(t) = θ∗ is unstable. All other solutions diverge as θ(t)→∞ exponentially with rate ρ.
• ρ < 1 =⇒ trivial solution is asymptotically stable, all other solutions converge as θ(t)→ θ∗ exponentially with rate ρ.
Proof. The proof of this theorem relies heavily on the well-established mathematics of Floquet theory and the stability result is contained in Theorem 1.9 and Theorem 1.10 of (Halanay, 1966). Theorem 1.10 states that for a system of differential equations of the form
ξ̇ = A(t)ξ, A(t+ T ) = A(t) (14)
where A(t) is piecewise continuous, the stability of the trivial solution ξ(t) ≡ 0 is determined by the spectral radius of the system’s monodromy matrix M defined below

heading: M = ψ−1(0)ψ(T ) where ψ̇(t) = A(t)ψ(t)
text: (Monodromy Matrix)
Matrix-valued functions ψ(t) are the system’s fundamental solution matrix, and elementary existence results for linear systems allow one to choose a ψ(t) such that ψ(0) = I so that the monodromy matrix simplifies as
M = ψ−1(0)ψ(T )
= I−1ψ(T )
= Iψ(T )
M = ψ(T )
Below we denote the spectral radius of ψ(T ) (and hence of M ) as ρ.
(First Order Linear Form, Stability via Floquet)
Below we show that equation 3 can be transformed to the form equation 14 with A(t) continuous such that trivial solution stability of equation 14 is equivalent to stability of the solution θ(t) = θ∗.
Let ξ(t) be the phase space transformation of θ(t), similarly to the proof of Proposition 2
ξ(t) := [ θ(t)− θ∗ θ̇(t) ] so that for each t, θ(t) ∈ Rd and ξ(t) ∈ R2d. This means equation 3 (restated below)
θ̈(t) + 1− µ √ η θ̇(t) +B(t)(θ(t)− θ∗) = 0
is equivalent to the following first order linear form ξ̇(t) = A(t)ξ(t) where A(t) = [
0d×d Id×d −B(t) − 1−µ√η Id×d ] Since B(t) is periodic and continuous, which is stronger than the requisite piecewise continuity. Since all other submatrices of A(t) are constant, we have that A(t) is also periodic and (piecewise) continuous. Now Theorem 1.10 in (Halanay, 1966) immediately provides the following
• ρ > 1 =⇒ trivial solution ξ(t) = 0 is unstable. All other solutions diverge as ξ(t)→∞ exponentially with rate ρ.
• ρ < 1 =⇒ trivial solution is asymptotically stable, all other solutions converge as ξ(t)→ 0 exponentially with rate ρ.
Since ξ(t) = [ θ(t)− θ∗ θ̇(t) ] :
• The trivial solution ξ(t) = 0 is equivalent to θ(t) = θ∗
• ξ(t)→ 0 is equivalent to θ(t)→ θ∗
• ξ(t)→∞ is equivalent to θ(t)→∞
The theorem is proved.

heading: A.4 DETAILED DESCRIPTIONS OF EXPERIMENTS
text: Below we provide details of the experiments from section 4 in the main body. The tables below contain details about the data generating process, learning algorithm and optimizer, and hyperparameter choices.

Review summaries:
This paper studies the convergence of SGD with momentum under linear regression with covariate shift. Data comes from $Y_i = <\theta^*, X_i> + \epsilon_i$ for a fixed $\theta^*$ and iid zero-mean noise $\epsilon_i$, but the distribution of $X_i$ varies over time ("covariate shift"), leading to non-iid samples. In this setting, the paper shows that the expected progress of SGD with momentum is a discretization of a second-order ODE. Then, the paper proceeds to show that if the covariate shift is periodic, the convergence/divergence of the ODE can be determined by looking at the spectral radius of a matrix called the monodromy matrix. The paper moves on to show experiments that test the theoretical characterizations and also investigate if different settings under relaxed assumptions show similar "resonance" behavior.  The keyword "resonance-driven divergence" in this paper can be understood as follows: divergence of SGD with momentum happens only at some specific frequency of covariate shift. Metaphorically, this is similar to breaking a wine glass with a sound wave tuned to the right frequency. 

This paper analyses the phenomenon of resonance in momentum SGD (SGDm) under a time-dependent covariate shift. This setting is useful to depart from iid-sampling-based theory, and could have useful implications for continual learning and reinforcement learning. In the paper, it is shown that such a setting corresponds to a _parametric_ oscillator of the parameters, which explains instability and divergence in a non-iid SGDm setting. The authors then test their hypothesis empirically, first in a simple setting that respects all the assumptions of the theorem, and then progressively getting rid of more and more assumptions, which ends up in finding a similar (albeit heavily dampened) phenomenon in non-linear ReLU-MLPs using Adam. 

This paper studies the behavior of SGD with momentum when the training data are not iid. The authors consider a linear regression problem with covariance drift. In this case, they show that SGDm is the Euler's method for a second order ODE with a changing coefficient. Under the assumption of periodic covariance of the data, the authors identify a resonance phenomenon between the oscillation of data covariance and the oscillation of SGDm iterator. SGDm is proven to diverge when resonance happens. Though, when there is no resonance the authors are unable to prove the convergence of SGDm.   Beside the theory, numerical experiments are conducted for cases covered and not covered by the theory. Numerical results show a dependence of the optimizer performance and the covariance pattern in much broader settings than linear regression, even including neural networks. 

The paper discusses conditions under which SGD with momentum (SGDm) would diverge. Specifically, training samples are assumed to be non-iid, connected through oscillations in first or second moments. For the specific Gaussian setting of linear regression with oscillatory covariate shift, a driven (with time-dependent parameters) linear oscillator ODE is formulated for the parameters under SGDm. This allows to derive conditions, related to the learning and momentum rates, under which parameters would con-/diverge. Theoretical predictions are tested empirically, and it is shown that the basic phenomenon (resonance or at least suboptimal convergence) is also present in setups with noise, non-harmonic oscillations, other optimizers (Adam), and nonlinear regression models (NN).