# Path differences between quasistatic and fatigue cracks in anisotropic media

Xinyuan Zhai  IMSIA, CNRS, EDF, ENSTA Paris, Institut Polytechnique de Paris,
91120 Palaiseau, France  Thomas Corre  IMSIA, CNRS, EDF, ENSTA Paris, Institut
Polytechnique de Paris, 91120 Palaiseau, France  Nantes Université, Ecole
Centrale Nantes, CNRS, GeM, UMR 6183, F-44000 Nantes, France  Ataollah
Mesgarnejad  Physics Department and Center for Interdisciplinary Research on
Complex Systems, Northeastern University, Boston, Massachusetts 02115, USA
Alain Karma  Physics Department and Center for Interdisciplinary Research on
Complex Systems, Northeastern University, Boston, Massachusetts 02115, USA
Véronique Lazarus<EMAIL_ADDRESS>IMSIA, CNRS, EDF, ENSTA
Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France

###### Abstract

The propagation path of quasistatic cracks under monotonic loading is known to
be strongly influenced by the anisotropy of the fracture energy in crystalline
solids or engineered materials with a regular microstructure. Such cracks
generally follow directions close to minima of the fracture energy. Here we
demonstrate both experimentally and computationally that fatigue cracks under
cyclic loading follow dramatically different paths that are predominantly
dictated by the symmetry of the loading with the microstructure playing a
negligible or subdominant role.

††preprint: APS/123-QED

Over and above the fascination with the shapes it creates [1], understanding
how and when cracks form is of broad importance [2], for instance in the
design of manufacturing processes [3, 4, 5], in geophysics for a better
understanding of earthquakes [6, 7], in mechanical engineering with obvious
societal issues of safety [8, 9]. Cracks or sharp defects are known to
particularly affect the resistance to fracture of any material as they lead to
stress concentration in their vicinity [10]. The theoretical framework [11,
12, 13] and computational tools to predict brittle crack paths in diverse
settings, including complex three-dimensional geometries [14, 15], are by now
well-developed. However, what criterion should be used to predict crack paths
is still not fundamentally understood in the practical case of fatigue crack
growth occurring progressively under a repetitive cyclic load. Such kind of
solicitations are common for industrial components (e.g. turbines, rotors,
bridges), as well as in earth and life sciences (e.g. waves, tidal ebb and
flow, weather or climate cycles, walking/running steps).

To date, it has been commonly assumed that fatigue cracks follow the same path
as quasistatic cracks that propagate infinitesimally slowly (_i.e.,_ at a
velocity much smaller than the sound speed) under a monotonic, as opposed to
fatigue cyclic load [16, 15]. This assumption is based on the qualitative
observation that fatigue and quasistatic cracks follow a similar path in
isotropic materials, _i.e.,_ materials with both isotropic elastic properties
and a fracture energy independent of the direction of the crack path. In such
materials, the principal of local symmetry (PLS) has been traditionally used
to predict quasistatic crack paths [11] for 2D plane strain or plane stress
within the traditional framework of LEFM; 3D out of plane mode III loading
causes crack front segmentation, but the propagation of segmented fronts
remains consistent with the PLS [15]. The PLS assumes that the stress
distribution remains symmetrical about the instantaneous crack axis during
propagation and corresponds to a vanishing mode II stress intensity factor
($K_{II}=0$). It stems from the fact that cracks with non-vanishing $K_{II}$
kink by a discontinuous change of propagation direction. Therefore, cracks
following curvilinear paths, with a continuous change of direction, should
propagate with $K_{II}=0$ [12]. Such crack paths also maximize the energy
release rate, _i.e.,_ , the rate at which the stored elastic energy in the
material is released taking into account the energy cost $\Gamma$ (per area)
of creating new fracture surfaces. The PLS and Maximum Energy Release Rate
(MERR) criteria predict slightly different kink angles [17] for a crack
initially loaded with a finite $K_{II}$, but the same curvilinear paths with
$K_{II}=0$.

For anisotropic materials, the PLS is not applicable but extensions of the
MERR criterion have been used to predict crack paths in various contexts [18,
19, 13, 20, 21]. A MERR criterion for anisotropic media was formally derived
within the diffuse interface phase-field approach to brittle fracture [19, 13]
by taking into account the anisotropy of the fracture energy quantified by the
function, $\Gamma(\theta)$, where $\theta$ is the angle between the
instantaneous crack propagation direction and some fixed reference axis. The
selected crack path was shown to be the one that maximizes with respect to
$\theta$ the energy released $G(\theta)-\Gamma(\theta)$, where $G(\theta)$
quantifies the elastic strain energy available by unit surface to propagate
the crack. Consequently, the crack is predicted to propagate with a finite
$K_{II}$ proportional to $d\Gamma(\theta)/d\theta$ in anisotropic media, and
thus $K_{II}=0$ in isotropic media consistent with the PLS.

The MERR criterion has been validated for quasistatic brittle fracture by
phase-field simulations [13, 21] as well as experiments using materials with a
composite microstructure engineered to have nearly isotropic elastic
properties but a strongly anisotropic fracture energy [21]. However, it is not
expected to apply to fatigue cracks that propagate under cyclic loads that
peak at $G(\theta)$ magnitudes well below the Griffith threshold
$\Gamma(\theta)$; fatigue cracks do not obey the energy balance condition
$G(\theta)=\Gamma(\theta)$, which is valid at each instant of propagation of a
quasistatic crack. Fatigue crack growth is more commonly described by a
phenomenological Paris-type law $da/dN=C\Delta K_{I}^{n}$, which relates the
crack advance per cycle to the variation $\Delta K_{I}$ of mode I stress
intensity factor during each cycle. For fatigue crack growth in isotropic
media, the PLS can be rationalized in the same way as for quasistatic fracture
by remarking that cracks with $K_{II}\neq 0$ kink, and that therefore smooth
cracks should propagate with $K_{II}=0$. However, there is no reason to expect
that the MERR criterion with $K_{II}\sim d\Gamma(\theta)/d\theta$ is
applicable to fatigue crack paths. Hence, it remains largely unknown how
fatigue and quasistatic crack paths differ from each other in anisotropic
media and, in the absence of a MERR criterion, what physically controls the
former.

In this Letter, we address those questions by comparing quasistatic and
fatigue crack paths both experimentally and computationally in two model
materials where the elastic properties are isotropic but the fracture energy
is strongly anisotropic. We demonstrate that in both cases quasistatic and
fatigue cracks follow strikingly different paths. The experiments use
3D-printed samples with a criss-cross microstructure yielding a fracture
energy $\Gamma(\theta)$ with a strong four-fold anisotropy. As illustrated in
Fig. 1, in samples with the same printed structure and loading in pure tensile
mode I or with the addition of a shear mode II component, quasistatic cracks
follow non-vanishing $K_{II}$ directions corresponding to sharp local minima
of the fracture energy, as expected from MERR. In contrast, somewhat
surprisingly, fatigue cracks follow vanishing $K_{II}$ directions despite the
strong anisotropy of the underlying microstructure. The phase-field
computations use a different two-fold anisotropic form of $\Gamma(\theta)$,
which was previously used to validate the MERR criterion for quasistatic crack
paths in anisotropic media [13, 21]. They yield similar results as the
experiments with the difference that fatigue cracks follow paths where
$K_{II}$ is much smaller than along quasistatic crack paths, but non-
vanishing. We describe the experiments and computations in what follows and
further interpret our findings.

Figure 1:  Different crack propagation paths observed under monotonic load
(left column) and cyclic load (right column) in polycarbonate CT (c-f) or CTS
(g, h) samples printed by FDM. Under increasing imposed displacement (a), the
crack propagates quasistatically in the weak plane of the microstructure (c,
e, g) until brutal failure in another direction. Under cyclic load (b),
propagation occurs along a different path (d, f, h) during all the fatigue
cycles. Bifurcation in the weak plane occurs an increasing load is applied to
completely break the specimen. (Samples geometry: 50 $\times$ 48 $\times$ 8
mm3. c,d: $\varphi=30\degree$; e, f: $\varphi=45\degree$;g, h:
$\varphi=45\degree$, $\beta=15\degree$)

To address the question experimentally, polycarbonate Compact Tensile (CT) and
Compact Tensile Shear samples (CTS) are printed by Fused Deposition Modeling
(FDM) following a 100% filling criss-cross pattern: each successive layer is
oriented at 90$\degree$ to the previous one, with a $\varphi$ inclination to
the initial crack Fig. 1). Tensile tests on samples with different filling
angles, showed that the material is nearly isotropic with Young modulus $E\sim
2000$ MPa and Poisson’ ratio $\nu=0.34$ [22]. Once printed, a sharp crack is
introduced by fatigue [23] and black-and-white speckle is applied on the
samples prior to testing for Digital Image Correlation (DIC). To propagate the
crack, the samples are loaded in a tensile machine either by increasing
monotonically the displacement at a constant rate for the quasistatic
propagation case or by applying a cyclic force load for fatigue (Fig. 1a,b).
While the load is mirror-symmetric to the crack plane for CT samples
corresponding to pure mode I, a non-zero $\beta$ tilt for CTS samples induces
non-symmetrical mode I+II loading.

Figure 1 highlights the difference between quasistatic and fatigue propagation
on a set of representative cases. Applying increasing displacement, the force-
displacement curve has a mountain shape (Fig. 1a). Before the peak, no crack
propagation occurs and the force displacement curve is quasi-linear, which is
consistent with linear elasticity. After the peak, the crack propagates quasi-
statically and the force drops slowly. Fracture occurs in the $\varphi$
direction (Fig. 1c, e, g) for both CT and CTS. At the microstructure scale, it
corresponds to separating two parallel threads every other layer, and to the
orthogonal breaking of the threads on the remaining layers. Applying cyclic
loading (Fig. 1b), the crack progressively advances with each loading cycle.
The crack advances straight ahead for CT samples (Fig. 1d, f), while a kink is
observed for CTS (Fig. 1h) but in a different direction to that obtained under
monotonic loading (Fig. 1g). At the microstructure scale, the crack does not
generally follow the orientation of the threads. We have performed experiments
for different values of $\beta$ and $\varphi$, with good reproducibility and
similar differences in the propagation paths. For each case, the crack is
oriented by the microstructure for monotonic loading and by the loading
orientation for fatigue.

Figure 2:  Evolution of the stress intensity factors during crack propagation
measured by DIC. They are measured along the path indicated by a white segment
in the insets. $K_{II}\neq 0$ during quasistatic propagation (a, b).
$K_{II}=0$ in both CT and CTS during fatigue propagation (c, d)

To find out more about the underlying crack propagation rules, and in
particular if the PLS is retrieved in fatigue, the stress intensity factors
($K_{I}$, $K_{II}$) are measured by DIC [24, 22] using the open-source
software UFreckels [25, 26]. Figure 2 shows their evolution during crack
propagation for the experiments of Fig. 1 (c, d, g, h). For monotonic load and
CT samples (Fig. 2a), only mode I builds up before propagation while the
applied load is increased; for CTS samples (Fig. 2b), mode I and II
contributions increase. This is due to the symmetry and non-symmetry of the CT
and CTS setups respectively. At the onset of propagation, the quasistatic
crack changes direction and a jump in the value of $K_{I}$ and $K_{II}$ is
observed. This jump is consistent with LEFM [17] as demonstrated in [22]. As
$K_{II}\neq 0$, the PLS is not verified but the crack propagation direction
can be predicted by a MERR criterion as in previous studies of quasistatic
crack propagation in anisotropic media [18, 19, 13, 20, 21, 23, 27]. For
cyclic load (Fig. 2c, d), the crack propagates as soon as cycling begins. The
stress intensity factors evolution is the same as that expected for an
isotropic material: (i) $K_{I}$ increases for fixed applied force from the
increase of the crack length [28], (ii) $K_{II}=0$ as soon as the crack
propagates in agreement with the Principle of Local Symmetry. In brief,
contrary to the quasistatic case, fatigue propagation follows the propagation
rules of isotropic material and is not affected by the anisotropy of the
microstructure.

