# Inverse Compton emission from heavy WIMP annihilations in the Galactic
Centre

Julia I. Djuvsland<EMAIL_ADDRESS>Jim Hinton  Brian Reville

###### Abstract

A thermal relic WIMP remains a prime candidate for the nature of Dark Matter,
particularly for the more poorly constrained case of a heavy ($\gtrsim$ 1 TeV)
WIMP. The highest fluxes from WIMP annihilations are expected in the region of
the Galactic Centre (GC) where current and near future gamma-ray observatories
can be exploited to place tight limits on the WIMP paradigm. It is regularly
noted that the annihilation flux of gammas will be accompanied by charged
secondary particles which can produce ‘delayed’ inverse Compton (IC) gamma-ray
emission, but this component is often neglected in indirect Dark Matter
searches. In this work the inverse Compton emission is studied for the
specific conditions of heavy WIMP annihilation in the GC. Using models for the
magnetic and radiation fields of the region, and taking into consideration the
transport of secondary particles, we find that for TeV WIMPs the IC component
cannot be neglected in the GC, with the particles produced cooling within the
region rather than propagating out in to the Galaxy. This effect changes the
predicted spectral shape substantially and thus boosts the detection prospects
for heavy WIMPs.

###### keywords:

Dark Matter , Galactic Centre , Inverse Compton emission

††journal: Physics of the Dark Universe

[MPIK]organization=Max-Planck-Institut für Kernphysik,
addressline=Saupfercheckweg 1, city=Heidelberg, postcode=69117,
country=Germany [IFT]organization=Department for Physics and Technology,
University of Bergen, addressline=Allegaten 55, city=Bergen, postcode=5020,
country=Norway

##  1 Introduction

A weakly interacting massive particle (WIMP), left over as a _thermal relic_
[1] of the Big Bang, remains a compelling candidate for the Dark Matter (DM)
that makes up 27% of the energy density in our universe [2]. While there is
ample evidence for the existence of DM through astrophysical observations like
gravitational lensing [3], galaxy rotation curves [4] and the study of the
Cosmic Microwave Background (CMB) [2], the particle nature of DM is still
poorly understood. To address this, many efforts are currently underway to
search for a DM particle candidate at particle colliders, in underground
laboratories and by exploring the gamma-ray sky – to name just a few.

A popular model explored in this context is that of annihilating WIMPs, whose
parameter space is narrowed by ongoing searches. The DM searches at CERN’s
Large Hadron Collider (LHC) explore the GeV WIMP mass region and seem to
advocate heavy candidates, e.g. [5, 6]. This mass range can also be explored
via indirect DM searches that measure gamma-rays arriving at Earth using
detectors such as the Fermi-LAT satellite, ground-based Imaging Atmospheric
Cherenkov Telescopes (IACTs) or high altitude air shower observatories.

While the current and future IACTs - H.E.S.S. [7], MAGIC [8], VERITAS [9] and
CTA [10] – are most sensitive to WIMPs with TeV masses, high altitude air
shower observatories such as HAWC, LHAASO and SWGO can constrain even higher
WIMP masses ($O$(100 TeV)) [11, 12, 13]. Fermi-LAT currently quotes its most
stringent exclusion limits in the GeV range with little sensitivity for TeV-
WIMPs [14]. These limits are however obtained by only taking the _direct_
photon component into account. With direct photons, we refer to photons that
are directly produced during the WIMP annihilation and the consecutive decay
of the immediate annihilation products, i.e. during final-state radiation and
hadronisation processes. Yet, gamma-rays can also be produced indirectly, in
particular by the _direct_ electrons111For simplicity the term electron is
used to describe both direct electrons and direct positrons arising from WIMP
annihilation throughout this paper. that arise from the WIMP annihilation.
They give rise to energetic photons via inverse Compton (IC) scattering of
ambient photons, and also synchrotron and Bremsstrahlung emission. These
_indirect_ photons, from the IC process in particular, contribute to the
expected gamma-ray signal from WIMPs and have typically lower energies than
the signal component from the direct photons.

The shape and intensity of the indirect photon component of the Galactic WIMP
gamma-ray spectrum depends highly on the transport properties of charged
particles near the GC. This is because the GC conditions impact the lifetime
of the radiating electrons. Cosmic-ray (CR) transport around the Galaxy
remains however an unsolved problem. The majority of observational constraints
are inferred from local measurements, which alone can not provide an accurate
picture of the relevant physical parameters in the GC region. Gamma-ray
observations provide an indirect probe of the CR distribution, and provide
compelling evidence for greatly reduced diffusion coefficients with respect to
the Galactic average [15] in this region. The complex magnetic field topology
[16, 17] and fluid motions (in particular outflows [18, 19]) in the GC region
influence the emission profiles of both the astrophysical foregrounds and
crucially, any delayed emission probed in indirect WIMP searches. The
transport of the emitting particles in this extreme environment are not
expected to match those of the Galactic average: using local CR observations
one would infer a scattering mean free path of approximately 10 pc for
particles with TeV energies. In this work we highlight the physical motivation
for reduced scattering lengths in the inner kpc region. The results that
follow make the further simplifying assumption of isotropic diffusion.

The importance of taking the Bremsstrahlung and IC-photons into account in
indirect WIMP searches was discussed previously, e.g. [20, 21], especially for
light WIMP masses [22] and in the context of the Fermi GeV excess [23]. Also
the importance of synchrotron radiation effects were studied [24, 21]. In [25]
the potential of a multiwavelength analysis combining the gamma ray and
microwave observations was explored.

Still, the indirect component of the gamma-ray signal is frequently ignored in
indirect DM searches. Therefore this work highlights the fact that the
indirect component is not negligible for TeV mass WIMP annihilation in the
Galactic Centre. As the indirect component leads to additional photons (with
lower energies) than expected from the direct photon contribution only,
ignoring this contribution would result in searches underestimating their
sensitivity and/or misinterpretation of a detected signal.

##  2 Cosmic ray transport and energy losses in the GC region

The central few hundred pc of our Galaxy is a complex and unique environment,
dominated observationally by the Central Molecular Zone (CMZ), with $\sim
5\times 10^{7}$ solar masses of gas in a thin ($\sim$30 pc) highly structured
disc. Densities in the CMZ are orders of magnitude larger than is typical in
the disc of the Milky Way (MW) and strong magnetic fields have been detected
[16]. Either the star formation of the CMZ, or the central supermassive black
hole Sgr A⋆ is responsible for blowing the Fermi/eROSITA bubbles [26, 27],
which extend for 55∘/80∘ degrees above and below the GC, and which appear to
be connected to the CMZ by _Chimneys_ recently detected in the X-ray band
[19]. The heart of the Galaxy is also naturally where the most intense large
scale radiation fields can be found [28]. Taken together it is clear that the
GC is a unique environment, in no way resembling typical MW conditions.

Whilst emission from the entire GC region is expected from the annihilation of
a hypothetical thermal relic WIMP, the emission within the CMZ itself would be
very difficult to disentangle from the foreground astrophysical emission [15,
29, 30, 31]. Above and below the CMZ, astrophysical foregrounds are greatly
reduced, and conditions are expected to be much less extreme than within it,
both in terms of densities and magnetic fields. These regions are nonetheless
expected to be atypical, with strong radiation fields and magnetic turbulence,
as well as fast outflows. In addition, the presence of a wind from the CMZ is
now well established [18, 32], and advection may therefore dominate the
transport of relativistic charged particles over a wide energy range in the
GC.

The radiation field model used in this work is a self-consistent model of the
broad-band continuum emission of our Galaxy and was derived from modelling
maps of the all-sky emission in the infrared and submillimetre regime [28].
The magnetic field employed here follows the prescription of Jansson & Farrar,
including large-scale regular fields, striated fields and small-scale random
fields [33, 34]. The regular field consists of a disk field and an extended
halo field including a large, out-of-plane component and the orientation of
the striated component is aligned with the regular field. The magnetic field
topology is undoubtedly more complex in the inner $\sim 100$ pc region [17],
though such complications are beyond the scope of the present work.

The diffusion coefficient in the region above the CMZ is poorly constrained at
present, though there are strong physical arguments that indicate it is
substantially smaller than the galactic average. TeV gamma-ray observations of
the diffuse emission in the central region favour an hadronic origin, as
evidenced by its close correlation with the gas maps [29]. The resulting
cosmic-ray distribution inferred from the gas maps is consistent with a $1/r$
radial profile, which constrains the diffusion coefficient to be substantially
smaller than the galactic average, at least in the CMZ [35].

One can estimate a diffusion coefficient within the standard quasi-linear
framework. The scattering rate of a magnetised particle, i.e. a particle
undergoing helical motion about a mean field $B_{0}$, is (in cgs units) [e.g.
36],

$\nu=\left.\frac{\pi}{4}\left(\frac{k\mathcal{E}_{k}}{B_{0}^{2}/8\pi}\right)\frac{c}{r_{\rm
g}}\right|_{k\mu r_{\rm g}=1},$

where $r_{\rm g}=E/eB_{0}$ is the particle Larmor radius, $\mu$ its pitch
angle with respect to the local mean magnetic field. $\mathcal{E}_{k}dk$ is
the energy density in waves, with wavelengths $\lambda=2\pi/k$, that are
_gyro-resonant_ with particles of energy $E$ and pitch $\mu$. Assuming a
turbulent power-law spectrum of scattering modes, $\mathcal{E}_{k}\propto
k^{-\alpha}$, sharply cut-off above $\lambda>L$, an approximate expression for
the average diffusion coefficient along $B_{0}$ is