To illustrate computationally the difference between quasistatic and fatigue
crack paths in anisotropic media, we used the two-fold anisotropy [13, 21]
$\Gamma(\theta)=\Gamma_{max}\sqrt{{\cal
A}^{-2}\cos^{2}(\theta-\phi)+\sin^{2}(\theta-\phi)}$ where ${\cal
A}=\Gamma_{max}/\Gamma_{0}$ measures the toughness ratio for propagation
perpendicular to ($\Gamma_{max}$) and parallel to ($\Gamma_{0}$) a weak plane
analogous to one of the thread axis in the experiments. This form causes
quasistatic cracks to follow directions close to weak planes for sufficiently
strong anisotropy [13, 21]. Quasistatic propagation was modeled using the
well-established alternate minimization scheme applied to anisotropic media in
[21]. Fatigue crack growth was investigated using a recently developed phase-
field model of fatigue [29] that enables crack growth below the Griffith
threshold ($G<\Gamma$) by degradation of the fracture toughness [30]. This
formulation has reproduced salient features of fatigue crack growth, such as
the Paris law and non-trivial 2D and 3D crack paths in isotropic media [29].
We simulate compact specimen of width $W$ [31]. Details of the numerical
methods are given in the supplementary material [27]. The loading is applied
via a vertical displacement $u_{y}$ on the top half of the upper hole and by
setting $u_{x}=u_{y}=0$ around the lower hole. In this geometry, the energy
release rate decreases for displacement-controlled loading as the crack
propagates. Therefore, in both quasi-static and fatigue simulations, we set
the displacement $u_{y}=k\ell$ where
$\ell=2\sqrt{\Gamma_{max}W/3E}={3E}K_{IC}^{max}\sqrt{3W}/2$ and slowly
increase $k$ such that quasi-static crack propagation remains energetically
stable and fatigue cracks are able to propagate an extended distance. For
quasi-static propagation, we vary $k$ from 0 to 4. The onset of crack
propagation occurs at $k\approx 1.6$ and $k\approx 2.6$ for $\phi=30^{\circ}$
and $\phi=90^{\circ}$, respectively, followed by quasi-static crack
propagation when $k$ is increased above these thresholds. For fatigue, we set
$k=1.5$ for $\phi=90^{\circ}$ below its quasi-static limit. We increase $k$
from 1 to $2.25$ for $\phi=30^{\circ}$, which ensures that the energy release
rate at the crack tip remains at all times above the minimum threshold
$G_{min}$ for fatigue crack growth but below the Griffith threshold $G=\Gamma$
for quasi-static propagation. Further details of the phase-field methods can
be found in [21] for quasi-static and [32] for fatigue crack propagation.

The results are shown in Fig. 3 which compares quasistatic and fatigue crack
paths for a large anisotropy $\mathcal{A}=4$ and two different weak plane
orientations. When the weak plane is perpendicular to the horizontal axis
(left panel of Fig. 3 for $\phi=90^{\circ}$), the quasistatic crack kinks at a
large angle to follow a direction close to the weak plane with $K_{II}\neq 0$.
In contrast, the fatigue crack propagates straight with $K_{II}=0$. When the
weak plane is at a $30^{\circ}$ angle with respect to the horizontal axis
(right panel of Fig. 3 for $\phi=30^{\circ}$), the quasistatic crack kinks to
follow again a direction close to the weak plane, while the fatigue crack is
much more weakly defected to propagate in a direction with $K_{II}$ much
smaller than the quasistatic crack, but non-vanishing.

From a general symmetry viewpoint, we would expect fatigue cracks to follow a
direction that is influenced both by the symmetry of the loading, which
dictates a preferential axis for the cyclic damage process (_e.g.,_ the
principal stress axis), and the symmetry of the underlying microstructure of
the material, which offers the possibility for damage to occur preferentially
in certain directions (_e.g.,_ parallel to the thread axis). The experiments
demonstrate that the symmetry of the loading plays a dominant role for the
printed materials investigated here. Under pure tensile loading, the stress
fields are mirror-symmetric about the horizontal propagation axis and fatigue
cracks propagate along this axis. They follow the PLS irrespective of the
orientation of the microstructure (_i.e.,_ $\varphi=45^{\circ}$ and
$\varphi=30^{\circ}$ yield identical fatigue crack paths with $K_{II}\approx
0$, Fig. 1d, f ). In contrast, in the computations, the symmetry of the
microstructure plays a subdominant, albeit not entirely negligible, role. For
$\phi=90^{\circ}$, both the loading and the microstructure (_i.e.,_
$\Gamma(\theta)$) are mirror-symmetric about the horizontal axis and the
fatigue crack propagates straight with $K_{II}=0$. For $\phi=30^{\circ}$, the
microstructure is no longer mirror-symmetric about the horizontal axis,
thereby causing the fatigue crack to propagate at a small angle. Further
studies are warranted to elucidate the relative importance of loading and
microstructure in dictating fatigue crack paths in different anisotropic
materials. The present experiments and computations clearly demonstrate that
loading plays a dominant role in two model cases where quasistatic crack paths
are largely dictated by the microstructure.

Figure 3:  Comparison of computed quasi-static (gray) and fatigue (red) crack
paths with the weak plane at (a) $\phi=90^{\circ}$ (b) $\phi=30^{\circ}$ with
respect to the horizontal propagation axis.

The research at Northeastern University was supported by the United States
NSF-CMMT grant 1827343. The research at ENSTA Paris was partially supported by
Agence de l’Innovation de Défense – AID - via Centre Interdisciplinaire
d’Etudes pour la Défense et la Sécurité – CIEDS (AID-2018 60 0071 00 470 75 01
and FracAddi project 2022) and by IP Paris Doctoral Allocation. The authors
thanks N. Thurieau, L. Cherfa for their help in the experiments.

For the purpose of Open Access, a CC-BY public copyright licence
(https://creativecommons.org/licenses/by/4.0/) has been applied by the authors
to the present document and will be applied to all subsequent versions up to
the Author Accepted Manuscript arising from this submission.

## References

  * Hull [1999] D. Hull,  _Fractography_ , edited by C. U. Press (Press syndicate of the university of Cambridge, UK, 1999).
  * Cotterell [2010] B. Cotterell,  _Fracture And Life_ (Imperial College Press, 2010).
  * Nam  _et al._ [2012] K. H. Nam, I. H. Park, and S. H. Ko, Nature 485, 221 (2012).
  * Lazarus [2017] V. Lazarus, EPL (Europhysics Letters) 117, 24002 (2017).
  * Ngo  _et al._ [2018] T. D. Ngo, A. Kashani, G. Imbalzano, K. T. Nguyen, and D. Hui, Composites Part B: Engineering 143, 172 (2018).
  * Fliss  _et al._ [2005] S. Fliss, H. S. Bhat, R. Dmowska, and J. R. Rice, Journal of Geophysical Research: Solid Earth 110, 10.1029/2004JB003368 (2005).
  * Rubino  _et al._ [2022] V. Rubino, N. Lapusta, and A. Rosakis, Nature 606, 922 (2022).
  * Anderson [2005] T. L. Anderson,  _Fracture Mechanics_ (CRC Press, Boca Raton, FL, 2005).
  * Besson  _et al._ [2023] J. Besson, J. Bleyer, S. Feld-Payet, A.-F. Gourgues-Lorenzon, F. Hannard, T. Helfer, J. Hure, D. Kondo, V. Lazarus, C. Le Bourlot, H. Maitournam, C. Maurini, N. Moes, T. Morgeneyer, L. Morin, T. Petit, and A. Simar,  _MEALOR II Damage Mechanics and Local Approach to Fracture_ (Zenodo, 2023).
  * Williams [1952] M. L. Williams, Journal of Applied Mechanics 19, 526 (1952).
  * Goldstein and Salganik [1974] R. V. Goldstein and R. L. Salganik, International Journal of Fracture 10, 507 (1974).
  * Cotterell and Rice [1980] B. Cotterell and J. R. Rice, International Journal of Fracture 16, 155 (1980).
  * Hakim and Karma [2009] V. Hakim and A. Karma, Journal of the Mechanics and Physics of Solids 57, 342 (2009).
  * Bourdin  _et al._ [2014] B. Bourdin, J.-J. Marigo, C. Maurini, and P. Sicsic, Physical Review Letters 112, 014301 (2014).
  * Chen  _et al._ [2015] C.-H. Chen, T. Cambonie, V. Lazarus, M. Nicoli, A. J. Pons, and A. Karma, Physical Review Letters 115, 265503 (2015).
  * Buchholz  _et al._ [2004] F.-G. Buchholz, A. Chergui, and H. A. Richard, Engineering Fracture Mechanics 71, 455 (2004).
  * Amestoy and Leblond [1992] M. Amestoy and J.-B. Leblond, International Journal of Solids and Structures 29, 465 (1992).
  * Marder [2004] M. Marder, Europhysics Letters 66, 364 (2004), publisher: IOP Publishing.
  * Hakim and Karma [2005] V. Hakim and A. Karma, Physical Review Letters 95, 235501 (2005), arXiv: cond-mat/0412775 Publisher: American Physical Society.
  * Takei  _et al._ [2013] A. Takei, B. Roman, J. Bico, E. Hamm, and F. Melo, Phys. Rev. Lett. 110, 144301 (2013).
  * Mesgarnejad  _et al._ [2020] A. Mesgarnejad, C. Pan, R. M. Erb, S. J. Shefelbine, and A. Karma, Physical Review E 102, 013004 (2020).
  * Corre and Lazarus [2021] T. Corre and V. Lazarus, International Journal of Fracture 230, 19 (2021).
  * Zhai [2023] X. Zhai,  _Crack propagation in elastic media with anisotropic fracture toughness : experiments and numerical modeling_ , Ph.D. thesis, Institut Polytechnique de Paris (2023).
  * Hamam  _et al._ [2007] R. Hamam, F. Hild, and S. Roux, Strain 43, 181 (2007).
  * Réthoré [2015] J. Réthoré, International Journal for Numerical Methods in Engineering 103, 516 (2015).
  * Réthoré [2018] J. Réthoré, UFreckles v2.0 , 10.5281/zenodo.1433776 (2018).
  * [27] See Supplemental Material at [URL will be inserted by publisher] for Complements about the MERR criterion and Details of the numerical methods.
  * Tada  _et al._ [2000] H. Tada, P. C. Paris, and G. R. Irwin,  _The Stress Analysis of Cracks Handbook_ , edited by T. A. S. of Mechanical Engineers (Professional Engineering Publishing, 2000).
  * Grossman-Ponemon  _et al._ [2022a] B. E. Grossman-Ponemon, A. Mesgarnejad, and A. Karma, Engineering Fracture Mechanics 264, 108255 (2022a).
  * Mesgarnejad  _et al._ [2019] A. Mesgarnejad, A. Imanian, and A. Karma, Theoretical and Applied Fracture Mechanics 103, 102282 (2019).
  * [31] ASTM,  _Standard Test Method for 1 Measurement of Fracture Toughness_ , Tech. Rep. E 1820-05 (ASTM).
  * Grossman-Ponemon  _et al._ [2022b] B. E. Grossman-Ponemon, A. Mesgarnejad, and A. Karma, Engineering Fracture Mechanics 264, 108255 (2022b).


# Computational EEG in Personalized Medicine:
A study in Parkinson’s Disease

Sebastian Mathias Keller  University of Basel  Maxim Samarin  University of
Basel  Antonia Meyer  Hospital of the University of Basel  Vitalii Kosak
(Cozak)  University of Basel  Hospital of the University of Basel  Ute
Gschwandtner  Hospital of the University of Basel  Peter Fuhr  Hospital of the
University of Basel  Volker Roth  University of Basel

###### Abstract

Recordings of electrical brain activity carry information about a person’s
cognitive health. For recording EEG signals, a very common setting is for a
subject to be at rest with its eyes closed. Analysis of these recordings often
involve a dimensionality reduction step in which electrodes are grouped into
10 or more regions (depending on the number of electrodes available). Then an
average over each group is taken which serves as a feature in subsequent
evaluation. Currently, the most prominent features used in clinical practice
are based on spectral power densities. In our work we consider a simplified
grouping of electrodes into two regions only. In addition to spectral features
we introduce a secondary, non-redundant view on brain activity through the
lens of Tsallis Entropy $S_{q=2}$. We further take EEG measurements not only
in an eyes closed (ec) but also in an eyes open (eo) state. For our cohort of
healthy controls (HC) and individuals suffering from Parkinson’s disease (PD),
the question we are asking is the following: How well can one discriminate
between HC and PD within this simplified, binary grouping? This question is
motivated by the commercial availability of inexpensive and easy to use
portable EEG devices. If enough information is retained in this binary
grouping, then such simple devices could potentially be used as personal
monitoring tools, as standard screening tools by general practitioners or as
digital biomarkers for easy long term monitoring during neurological studies.

##  1 Introduction

The recording of the first human EEG was performed in 1924 by German physician
Hans Berger. Since then relative band power has become an established measure
to quantify deviations from normal oscillatory brain activity: Several minutes
of EEG signal, usually recorded under resting state condition and in an ‘eyes
closed’ (ec) setting, are filtered into four to seven non-overlapping
frequency bands covering a range from 0.5Hz up to 70Hz. Then, for each of
those bands, the relative signal power is calculated. Based on this signal
power, multiple studies have established a statistical connection between the
neurological condition of their subjects (as measured e.g. by psychological
testing) and the distribution of spectral power within these bands [1], [2],
[3].
High density EEG machines available today provide up to 256 individual
electrodes recording brain activity at a frequency of 1000Hz or higher. If one
is not interested in performing source reconstruction or spatial filtering
such a high number of electrodes is usually not necessary, especially as
electrodes in proximity to each other are very highly correlated. It is thus
expected that features based on a grouping of electrodes into 10 or more
regions as shown in fig. 4(b) are only slightly more informative as a grouping
into only 2 such regions as shown in fig. 4(a). In fact a color coded
representation of the mean Tsallis entropy $\langle\hat{S}_{q=2}\rangle$ per
electrode, as shown in fig. 1, indicates that the simplest viable partitioning
is obtained by grouping electrodes into a frontal and an occipital region. If
these two regions still carry enough information about a persons cognitive
health it is potentially possible that two electrodes located at the front and
back of the head would capture enough information to present a simplified
alternative to a much denser electrode setting. Consequentially, as portable,
longterm EEG devices with up to five electrodes are already available, these
devices might prove useful as future long term monitoring tools analog to how
ECG devices are used today to monitor heart rate variability.

Figure 1: Estimated Tsallis entropy during eo state for the 4-8Hz band
(averaged over several minutes) is shown for every electrode. Most but not all
subjects of our cohort fall into this “front versus back” pattern where
entropy is higher in the frontal brain regions. As subjects transition into ec
state, frontal entropy is significantly reduced (not shown). CT068 and CT086
are healthy controls whereas CT016 and CT040 have been diagnosed with
Parkinson’s disease.

###  1.1 Parkinson’s Disease and Dementia

According to the American Psychiatric Association [4], Dementia is a syndrome
that consists of a decline in cognitive and intellectual abilities occurring
in an awake and alert patient. As a symptom complex it is known to occur in
over 70 disorders [5], among which Parkinson’s Disease (PD) is the second most
common age-related neurodegenerative disorder after Alzheimer’s disease. It is
estimated to affect nearly 2 percent of those over age 65. According to the
National Parkinson Foundation (NPF), over the entire course of their illness,
50 to 80 percent of those suffering from Parkinson’s disease eventually
develop dementia. The average time from onset of movement problems to the
development of dementia is about 10 years. Accurately diagnosing Parkinson’s
Disease – especially in its early stages – requires experienced practitioners.
Thus providing tools for detecting early changes in brain activity that are as
easy to use as taking ones own blood pressure is an important first step.

Figure 2: A significant change in brain oscillation can be observed when
transitioning from an eyes open to an eyes closed state and vice versa. This
change is known as the “Berger effect” or “$\alpha$-blockade” as the amplitude
of EEG oscillations in the $\alpha$-band (8-13 Hz) decreases significantly
when subjects open their eyes. Visible are also the three eye blinks during
eyes open state.

##  2 Materials and Methods

In this study, we focus on how to complement and improve existing signal power
based classification methods by (i) introducing additional entropy based
features, (ii) including eo and ec states into the measurement process and
(iii) using the Berger effect (see fig. 2) as a potential feature indicative
of abnormal brain activity. At the same time we simplify the grouping of
electrodes and investigate whether a measurement of frontal and occipital
brain activity might potentially suffice to indicate a deviation from healthy
brain oscillations.

###  2.1 Sample Description

We recruited participants for this study from the outpatient clinic for
movement disorders of the Department of Neurology and Neurophysiology of the
University Hospital of Basel (Universitätsspital Basel), Switzerland in the
period from 2011 to 2012. We selected 46 persons with Parkinson’s disease and
labelled them as “PD”. Parkinson’s disease was diagnosed according to the
United Kingdom Parkinson’s Disease Society Brain Bank criteria. The patients
who had dementia (Diagnostic and Statistical Manual of Mental Disorders, 4th
Edition), history of stroke, epilepsy, multiple sclerosis and surgical
interventions to the brain, or/and insufficient knowledge of German language,
were excluded. We selected 16 persons without Parkinson’s disease, matched to
patients according to sex, age, and education, and labelled them as “HC”. For
controls, inclusion criteria were a subjective report of good health and a
neuropsychological examination within normal limits. Exclusion criteria were
identical to those for patients. The local ethics committee approved the study
protocol. All participants provided written informed consent.

###  2.2 Tsallis Entropy $\boldsymbol{S_{q=2}}$

The q-entropy, also known as the Tsallis entropy, was proposed in [6] and
takes the following form

$S_{q}(p_{i})=\frac{1-\sum_{i=1}^{W}p_{i}^{q}}{q-1}$ (1)

with probabilities $p_{i}$ and $q\in\mathbb{R}$. In the limit $q\rightarrow
1$, eq. 1 recovers the standard Shannon entropy. In analogy to the Normal
distribution, which is the maximum entropy distribution of the Shannon entropy
within the class of distributions with $E(x)=\mu$ and
$E((x-\mu)^{2})=\sigma^{2}$, Tsallis entropy $S_{q}$ is maximized by a
q-Gaussian distribution. The pdf $p_{qg}$ of a q-Gaussian is given by

$p_{qg}(x)=p_{0}\left[1-(1-q)\frac{x}{x_{0}}^{2}\right]^{1/(1-q)},$ (2)

for details see section 3 of [7]. Fig. 5(b) shows the shape of such
q-Gaussians for four different values of q where $q=1$ is the Normal
distribution. Higher values of q lead to an increased probability in the tails
of the distribution. Random samplings from a heavy tailed Gaussian with $q=2$
are shown in fig. 5(a) and compared to samples obtained for $q=1$. When band
filtering a Gaussian white noise signal and a q-Gaussian ($q=2$) white noise
signal we find that the signal obtained from the q-Gaussians has more
similarities to a band filtered EEG signal for the same band in terms of
amplitude variance than its Gaussian counterpart. This leads to the conclusion
that Tsallis entropy might be a more appropriate entropy measure for EEG
signals. This view is supported by several studies achieving improved
classification results in similar settings based on Tsallis entropy, e.g. [8],
[9], [10].

For estimating q-Tsallis entropy for $q=2$ an estimator $\hat{S}_{q=2}$ has
been proposed in [11]. It is parameter-free, easy to calculate, designed to
estimate the entropy of time series and given by

$\hat{S}_{q=2}=1-\frac{\frac{1}{N}\sum s_{i}^{2}}{\sigma^{2}}$ (3)

with $s_{i}^{2}=\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\mu)^{2}$. To estimate the
entropy of a time series it is first divided into $N$ bins where each extreme
point of the signal marks the beginning, resp. the end, of a bin. The variance
within each bin is given by $s_{i}^{2}$ whereas $\sigma^{2}$ is the variance
of the entire time series.

(a) HC: ec $\rightarrow$ eo: for HC the occipital region has a higher increase
in entropy in the eo state than for PD

(b) PD: ec $\rightarrow$ eo: for PD the increase in occipital entropy content
is less pronounced

(c) HC and PD group: the front to back ratio of entropy in the 4Hz-8Hz band is
higher for healthy controls

Figure 3:  4Hz-8Hz band: HC vs. PD in ec and eo state. In each panel the stars
mark the means of the resp. groups. In panel (c) the crosses mark the extremes
for HC and PD (during eo state) calculated using Archetype analysis [12]. They
turn out to be the two most important features in table 1, row 4. Proximity to
either extreme is an indicator of a person’s neurological status, i.e. PD or
HC.

##  3 Results and Discussion

Table 1: Binary classification with repeated cross validation (80/20) using gradient boosting machines (gbm) on a 2-region frontal vs. occipital setting with target PD vs. HC. BP is band power, TS is Tsallis Entropy. Measure is area under the ROC curve. Due to the imbalance between HC and PD group weighted classification and SMOTE were used but no significant difference in ROC was obtained. Classification (gbm) |  |  |
---|---|---|---
Feature Set | Ratios | Modality | five most important Features (descending order) | ROC
BP | no | EC | bp.b.ec.4-8, bp.b.ec.13-45, bp.b.ec.1-4, bp.f.ec.1-4, bp.f.ec.4-8 | $64.0\%$
BP | no | EO | bp.b.eo.1-4, bp.b.eo.4-8, bp.b.eo.13-45, bp.f.eo.13-45, bp.f.eo.4-8 | $74.8\%$
TS | no | EC | ts.f.ec.1-4, ts.f.ec.4-8, ts.b.ec.8-13, ts.b.ec.13-45, ts.b.ec.4-8 | $65.5\%$
TS | no | EO | ts.b.eo.4-8, ts.f.eo.4-8, ts.b.eo.13-45, ts.f.eo.8-13, ts.b.eo.8-13 | $71.8\%$
BP + TS | no | EC + EO | bp.b.eo.1-4, bp.f.ec.8-13, ts.f.eo.4-8, ts.b.ec.8-13, ts.b.eo.1-4 | $\mathbf{75.7\%}$
Ratios(BP, TS) | yes | EC + EO | bp.feo.vs.beo.1-4, ts.fec.vs.bec.8-13, bp.feo.vs.beo.4-8, | $67.0\%$
|  |  | ts.bec.vs.beo.13-45, bp.fec.bs.bec.1-4 |
BP + TS | yes | EC + EO | bp.b.eo.1-4, bp.feo.vs.beo.1-4, bp.feo.vs.beo.4-8, | $71.5\%$
|  |  | ts.fec.vs.bec.8-13, ts.f.eo.4-8 |

EEG measurements in ec and eo states were filtered into four frequency bands
(1-4Hz, 4-8Hz, 8-13Hz and 13-45Hz). For each band Tsallis entropy (ts) for
$q=2$ and relative band power (bp) was calculated for the frontal (f) and
occipital (b) region. PD vs. HC classification was performed with repeated
cross validation on a train/test split of 80/20 using gradient boosting
machines (GBM). Due to class imbalance the loss function was weighted
accordingly and classification was measured using ROC. Our best result in the
binary grouping setting is achieved with a combination of ec and eo
measurements using entropy and band power features. Using only the ratios of
band power and entropy in front and back gives only a moderate classification
of ROC=67%. Thus the absolute value of measured quantities in back and front
are meaningful – not only individual ratios. Combining ratios and absolute
values of power and entropy does not improve overall classification: the
possible benefit of additional features is outweighed by the larger search
space. The largest increase in ROC was obtained when switching from “eyes
closed” to “eyes open” measurements. We thus recommend to make it a standard
procedure to record EEG measurements during both ec and especially eo state.
Our best classification result compares favourably – given the reduced
dimensionality – with results based on 10 region groupings presented in [13],
where ROC=80% is achieved on a similar patient cohort also from the University
Hospital Basel.