$D\lesssim 10^{27}\eta^{-1}\left\\{\begin{array}[]{cl}E_{\rm
TeV}^{1/3}B_{10\mu{\rm G}}^{-1/3}L_{\rm pc}^{2/3}&\text{Kolmogorov}\\\ E_{\rm
TeV}^{1/2}B_{10\mu{\rm G}}^{-1/2}L_{\rm
pc}^{1/2}&\text{Kraichnen}\end{array}\right.$

Here, Kolmogorov and Kraichnen turbulent spectra correspond to $\alpha=5/3$
and $3/2$ respectively. The quantity $\eta$ (assumed to be $\leq 1$) is the
ratio of total integrated turbulent magnetic field energy density to that in
the mean field, while $L_{\rm pc}$ is $L$ in units of parsecs, $E_{\rm Tev}$
the particle energy in TeV. At large galacto-centric radii $L_{\rm pc}\approx
100$ [e.g. 37]. If we assume the turbulence in the inner regions is driven by
stellar activity, $L$ should decrease as one moves closer to the CMZ. In the
following, we set $\eta=1$ and take the outer scale of the turbulence $L$ to
be $10\%$ of its distance from the Galactic Centre (i.e. approaching the
average ISM value at 1 kpc). To simplify the calculations that follow, we take
$D$ as a proxy for the radial diffusion coefficient and focus exclusively on
Kolmogorov turbulence.

With these assumptions about the conditions in the GC region, we use the open-
source GAMERA package [38], as employed e.g. in [39], to calculate the
timescales of the electron cooling. In Figure 1 we compare them to the
timescales of electron propagation in the GC. The timescales are shown for
electrons with two exemplary energies, namely 10 GeV (blue lines) and 100 GeV
(orange lines). The total cooling timescales are shown by the solid lines,
while the timescales of the inverse Compton emission are shown by the dashed
lines. The dotted lines show the diffusion timescales for the different
electron energies and the green dash-dotted line shows the advection
timescales for a constant wind speed of 200 km/s.

Also shown are $r_{P}$, the radii with maximal signal significance222Where the
signal significance corresponds to the ratio of the number of signal events
(direct photon counts) and the square root of the number of background events
for a constant background distribution within a sphere., for a dark matter
distribution according to a Navarro, Frenk and White (NFW) [40, 41] (grey
solid line) and an Einasto [42] (grey dashed line) profile. These radii
correspond to the typical distances that are employed in galactic indirect DM
searches as can be seen from the grey area that indicates the region of
interest (ROI) of the galactic DM search of the H.E.S.S. Collaboration [43].

The shape of the WIMP distribution in the MW is however poorly constrained.
While the cusp NFW and Einasto profiles have been widely employed in indirect
DM searches, there is experimental evidence that prefers shallower cusps or
even cored density profiles [44]. The latter profiles have larger values for
$r_{P}$, which means that the expected signal strength is weaker in the GC
than for the cusp profiles and that the shallow cusp or cored morphology is
disadvantageous for WIMP searches at the GC.

From Figure 1 it is evident that for the GC distances that are explored by
gamma ray telescopes the IC emission dominates the cooling. At larger
distances synchrotron effects become more important while at smaller distances
both synchroton emission and Bremsstrahlung contribute significantly to the
electron cooling (not shown in Figure 1). At about 20 pc the synchroton
emission and Bremsstrahlung lead to a sudden step in the cooling timescales
and at even smaller distances the cooling times are only tentative as the
magnetic field model breaks down inside the CMZ.

Figure 1: Timescales of electron cooling (solid lines), inverse Compton
emission (dashed lines) and diffusion (dotted lines) as a function of the
distance from the Galactic Centre. The timescales are shown for the two
different electron energies 10 GeV and 100 GeV. The dash-dotted green line
shows the timescale of advection with a velocity of 200 km/s. The vertical
lines show the radius at which the signal significance peaks for a DM
distribution according to a NFW (solid line) and an Einasto (dashed line)
profile. The grey area indicates a typical region of interest for galactic
indirect DM searches, here the one from [43].

All in all, we can see from Figure 1 that the cooling timescales are shorter
than the propagation timescales for energetic electrons ($>$ 10 GeV) meaning
that the electrons are more likely to cool locally/in-situ than to leave the
GC region.

##  3 Inverse Compton Emission

Electrons that are produced via the annihilation of heavy WIMPs in the GC
mostly cool locally/in-situ as shown in Figure 1. This energy loss is likely
dominated by inverse Compton scattering especially in the region that is
investigated by DM searches with IACTs. During the IC process electrons
scatter off photons from the CMB, starlight, UV and IR light and transfer
their energy such that energetic photons are produced. Consequently, the
electrons from WIMP annihilation in the GC give rise to an additional
component in the photon spectrum apart from the photons originating directly
from the annihilation process, from the same spatial region as the direct
signal.

For our estimations, we use the spectra of electrons and positrons from WIMP
annihilation described in the widely used ’cookbook’ [21] and available in
numerical form [45, 46]. For the calculation of the resulting photon spectra,
we employ the GAMERA package. We assume a constant injected power over a given
timescale, resulting in an equilibrium situation above a certain energy. The
WIMP density profile governs at which position the IC emitting particles are
produced but the results presented here are largely independent of the density
profile, at least in the energy range where in-situ/local cooling can safely
be assumed.

(a) Total photon spectrum and its components

(b) Comparison of expected spectra with data

Figure 2: Total photon flux from WIMP annihilation to $\tau$ leptons in the
GC. For the normalisation, the J-factor is set to $1.53\times 10^{22}$
GeV2cm-5 and the thermal cross section of $\langle\sigma v\rangle=3\times
10^{-26}$ cm3/s is used. Left panel: The gray solid line shows the total
photon spectrum that is composed of the direct photon component (solid black
line) and the indirect photon components: Inverse Compton emission,
synchrotron radiation and Bremsstrahlung (dashed, dotted and dash-dotted line
respectively). The yellow line shows the spectrum of the electrons and
positrons that are directly emitted from the WIMP annihilation and that give
rise to the indirect photon contributions. The WIMP mass is set to 100 GeV and
the model was evaluated at a distance of 100 pc from the GC. Right panel:
Comparison of the total photon spectra for 3 different WIMP masses with
experimental data. The WMAP haze measured with Planck [47] (green), the Fermi
GC excess [14] (blue) and the extrapolated flux measured with H.E.S.S. [43]
(red).

Figure 2 shows the differential electron and photon flux times squared energy
(i.e. total power output per logarithmic energy interval) expected for WIMP
annihilation to $\tau$ leptons. The direct photon spectrum (black solid line
in the left panel) is taken from [21, 45, 46] and corresponds to a model-
independent baseline photon contribution, i.e. three-body final states of WIMP
annihilation are omitted – even though they might be viable in some models.
The component denoted “direct $e^{\pm}$” (yellow line) corresponds to twice
the tabulated positron flux from [21, 45, 46] to account for electron as well
as positron production. This direct electron-positron spectrum is used as
input to our model and evaluated in this case at a distance of 100 pc from the
GC as this is the typical region of interest for gamma-ray WIMP searches (as
illustrated in Figure 1). The model is evolved for $10^{6}$ years and the
escape of low energetic electrons due to advection is implemented in GAMERA in
order to obtain the indirect components of the photon spectrum (broken lines).
The indirect photon components add a low energy tail to the direct component.
The former is dominated by the IC emission for $>$GeV-mass WIMPs which is why
we concentrate on this contribution in this work. Bremsstrahlung is not
expected to contribute significantly in the low density region outside of the
Galactic disc / CMZ, synchrotron radiation contributes at much lower energies.

The total photon spectrum is given by the sum of all broken lines in the left
panel of the figure and illustrated as a gray solid line. The gray curves in
the left and right panel are equivalent and show the expected photon spectrum
for a thermal relic WIMP with $\langle\sigma v\rangle=3\times 10^{-26}$ cm3/s
and a mass of 100 GeV. The branching ratio to $\tau$ leptons is set to 1. For
the normalisation, a $J$-factor of $1.53\times 10^{22}$ GeV2cm-5 is used and
corresponds to the value computed in [14] for a NFW profile.

In the right panel of Figure 2, the expected photon spectra are compared to
experimental data to illustrate the respective constraining power. The green
curve corresponds to the WMAP haze measured with Planck [47] in a field of
view of 10° and constrains the synchrotron contribution of the WIMP
annihilation spectra. In blue the data points of the so-called Fermi GC excess
[14] are shown. These data constrain thermal relic WIMPs up to tens of GeV in
mass and is particularly sensitive to the indirect component of the photon
spectrum. The red line indicates the flux that should have been measured with
H.E.S.S. if its ROI corresponded to the one of Fermi. To extrapolate this
value, the expected direct photon spectrum of [21] for a 1 TeV WIMP is used
together with the excluded $\langle\sigma v\rangle$ value from [43]. To
account for the difference in size between the Fermi and the H.E.S.S. ROI, the
computed H.E.S.S. flux is scaled by the ratio of the J-factors for NFW
profiles from [14] and [43]. The H.E.S.S. data excludes thermal relic WIMPs
with TeV mass for this annihilation channel and assumed profile. This figure
illustrates that Galactic WIMPs can be probed in multiple wavelength ranges if
all emission components are taken into account and that the direct photon
component covers only a comparatively narrow energy band.

Having established the nature and importance of the IC component of the photon
spectrum, we are investigating its robustness in the remaining part of this
section. In Figure 3 the evolution of the IC component of the photon spectrum
from a WIMP signal is shown as a function of time (left panel) and as a
function of distance (right panel). The spectra are normalised to the total
luminosity of the direct photons, $L_{\gamma}$.

From Figure 1 it is clear that the relevant timescales for residence of
injected particles in the Galactic Centre region ranges from $10^{5}$ –
$10^{7}$ years. Figure 3(a) shows the impact of this confinement timescale on
the resulting spectral energy distribution of IC emission. The coloured
components show the IC emission of the electrons after evolving our model for
a given time. The black dashed line illustrates the nominal implementation of
our model and corresponds to an evolution of the system for $10^{6}$ years
including the escape of low energetic electrons due to advection. By
accounting for the escape of electrons, the timescale of the evolution of the
model becomes irrelevant as long as it is larger than $10^{6}$ years. On very
short timescales the emission is dominated by the highest energy electrons,
and the illustrated case of a 1 TeV WIMP ($W\overline{W}$ channel),
significantly suppressed by the Klein-Nishina (KN) effect. On the timescales
relevant for electrons in the GC ( $>10^{5}$ years) the IC emission
constitutes a significant component of the total photon spectrum.

The flux of the IC component of the photon spectrum depends on the local ratio
of magnetic to radiation energy density. The spectral shape of the IC
component depends on the local radiation spectral energy distribution (SED).
Figure 3(b) shows the dependence of the resulting IC SED on the distance from
the GC within the framework of our model, for the same injected spectrum as in
Figure 3(a). Overall the impact of changing conditions is relatively modest,
provided the distance of the emission region from the GC is $<$1 kpc – as
expected for the bulk of the annihilation signal for almost any reasonable
density profile.

(a) Time dependence of IC emission

(b) Distance dependence of IC emission

Figure 3: Spectrum of direct and inverse Compton photons from WIMP
annihilation to $W$ bosons in the GC, normalised to the total luminosity of
the direct photons. The solid line shows the spectrum for direct photons while
the dashed lines show the spectra of IC photons. The latter are shown for
different timescales of ongoing DM annihilation in the left panel (at 100 pc
distance from the GC) and for different distances from the GC in the right
panel (for a fixed timescale of $10^{6}$ years). The WIMP mass is set to 1 TeV
in this example.

In Figure 4 the direct and IC components of the photon spectra are compared
for different WIMP masses and decay channels. Again, the spectra are
normalised to the total luminosity of the direct photons. The model is
evaluated at 100 pc from the GC after an evolution time of $10^{6}$ years and
includes electron escape with the nominal outflow velocity as illustrated in
Figure 1.

Changing the mass of the annihilating WIMP modifies both the spectrum of
injected electrons and the impact of KN suppression on the resulting IC
emission. With electroweak corrections, the $W\overline{W}$ channel from [21]
produces harder photon and electron spectra at higher masses, at the same time
the emission from injected high-energy electrons is increasingly suppressed by
the KN effect. As shown in Figure 4(a), the result of these two competing
effects is a relatively stable IC spectrum over a wide mass range (0.1–10
TeV).

Considering different production channels for electrons and photons, however,
has a dramatic effect on the resulting IC emission (see Figure 4(b)). The very
hard photon and electron production associated to the $\tau\overline{\tau}$
channel results in a high peak energy and reduction in low energy IC emission
in this case. The soft spectra produced in the $b\overline{b}$ channel result
in correspondingly soft IC emission at high energies. In all cases however,
the IC emission represents a significant low energy shoulder to the total
gamma-ray emission of the GC region due to WIMP annihilation.

(a) WIMP mass dependence of IC emission

(b) Decay channel dependence of IC emission

Figure 4: Spectrum of direct and inverse Compton photons from WIMP
annihilation in the GC normalised to the total luminosity of the direct
photons. The solid lines show the spectrum for direct photons while the dashed
lines show the spectra of IC photons. In the left panel the spectra are shown
for different WIMP masses in the $W$ annihilation channel, while the right
panel shows the spectra for different annihilation channels and a WIMP mass of
1 TeV.

The stability of the IC component under a change in the model is illustrated
in Figure 5. The different colours correspond to different energy densities of
the magnetic field and the radiation field model. They are varied by a factor
of 2 from the nominal value to illustrate the uncertainty associated with
these components.

Figure 5: Photon spectrum from WIMP annihilation to $W$ bosons in the GC
normalised to the total luminosity of the direct photon component (solid
line). The dashed lines show the IC component of the spectrum for a variation
of the magnetic field strength and the radiation field of the model. The WIMP
mass is set to 1 TeV and the model was evaluated at a distance of 100 pc from
the GC and evolved for $10^{6}$ years.

##  4 Discussion

Based on the results presented here, the presence of significant additional
low energy photons associated to IC emission in the GC from annihilating Dark
Matter appears to be a robust prediction for $\sim$TeV mass WIMPs. As such it
should be incorporated in future experimental searches. A particular concern
is that searches increasingly make use of a template spectrum and a ‘3D’
likelihood fitting approach, and in neglecting the IC component may reach
misleading conclusions.

As mentioned above, there is a competition between hardening production
spectra and increasing KN suppression as WIMP masses get heavier. Figure 6(a)
shows the overall power input in to electrons (solid lines) and subsequently
in to IC emission (dashed lines), relative to the power of the direct photon
component. The three different decay channels $W\overline{W}$,
$\tau\overline{\tau}$ and $b\overline{b}$ are shown. The figure shows that not
all the power of the electrons are put into the IC emission and that for
extremely heavy ($\sim$30 TeV) WIMPs, the IC power starts to be quite strongly
suppressed. Still, the effect of IC emission is certainly relevant over a very
wide mass range in all production channels.

The impact that the IC emission can have on detectability is strongly
dependent on the energy-dependent sensitivity of the gamma-ray detector
considered. For Fermi-LAT the energy-flux sensitivity reaches a maximum in the
GeV region, which is why the IC component does have a considerable impact on
the expected signal counts for heavy WIMPs. Whilst a complete DM search with
Fermi-LAT, including the evaluation of background effects, is beyond the scope
of this work, we estimate the increase in the signal size seen in Fermi-LAT
with and without the IC component, in the case of heavy WIMPs. To this end, we
parameterise the acceptance of the Fermi-LAT Pass 8 data Release 3
(UltraCleanVeto) Version 3 [48] and estimate the ratio of expected photon
counts with and without taking the IC component into account. This ratio is
shown as a function of the WIMP mass in Figure 6(b). The $W\overline{W}$,
$\tau\overline{\tau}$ and $b\overline{b}$ annihilation channels are shown in
black, purple and violet, respectively. The solid lines show the expected
count ratio for an energy threshold of 1 GeV and the dotted lines correspond
to an energy threshold of 10 GeV. The figure shows that including the IC
component of the WIMP signal, increases the expected signal yield for Fermi-
LAT significantly – in the case of heavy WIMPs decaying to $\tau$ leptons by
more than a factor 100. Therefore, we conclude firstly that the current Fermi-
LAT limits that disregard the IC component are too pessimistic and would
improve when adding the IC component to the signal model. Secondly, due to the
existence of the IC component, the Fermi-LAT data could confirm or challenge a
potential claim of heavy WIMP observation from future ground-based gamma-ray
telescopes, making the data an important asset for DM searches for many years
to come.

(a) Relative production efficiency

(b) Signal count ratio with/without IC photons

Figure 6: Left: Overall power input to electrons from WIMP annihilation
relative to the overall power input to direct photons (solid lines). The
dashed lines show the same quantity for the resulting IC emission. Right:
Ratio of the photon counts of a DM signal with and without taking the IC
component into consideration as a function of the WIMP mass. The solid
(dotted) lines correspond to an energy threshold of the experiment of 1 GeV
(10 GeV). Three different WIMP annihilation channels $W\overline{W}$,
$\tau\overline{\tau}$ and $b\overline{b}$ are shown in black, purple and
violet, respectively.

The impact of the IC component of the WIMP signal on ground based gamma-ray
observatories is expected to be smaller than the effect it has for Fermi-LAT,
except in the case of extremely high mass WIMPs.

This work provides a first estimate of the IC component of the gamma-ray
signal from WIMP annihilation in the GC. More sophisticated models of the GC
environment should be considered in the future, as well as the intermediate
regime ($\sim$1-10 GeV electrons) where particles are transported significant
distances before cooling, with an impact on the measured emission profile.

Finally we note that any claim of a gamma-ray WIMP annihilation signal must be
tested against radio-mm limits on the accompanying synchrotron emission.
Although, as Figure 2(b) shows, synchrotron constraints for a thermal relic
are expected to be relatively weak unless the GC magnetic field strength is
much higher than the model adopted here.

##  5 Conclusions

The gamma ray spectrum of WIMP annihilation in the centre of our galaxy
consists of two components originating from direct and indirect photons
respectively. We have shown that the latter component is dominated by IC
emission for heavy WIMPs as the electrons produced by the annihilation radiate
locally, at least at sufficiently high energies. Whilst uncertainties related
to transport and cooling are significant, they do not exceed the fundamental
uncertainty due to the currently poorly constrained DM halo profile of the
Milky Way. We evaluated the robustness of the IC component for different
timescales, distances from the GC, WIMP masses and production channels. All in
all, we found that the IC component is significant and should not be neglected
in WIMP searches.

In addition, we evaluated the effect the IC component has on the detectability
of heavy ($\sim$TeV) WIMPs with Fermi-LAT. Including the IC component
significantly strengthens the WIMP signal and revised limits from the Fermi
collaboration incorporating IC emission are desirable.

In order to solve the long standing mystery of the particle nature of DM all
aspects of the signal should be taken into account. We therefore urge that
also the indirect component of the gamma-ray spectrum should be considered in
the next generation of Galactic WIMP searches with gamma-ray observatories.
This is especially true in case a WIMP signal is observed, as the shape of the
spectrum will give crucial hints on the nature of the particle.

## Acknowledgements

J.D., is funded by the Research Council of Norway, project number 301718 and
grateful for the hospitality of the Non-Thermal Astrophysics division at MPIK.

## References

  * [1] G. Arcadi, M. Dutra, P. Ghosh, et al., The waning of the WIMP? a review of models, searches, and constraints, The European Physical Journal C 78 (3) (Mar 2018).  doi:10.1140/epjc/s10052-018-5662-y.
  * [2] P. Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, et al., Planck 2018 results. VI. Cosmological parameters, A&A641 (2020) A6.  arXiv:1807.06209, doi:10.1051/0004-6361/201833910.
  * [3] R. Massey, J. Rhodes, R. Ellis, et al., Dark matter maps reveal cosmic scaffolding, Nature 445 (7125) (2007) 286–290.  doi:10.1038/nature05497.
  * [4] Y. Sofue, V. Rubin, Rotation curves of spiral galaxies, Annual Review of Astronomy and Astrophysics 39 (1) (2001) 137–174.  doi:10.1146/annurev.astro.39.1.137.
  * [5] ATLAS Collaboration, Search for associated production of a z boson with an invisibly decaying higgs boson or dark matter candidates at $sqrt{s}=13\,$tev with the ATLAS detector, Physics Letters B 829 (2022) 137066.  doi:10.1016/j.physletb.2022.137066.
  * [6] CMS Collaboration, Search for dark matter produced in association with a leptonically decaying $z$ boson in proton-proton collisions at $sqrt{s}=13\,$tev, Eur. Phys. J. C (Jan 2021).  doi:10.1140/epjc/s10052-020-08739-5.
  * [7] H. E. S. S. Collaboration, H. Abdalla, F. Aharonian, F. A. Benkhali, et al., Search for dark matter annihilation signals in the h.e.s.s. inner galaxy survey (2022).  doi:10.48550/ARXIV.2207.10471.
  * [8] V. Acciari, S. Ansoldi, L. Antonelli, et al., Combined searches for dark matter in dwarf spheroidal galaxies observed with the magic telescopes, including new data from coma berenices and draco, Physics of the Dark Universe 35 (2022) 100912.  doi:https://doi.org/10.1016/j.dark.2021.100912.
  * [9] S. Archambault, A. Archer, W. Benbow, et al., Dark matter constraints from a joint analysis of dwarf spheroidal galaxy observations with VERITAS, Physical Review D 95 (8) (apr 2017).  doi:10.1103/physrevd.95.082001.
  * [10] A. Acharyya, R. Adam, C. Adams, et al., Sensitivity of the cherenkov telescope array to a dark matter signal from the galactic centre, Journal of Cosmology and Astroparticle Physics 2021 (01) (2021) 057–057.  doi:10.1088/1475-7516/2021/01/057.
  * [11] A. Abeysekara, A. Albert, R. Alfaro, et al., A search for dark matter in the galactic halo with HAWC, Journal of Cosmology and Astroparticle Physics 2018 (02) (2018) 049–049.  doi:10.1088/1475-7516/2018/02/049.
  * [12] D.-Z. He, X.-J. Bi, S.-J. Lin, et al., Prospect for dark matter annihilation signatures from gamma-ray observation of dwarf galaxies by LHAASO, Physical Review D 100 (8) (oct 2019).  doi:10.1103/physrevd.100.083003.
  * [13] A. Viana, H. Schoorlemmer, A. Albert, et al., Searching for dark matter in the Galactic halo with a wide field of view TeV gamma-ray observatory in the Southern Hemisphere, jcap 2019 (12) (2019) 061.  arXiv:1906.03353, doi:10.1088/1475-7516/2019/12/061.
  * [14] M. Ackermann, M. Ajello, A. Albert, et al., The fermi galactic center GeV excess and implications for dark matter, The Astrophysical Journal 840 (1) (2017) 43.  doi:10.3847/1538-4357/aa6cab.
  * [15] HESS Collaboration, A. Abramowski, F. Aharonian, F. A. o. Benkhali, Acceleration of petaelectronvolt protons in the Galactic Centre, Nature531 (7595) (2016) 476–479.  arXiv:1603.07730, doi:10.1038/nature17147.
  * [16] K. Ferrière, Interstellar magnetic fields in the Galactic center region, A&A505 (3) (2009) 1183–1198.  arXiv:0908.2037, doi:10.1051/0004-6361/200912617.
  * [17] M. Guenduez, J. Becker Tjus, K. Ferrière, et al., A novel analytical model of the magnetic field configuration in the Galactic center, A&A644 (2020) A71.  arXiv:1906.05211, doi:10.1051/0004-6361/201936081.
  * [18] J. Bland-Hawthorn, M. Cohen, The Large-Scale Bipolar Wind in the Galactic Center, ApJ582 (1) (2003) 246–256.  arXiv:astro-ph/0208553, doi:10.1086/344573.
  * [19] G. Ponti, F. Hofmann, E. Churazov, et al., An X-ray chimney extending hundreds of parsecs above and below the Galactic Centre, Nature567 (7748) (2019) 347–350.  arXiv:1904.05969, doi:10.1038/s41586-019-1009-6.
  * [20] I. Cholis, G. Dobler, D. P. Finkbeiner, et al., The case for a 700+ gev wimp: Cosmic ray spectra from PAMELA, fermi, and ATIC, Physical Review D 80 (12) (dec 2009).  doi:10.1103/physrevd.80.123518.
  * [21] M. Cirelli, G. Corcella, A. Hektor, et al., PPPC 4 DM ID: a poor particle physicist cookbook for dark matter indirect detection, Journal of Cosmology and Astroparticle Physics 2011 (03) (2011) 051–051.  doi:10.1088/1475-7516/2011/03/051.
  * [22] M. Cirelli, P. D. Serpico, G. Zaharijas, Bremsstrahlung gamma rays from light dark matter, Journal of Cosmology and Astroparticle Physics 2013 (11) (2013) 035–035.  doi:10.1088/1475-7516/2013/11/035.
  * [23] T. Lacroix, C. Bœhm, J. Silk, Fitting the fermi-LAT GeV excess: On the importance of including the propagation of electrons from dark matter, Physical Review D 90 (4) (aug 2014).  doi:10.1103/physrevd.90.043508.
  * [24] J. Zhang, X.-J. Bi, J. Liu, et al., Discriminating different scenarios to account for the cosmic e± excess by synchrotron and inverse compton radiation, Physical Review D 80 (2) (jul 2009).  doi:10.1103/physrevd.80.023007.
  * [25] A. E. Egorov, J. M. Gaskins, E. Pierpaoli, et al., Dark matter implications of the wmap-planck haze, Journal of Cosmology and Astroparticle Physics 2016 (03) (2016) 060.  doi:10.1088/1475-7516/2016/03/060.
  * [26] M. Ackermann, A. Albert, W. B. Atwood, et al., The Spectrum and Morphology of the Fermi Bubbles, ApJ793 (1) (2014) 64.  arXiv:1407.7905, doi:10.1088/0004-637X/793/1/64.
  * [27] P. Predehl, R. A. Sunyaev, W. Becker, et al., Detection of large-scale X-ray bubbles in the Milky Way halo, Nature588 (7837) (2020) 227–231.  arXiv:2012.05840, doi:10.1038/s41586-020-2979-0.
  * [28] C. C. Popescu, R. Yang, R. J. Tuffs, et al., A radiation transfer model for the Milky Way: I. Radiation fields and application to high-energy astrophysics, Monthly Notices of the Royal Astronomical Society 470 (3) (2017) 2539–2558.  arXiv:https://academic.oup.com/mnras/article-pdf/470/3/2539/18245516/stx1282.pdf, doi:10.1093/mnras/stx1282.
  * [29] H. E. S. S. Collaboration, H. Abdalla, A. Abramowski, F. Aharonian, et al., Characterising the VHE diffuse emission in the central 200 parsecs of our Galaxy with H.E.S.S., A&A612 (2018) A9.  arXiv:1706.04535, doi:10.1051/0004-6361/201730824.
  * [30] MAGIC Collaboration, V. A. Acciari, S. Ansoldi, L. A. Antonelli, et al., MAGIC observations of the diffuse $\gamma$-ray emission in the vicinity of the Galactic center, A&A642 (2020) A190.  arXiv:2006.00623, doi:10.1051/0004-6361/201936896.
  * [31] C. B. Adams, W. Benbow, A. Brill, et al., VERITAS Observations of the Galactic Center Region at Multi-TeV Gamma-Ray Energies, ApJ913 (2) (2021) 115\.  arXiv:2104.12735, doi:10.3847/1538-4357/abf926.
  * [32] E. M. Di Teodoro, N. M. McClure-Griffiths, F. J. Lockman, et al., Blowing in the Milky Way Wind: Neutral Hydrogen Clouds Tracing the Galactic Nuclear Outflow, ApJ855 (1) (2018) 33.  arXiv:1802.02152, doi:10.3847/1538-4357/aaad6a.
  * [33] R. Jansson, G. R. Farrar, A NEW MODEL OF THE GALACTIC MAGNETIC FIELD, The Astrophysical Journal 757 (1) (2012) 14.  doi:10.1088/0004-637x/757/1/14.
  * [34] R. Jansson, G. R. Farrar, THE GALACTIC MAGNETIC FIELD, The Astrophysical Journal 761 (1) (2012) L11.  doi:10.1088/2041-8205/761/1/l11.
  * [35] A. Scherer, J. Cuadra, F. E. Bauer, Galactic center gamma-ray production by cosmic rays from stellar winds and Sgr A East, A&A659 (2022) A105.  arXiv:2202.05429, doi:10.1051/0004-6361/202142401.
  * [36] R. M. Kulsrud, Plasma Physics for Astrophysics, 2005.
URL https://ui.adsabs.harvard.edu/abs/2005ppa..book.....K

  * [37] M. C. Beck, A. M. Beck, R. Beck, et al., New constraints on modelling the random magnetic field of the MW, J. Cosmology Astropart. Phys.2016 (5) (2016) 056.  arXiv:1409.5120, doi:10.1088/1475-7516/2016/05/056.
  * [38] J. Hahn, GAMERA \- a new modeling package for non-thermal spectral modeling, in: 34th International Cosmic Ray Conference (ICRC2015), Vol. 34 of International Cosmic Ray Conference, 2015, p. 917.
URL https://ui.adsabs.harvard.edu/abs/2015ICRC...34..917H

  * [39] M. Breuhaus, B. Reville, J. A. Hinton, Pulsar wind nebula origin of the LHAASO-detected ultra-high energy $\gamma$-ray sources, Astronomy & Astrophysics 660 (2022) A8.  doi:10.1051/0004-6361/202142097.
  * [40] J. F. Navarro, C. S. Frenk, S. D. M. White, The structure of cold dark matter halos, The Astrophysical Journal 462 (1996) 563.  doi:10.1086/177173.
  * [41] J. F. Navarro, C. S. Frenk, S. D. M. White, A universal density profile from hierarchical clustering, The Astrophysical Journal 490 (2) (1997) 493–508.  doi:10.1086/304888.
  * [42] J. Einasto, On the Construction of a Composite Model for the Galaxy and on the Determination of the System of Galactic Parameters, Trudy Astrofizicheskogo Instituta Alma-Ata 5 (1965) 87–100.
URL https://ui.adsabs.harvard.edu/abs/1965TrAlm...5...87E

  * [43] H. Abdallah, A. Abramowski, F. Aharonian, et al., Search for dark matter annihilations towards the inner galactic halo from 10 years of observations with h.e.s.s., Physical Review Letters 117 (11) (sep 2016).  doi:10.1103/physrevlett.117.111301.
  * [44] M. Portail, O. Gerhard, C. Wegg, et al., Dynamical modelling of the galactic bulge and bar: the milky way's pattern speed, stellar and dark matter mass distribution, Monthly Notices of the Royal Astronomical Society 465 (2) (2016) 1621–1644.  doi:10.1093/mnras/stw2819.
  * [45] Website: Pppc 4 dm id - a poor particle physicist cookbook for dark matter indirect detection, http://www.marcocirelli.net/PPPC4DMID.html, accessed: 2022-07-07.
  * [46] P. Ciafaloni, D. Comelli, A. Riotto, et al., Weak corrections are relevant for dark matter indirect detection, Journal of Cosmology and Astroparticle Physics 2011 (03) (2011) 019–019.  doi:10.1088/1475-7516/2011/03/019.
  * [47] P. A. R. Ade, N. Aghanim, M. Arnaud, et al., Planck intermediate results, Astronomy & Astrophysics 554 (2013) A139.  doi:10.1051/0004-6361/201220271.
  * [48] Website: Fermi lat performance, https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm, accessed: 2022-07-08.


# Quadratic minimization: from conjugate gradient to an adaptive Heavy-ball
method with Polyak step-sizes

Baptiste Goujaud111CMAP, École Polytechnique, Institut Polytechnique de Paris.
<EMAIL_ADDRESS>Adrien Taylor222INRIA, École Normale
Supérieure, CNRS, PSL Research University, Paris<EMAIL_ADDRESS>Aymeric Dieuleveut333CMAP, École Polytechnique, Institut Polytechnique de
Paris<EMAIL_ADDRESS>
###### Abstract

In this work, we propose an adaptive variation on the classical Heavy-ball
method for convex quadratic minimization. The adaptivity crucially relies on
so-called “Polyak step-sizes”, which consists in using the knowledge of the
optimal value of the optimization problem at hand instead of problem
parameters such as a few eigenvalues of the Hessian of the problem. This
method happens to also be equivalent to a variation of the classical conjugate
gradient method, and thereby inherits many of its attractive features,
including its finite-time convergence, instance optimality, and its worst-case
convergence rates.

The classical gradient method with Polyak step-sizes is known to behave very
well in situations in which it can be used, and the question of whether
incorporating momentum in this method is possible and can improve the method
itself appeared to be open. We provide a definitive answer to this question
for minimizing convex quadratic functions, a arguably necessary first step for
developing such methods in more general setups.

##  1 Introduction

Consider the convex quadratic minimization problem in the form

$\min_{x\in\mathbb{R}^{d}}\left\\{f(x)\triangleq\frac{1}{2}\langle
x,\,Hx\rangle+\langle h,\,x\rangle\triangleq\frac{1}{2}\langle
x-x_{\star},\,H(x-x_{\star})\rangle+f_{\star}\right\\}$ (1)

where $H\succcurlyeq 0$ is a symmetric positive semi-definite matrix, and we
denote $f_{\star}$ the minimum value of $f$ (a few instances of such problems
are presented in Section 3). In the context of large-scale optimization (i.e.
$d\gg 1$), we are often interested in using first-order iterative methods for
solving eq. 1. There are many known and celebrated iterative methods for
solving such quadratic problems, including conjugate gradient, Heavy-ball
methods (a.k.a., Polyak momentum), Chebyshev methods, and gradient descent.
Each of those methods having different specifications, the choice of the
method largely depends on the application at hand. In particular, a typical
drawback of momentum-based methods is that they generally require the
knowledge of some problem parameters (such as extreme values of the spectrum
of $H$). This problem is typically not as critical for simpler gradient
descent schemes with no momentum, although it generally still requires some
knowledge on problem parameters if we want to avoid using linesearch-based
strategies. This limitation motivates the search for adaptive strategies,
fixing step-size using past observations about the problem at hand. In the
context of (sub)gradient descent, a famous adaptive strategy is the so-called
Polyak step-size, which relies on the knowledge of the optimal value
$f_{\star}$:

$x_{t+1}=x_{t}-\gamma_{t}\nabla f(x_{t}),\quad\text{with}\
\gamma_{t}=\frac{f(x_{t})-f_{\star}}{\|\nabla f(x_{t})\|^{2}}.$ (2)

Polyak steps were originally proposed in Polyak (1987) for nonsmooth convex
minimization; it is also discussed in Boyd et al. (2003) and a few variants
are proposed by, e.g., Barré et al. (2020); Loizou et al. (2021); Gower et al.
(2022) including for stochastic minimization. In terms of speed, this strategy
(and variants) enjoy similar theoretical convergence properties as those for
gradient descent. This methods appears to perform quite well in applications
where $f_{\star}$ can be efficiently estimated—see, e.g., Hazan and Kakade
(2019) for an adaptation of the method for estimating it online (although
Hazan and Kakade (2019) contains a number of mistakes, the approach can be
corrected to achieve the claimed target). Therefore, a remaining open question
in this context is whether the performances of this method can be improved by
incorporating momentum in it. A first answer to this question was provided by
Barré et al. (2020), although it is not clear that it can match the same
convergence properties as optimal first-order methods.

In this work, we answer this question for the class of quadratic problems. In
short, it turns out that the following conjugate gradient-like iterative
procedure

$\displaystyle x_{t+1}=\arg\\!\min_{x}\left\\{\|x-x_{\star}\|^{2}\text{ s.t.
}x\in x_{0}+\mathrm{span}\\{\nabla f(x_{0}),\nabla f(x_{1}),\ldots,\nabla
f(x_{t})\\}\right\\},$ (3)

can be rewritten exactly as an Heavy-ball type method whose parameters are
chosen adaptively using the value of $f_{\star}$. This might come as a
surprise, as the iteration eq. 3 might seem impractical due to its formulation
relying on the knowledge of $x_{\star}$. More precisely, eq. 3 is exactly
equivalent to:

$\displaystyle x_{t+1}$ $\displaystyle=x_{t}-(1+m_{t})\times h_{t}\nabla
f(x_{t})+m_{t}(x_{t}-x_{t-1}),$ (4)

with parameters

$\displaystyle\forall t\geqslant 0,\quad h_{t}$ $\displaystyle\triangleq$
$\displaystyle\frac{2(f(x_{t})-f_{\star})}{\|\nabla f(x_{t})\|^{2}},\quad$ (5)
$\displaystyle m_{0}\triangleq 0\ \ \text{and}\ \ \forall t\geqslant 1,\quad
m_{t}$ $\displaystyle\triangleq$
$\displaystyle\frac{-(f(x_{t})-f_{\star})\langle\nabla f(x_{t}),\nabla
f(x_{t-1})\rangle}{(f(x_{t-1})-f_{\star})\|\nabla
f(x_{t})\|^{2}+(f(x_{t})-f_{\star})\langle\nabla f(x_{t}),\nabla
f(x_{t-1})\rangle}\ .$ (6)

In eq. 4, $m_{t}$ corresponds to the momentum coefficient and $h_{t}$ to a
step-size. With the tuning of eq. 5, this step-size is twice the Polyak step-
size in eq. 2. This Heavy-ball momentum method with Polyak step-sizes is
summarized in Algorithm 1 and illustrated in Figure 1. Due to its equivalence
with eq. 3, the Heavy-ball method eq. 4 inherits nice advantageous properties
of conjugate gradient-type methods, including:

  1. (i)

finite-time convergence: the problem eq. 1 is solved exactly after at most $d$
iterations,

  2. (ii)

instance optimality: for all $H\succcurlyeq 0$, no first-order method
satisfying $x_{t+1}\in x_{0}+\mathrm{span}\\{\nabla f(x_{0}),\ldots,\nabla
f(x_{t})\\}$ results in a smaller $\|x_{t}-x_{\star}\|$,

  3. (iii)

it inherits optimal worst-case convergence rates on quadratic functions.

Of course, a few of those points needs to be nuanced in practice due to finite
precision arithmetic. The equivalence between eq. 3 and eq. 4 is formally
stated in the following theorem.

###### Theorem 1.1

Let $(x_{t})_{t\in\mathbb{N}}$ be the sequence defined by the recursion eq. 3,
namely such that for any $t$, $x_{t+1}$ is the Euclidean projection of
$x_{\star}$ onto the affine subspace $x_{0}+\mathrm{span}\\{\nabla
f(x_{0}),\nabla f(x_{1}),\ldots,\nabla f(x_{t})\\}$. Then
$(x_{t})_{t\in\mathbb{N}}$ is the sequence generated by Algorithm 1.

Algorithm 1 Adaptive Heavy-ball algorithm

Input $T$ and $f:x\mapsto f(x)\triangleq\frac{1}{2}\langle
x-x_{\star},\,H(x-x_{\star})\rangle+f_{\star}$
Initialize $x_{0}\in\mathbb{R}^{d}$, $m_{0}=0$
for  _$t=0\cdots T-1$_ do

$h_{t}=\frac{2(f(x_{t})-f_{\star})}{\|\nabla f(x_{t})\|^{2}}$
$x_{t+1}=x_{t}-(1+m_{t})\times h_{t}\nabla f(x_{t})+m_{t}(x_{t}-x_{t-1})$
$m_{t+1}=\frac{-(f(x_{t+1})-f_{\star})\langle\nabla f(x_{t+1}),\nabla
f(x_{t})\rangle}{(f(x_{t})-f_{\star})\|\nabla
f(x_{t+1})\|^{2}+(f(x_{t+1})-f_{\star})\langle\nabla f(x_{t+1}),\nabla
f(x_{t})\rangle}$

end for

Result: $x_{T}$

Figure 1: Comparison in semi-log scale over 50 iterations of different first-
order methods applied on a 25-dimensional quadratic objective with condition
number $10$. GD with constant step-size, GD with Polyak step-size and GD with
variant of Polyak step-size refer to the GD method tuned respectively with the
step-size $\gamma=2/(L+\mu)$, $\gamma_{t}=(f(x_{t})-f_{\star})/\|\nabla
f(x_{t})\|^{2}$ and $\gamma_{t}=2(f(x_{t})-f_{\star})/\|\nabla
f(x_{t})\|^{2}$. HB with constant tuning is the HB method tuned with constant
parameters $\gamma_{t}=(2/(\sqrt{L}+\sqrt{\mu}))^{2}$ and
$m_{t}=((\sqrt{L}-\sqrt{\mu})(\sqrt{L}+\sqrt{\mu}))^{2}$ while HB with Polyak
step-size based tuning refers to Algorithm 1.

1.1 turns out to be a particular case of a more general result stating that
the iterates of any conjugate gradient-type method described with a polynomial
$Q$ as

$x_{t+1}=\mathrm{argmin}_{x}\left\\{\langle
x-x_{\star},\,Q(H)(x-x_{\star})\rangle\text{ s.t. }x\in
x_{0}+\mathrm{span}\\{\nabla f(x_{0}),\ldots,\nabla f(x_{t})\\}\right\\},$
(Q-minimization)

are equivalently written in terms of an adaptive Heavy-ball iteration. In
particular, eq. 3 corresponds to eq. Q-minimization with $Q(x)=1$. Similarly,
classical conjugate gradient method corresponds to eq. Q-minimization with
$Q(x)=x$ (this fact is quite famous, see, e.g., Nocedal and Wright (1999)). We
were surprised not to find this general result written as is in the
literature, and we therefore provide it in Section 2. The key point of this
work is that the equivalent Heavy-ball reformulation of eq. 3 can be written
in terms of $f_{\star}$, thereby obtaining a momentum-based Polyak step-size.

#### Notations.

We denote $\preccurlyeq$ the order between symmetric matrices;
$\operatorname{Sp}{H}$ the spectrum of the matrix $H$, namely its set of
eigenvalues; $\mathbb{R}_{d}[X]$ the set of polynomials with degree at most
$d$.

###  1.1 Preliminary material

#### Worst-case optimality.

Solving eq. 1 is a very important problem and several methods have been
proposed to achieve this goal. They are compared with each other through
notions of performance. This consists in evaluating the precision of an
algorithm over the functions of a given class after a given number $T$ of
iterations. The main framework is _worst-case analysis_ and the precision is
the value of a given metric, e.g. the distance of the last iterate to optimum
$\|x_{T}-x_{\star}\|$, the function value of the last iterate
$f(x_{T})-f(x_{\star})$, or its gradient norm $\|\nabla f(x_{T})\|$. The
_worst-case analysis_ framework consists in finding the guarantees of a method
that hold for each and every functions of a given class, as for instance the
class of $L$-smooth $\mu$-strongly convex quadratic functions described as
quadratic functions verifying $\mu I\preccurlyeq H\preccurlyeq LI$ for given
$0<\mu\leqslant L$. The gradient descent (GD) method characterized by the
update

$x_{t+1}=x_{t}-\gamma_{t}\nabla f(x_{t})$ (7)

therefore verifies $\|x_{t}-x_{\star}\|=O((\frac{L-\mu}{L+\mu})^{t})$ on all
such functions for $\gamma_{t}=\frac{2}{L+\mu}$. Thanks to a relationship with
polynomial analysis, Golub and Varga (1961) proved that the Chebyshev method,
described as

$x_{t+1}=x_{t}-\gamma_{t}\nabla f(x_{t})+m_{t}(x_{t}-x_{t-1}),$ (8)

for a well chosen tuning of the parameters $\gamma_{t}$ and $m_{t}$
($m_{t}=(\frac{\sqrt{L}-\sqrt{\mu}}{\sqrt{L}+\sqrt{\mu}})^{2}\frac{1+((\sqrt{L}-\sqrt{\mu})/(\sqrt{L}+\sqrt{\mu}))^{2(t-1)}}{1+((\sqrt{L}-\sqrt{\mu})/(\sqrt{L}+\sqrt{\mu}))^{2(t+1)}}$,
$\gamma_{t}=\frac{2}{L+\mu}(1+m_{t})$), is _worst-case optimal_ on this class
of function, achieving the guarantee
$\|x_{t}-x_{\star}\|=O((\frac{\sqrt{L}-\sqrt{\mu}}{\sqrt{L}+\sqrt{\mu}})^{t})$
(often referred to as “acceleration”). Methods based on this two-term
recursion are called “Heavy-ball” or “Polyak momentum” (Polyak, 1964). In
particular, the stationary regime of the Chebyshev method is the Heavy-ball
(HB) method tuned with
$m_{t}=(\frac{\sqrt{L}-\sqrt{\mu}}{\sqrt{L}+\sqrt{\mu}})^{2}$ and
$\gamma_{t}=\frac{2}{L+\mu}(1+m_{t})=(\frac{2}{\sqrt{L}+\sqrt{\mu}})^{2}$ and
achieves the worst-case guarantee
$\|x_{t}-x_{\star}\|=O(t(\frac{1-\sqrt{\kappa}}{1+\sqrt{\kappa}})^{t})$, close
to the optimal one achieved by the Chebyshev method. Note that eq. 8 is
another formulation of eq. 4 where $\gamma_{t}=(1+m_{t})\times h_{t}$. In all
the aforementioned tuning, $h_{t}$ has the same value: $h_{t}=\frac{2}{L+\mu}$
(see Section 4 for more detailed discussion on this).

#### Span of gradients and Krylov subspaces.

All methods described above can be defined using a recursion:

$x_{t}=x_{0}-\sum_{i=0}^{t-1}\gamma_{i}^{(t)}\nabla f(x_{i})$ (9)

for some sequence $(\gamma_{i}^{(t)})_{i\in\llbracket 0,t-1\rrbracket}$ as
suggested in (Nemirovski, 1992, 1994). Note that the recursion eq. 9 can also
be explicitly written as $x_{t}=x_{0}-H\sum_{i=0}^{t-1}\gamma_{i}^{(t)}x_{i}$,
and therefore, $x_{t}-x_{0}\in H\mathrm{span}(\\{x_{i}\\}_{i\in\llbracket
0,t-1\rrbracket})$. We deduce by recursion that $x_{t}-x_{0}\in
H\mathcal{K}_{t}(H,x_{0})$ where
$\mathcal{K}_{t}(H,x_{0})\triangleq\mathrm{span}(\\{H^{i}x_{0}\\}_{i\in\llbracket
0,t-1\rrbracket})$ is called _order-t Krylov subspace generated by $H$ and
$x_{0}$_. This creates a link between first-order algorithms and polynomials,
summarized in the following lemma (which is implicitely used in Golub and
Varga (1961) and formally stated, e.g., in (Goujaud et al., 2022, Proposition
4.1)).

###### Lemma 1.2

. Let $f\in\mathcal{C}_{\Lambda}$. The iterates $x_{t}$ satisfy

$x_{t}\in x_{0}+\mathrm{span}\\{\nabla f(x_{0}),\ldots,\nabla
f(x_{t-1})\\}\,,$ (10)

where $x_{0}$ is the initial approximation of $x_{\star}$, if and only if
there exists a sequence of polynomials $(P_{t})_{t\in\mathbb{N}}$, each of
degree at most 1 more than the highest degree of all previous polynomials and
$P_{0}$ of degree 0 (hence the degree of $P_{t}$ is at most $t$), such that

$\vphantom{\sum_{i}}\forall\;t,\quad
x_{t}-x_{\star}=P_{t}(H)(x_{0}-x_{\star}),\quad P_{t}(0)=1\,.$ (11)

Similar to the way we use this technique below, this lemma has already been
extensively used to design methods; see, e.g., d’Aspremont et al. (2021,
Chapter 1) or the blog post by Pedregosa (2020) for gentle introductions to
this technique. For instance, we can use this technique for optimizing the
step-size of the gradient method, or to derive the Chebyshev method, which
optimizes the worst-case on the class of smooth and strongly convex quadratic
functions (see Fischer, 2011). More recently, Goujaud et al. (2022) used it to
derive a method which can take advantage of a possible gap in the spectrum of
$H$. This approach has also been used for other applications such as
accelerated gossip algorithms (Berthier et al., 2020).

#### Adaptive methods.

In 1.2, while $P_{t}(H)$ is a polynomial evaluated on the matrix $H$, its
scalar coefficients might or might not depend on $H$. If they depend on $H$,
we say that the associated method is adaptive. Non-adaptive methods suffer
from two main drawbacks: (i) they use the same parameters for all the
functions within the class of problems, not taking advantage of the observed
quantities; (ii) the underlying parameters must scale with the function class
parameters, and therefore depends on the values of $L$ and $\mu$, which are
generally difficult to estimate (and actually do not correspond to first-order
information, as they rely on the Hessian of the function at hand). Ultimately,
adaptive methods aim at solving those issues by choosing parameters (step-
size, momentum, etc.) on the fly.

#### Polyak steps.

It is straightforward that a gradient descent update verifies on any convex
function $f$ that
$\|x_{t+1}-x_{\star}\|^{2}\leqslant\|x_{t}-x_{\star}\|^{2}-2\gamma_{t}(f-f_{\star})+\gamma_{t}^{2}\|\nabla
f(x_{t})\|^{2}$. Polyak (1987) argues that, based on this inequality, the best
guaranteed progression is then achieved for
$\gamma_{t}=\frac{f(x_{t})-f_{\star}}{\|\nabla f(x_{t})\|^{2}}$. This choice
is called “Polyak step-size” and has been studied intensively even recently
Loizou et al. (2021); Gower et al. (2022).

Other variants of the latter have been proposed. For instance (Barré et al.,
2020, Variant 1) suggested the step-size
$\gamma_{t}=2\frac{f(x_{t})-f_{\star}}{\|\nabla f(x_{t})\|^{2}}$. This also
optimizes the exact progress of one gradient descent update over quadratics
realizing a projection of $x_{t}-x_{\star}$ over the orthogonal subspace of
$\nabla f(x_{t})$.

Therefore, the Polyak step-size strategy applied to the gradient descent
method achieves the same worst-case guarantee of the well tuned fixed step-
size gradient descent method, while not relying on Hessian information.
Moreover, due to its adaptivity to each function, and since generic functions
do not look like worst-cases, the Polyak step-size strategy applied to the
gradient descent method perform very well in practice (See Figure 1 and (Barré
et al., 2020, Figure 1)), sometimes even beating the well tuned non-adaptive
Heavy-ball method even if the worst-case guarantees are sorted in a different
order.

#### Instance-optimality.

While optimal worst-case method of the form of eq. 9 have been found with
predetermined parameters, it would be better to find a method under the form
of eq. 9 that is optimal (for some performance metric), not only on worst-case
analysis, but on each function individually, taking advantage of the
adaptivity of the parameters. The well-known conjugate gradient method
achieves this goal when the performance metric is the function value of the
last iterate. The MinRes method attacks the problem minimizing the gradient
norm of the last iterate.

#### Contributions.

In this work we derive iterative methods of the form of eq. 9 (which iterates
lie in the span of previously observed gradients) that are instance-optimal
for a variety of performance metrics. All those methods updates are variations
of the Heavy-ball two-term recursion eq. 8 with only parameters $\gamma_{t}$
and $m_{t}$ changing from one method to another. Finally, we show (see 1.1)
that for a well chosen yet classical performance metric, this associated
method Algorithm 1 is not relying on second-order information at all (not even
$L$ and $\mu$). Instead it is using a classical variant of the Polyak step-
size $\frac{2(f(x_{t})-f_{\star})}{\|\nabla f(x_{t})\|^{2}}$, providing an
answer to the question “can we accelerate methods with Polyak step-size?”.

###  1.2 Related works

Polyak step-sizes were proposed in (Polyak, 1987). Despite the dependency on
$f_{\star}$, the Polyak step-size is more studied theoretically and used in
practice due to its efficiency when applied to real world problems. Recent
works (e.g. Loizou et al., 2021; D’Orazio et al., 2021) argue that this
dependency is not a practical issue for many problems which we can assume
verify $f_{\star}=0$ (see Section 3). A few variants of the Polyak step-size
strategy were proposed by Barré et al. (2020), including a version
incorporating a Nesterov-type momentum Nesterov (1983), achieving a worst-case
guarantee of $\|x_{t}-x_{\star}\|^{2}=O((1-2(\mu/L)^{2/3})^{2t})$ over the
class of (non-necessarily quadratic) $L$-smooth $\mu$-strongly convex
functions, thereby improving over previous works on adaptive first-order
methods. However, the proposed method does not allow to remove the dependency
on $L$, and does not achieve the black-box complexity of smooth strongly
convex minimization Nesterov (2004). In (Loizou et al., 2021), the authors
study the stochastic Polyak step-size, whereas (D’Orazio et al., 2021) applies
it to Mirror descent.

Many alternative adaptive methods have been proposed in the past. Among them,
let us mention (Barzilai and Borwein, 1988) which introduced the so-called
Barzilai-Borwein step-size rule, and the more recent (Malitsky and Mishchenko,
2020) which developed a step-size policy that adapts to the local geometry
with convergence guarantees beyond quadratic minimization.

##  2 Main theorem

This section states and proves 1.1. In short, given a certain function $f$
(characterized by $H$ and $x_{\star}$ here) and a starting point $x_{0}$, we
search for an iterative procedure, possibly adaptive, verifying the
polynomial-based expression eq. 11 such that $x_{t}$ converges as fast as
possible to $x_{\star}$ for some predefined performance metric. Most classical
ways to measure the performance of such optimization schemes include the
distance to optimum $\|x_{t}-x_{\star}\|^{2}$, the function accuracy gap
$f(x_{t})-f_{\star}$, the squared gradient norm $\|\nabla f(x_{t})\|^{2}$, and
linear combinations of the former. Let us abstract those notions by denoting
the performance measure of choice by $\langle
x_{t}-x_{\star},\,Q(H)(x_{t}-x_{\star})\rangle$ (with $Q$ a predefined
polynomial that is positive on $\mathbb{R}_{>0}$). Then, we consider the
iterative scheme given by (Q-minimization).

$x_{t+1}=\mathrm{argmin}_{x}\left\\{\langle
x-x_{\star},\,Q(H)(x-x_{\star})\rangle\text{ s.t. }x\in
x_{0}+\mathrm{span}\\{\nabla f(x_{0}),\ldots,\nabla f(x_{t})\\}\right\\},$
(Q-minimization)

The next theorem provides an explicit instance-optimal method to solve eq.
Q-minimization.

######  Theorem 2.1 (Main)

The unique solution to eq. Q-minimization is given by the Heavy-ball procedure

$x_{t+1}=x_{t}-(1+m_{t})\times h_{t}\nabla f(x_{t})+m_{t}(x_{t}-x_{t-1})$ (12)

where

$\left\\{\begin{array}[]{ll}h_{t}&=\frac{\langle
x_{t}-x_{\star},\,HQ(H)(x_{t}-x_{\star})\rangle}{\langle
x_{t}-x_{\star},\,H^{2}Q(H)(x_{t}-x_{\star})\rangle};\\\
m_{t}&=\frac{-b_{t}h_{t}}{1+b_{t}h_{t}},\quad\mathrm{\leavevmode\nobreak\
with\leavevmode\nobreak\ }b_{t}=\frac{\langle
x_{t}-x_{\star},\,H^{2}Q(H)(x_{t-1}-x_{\star})\rangle}{\langle
x_{t-1}-x_{\star},\,HQ(H)(x_{t-1}-x_{\star})\rangle}.\end{array}\right.$ (13)

Remark that setting $Q(X)=X$ leads to a nice expression of the conjugate
gradient method. Indeed, setting $Q(X)=X$ corresponds to optimally minimizing
the excess loss $f(x_{t})-f_{\star}$.

As already known, the conjugate gradient method requires the knowledge of H
(or an Hessian vector product) to proceed. This is also a priori the case for
all other choices of $Q(\cdot)$. In the case of $Q(X)=1$, which corresponds to
minimizing the distance to the optimum (see eq. 3), we can use an alternate
writing making use of $f_{\star}$:

$\left\\{\begin{array}[]{ll}h_{t}&=\frac{2(f(x_{t})-f_{\star})}{\|\nabla
f(x_{t})\|^{2}}\\\
m_{t}&=\frac{-b_{t}h_{t}}{1+b_{t}h_{t}},\quad\mathrm{\leavevmode\nobreak\
with\leavevmode\nobreak\ }b_{t}=\frac{\left<\nabla f(x_{t}),\nabla
f(x_{t-1})\right>}{2(f(x_{t-1})-f_{\star})}.\end{array}\right.$ (14)

Proof.

#### Designing methods from the polynomial point of view.

As suggested by 1.2, we look for an iterative method which can be expressed in
the form $x_{t}-x_{\star}=P_{t}(H)(x_{0}-x_{\star})$, where $P_{t}$ is a
$t^{th}$ degree polynomial with $P_{t}(0)=1$. Furthermore, as we look for an
instance-optimal method, the latter polynomial must be instance-specific, and
the coefficients of $P_{t}$ should depend on $H$ (and should describe the
iterative procedure (Q-minimization)).

Recalling that $H$ is real symmetric matrix, we denote by
$\lambda\in\operatorname{Sp}(H)$ its eigenvalues and by $v_{\lambda}$ th
associated orthonormal basis of eigenvectors, leading to
$H=\sum_{\lambda\in\operatorname{Sp}(H)}\lambda v_{\lambda}v_{\lambda}^{T}$.
The quantity to be minimized can now be written as:

$\displaystyle\langle x_{t}-x_{\star},\,Q(H)(x_{t}-x_{\star})\rangle$
$\displaystyle\leavevmode\nobreak\ =\langle
x_{0}-x_{\star},\,P_{t}(H)^{T}Q(H)P_{t}(H)(x_{0}-x_{\star})\rangle$ (15)
$\displaystyle\leavevmode\nobreak\
=\sum_{\lambda\in\operatorname{Sp}(H)}Q(\lambda)P_{t}(\lambda)^{2}\langle
x_{0}-x_{\star},v_{\lambda}\rangle^{2}$ (16)
$\displaystyle\leavevmode\nobreak\
=\displaystyle\int_{\lambda\in\mathbb{R}^{+}}P_{t}(\lambda)^{2}\leavevmode\nobreak\
\mathrm{d\lambda_{Q}(\lambda)}$ (17)

with $\mathrm{\lambda_{Q}}$ the discrete measure
$\sum_{\lambda\in\operatorname{Sp}(H)}Q(\lambda)\langle
x_{0}-x_{\star},v_{\lambda}\rangle^{2}\leavevmode\nobreak\ \delta_{\lambda}$
(we sometimes use the shorthand notation $\int P_{t}^{2}\leavevmode\nobreak\
\mathrm{d\lambda_{Q}}$ for eq. 17 in what follows). It is clear that eq. 17 is
$0$ if and only if $P_{t}(\lambda)=0$ for all
$\lambda\in\operatorname{Sp}(H)$. As a consequence, we conclude that (i)
choosing the right sequence of polynomials leads to convergence in exactly
$|\operatorname{Sp}(H)|$ iterations, and (ii) $\langle
P^{(1)},P^{(2)}\rangle_{Q}\triangleq\int P^{(1)}P^{(2)}\leavevmode\nobreak\
\mathrm{d\lambda_{Q}}$ is an inner product on
$\mathbb{R}_{|\operatorname{Sp}(H)|-1}[X]$. We therefore want to solve

$\left\\{\begin{array}[]{ll}\underset{P_{t}\in\mathbb{R}_{t}[X]}{\mathrm{minimize\leavevmode\nobreak\
}}&\|P_{t}\|_{Q}^{2}\\\ \mathrm{subject\leavevmode\nobreak\
to\leavevmode\nobreak\ }&P_{t}(0)=1\end{array}\right.$ (18)

for any $t\leqslant|\operatorname{Sp}(H)|-1$ where
$\|P\|_{Q}^{2}\triangleq\langle P,P\rangle_{Q}=\int P^{2}\leavevmode\nobreak\
\mathrm{d\lambda_{Q}}$ denotes the underlying norm of the inner product
$\langle\cdot,\cdot\rangle_{Q}$. For $t\geqslant|\operatorname{Sp}(H)|$, we
consider instead $P_{t}$ as a multiple of the polynomial
$\prod_{\lambda\in\operatorname{Sp}(H)}(X-\lambda)$ in $X$. The next steps are
somewhat standard and follow a classical pattern for solving eq. 18 (see, e.g.
Berthier et al. (2020) and the references therein).

#### From minimal norm to orthogonality.

The solution to eq. 18 is the projection of the polynomial $0$ over the affine
space $\left\\{P\in\mathbb{R}_{t}[X]\mid P(0)=1\right\\}$ with respect to the
inner product $\langle\cdot,\cdot\rangle_{Q}$. A necessary and sufficient
condition for $P$ to be the solution of problem eq. 18 is therefore to verify
$\langle 0-P,\Delta P\rangle_{Q}=0$ for any $\Delta P$ in the vectorial
subspace $\left\\{P\in\mathbb{R}_{t}[X]\mid
P(0)=0\right\\}=X\mathbb{R}_{t-1}[X]$. Hence $P_{t}$ solves problem eq. 18 iff

$\langle P_{t},XR\rangle_{Q}=0,\forall R\in\mathbb{R}_{t-1}[X].$ (19)

Note however, that for any $(P,R)\in\mathbb{R}[X]^{2}$,

$\displaystyle\langle P,XR\rangle_{Q}$
$\displaystyle=\displaystyle\int_{\lambda\in\mathbb{R}^{+}}P(\lambda)\times\lambda
R(\lambda)\leavevmode\nobreak\ \mathrm{d\lambda_{Q}(\lambda)}$
$\displaystyle=\displaystyle\int_{\lambda\in\mathbb{R}^{+}}P(\lambda)\times
R(\lambda)\leavevmode\nobreak\ \mathrm{d\lambda_{XQ}(\lambda)}$
$\displaystyle\triangleq\langle P,R\rangle_{XQ}$

with $\mathrm{d\lambda_{XQ}(\lambda)}\triangleq\mathrm{\lambda
d\lambda_{Q}(\lambda)}=\displaystyle\sum_{\lambda\in\operatorname{Sp}(H)}\lambda
Q(\lambda)\times\langle
x_{0}-x_{\star},v_{\lambda}\rangle^{2}\leavevmode\nobreak\ \delta_{\lambda}$.
Using the latter inner product, the condition for $P_{t}$ to be the solution
to problem eq. 18 becomes:

$P_{t}\in\mathbb{R}_{t-1}[X]^{\perp_{XQ}}$ (20)

Hence, $(P_{t})_{t\in\mathbb{N}}$ is a family of orthogonal polynomials for
the inner product $\langle\cdot,\cdot\rangle_{XQ}$.

#### From orthogonality to recursion.

We now focus on finding an explicit expression for the polynomials $P_{t}$. As
for all families of orthogonal polynomials, $(P_{t})_{t\in\mathbb{N}}$ can be
obtained through a two-term recursion of the form:

$P_{t+1}(X)=(a_{t}X+b_{t})P_{t}(X)+c_{t}P_{t-1}(X),\leavevmode\nobreak\ \text{
for some }(a_{t},b_{t},c_{t})\in\mathbb{R}^{3},$ (21)

which is easy to verify by induction. Our goal is to find $a_{t}$, $b_{t}$ and
$c_{t}$. First, notice that $(a_{t}X+b_{t})P_{t}(X)+c_{t}P_{t-1}(X)$ is
orthogonal to $\mathbb{R}_{t-2}[X]$ independently of the values of $a_{t}$,
$b_{t}$ and $c_{t}$. Those three coefficients can be found via the following
three conditions: (i) $\langle P_{t+1},P_{t}\rangle_{XQ}=0$, (ii) $\langle
P_{t+1},P_{t-1}\rangle_{XQ}=0$, and (iii) $P_{t+1}(0)=1$.

More precisely, it is clear that $a_{t}\neq 0$ for $P_{t+1}$ to be of degree
$t+1$. Therefore, one can factorize by $a_{t}$. Reparametrizing eq. 21, one
can write

$P_{t+1}(X)=\frac{(\tilde{a}_{t}-X)P_{t}(X)+\tilde{b}_{t}P_{t-1}(X)}{\tilde{c}_{t}},\leavevmode\nobreak\
\text{with }(\tilde{a}_{t},\tilde{b}_{t},\tilde{c}_{t})\in\mathbb{R}^{3}.$

Moreover, evaluation at $X=0$ gives
$\frac{\tilde{a}_{t}+\tilde{b}_{t}}{\tilde{c}_{t}}=1$, thereby enforcing
$\tilde{c}_{t}=\tilde{a}_{t}+\tilde{b}_{t}$. It remains to verify the two
orthogonality conditions (independent of $\tilde{c}_{t}$):

$\begin{array}[]{lll}\tilde{a}_{t}\langle
P_{t},P_{t}\rangle_{XQ}&+\tilde{b}_{t}\langle
P_{t-1},P_{t}\rangle_{XQ}&=\langle XP_{t}(X),P_{t}(X)\rangle_{XQ},\\\
\tilde{a}_{t}\langle P_{t},P_{t-1}\rangle_{XQ}&+\tilde{b}_{t}\langle
P_{t-1},P_{t-1}\rangle_{XQ}&=\langle
XP_{t}(X),P_{t-1}(X)\rangle_{XQ}.\end{array}$

Note that this system of equations is decoupled since $\langle
P_{t-1},P_{t}\rangle_{XQ}=0$, and we finally arrive to

$P_{t+1}(X)=\frac{(\tilde{a}_{t}-X)P_{t}(X)+\tilde{b}_{t}P_{t-1}(X)}{\tilde{a}_{t}+\tilde{b}_{t}},$
(22)

with

$\left\\{\begin{array}[]{ll}\tilde{a}_{t}&=\frac{\langle
XP_{t}(X),P_{t}(X)\rangle_{XQ}}{\langle P_{t},P_{t}\rangle_{XQ}},\\\
\tilde{b}_{t}&=\frac{\langle XP_{t}(X),P_{t-1}(X)\rangle_{XQ}}{\langle
P_{t-1},P_{t-1}\rangle_{XQ}}.\end{array}\right.$ (23)

#### From a polynomial recursion to an iterative optimization method.

For reaching the final desired result, we simply multiply eq. 22 (evaluated in
$H$) by $x_{0}-x_{\star}$:

$\displaystyle x_{t+1}-x_{\star}$ $\displaystyle\leavevmode\nobreak\
=\frac{\tilde{a}_{t}(x_{t}-x_{\star})-H(x_{t}-x_{\star})+\tilde{b}_{t}(x_{t-1}-x_{\star})}{\tilde{a}_{t}+\tilde{b}_{t}},$
$\displaystyle\leavevmode\nobreak\
=x_{t}-x_{\star}-\frac{1}{\tilde{a}_{t}+\tilde{b}_{t}}\nabla
f(x_{t})+\frac{-\tilde{b}_{t}}{\tilde{a}_{t}+\tilde{b}_{t}}(x_{t}-x_{t-1}),$

thereby reaching the desired

$x_{t+1}\leavevmode\nobreak\ =x_{t}-h_{t}\nabla f(x_{t})+m_{t}(x_{t}-x_{t-1})$
(24)

with

$h_{t}=\frac{1}{\tilde{a}_{t}+\tilde{b}_{t}},\quad\text{and}\quad
m_{t}=\frac{-\tilde{b}_{t}}{\tilde{a}_{t}+\tilde{b}_{t}}.$ (25)

From eq. 24, we recognize an Heavy-ball method with some variable step-size
$h_{t}$ and momentum term $m_{t}$, thereby concluding the proof.
$\hfill\blacksquare$

######  Remark 2.2 (Step-size parametrization.)

While the $\gamma_{t}$ plays a different role in eq. 8 and eq. 7, they both
usually are called “step-size” by default. But we noticed that both in
Chebyshev method and the Heavy-ball method (optimally tuned),
$h_{t}=\frac{\gamma_{t}}{1+m_{t}}$ is exactly $\frac{2}{L+\mu}$, value of the
optimal step-size for gradient descent. In eq. 5, we notice again that the
value of $h_{t}$ is the optimal step-size for a single step of gradient
descent. For this reason, we believe that the natural parametrization of the
Heavy-ball methods should be $x_{t+1}=x_{t}-(1+m_{t})\times h_{t}\nabla
f(x_{t})+m_{t}(x_{t}-x_{t-1})$ and that $h_{t}$ should be referred to as the
“natural” step-size. Indeed, when one thinks of the Heavy-ball method with
Polyak step-sizes, they would set $\gamma_{t}$ to the Polyak step-size, not
$h_{t}=\frac{\gamma_{t}}{1+m_{t}}$. We therefore provide a novel view on what
should be tested.

##  3 Numerical experiments

In this section, we compare gradient descent, Heavy-ball and conjugate
gradient method in an adaptive setting or not. Figure 2 shows the performance
of all these methods on a quadratic objective with known minimal value
$f_{\star}$. The hessian of this quadratic objective has been generating from
a sequence of eigenvalues with geometric increase, and a random orthogonal
transformation. The difference between Figure 1 and Figure 2 is the dimension
of the problem as well as the condition number of the objective function. Due
to finite precision arithmetic, the finite-time convergence is not visible
when the condition number is too large. However both figures show that our
method and the conjugate gradient algorithm behave similarly and faster than
the other methods. The code can be found on this GitHub repository.

(a) Comparison of distances to optimum

(b) Comparison of excess losses

Figure 2: Comparison in semi-log scale over 2000 iterations of different
first-order methods applied on a 1000-dimensional quadratic objective with
condition number $10^{5}$. GD with constant step-size, GD with Polyak step-
size and GD with variant of Polyak step-size refer to the GD method tuned
respectively with the step-size $\gamma=2/(L+\mu)$,
$\gamma_{t}=(f(x_{t})-f_{\star})/\|\nabla f(x_{t})\|^{2}$ and
$\gamma_{t}=2(f(x_{t})-f_{\star})/\|\nabla f(x_{t})\|^{2}$. HB with constant
tuning is the HB method tuned with constant parameters
$\gamma_{t}=(2/(\sqrt{L}+\sqrt{\mu}))^{2}$ and
$m_{t}=((\sqrt{L}-\sqrt{\mu})(\sqrt{L}+\sqrt{\mu}))^{2}$ while HB with Polyak
step-size based tuning refers to Algorithm 1.

##  4 Concluding remarks and discussion

Polyak step-sizes are known for their general good working performances when
the optimal value to the optimization problem at hand is known. The question
of whether Polyak step-sizes can be used together with momentum for obtaining
accelerated first-order methods appears to be an open question (Barré et al.,
2020), which we answer in the simpler case of convex quadratic minimization.
In this context, we argue that not only this tuning works well, but also it
pops up naturally when investigating instance-optimal first-order iterative
methods. Furthermore, we believe it is a necessary step for being able to
understand more general optimization settings beyond quadratics. As our method
does not seem to work well beyond quadratics, we leave further investigations
on this topic for future work.

Among our competitors, we note that the celebrated conjugate gradient (CG)
method is another instance-optimal algorithm for quadratics. Whereas our
method minimizes the distance to the solution at each iteration, CG is
instance-optimal for minimizing function values at each iteration. Perhaps
interestingly, the two methods appeared to behave similarly in our numerical
experiments. That being said, the main practical differences between the two
methods are that CG Heavy-ball-like formulation naturally relies on higher
order information while Polyak step-sizes do require knowledge of $f_{\star}$.
In typical optimization problems, this value is not known. However, there are
a few settings where this value is actually well-known, typically when
$f_{\star}=0$ generically (in machine learning, this setting is known as the
“interpolation” regime; an alternative could be to use Polyak-steps as a
competitor to MinRes). Finally, let us mention that although a few
generalization of CG, often referred to as nonlinear conjugate gradient, were
studied in the literature (see, e.g., (Bonnans et al., 2006; Nocedal and
Wright, 1999; Hager and Zhang, 2006)).

## Acknowledgement

We would to thank Raphael Berthier for his insightful feedbacks. The work of
B. Goujaud and A. Dieuleveut is partially supported by
ANR-19-CHIA-0002-01/chaire SCAI, and Hi!Paris. A. Taylor acknowledges support
from the European Research Council (grant SEQUOIA 724063). This work was
partly funded by the French government under management of Agence Nationale de
la Recherche as part of the “Investissements d’avenir” program, reference
ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

## References

  * Barré et al. [2020] Mathieu Barré, Adrien Taylor, and Alexandre d’Aspremont.  Complexity guarantees for Polyak steps with momentum.  In _Conference on Learning Theory (COLT)_ , pages 452–478. PMLR, 2020\.
  * Barzilai and Borwein [1988] Jonathan Barzilai and Jonathan M. Borwein.  Two-point step size gradient methods.  _IMA journal of numerical analysis_ , 8(1):141–148, 1988.
  * Berthier et al. [2020] Raphaël Berthier, Francis Bach, and Pierre Gaillard.  Accelerated gossip in networks of given dimension using jacobi polynomial iterations.  _SIAM Journal on Mathematics of Data Science_ , 2(1):24–47, 2020.
  * Bonnans et al. [2006] Joseph-Frédéric Bonnans, Jean-Charles Gilbert, Claude Lemaréchal, and Claudia A. Sagastizábal.  _Numerical optimization: theoretical and practical aspects_.  Springer Science & Business Media, 2006.
  * Boyd et al. [2003] Stephen Boyd, Lin Xiao, and Almir Mutapcic.  Subgradient methods.  _lecture notes of EE392o, Stanford University, Autumn Quarter_ , 2003\.
  * D’Orazio et al. [2021] Ryan D’Orazio, Nicolas Loizou, Issam Laradji, and Ioannis Mitliagkas.  Stochastic mirror descent: Convergence analysis and adaptive variants via the mirror stochastic Polyak stepsize.  _arXiv preprint arXiv:2110.15412_ , 2021.
  * d’Aspremont et al. [2021] Alexandre d’Aspremont, Damien Scieur, and Adrien Taylor.  Acceleration methods.  _Foundations and Trends® in Optimization_ , 5(1-2):1–245, 2021.
  * Fischer [2011] Bernd Fischer.  _Polynomial based iteration methods for symmetric linear systems_.  SIAM, 2011.
  * Golub and Varga [1961] Gene H. Golub and Richard S. Varga.  Chebyshev semi-iterative methods, successive overrelaxation iterative methods, and second order richardson iterative methods.  _Numerische Mathematik_ , 3(1):157–168, 1961\.
  * Goujaud et al. [2022] Baptiste Goujaud, Damien Scieur, Aymeric Dieuleveut, Adrien B. Taylor, and Fabian Pedregosa.  Super-acceleration with cyclical step-sizes.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 3028–3065. PMLR, 2022.
  * Gower et al. [2022] Robert M. Gower, Mathieu Blondel, Nidham Gazagnadou, and Fabian Pedregosa.  Cutting some slack for SGD with Adaptive Polyak Stepsizes.  _arXiv preprint arXiv:2202.12328_ , 2022.
  * Hager and Zhang [2006] William W. Hager and Hongchao Zhang.  A survey of nonlinear conjugate gradient methods.  _Pacific journal of Optimization_ , 2(1):35–58, 2006.
  * Hazan and Kakade [2019] Elad Hazan and Sham Kakade.  Revisiting the Polyak step size.  _arXiv preprint arXiv:1905.00313_ , 2019.
  * Loizou et al. [2021] Nicolas Loizou, Sharan Vaswani, Issam Hadj Laradji, and Simon Lacoste-Julien.  Stochastic Polyak step-size for SGD: An adaptive learning rate for fast convergence.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 1306–1314. PMLR, 2021.
  * Malitsky and Mishchenko [2020] Yura Malitsky and Konstantin Mishchenko.  Adaptive gradient descent without descent.  In _International Conference on Machine Learning (ICML)_ , pages 6702–6712. PMLR, 2020.
  * Nemirovski [1992] Arkadi S. Nemirovski.  Information-based complexity of linear operator equations.  _Journal of Complexity_ , 8(2):153–175, 1992\.
  * Nemirovski [1994] Arkadi S. Nemirovski.  Information-based complexity of convex programming.  _Lecture notes,http://www2.isye.gatech.edu/~nemirovs/Lec_EMCO.pdf_, 1994.
  * Nesterov [1983] Yurii Nesterov.  A method of solving a convex programming problem with convergence rate ${O}(1/k^{2})$.  _Soviet Mathematics Doklady_ , 27(2):372–376, 1983.
  * Nesterov [2004] Yurii Nesterov.  _Introductory lectures on convex optimization: a basic course_.  Applied optimization. Kluwer Academic Publishers, 2004.
  * Nocedal and Wright [1999] Jorge Nocedal and Stephen J. Wright.  _Numerical optimization_.  Springer, 1999.
  * Pedregosa [2020] Fabian Pedregosa.  Residual polynomials and the chebyshev method, 2020.  URL http://fa.bianp.net/blog/2020/polyopt/.
  * Polyak [1964] Boris T. Polyak.  Some methods of speeding up the convergence of iteration methods.  _USSR computational mathematics and mathematical physics_ , 4(5):1–17, 1964.
  * Polyak [1987] Boris T. Polyak.  Introduction to optimization. 1987.  _Optimization Software, Inc, New York_ , 1987.