#### Acknowledgment

This project is supported by the the Swiss National Science Foundation project
CR32I2-159682

## References

  * Dauwels et al. [2011] Justin Dauwels, François-Benoît Vialatte, and Andrzej Cichocki.  On the early diagnosis of alzheimer’s disease from eeg signals: A mini-review, 2011.
  * Bonanni et al. [2008] Laura Bonanni, Astrid Thomas, Pietro Tiraboschi, Bernardo Perfetti, Sara Varanese, and Marco Onofrj.  Eeg comparisons in early alzheimer’s disease, dementia with lewy bodies and parkinson’s disease with dementia patients with a 2-year follow-up.  _Brain : a journal of neurology_ , 131(Pt 3):690–705, Mar 2008.  doi: 10.1093/brain/awm322.
  * Hardmeier et al. [2014] Martin Hardmeier, Florian Hatz, Habib Bousleiman, Christian Schindler, Cornelis Jan Stam, and Peter Fuhr.  Reproducibility of functional connectivity and graph measures based on the phase lag index (pli) and weighted phase lag index (wpli) derived from high resolution eeg.  _PloS one_ , 9(10):e108648, 2014.
  * Ame [1994] _Diagnostic and Statistical Manual of Mental Disorders: DSM-IV_.  American Psychiatric Association, Washington, DC, 4 edition, 1994.  Committee on Nomenclature and Statistics.
  * D.S. and P.J [1997] Geldmacher D.S. and Whitehouse P.J.  Differential diagnosis of alzheimer’s disease.  _Neurology 48(suppl 6):S2-S9_ , pages 1125–1138, 1997.
  * Tsallis [1988] Constantino Tsallis.  Possible generalization of boltzmann-gibbs statistics.  _Journal of Statistical Physics_ , 52(1):479–487, Jul 1988.  ISSN 1572-9613.  doi: 10.1007/BF01016429.
  * Picoli et al. [2009] Jr. S. Picoli, R. S. Mendes, L. C. Malacarne, and R. P. B. Santos.  q-distributions in complex systems: a brief review.  _Brazilian Journal of Physics_ , 39(2a):468–474, 2009.  doi: 10.1590/S0103-97332009000400023.
  * McBride et al. [2013] J. McBride, X. Zhao, T. Nichols, V. Vagnini, N. Munro, D. Berry, and Y. Jiang.  Scalp eeg-based discrimination of cognitive deficits after traumatic brain injury using event-related tsallis entropy analysis.  _IEEE Transactions on Biomedical Engineering_ , 60(1):90–96, Jan 2013.  ISSN 0018-9294.  doi: 10.1109/TBME.2012.2223698.
  * Bock et al. [2010] T. J. De Bock, S. Das, M. Mohsin, N. B. Munro, L. M. Hively, Y. Jiang, C. D. Smith, D. R. Wekstein, G. A. Jicha, A. Lawson, J. Lianekhammy, E. Walsh, S. Kiser, and C. Black.  Early detection of alzheimer’s disease using nonlinear analysis of eeg via tsallis entropy.  In _2010 Biomedical Sciences and Engineering Conference_ , pages 1–4, May 2010.  doi: 10.1109/BSEC.2010.5510813.
  * Capurro et al. [1999] Alberto Capurro, Luis Diambra, Daisy Lorenzo, Macadar O, M.T. Martin, Mostaccio C, A Plastino, Perez J, Rozman E, Maria Eugenia Torres, and Velluti J.  Human brain dynamics: The analysis of eeg signals with tsallis information measure.  _Physica A Statistical and Theoretical Physics_ , 265, 1999.  doi: 10.1016/S0378-4371.
  * Sneddon [2007] Robert Sneddon.  The tsallis entropy of natural information.  _Physica A: Statistical Mechanics and its Applications_ , 386(1):101–118, 2007.  URL https://EconPapers.repec.org/RePEc:eee:phsmap:v:386:y:2007:i:1:p:101-118.
  * Cutler and Breiman [1994] Adele Cutler and Leo Breiman.  Archetypal analysis.  _Technometrics_ , 36(4):338–347, 1994.
  * Chaturvedi et al. [2017] M. Chaturvedi, F. Hatz, and U. Gschwandtner.  Quantitative eeg (qeeg) measures differentiate parkinson’s disease (pd) patients from healthy controls (hc).  _Frontiers in aging neuroscience_ , 9(3), 2017.  doi: 10.3389/fnagi.2017.00003.

##  Appendix A Appendices

###  A.1 Simplified Partitioning: 10 regions vs. 2 regions

(a) Electrodes partitioned into frontal and occipital
region

(b) A standard partitioning into 10 brain regions as
used in a clinical setting

Figure 4: The observed division into frontal and occipital region shown in
fig. 1 motivates the simplification of the standard partitioning with its 10
regions to a binary “front vs. back” partitioning of electrodes.

###  A.2 Family of q-Gaussians

(a) Random samples from a Normal distri-
bution and from a heavy tailed q-Gaussian.

(b) Within the family of q-Gaussian, the Normal distribution is recovered for
$q=1$. Any q-Gaussian with $q>1$ has tails with higher probabilities.

Figure 5: In the family of q-Gaussians a distribution with lighter or heavier
tails than the Normal distribution ($q=1$) is obtained for $q<1$ or $q>1$.


# Stochastic Learning Rate Optimization in the Stochastic Approximation and
Online Learning Settings

Theodoros Mamalis
Department of Electrical and Computer Engineering
University of Illinois at Urbana-Champaign
306 N Wright St, Urbana, IL 61801, USA
<EMAIL_ADDRESS>
Dušan Stipanović
Coordinated Science Laboratory
University of Illinois at Urbana-Champaign
1308 W Main St, Urbana, IL 61801, USA
<EMAIL_ADDRESS>
Petros Voulgaris
Department of Mechanical Engineering
University of Nevada
Reno, NV 89557, USA
<EMAIL_ADDRESS>

###### Abstract

In this work, multiplicative stochasticity is applied to the learning rate of
stochastic optimization algorithms, giving rise to stochastic learning-rate
schemes. In-expectation theoretical convergence results of Stochastic Gradient
Descent equipped with this novel stochastic learning rate scheme under the
stochastic setting, as well as convergence results under the online
optimization settings are provided. Empirical results consider the case of an
adaptively uniformly distributed multiplicative stochasticity and include not
only Stochastic Gradient Descent, but also other popular algorithms equipped
with a stochastic learning rate. They demonstrate noticeable optimization
performance gains, with respect to their deterministic-learning-rate versions.

00footnotetext: This work is supported by the grant from the National Robotics
Initiative grant titled “NRI: FND:COLLAB: Multi-Vehicle Systems for Collecting
Shadow-Free Imagery in Precision Agriculture” (grant no.2019-04791/project
accession no. 1020285) from the USDA National Institute of Food and
Agriculture.

##  1 Introduction

Machine learning models today range from simple linear regression to deep
neural networks. These models try to capture the underlying behavior of
systems which produce outputs $y$ when given inputs $x$. The performance of
these models is usually measured via a loss function $\ell(h(x;\theta),y)$
where $h(.;.)$ is the predicted output of the model for an input $x$ for given
model parameters $\theta\in\mathbb{R}^{d}$. Then, if the joint distribution
$P$ of $x$ and $y$ is known, the estimates $\theta$ are the minimizers of the
expected loss ${\mathbb{E}_{x,y\sim P}[\ell(h(x;\theta),y)]}$. Nonetheless,
complete knowledge of $P$ is absent in most practical scenarios, and a limited
set of training data is available instead. In that case, an estimate of the
expected loss is minimized, known as the empirical loss function, giving rise
to the empirical loss minimization (ERM) problem:

$\mathop{\min}\limits_{\theta\in{R^{d}}}\frac{1}{n}\sum\nolimits_{i=1}^{n}{\ell(h(x_{i};\theta),y_{i})},$
(1)

where the sequence $\\{(x_{i},y_{i})\\}_{i=1}^{n}$ denotes the available
training data.

Then, if $\xi$ represents a random sample from the the available training
data, let the loss function as a composition of $\ell(.,.)$ and $h(.;.)$ be
denoted by $L(\theta;\xi)$. Also let
$f(\theta):={\mathbb{E}_{\xi}[L(\theta,\xi)]}$. This yields a stochastic
optimization setting where the computed gradient $g(\theta,\xi)$ is a function
of the true gradient $\nabla f(\theta)$. If the former is assumed to be an
unbiased estimator of the latter as in this work, then:

$\mathbb{E}_{\xi_{t}}[g({\theta_{t}},{\xi_{t}})]=\nabla f(\theta_{t}),$ (2)

where $\xi_{t}$ and $\theta_{t}$ are independent. Due to its simplicity, cost-
effectiveness and well studied properties, one of the most commonly used
algorithms for ERM is stochastic gradient descent (SGD) originally developed
in Robbins and Monro (1951), which iteratively gives an estimate of $\theta$
in (1) via the rule:

${\theta_{t+1}}={\theta_{t}}-{a_{t}}g_{t},$ (3)

for $t=1,2,\cdots$ and where $g_{t}$ is short for $g({\theta_{t}},{\xi_{t}})$.
This rule will be referred to as original SGD or simply SGD throughout the
paper.

Even though commonly the learning rate $a_{t}$ is taken to be deterministic,
the novelty in this work is that it introduces the setting where the learning
rate in SGD becomes stochastic, by equipping it with multiplicative
stochasticity. That is, the learning rate now becomes $a_{t}=\eta_{t}u_{t}$
(where e.g., $\eta_{t}=a/\sqrt{t}$ or $\eta_{t}=a$, with $a$ a positive
constant; with some abuse of terminology, $\eta_{t}$ will be referred to as
step size) and $u_{t}$ is a random variable, referred to as stochastic factor
(SF) in the rest of the paper. A Multiplicative-Stochastic-Learning-Rate is
referred to as MSLR.

In general, it has been observed that minimization performance can be
noticeably improved under appropriate learning schemes e.g., under
deterministic adaptive learning rate schemes. Examples include ADAM Kingma and
Ba (2015) and some variants (e.g., AMSGrad Reddi et al. (2018) or ADAMW
Loshchilov and Hutter (2019)) or precursors (e.g. McMahan and Streeter (2010),
Duchi et al. (2011)). This work demonstrates that performance can be
significantly improved if instead of just adaptive, the learning rate becomes
both adaptive and stochastic.

Convergence analysis in expectation but for deterministic learning rates has
been studied in Bertsekas and Tsitsiklis (2000); Nemirovski et al. (2009);
Moulines and Bach (2011); Ghadimi and Lan (2013); Bottou et al. (2018). This
work provides convergence results in-expectation for SGD using stochastic
learning rate as well as results under the online learning framework
(introduced in Zinkevich (2003)). It shows that with stochastic-learning-rates
optimization performance is significantly improved, while providing similar
technical results to the deterministic-learning-rate case. Convergence
analysis in the almost-surely setting using stochastic learning rate schemes,
and in specific MSLR, for algorithms including SGD has been done in Mamalis et
al. (2021).

As discussed, both of the stochastic optimization and online learning settings
will be investigated. This is because the main theme of the former is
smoothness of the function $f(\theta)$, and not necessarily convexity, and the
main theme of the latter is convexity and not necessarily smoothness. This
provides results which characterize the convergence, and behaviour in general,
of the algorithm under the assumptions of each setting.

Moreover, the analysis will allow for any SF distribution with bounded
support. This includes resetting SF distributions that yield resetting
learning rate schemes, which have been shown to improve optimization
performance (O’Donoghue and Candes (2012); Loshchilov and Hutter (2017)). A
Resetting-MSLR scheme will be referred to as RMSLR. If, moreover, its SF is
uniformly distributed, the resetting learning rate scheme will be referred to
as RUMSLR (e.g. as in Fig. 1).

The remainder of the paper is organized as follows. Section 2 provides the
stochastic optimization setting of this work and provides the in-expectation
convergence results. The online learning setting along with the respective
convergence results is given in Section 3. Section 4 provides a discussion on
RUMSLR and its particular SF along with a discussion on selecting the SF
hyperparameters. In Section 5, the experimental results for various
optimization algorithms using this RUMSLR scheme are presented and compared to
the deterministic-learning rate case. Section 6 concludes the paper with paths
for future work.

Figure 1: Plot that illustrates the stochastic factor for an RUMSLR scheme
with SF
$u_{t}\sim\mathcal{U}(\sqrt[\varepsilon_{t}]{{{c_{1}}}},\sqrt[\varepsilon_{t}]{{{c_{2}}}})$
where $\varepsilon_{t}=t\,\,(\bmod\beta)+1$, and $c_{1}=0.8,c_{2}=1.2$. In
this case, the value of the resetting parameter is $\beta=100$, i.e. the
distribution of $u_{t}$ resets every $100$ iterations to
$\mathcal{U}(0.8,1.2)$.

##  2 Stochastic Optimization Setting

This section will provide convergence results for stochastic-learning-rate SGD
under the stochastic optimization setting, briefly discussed in Section 1.
These results hold for both MSLR and RMSLR schemes. The following assumptions
will be used in the theorems. Firstly, function $f(\theta)$ is $L$-smooth,
i.e.:

$f(\theta)\leq f(\theta^{\prime})+\nabla
f{(\theta^{\prime})^{T}}(\theta-\theta^{\prime})+\frac{L}{2}{\left\|{\theta-\theta^{\prime}}\right\|^{2}},$
(4)

for all $\theta^{\prime},\theta\in\mathbb{R}^{d}$. Furthermore, for all
$\theta\in\mathbb{R}^{d}$, it satisfies:

${\left\|{\nabla f(\theta)}\right\|^{2}}\geq 2c(f(\theta)-{f_{min}}),$ (5)

where ${f_{min}}$ is the global minimum of $f(\theta)$, and $0<c\leq L$. It is
noted that this assumption only requires the function to attain a global
minimum (more detailed discussion in Karimi et al. (2016)). Also:

${\mathbb{E}_{{\xi_{t}}}}[{\left\|{g({\theta_{t}},{\xi_{t}})}\right\|^{2}}]\leq{M}+{M_{G}}{\left\|{\nabla
f({\theta_{t}})}\right\|^{2}}\,\,\text{ for }\,\,t\in\mathbb{N},$ (6)

for some nonnegative scalars $M$ and $M_{G}\geq 1$.

The theorem that follows provides a bound on the optimality gap of the SGD
algorithm when using a stochastic-learning-rate with a constant step size:

######  Theorem 1.

Assume (4-6) are satisfied. Assume a bounded stochasticity factor $c_{1}\leq
u_{t}\leq c_{2}$. Furthermore, assume a positive learning rate of:

$a_{t}=\frac{c_{1}}{c_{2}^{2}}au_{t},$ (7)

and:

$a<\frac{c_{2}}{c_{1}}\frac{1}{LM_{G}},$ (8)

then the iterates of the SGD algorithm in (3) satisfy:

$\mathop{\lim}\limits_{t\to\infty}{\mathbb{E}}[f({\theta_{t+1}})-{f_{min}}]\leq\frac{{aLM}}{{2c}}.$
(9)

The guarantee in (9) is the same as when SGD uses a deterministic-learning-
rate, e.g., as given in Bottou et al. (2018). The step size is
$\frac{c_{1}}{c_{2}^{2}}$ times smaller than in the deterministic-learning-
rate case. This is pragmatically nonrestrictive for convergence speed as
supported by experimental results, since $c_{1}$ and $c_{2}$ are adequately
close in practice (see discussion in Section 4).

###### Proof.

Taking expectations with respect to $\xi_{t}$ in (4), and from (3)
substituting $\theta_{t+1}$ for $\theta$ and $\theta_{t}$ for
$\theta^{\prime}$, and in (7) letting
$\eta:=\frac{\sqrt{c_{1}}}{c_{2}\sqrt{c_{2}}}a$ (thus $a_{t}=\eta u_{t}$)
yields:

$\displaystyle{\mathbb{E}_{{\xi_{t}}}}[f({\theta_{t+1}})]$ $\displaystyle\leq
f({\theta_{t}})-{a_{t}}{\nabla^{T}}f({\theta_{t}})\mathbb{E}_{{\xi_{t}}}[g({\theta_{t}},{\xi_{t}})]+a_{t}^{2}\frac{L}{2}{\mathbb{E}_{{\xi_{t}}}}\left[{{{\left\|{g({\theta_{t}},{\xi_{t}})}\right\|}^{2}}}\right]$
(10) $\displaystyle\leq f({\theta_{t}})-{a_{t}}{\left\|{\nabla
f({\theta_{t}})}\right\|^{2}}+a_{t}^{2}\frac{L}{2}{\mathbb{E}_{{\xi_{t}}}}\left[{{{\left\|{g({\theta_{t}},{\xi_{t}})}\right\|}^{2}}}\right]$
$\displaystyle\leq
f({\theta_{t}})+{a_{t}}({a_{t}}\frac{{{M_{G}}L}}{2}-1){\left\|{\nabla
f({\theta_{t}})}\right\|^{2}}+a_{t}^{2}\frac{L}{2}M$ $\displaystyle\leq
f({\theta_{t}})+\frac{1}{2}{a_{t}}({a_{t}}{M_{G}}L-2){\left\|{\nabla
f({\theta_{t}})}\right\|^{2}}+a_{t}^{2}\frac{L}{2}M$ $\displaystyle\leq
f({\theta_{t}})+\frac{1}{2}{a_{t}}({\eta}c_{2}{M_{G}}L-2){\left\|{\nabla
f({\theta_{t}})}\right\|^{2}}+a_{t}^{2}\frac{L}{2}M$ $\displaystyle\leq
f({\theta_{t}})-\frac{1}{2}{a_{t}}{\left\|{\nabla
f({\theta_{t}})}\right\|^{2}}+a_{t}^{2}\frac{L}{2}M$ $\displaystyle\leq
f({\theta_{t}})-c{a_{t}}(f({\theta_{t}})-{f_{min}})+a_{t}^{2}\frac{L}{2}M.$

The second inequality followed from (2), the third from (6), the fifth from
$a_{t}=\eta u_{t}\leq\eta c_{2}$, the sixth from (8) which gives that
$ac_{2}M_{G}L-2\leq-1$, and the seventh from (5).

Then, using $\eta c_{1}\leq a_{t}\leq\eta c_{2}$ in (10), taking expectations
with respect to all $\xi_{i}$ up to time $t$ and denoting the result as
$\mathbb{E}_{\xi}[\cdot]$:

$\displaystyle{\mathbb{E}_{\xi}}[f({\theta_{t+1}})]\leq{\mathbb{E}_{\xi}}[f({\theta_{t}})]-c\eta
c_{1}({\mathbb{E}_{\xi}}[f({\theta_{t}})]-{f_{min}})+c_{2}^{2}\eta^{2}\frac{{LM}}{2},$
(11)

or equivalently:

$\displaystyle{\mathbb{E}_{\xi}}[f({\theta_{t+1}})]$
$\displaystyle\leq{\mathbb{E}_{\xi}}[f({\theta_{t}})]-c\eta
c_{1}({\mathbb{E}_{\xi}}[f({\theta_{t}})]-{f_{min}})$ (12)
$\displaystyle+\eta^{2}c_{2}^{2}\frac{{LM}}{2}+\eta
c_{2}^{2}\frac{{LM}}{{2cc_{1}}}-\eta
c_{2}^{2}\frac{{LM}}{{2cc_{1}}}+{f_{min}}-{f_{min}}.$

This further results in:

$\displaystyle{\mathbb{E}}[f({\theta_{t+1}})-{f_{min}}]-\eta\frac{{LMc_{2}^{2}}}{{2cc_{1}}}\leq(1-c{{\eta}}c_{1})({\mathbb{E}}[f({\theta_{t}})-{f_{min}}]-{{\eta}}\frac{{LMc_{2}^{2}}}{{2cc_{1}}}).$
(13)

By unrolling the inequality:

$\displaystyle{\mathbb{E}}[f({\theta_{t+1}})-{f_{min}}]-\eta\frac{{LMc_{2}^{2}}}{{2cc_{1}}}\leq(1-c{{\eta}}c_{1})^{t}({\mathbb{E}}[f({\theta_{1}})-{f_{min}}]-{{\eta}}\frac{{LMc_{2}^{2}}}{{2cc_{1}}}).$
(14)

Furthermore, from (8) it follows that:

$a<\sqrt{\frac{c_{2}}{c_{1}}}\frac{1}{LM_{G}}\Leftrightarrow\eta=\frac{a\sqrt{c_{1}}}{c_{2}\sqrt{c_{2}}}<\frac{1}{LM_{G}c_{2}},$
(15)

which for $M_{G}\geq 1$, $c\leq L$ yields:

$c\eta c_{1}\leq ca_{t}\leq c\eta
c_{2}<\frac{cc_{2}}{{c_{2}{M_{G}}L}}\leq\frac{c}{{L}}\leq 1.$ (16)

Therefore from (14) it is that:

$\mathop{\lim}\limits_{t\to\infty}{\mathbb{E}}[f({\theta_{t+1}})-{f_{min}}]-\eta\frac{{LMc_{2}^{2}}}{{2cc_{1}}}\leq
0,$ (17)

or equivalently:

$\mathop{\lim}\limits_{t\to\infty}{\mathbb{E}}[f({\theta_{t+1}})-{f_{min}}]\leq\eta\frac{{LMc_{2}^{2}}}{{2cc_{1}}}.$
(18)

Substituting back $\frac{c_{1}}{c_{2}^{2}}a$ for $\eta$ proves the theorem. ∎

Next, a convergence theorem for diminishing step size is provided.

######  Theorem 2.

Assume (4-6) are satisfied. Assume a bounded stochasticity factor $c_{1}\leq
u_{t}\leq c_{2}$. Furthermore, assume a positive learning rate of
$a_{t}=\eta_{t}u_{t},$ Then, if $\eta_{t}$ is decreasing,
$\sum\nolimits_{t=1}^{\infty}{\eta_{t}^{2}}<\infty$ and:

$\eta_{1}<\frac{1}{LM_{G}c_{2}},$ (19)

the iterates of the SGD algorithm in (3) satisfy:

$\mathop{\lim}\limits_{t\to\infty}{\mathbb{E}}[f({\theta_{t+1}})]={f_{min}}.$
(20)

The guarantee in (20) provides the same convergence result as in the
deterministic-learning-rate case, it holds for arbitrary decreasing step sizes
and does not depend on the initial point unlike previous works (cf. Bottou et
al. (2018)).

###### Proof.

The proof for diminishing step size is identical to the proof of Theorem 1
until equation (10), except that the fifth inequality in (10) now follows from
$a_{t}=\eta_{t}u_{t}\leq\eta_{1}c_{2}$, and the sixth from (19) which gives
$\eta_{1}c_{2}M_{G}L-2\leq-1$ and therefore $\eta_{t}c_{2}M_{G}L-2\leq-1$,
since $\eta_{t}\leq\eta_{1}$ for $t\in\mathbb{N}$. Then, taking expectations
with respect to all $\xi_{i}$ up to time $t$:

$\displaystyle{\mathbb{E}_{\xi}}[f({\theta_{t+1}})]\leq{\mathbb{E}_{\xi}}[f({\theta_{t}})]-c{a_{t}}({\mathbb{E}_{\xi}}[f({\theta_{t}})]-{f_{min}})+a_{t}^{2}\frac{{LM}}{2},$
(21)

or:

$\displaystyle{\mathbb{E}_{\xi}}[f({\theta_{t+1}})]$
$\displaystyle\leq{\mathbb{E}_{\xi}}[f({\theta_{t}})]-c{a_{t}}({\mathbb{E}_{\xi}}[f({\theta_{t}})]-{f_{min}})+a_{t}^{2}\frac{{LM}}{2}+a_{t}\frac{{LM}}{{2c}}-a_{t}\frac{{LM}}{{2c}}$
(22)
$\displaystyle+a_{t+1}\frac{{LM}}{{2c}}-a_{t+1}\frac{{LM}}{{2c}}+{f_{min}}-{f_{min}}.$

Taking expectations with respect to all learning rate random variables
$a_{t}$, denoting the result $\mathbb{E}_{a}[\cdot]$ and also denoting
$\mathbb{E}_{a}[\mathbb{E}_{\xi}[\cdot]]$ as $\mathbb{E}[\cdot]$, gives:

$\displaystyle{\mathbb{E}}[f({\theta_{t+1}})]$
$\displaystyle\leq{\mathbb{E}}[f({\theta_{t}})]-c{\mathbb{E}_{a}[a_{t}]}({\mathbb{E}}[f({\theta_{t}})]-{f_{min}})+{\mathbb{E}_{a}[a_{t}^{2}]}\frac{{LM}}{2}+{\mathbb{E}_{a}[a_{t}]}\frac{{LM}}{{2c}}-{\mathbb{E}_{a}[a_{t}]}\frac{{LM}}{{2c}}$
(23)
$\displaystyle+{\mathbb{E}_{a}[a_{t+1}]}\frac{{LM}}{{2c}}-{\mathbb{E}_{a}[a_{t+1}]}\frac{{LM}}{{2c}}+{f_{min}}-{f_{min}}.$

Then using
$\mathbb{V}_{a}[a_{t}]={\mathbb{E}_{a}[a^{2}_{t}]}-{\mathbb{E}_{a}[a_{t}]}^{2}$
to replace ${\mathbb{E}_{a}[a^{2}_{t}]}$, and then rearranging terms:

$\displaystyle{\mathbb{E}}[f({\theta_{t+1}})-{f_{min}}]-{\mathbb{E}_{a}[a_{t+1}]}\frac{{LM}}{{2c}}$
$\displaystyle\leq(1-c{\mathbb{E}_{a}[{a_{t}}]})({\mathbb{E}}[f({\theta_{t}})-{f_{min}}]-{\mathbb{E}_{a}[{a_{t}}]}\frac{{LM}}{{2c}})$
(24)
$\displaystyle+({\mathbb{E}_{a}[a_{t}]}-{\mathbb{E}_{a}[a_{t+1}]})\frac{{LM}}{{2c}}+\mathbb{V}_{a}[{a_{t}}]\frac{{LM}}{2}.$

By unrolling the sequence it is:

$\displaystyle\mathbb{E}[f({\theta_{t+1}})-{f_{min}}]-{\mathbb{E}_{a}}[a_{t+1}]\frac{{LM}}{{2c}}$
$\displaystyle\leq\prod\limits_{j=1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{j}}])}({\mathbb{E}}[f({\theta_{1}})-{f_{min}}]-{\mathbb{E}_{a}[a_{1}]}\frac{{LM}}{{2c}})$
(25)
$\displaystyle+\frac{{LM}}{{2c}}\sum\limits_{j=1}^{t}{({\mathbb{E}_{a}[a_{j}]}-{\mathbb{E}_{a}[a_{j+1}]})}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}$
$\displaystyle+\frac{{LM}}{{2}}\sum\limits_{j=1}^{t}{\mathbb{V}_{a}[{a_{j}}]}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}.$

Furthermore, from (19), $M_{G}\geq 1$, and $c\leq L$:

$ca_{t}\leq ca_{1}\leq
c\eta_{1}c_{2}<\frac{cc_{2}}{{c_{2}{M_{G}}L}}\leq\frac{c}{{L}}\leq 1,$ (26)

it is $0\leq 1-c\mathbb{E}[a_{t}]<1$, since $c\mathbb{E}[a_{t}]$ is positive.
Moreover, $c\eta_{1}c_{1}\leq c\eta_{1}c_{2}<1$. Let
$\delta:=1-c\eta_{1}c_{1}$. Then, it is shown that the last two terms in (25)
converge. This is since:

$\displaystyle\sum\limits_{j=1}^{t}{({\mathbb{E}_{a}[a_{j}]}-{\mathbb{E}_{a}[a_{j+1}]})}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}\leq{(\eta_{1}c_{2}-\eta_{1}c_{1})}\sum\limits_{j=1}^{t}\prod\limits_{n=j+1}^{t}{(1-c\eta_{1}c_{1})}$
$\displaystyle=\eta_{1}(c_{2}-c_{1})\sum\limits_{j=1}^{t}(1-c\eta_{1}c_{1})^{t-j-1}=\eta_{1}(c_{2}-c_{1})\frac{{1-{\delta^{t}}}}{{(1-\delta)\delta}},$

which converges to $\frac{\eta_{1}(c_{2}-c_{1})}{(1-\delta)\delta}$ as
$t\to\infty$. For the other term:

$\displaystyle\sum\limits_{j=1}^{t}{\mathbb{V}_{a}[{a_{j}}]}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}$
(27)
$\displaystyle\leq\sum\limits_{j=1}^{t}\frac{{{\eta_{j}^{2}}}}{{4}}{{{(\sqrt[j]{{{c_{2}}}}-\sqrt[j]{{{c_{1}}}})}^{2}}}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}$
$\displaystyle\leq\frac{{{\eta_{1}^{2}{{{({c_{2}}-{c_{1}})}^{2}}}}}}{{4}}\sum\limits_{j=1}^{t}\prod\limits_{n=j+1}^{t}{(1-c\eta_{1}c_{1})}$
$\displaystyle\leq\frac{{{\eta_{1}^{2}{{{({c_{2}}-{c_{1}})}^{2}}}}}}{{4}}\sum\limits_{j=1}^{t}\delta^{t-j-1}$
$\displaystyle\leq\frac{{{\eta_{1}^{2}{{({c_{2}}-{c_{1}})}^{2}}}}}{{4}}\frac{1-\delta^{t}}{(1-\delta)\delta},$

which converges to
$\frac{{{\eta_{1}^{2}}{{({c_{2}}-{c_{1}})}^{2}}}}{{4(1-\delta)\delta}}$ as
$t\to\infty$, and where the second inequality followed from the variance-bound
of bounded random variables. Then let the algorithm run for long enough, so
that the last two terms in (14) are close enough to their limits so that for
their tails after some timestep $K$ the following holds:

$\displaystyle\sum\limits_{j=K+1}^{t}{({\mathbb{E}_{a}[a_{j}]}-{\mathbb{E}_{a}[a_{j+1}]})}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}$
$\displaystyle\approx 0$ (28)
$\displaystyle\sum\limits_{j=K+1}^{t}{\mathbb{V}_{a}[{a_{j}}]}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}$
$\displaystyle\approx 0.$

Then, partitioning the last two terms in (27) at $K$, and using (28), it is
that:

$\displaystyle\mathbb{E}[f({\theta_{t+1}})-{f_{min}}]-{\mathbb{E}_{a}}[a_{t+1}]\frac{{LM}}{{2c}}$
$\displaystyle\leq\prod\limits_{j=1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{j}}])}({\mathbb{E}}[f({\theta_{1}})-{f_{min}}]-{\mathbb{E}_{a}[a_{1}]}\frac{{LM}}{{2c}})$
(29)
$\displaystyle+\frac{{LM}}{{2c}}\sum\limits_{j=1}^{K}{({\mathbb{E}_{a}[a_{j}]}-{\mathbb{E}_{a}[a_{j+1}]})}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}$
$\displaystyle+\frac{{LM}}{{2}}\sum\limits_{j=1}^{K}{\mathbb{V}_{a}[{a_{j}}]}\prod\limits_{n=j+1}^{t}{(1-c{\mathbb{E}_{a}}[{a_{n}}])}.$

Therefore, from (29), and from
${\mathbb{E}_{a}[a_{t+1}]}=\eta_{t+1}{\mathbb{E}[u_{t+1}]}$ it follows that:

$\mathop{\lim}\limits_{t\to\infty}\left({\mathbb{E}}[f({\theta_{t+1}})-{f_{min}}]-\eta_{t+1}{\mathbb{E}[u_{t+1}]}\frac{{LM}}{{2c}}\right)\leq
0,$ (30)

which means:

$\mathop{\lim}\limits_{t\to\infty}{\mathbb{E}}[f({\theta_{t+1}})-{f_{min}}]=0,$
(31)

since
$\mathop{\lim}\limits_{t\to\infty}\eta_{t+1}{\mathbb{E}[u_{t+1}]}\frac{{LM}}{{2c}}=0$
from $c_{1}\leq{\mathbb{E}[u_{t+1}]}\leq c_{2}$ and $\eta_{t+1}\to 0$ as
$t\to\infty$. ∎

Concluding this section, from Theorems 1-2, it is thus observed that for the
stochastic optimization setting the optimality gap bounds for SGD using
stochastic-learning-rate do not differ from the ones when using a
deterministic-learning-rate.

##  3 Online Learning Setting

Now, an alternative analysis for stochastic-learning-rate SGD which holds for
both MSLR and RMSLR is given by considering the following optimization
problem, referred to as online optimization:

${R_{T}}=\sum\limits_{t=1}^{T}{{f_{t}}({\theta_{t}})}-\mathop{\min}\limits_{\theta\in
X}\sum\limits_{t=1}^{T}{{f_{t}}({\theta})},$ (32)

where ${\theta^{*}}=\mathop{\arg\min}_{\theta\in
X}\sum\nolimits_{t=1}^{T}{{f_{t}}(\theta)}$, functions $f_{1}$ through $f_{T}$
are convex, $g_{t}:=\nabla f_{t}(\theta_{t})$, and $X\subseteq{R^{d}}$ is a
convex and closed set. The function $R_{T}$ is called the regret up to
timestep $T$. This optimization framework was introduced in Zinkevich (2003).
Then, let the online algorithm:

$\displaystyle{\nu_{t}}={\theta_{t}}-{a_{t}}g_{t}$ (33)
$\displaystyle{\theta_{t+1}}=\Pi_{X}({\nu_{t}}),$

where ${\Pi_{X}}(y)=\mathop{\arg\min}_{\theta\in X}\left\|{x-y}\right\|$ is
the projection of $y$ on $X$. Studying an online algorithm with vanishing
average regrets amounts to studying a stochastic algorithm for ERM as
discussed in Cesa-Bianchi et al. (2004).

The following assumptions are used in the online minimization setting. The
first is that functions $f_{t}(\theta)$ are convex, i.e.:

$f_{t}(\theta)\geq f_{t}(\theta^{\prime})+\nabla
f_{t}{(\theta^{\prime})^{T}}(\theta-\theta^{\prime})\,\,\text{ for all
}\,\,\theta^{\prime},\theta\in\mathbb{R}^{d}.$ (34)

Furthermore, if $X$ is a set with bounded diameter $D$ then:

$\left\|{{\theta^{\prime}}-\theta}\right\|\leq D,$ (35)

for all ${{\theta^{\prime}},\theta}\in X$. Moreover, if $X\subseteq{R^{d}}$ is
a convex and closed set, then for any $\theta\in X$ and $y\in\mathbb{R}^{d}$,
it is that:

$\left\|{y-\theta}\right\|\geq\left\|{{\theta^{\prime}}-\theta}\right\|,$ (36)

for ${\theta^{\prime}}={\Pi_{X}}(y)=\mathop{\arg\min}_{\theta\in
X}\left\|{y-\theta}\right\|$. It is noted, that the regret bounds for the
deterministic-learning-rate version of SGD are obtained by following the
proofs of the next Theorems 3,4 but using $a_{t}=\frac{a}{\sqrt{t}}$ as the
learning rate.

Then, it will be shown that stochastic-learning-rate SGD has vanishing average
regret for diminishing step size $\eta_{t}=\frac{a}{\sqrt{t}}$. The first
theorem for the online learning setting follows:

######  Theorem 3.

Assume that $X$ is convex, closed, and has bounded diameter $D$, and that the
gradient of $f_{t}$ is bounded, i.e.
$\left\|{\nabla{f_{t}}(\theta)}\right\|_{\infty}\leq G$. Furthermore, assume a
positive learning rate of $a_{t}=\frac{a}{\sqrt{t}}u_{t},\,\,\,c_{1}\leq
u_{t}\leq c_{2}$ with $a=\frac{D}{{G\sqrt{{c_{2}}{c_{1}}}}}$. Then for all
$T\geq 1$ it is that:

${R_{T}}\leq\sqrt{\frac{{{c_{2}}}}{{{c_{1}}}}}3DG\sqrt{T},$ (37)

Therefore, $R_{T}/T\to 0$ for $T\to\infty$.

###### Proof.

Firstly, it is noted that:

$\displaystyle{\left\|{{\theta_{t+1}}-{\theta^{*}}}\right\|^{2}}$
$\displaystyle={\left\|{{\Pi_{X}}({\theta_{t}}-{a_{t}}{g_{t}})-{\theta^{*}}}\right\|^{2}}$
(38)
$\displaystyle\leq{\left\|{{\theta_{t}}-{a_{t}}{g_{t}}-{\theta^{*}}}\right\|^{2}}$
$\displaystyle={\left\|{{\theta_{t}}-{\theta^{*}}}\right\|^{2}}-2{a_{t}}g_{t}^{T}({\theta_{t}}-{\theta^{*}})+a_{t}^{2}{\left\|{{g_{t}}}\right\|^{2}},$

where the first inequality followed from (36) replacing $\theta$ by
$\theta^{*}$ and $y$ by ${\theta_{t}}-{a_{t}}{g_{t}}$. Rearranging with
respect to $g_{t}^{T}({\theta_{t}}-{\theta^{*}})$ and using (34), it is:

$\displaystyle{f_{t}}({\theta_{t}})-{f_{t}}({\theta^{*}})$ $\displaystyle\leq
g_{t}^{T}({\theta_{t}}-{\theta^{*}})$ (39)
$\displaystyle\leq\frac{{{{\left\|{{\theta_{t}}-{\theta^{*}}}\right\|}^{2}}-{{\left\|{{\theta_{t+1}}-{\theta^{*}}}\right\|}^{2}}}}{{2{a_{t}}}}+\frac{{{a_{t}}}}{2}{\left\|{{g_{t}}}\right\|^{2}}.$

Then, summing over all steps:

$\displaystyle\sum\limits_{t=1}^{T}{[{f_{t}}(\theta)-{f_{t}}({\theta^{*}})]}$
$\displaystyle\leq\sum\limits_{t=1}^{T}{\frac{1}{{2{a_{t}}}}\left({{{\left\|{{\theta_{t}}-{\theta^{*}}}\right\|}^{2}}-{{\left\|{{\theta_{t+1}}-{\theta^{*}}}\right\|}^{2}}}\right)}+\frac{1}{2}\sum\limits_{t=1}^{T}{{a_{t}}{{\left\|{{g_{t}}}\right\|}^{2}}}$
(40)
$\displaystyle=\frac{1}{{2{a_{1}}}}{\left\|{{\theta_{1}}-{\theta^{*}}}\right\|^{2}}+\sum\limits_{t=2}^{T}{{{\left\|{{\theta_{t}}-{\theta^{*}}}\right\|}^{2}}\left({\frac{1}{{2{a_{t}}}}-\frac{1}{{2{a_{t-1}}}}}\right)}+\frac{1}{2}\sum\limits_{t=1}^{T}{{a_{t}}{{\left\|{{g_{t}}}\right\|}^{2}}}$
$\displaystyle\leq\frac{{{D^{2}}}}{{2{a_{1}}}}+{D^{2}}\sum\limits_{t=2}^{T}{\left({\frac{1}{{2{a_{t}}}}-\frac{1}{{2{a_{t-1}}}}}\right)}+\frac{1}{2}\sum\limits_{t=1}^{T}{{a_{t}}{{\left\|{{g_{t}}}\right\|}^{2}}}$
$\displaystyle\leq\frac{{{D^{2}}}}{{2{a_{T}}}}+\frac{1}{2}\sum\limits_{t=1}^{T}{{a_{t}}{{\left\|{{g_{t}}}\right\|}^{2}}},$

where the third line followed by using the diameter $D$ of the convex set $X$,
and the fourth line by using the telescoping sum of the middle term. Using
$a_{t}=\frac{a}{\sqrt{t}}u_{t}$ and that $c_{1}\leq u_{t}\leq c_{2}$ it is:

$\displaystyle\sum\limits_{t=1}^{T}{[{f_{t}}(\theta)-{f_{t}}({\theta^{*}})]}$
$\displaystyle\leq\frac{{{D^{2}}\sqrt{T}}}{{2a{c_{1}}}}+\frac{{{ac_{2}}}}{2}\sum\limits_{t=1}^{T}{\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\sqrt{t}}}}$
(41)
$\displaystyle\leq\frac{{{D^{2}}\sqrt{T}}}{{2a{c_{1}}}}+\frac{{{a{G^{2}}c_{2}}}}{2}\sum\limits_{t=1}^{T}{\frac{1}{{\sqrt{t}}}}\leq\frac{{\sqrt{T}}}{{2{c_{1}}a}}{D^{2}}+\frac{{a{G^{2}}{c_{2}}}}{2}2\sqrt{T},$

where the second inequality followed by using the bound on the gradients, and
third inequality by using the fact that
$\sum\nolimits_{t=1}^{T}{\frac{1}{{\sqrt{t}}}}\leq 2\sqrt{T}-1\leq 2\sqrt{T}$.
Then, minimizing the right hand side of the third inequality yields
$a={D}/({{G\sqrt{{c_{2}}{c_{1}}}}})$ for which:

${R_{T}}\leq\sqrt{\frac{{{c_{2}}}}{{{c_{1}}}}}3DG\sqrt{T},$ (42)

thus concluding the proof. ∎

The next theorem connects SGD with stochastic-learning-rate schemes with the
data dimension $d$.

######  Theorem 4.

Assume that $X$ is convex, closed, and has bounded diameter D, and that the
gradient of $f_{t}$ is bounded, i.e.
$\left\|{\nabla{f_{t}}(\theta)}\right\|_{\infty}\leq G$. Furthermore, assume a
positive learning rate of $a_{t}=\frac{a}{\sqrt{t}}u_{t},\,\,\,c_{1}\leq
u_{t}\leq c_{2}$ where $\left\|{{g_{1:t}}}\right\|$ denotes
$\sqrt{\sum\nolimits_{t=1}^{t}{{{\left\|{{g_{t}}}\right\|}^{2}}}}$. Then, for:

$a=\frac{1}{{\sqrt{{c_{1}}{c_{2}}}}}\frac{D}{{\sqrt{2}\left\|{{g_{1:t}}}\right\|}}$
(43)

it is that, for all $T\geq 1$ :

${R_{T}}\leq\sqrt{\frac{{{c_{2}}}}{{{c_{1}}}}}\sqrt{2}DG\sqrt{dT}.$ (44)

###### Proof.

The proof is identical to the one of Theorem 3 until (40). Then, using
$a_{t}=\frac{a}{\sqrt{t}}u_{t}$ and (43) in (40) it is:

$\displaystyle\sum\limits_{t=1}^{T}{[{f_{t}}(\theta)-{f_{t}}({\theta^{*}})]}$
$\displaystyle\leq\frac{{{D^{2}}\sqrt{2{c_{1}}{c_{2}}}}}{{2D{c_{1}}}}\left\|{{g_{1:T}}}\right\|+\frac{{{c_{2}}D}}{{2\sqrt{2{c_{1}}{c_{2}}}}}\sum\limits_{t=1}^{T}{\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\left\|{{g_{1:t}}}\right\|}}}$
(45)
$\displaystyle\leq\sqrt{\frac{{{c_{2}}}}{{{c_{1}}}}}\frac{\sqrt{2}D}{2}\left\|{{g_{1:T}}}\right\|+\sqrt{\frac{{{c_{2}}}}{{{c_{1}}}}}\frac{D}{2\sqrt{2}}\sum\limits_{t=1}^{T}{\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\left\|{{g_{1:t}}}\right\|}}}$
$\displaystyle\leq\sqrt{\frac{{{c_{2}}}}{{{c_{1}}}}}\sqrt{2}D\left\|{{g_{1:T}}}\right\|.$

The third inequality followed from using
$\sum\nolimits_{t=1}^{T}{\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\left\|{{g_{1:t}}}\right\|}}}\leq
2\left\|{{g_{1:T}}}\right\|$. This can be proved by induction, since for $t=1$
it holds that $\left\|{{g_{1}}}\right\|\leq 2\left\|{{g_{1}}}\right\|$. Then,
assuming the inequality holds for $T-1$, for timestep $T$ it is:

$\displaystyle\sum\limits_{t=1}^{T}{\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\left\|{{g_{1:t}}}\right\|}}}=\sum\limits_{t=1}^{T-1}{\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\left\|{{g_{1:t}}}\right\|}}+\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\left\|{{g_{1:T}}}\right\|}}}$
$\displaystyle\leq
2\left\|{{g_{1:T-1}}}\right\|+\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\left\|{{g_{1:T}}}\right\|}}$
(46) $\displaystyle\leq
2\sqrt{{{\left\|{{g_{1:T}}}\right\|}^{2}}-g_{T}^{2}}+\frac{{{{\left\|{{g_{t}}}\right\|}^{2}}}}{{\left\|{{g_{1:T}}}\right\|}}\leq
2\left\|{{g_{1:T}}}\right\|.$

The last inequality follows from:

$\displaystyle 2\sqrt{{{\left\|{{g_{1:T}}}\right\|}^{2}}-g_{T}^{2}}$
$\displaystyle=2\sqrt{{{\left\|{{g_{1:T}}}\right\|}^{2}}-\frac{1}{2}\left\|{{g_{1:T}}}\right\|\frac{1}{2}\frac{{g_{T}^{2}}}{{\left\|{{g_{1:T}}}\right\|}}+\frac{{g_{T}^{4}}}{{4{{\left\|{{g_{1:T}}}\right\|}^{2}}}}-\frac{{g_{T}^{4}}}{{4{{\left\|{{g_{1:T}}}\right\|}^{2}}}}}$
(47)
$\displaystyle=2\sqrt{{{\left({{{\left\|{{g_{1:T}}}\right\|}^{2}}-\frac{{g_{T}^{2}}}{{2\left\|{{g_{1:T}}}\right\|}}}\right)}^{2}}-\frac{{g_{T}^{4}}}{{4{{\left\|{{g_{1:T}}}\right\|}^{2}}}}}$
$\displaystyle\leq
2\sqrt{{{\left({{{\left\|{{g_{1:T}}}\right\|}^{2}}-\frac{{g_{T}^{2}}}{{2\left\|{{g_{1:T}}}\right\|}}}\right)}^{2}}}$
$\displaystyle\leq
2{\left\|{{g_{1:T}}}\right\|^{2}}-\frac{{g_{T}^{2}}}{{\left\|{{g_{1:T}}}\right\|}}.$

Then, since $g_{t}:=\nabla f_{t}(\theta_{t})$ for $\theta_{t}\in{R^{d}}$, the
norm $\left\|{{g_{1:T}}}\right\|$ in (45) is written as
$\sqrt{\sum\nolimits_{t=1}^{T}{\sum\nolimits_{j=1}^{d}{{{\left|{{{g_{t,j}}}}\right|}^{2}}}}}$
where $g_{t,j}$ is the $j$th element of the gradient vector at timestep $t$.
Using the bound $G$ on the gradients it is:

$\sum\limits_{t=1}^{T}{[{f_{t}}(\theta)-{f_{t}}({\theta^{*}})]}\leq\sqrt{\frac{{{c_{2}}}}{{{c_{1}}}}}\sqrt{2}DG\sqrt{dT},$
(48)

thus proving the theorem. ∎

In conclusion, from Theorems 3,4 it is thus observed that for the online
optimization setting the regret bounds when using stochastic-learning-rates
differ by $\sqrt{\frac{c_{2}}{c_{1}}}$ from the bounds of the original SGD
algorithm. In specific, the guarantees in (37) and (44) are the same as in the
online learning setting for the original SGD, except for a multiplying factor
of $\sqrt{\frac{{{c_{2}}}}{{{c_{1}}}}}$ compared to previous works (e.g. Duchi
et al. (2011)).

##  4 Discussion on the Stochasticity Factor

In this work, the stochastic learning rate scheme considered for the empirical
results is a RUMSLR scheme and specifically
$u_{t}\sim\mathcal{U}(\sqrt[\varepsilon_{t}]{{{c_{1}}}},\sqrt[\varepsilon_{t}]{{{c_{2}}}})$
where $\varepsilon_{t}=t\,\,(\bmod\beta)+1$, and $\beta\geq 1$, shown in Fig.
1. In this case the distribution of $u_{t}$ resets every $\beta$ number of
iterations to $\mathcal{U}(c_{1},c_{2})$, named $\beta-$RUMSLR to distinguish
it from other RUMSLR schemes.

Furthermore, the bounds in the online learning setting in Theorems 3,4 as well
as the step size in the stochastic optimization setting in Theorem 1 present a
trade-off for the choice of $c_{1},c_{2}$. On one side they should be kept
close enough so that the stochastic-learning-rate SGD step size and online
bounds are close to that of the original SGD, since they differ by a factor of
$c_{1}/c_{2}^{2}$ and $\sqrt{c_{2}/c_{1}}$ respectively. For example
$c_{1}=0.8,c_{2}=1.2$ yield a modest 0.2 increase in the bound of RUMSLR SGD
in Theorems 3,4 compared to the original SGD, but the further away they are,
the larger the bound becomes, and the weaker the guarantees (or smaller the
step size in the case of Theorem 1) become. Indeed, it has been observed that
choosing low values for $c_{1}$ or large values for $c_{2}$ slows down the
algorithm significantly or destabilizes it. On the other side, they should be
kept adequately apart so that the bursts are significant enough. In this work,
the support of the distribution of $u_{t}$ is taken to be symmetric about 1
when it resets, even if the theoretical analysis holds for arbitrary positive
$c_{1}<c_{2}$. Moreover, it has been observed that continuous stochastic
bursts can increase convergence time significantly. Therefore, the
distribution of $u_{t}$ is chosen to be reset after it is adequately close to
unity and enough, $\beta$ timesteps have passed, e.g., as shown in Fig. 1,
where $\beta=100$, indicated by the bursts that occur every 100 iterations.

Experimentally, however, it will be demonstrated that $\beta-$RUMSLR SGD
yields noticeable improvements in minimization performance for datasets such
as CIFAR-10, and CIFAR-100 (Krizhevsky (2009)) when compared to SGD without
stochastic learning rate. This scheme was equipped to the learning rates of
other algorithms as well such as ADAM Kingma and Ba (2015), AMSGrad Reddi et
al. (2018), ADAMW Loshchilov and Hutter (2019), and two SGD with momentum
algorithms originally presented in Polyak (1964) and Nesterov (1983)
(subsequently referred to as SGD with momentum first and second versions
respectively), all of which, even though not theoretically studied in this
work, yielded similar improved performance when using $\beta-$RUMSLR.

##  5 Results

Experimental results using the $\beta-$RUMSLR schemes applied to SGD, SGD with
momentum versions first and second, ADAM, AMSGrad, and ADAMW are presented in
this section. For $\beta-$RUMSLR the plots show the average running loss every
20 mini-batches for CIFAR-10 and CIFAR-100, and every 200 mini-batches for
MNIST. The batch size was 128 in all experiments. The rest of the
hyperparameters, i.e., learning rates, momentums (when applicable) and SF
constants, are chosen via grid search. The datasets were standardized before
use. For CIFAR-10 and CIFAR-100, the network architecture was Pytorch’s
ResNet18. For MNIST, logistic regression was used. For the former two datasets
a constant stepsize was used, and a diminishing stepsize for the latter, for
all algorithms.

For the CIFAR-10 and CIFAR-100 datasets the train and test sets are composed
of 50000 train and 10000 test images respectively out of a total of 60000
images (an approximately 83,33%/16,67% train/test split). For the MNIST
dataset the train and test sets are composed of 60000 train and 10000 test
images respectively out of a total of 70000 images (an approximately
85,71%/14,29% train/test split).

The MNIST dataset was used to validate the online framework theoretical
guarantees for $\beta-$RUMSLR SGD using a convex loss function, in specific
logistic regression as aforementioned. Indeed, the results in Fig. 2 verify
the convergence of the MNIST $\beta-$RUMSLR scheme in a convex setting as
predicted by Theorems 3,4. Moreover, CIFAR-10 and CIFAR-100, the
$\beta-$RUMSLR scheme demonstrated significantly improved convergence than
when using a deterministic learning rate as shown in Figs. 3 and 4
respectively. As can be seen by the plots, $\beta-$RUMSLR algorithms
convergence much faster than their deterministic-learning-rate counterparts.

Figure 2: Plots that illustrate convergence for the $\beta$-RUMSLR SGD for
$\beta=100$, $c_{1}=0.7$ and $c_{2}=1.3$ in a convex setting.

Figure 3: Plots that illustrate the $\beta$-RUMSLR scheme (in orange) and the
deterministic learning rate (in blue) schemes for $\beta=100$, $c_{1}=0.8$ and
$c_{2}=1.2$. It is observed that for the CIFAR-10 dataset algorithms with
$\beta$-RUMSLR yield noticeable gains in minimization performance.

Figure 4: Plots that illustrate the $\beta$-RUMSLR (in orange) and the
deterministic learning rate (in blue) schemes for $\beta=100$, $c_{1}=0.9$ and
$c_{2}=1.1$. AMSGrad uses $c_{1}=0.85$ and $c_{2}=1.15$. For the CIFAR-100
dataset, algorithms with $\beta$-RUMSLR yield noticeable minimization
performance gains.

##  6 Conclusion and Future Work

This work introduced learning-rate schemes that make the learning rate
stochastic by multiplicatively equipping it with a random variable, coined
stochastic factor. Convergence of the SGD algorithm employing stochastic
learning rate schemes was theoretically analyzed and compared with known
results for the algorithm’s deterministic-learning-rate version. The
discussion that followed suggested how the hyperparameter values introduced by
the two stochastic-learning-rate schemes should be chosen. Empirical results
on popular algorithms demonstrated noticeable increase in optimization
performance, presenting stochastic-learning-rate schemes as a viable option
for enhancing performance. In-depth generalization performance studies,
convergence analysis of algorithms besides SGD, and investigating the effect
of various stochastic factor distributions and hyperparameters on algorithm
performance are some paths for future work.

## References

  * Bertsekas and Tsitsiklis [2000] D. Bertsekas and J. Tsitsiklis.  Gradient convergence in gradient methods with errors.  _SIAM Journal on Optimization_ , 10(3):627–642, 2000.  ISSN 1052-6234.  doi: 10.1137/S1052623497331063.
  * Bottou et al. [2018] L. Bottou et al.  Optimization methods for large-scale machine learning.  _SIAM Review_ , 60(2):223–311, jan 2018.  doi: 10.1137/16m1080173.
  * Cesa-Bianchi et al. [2004] N. Cesa-Bianchi et al.  On the generalization ability of on-line learning algorithms.  _IEEE Transactions on Information Theory_ , 50(9):2050–2057, 2004.  doi: 10.1109/TIT.2004.833339.
  * Duchi et al. [2011] J. Duchi et al.  Adaptive subgradient methods for online learning and stochastic optimization.  _J. Mach. Learn. Res._ , 12:2121–2159, July 2011.  ISSN 1532-4435.
  * Ghadimi and Lan [2013] S. Ghadimi and G. Lan.  Stochastic first- and zeroth-order methods for nonconvex stochastic programming.  _SIAM J. Optim._ , 23:2341–2368, 2013.
  * Karimi et al. [2016] H. Karimi et al.  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In P. Frasconi, N. Landwehr, G. Manco, and J. Vreeken, editors, _Machine Learning and Knowledge Discovery in Databases_ , pages 795–811, Cham, 2016. Springer International Publishing.  ISBN 978-3-319-46128-1.
  * Kingma and Ba [2015] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In _ICLR (Poster)_ , 2015.  URL http://arxiv.org/abs/1412.6980.
  * Krizhevsky [2009] A. Krizhevsky.  Learning multiple layers of features from tiny images.  Technical report, 2009.
  * Loshchilov and Hutter [2017] I. Loshchilov and F. Hutter.  Sgdr: Stochastic gradient descent with warm restarts, 2017.
  * Loshchilov and Hutter [2019] I. Loshchilov and F. Hutter.  Decoupled weight decay regularization.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=Bkg6RiCqY7.
  * Mamalis et al. [2021] T. Mamalis et al.  Accelerated almost-sure convergence rates for nonconvex stochastic gradient descent using stochastic learning rates, 2021.
  * McMahan and Streeter [2010] H. B. McMahan and M. J. Streeter.  Adaptive bound optimization for online convex optimization.  _CoRR_ , abs/1002.4908, 2010.  URL http://dblp.uni-trier.de/db/journals/corr/corr1002.html#abs-1002-4908.
  * Moulines and Bach [2011] E. Moulines and F. Bach.  Non-asymptotic analysis of stochastic approximation algorithms for machine learning.  In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 24. Curran Associates, Inc., 2011.  URL https://proceedings.neurips.cc/paper/2011/file/40008b9a5380fcacce3976bf7c08af5b-Paper.pdf.
  * Nemirovski et al. [2009] A. Nemirovski et al.  Robust stochastic approximation approach to stochastic programming.  _SIAM J. on Optimization_ , 19(4):1574–1609, Jan. 2009.  ISSN 1052-6234.  doi: 10.1137/070704277.  URL https://doi.org/10.1137/070704277.
  * Nesterov [1983] Y. Nesterov.  A method for solving the convex programming problem with convergence rate o(1/$k^{2}$).  _Proceedings of the USSR Academy of Sciences_ , 269:543–547, 1983.
  * O’Donoghue and Candes [2012] B. O’Donoghue and E. Candes.  Adaptive restart for accelerated gradient schemes, 2012.
  * Polyak [1964] B. T. Polyak.  Some methods of speeding up the convergence of iteration methods.  _USSR Computational Mathematics and Mathematical Physics_ , 4(5):1–17, jan 1964.  doi: 10.1016/0041-5553(64)90137-5.
  * Reddi et al. [2018] S. J. Reddi et al.  On the convergence of adam and beyond.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=ryQu7f-RZ.
  * Robbins and Monro [1951] H. Robbins and S. Monro.  A stochastic approximation method.  _Ann. Math. Statist._ , 22(3):400–407, 09 1951\.  doi: 10.1214/aoms/1177729586.  URL https://doi.org/10.1214/aoms/1177729586.
  * Zinkevich [2003] M. Zinkevich.  Online convex programming and generalized infinitesimal gradient ascent.  In _In ICML_ , pages 928–936, 2003.