[a,b]Diego Blas

# Probing ultralight and degenerate dark matter with galactic dynamics

###### Abstract

This short contributions summarizes a couple of recent results to test dark
matter properties with galactic dynamics. First, I will present the impact in
rotation curves from solitonic structures expected at the center of galaxies
for ultralight bosonic dark matter. As a result, one can claim that masses of
the order $m_{\rm DM}\lesssim 10^{-21}$eV are in tension with data. Second, I
will discuss how the dark matter medium properties change the way a ‘probe’
interacts with the halo. I will focus on dynamical friction and show how it is
modified in the case of degenerate fermions. This result may be used to
address the Fornax timing problem. I hope that this contribution represents an
inspiration to continue exploring other ideas in this direction of using
galactic dynamics to tell apart different dark matter models.

##  1 Introduction

The existence of a dark matter (DM) component plays a key role in the
properties of visible matter in galaxies. One of the most significant is the
modification of the gravitational potential felt by the visible matter
component, which eventually modifies properties such as the rotation curves of
visible matter at galactic scales [1]. However, when one considers dark matter
as a _medium_ in which the visible matter is moving, several new dynamical
effects may be present, e.g. tidal disruption, dynamical friction, dynamical
heating [2, 3]. These effects have been clarified and a significant amount of
dark matter phenomenology is known for the models where dark matter at
galactic scales can be understood as a collection of particles of small mass
and only interacting gravitationally. However, these _medium_ properties may
change if dark matter is made of compact objects, ultralight DM (ULDM) or if
other forces act in the dark sector. As a result, these other proper
properties of DM can be used to differentiate among dark matter models, and
may be even enhanced in some cases. A long-term goal suggest to study them for
the different dark matter models, and use _galactic dynamics for fundamental
physics_ , inspired by the related efforts related to stars [4].

In this short contribution, I will focus on two particular directions I have
worked on. First, I will discuss how _galactic rotation curves_ can be used to
test ultra-light bosonic111Due to lack of time, I won’t discuss in detail
similar bounds on the fermionic DM case, see e.g. [6]. DM [5]. Second, I will
show how _dynamical friction_ can also allow us to study the properties of
ultra-light fermionic DM. Before moving to the main body of this work, I want
to spend a few words on the landscape of DM models. A remarkable aspect of DM
is that, despite the large amount of data pointing toward its existence [1],
which furthermore comes from phenomena at very different time and length
scales, not even the most fundamental property of DM (its mass) is known. If
considered a particle, the constraints are roughly $m_{\rm DM}\gtrsim
10^{-22}$ eV for bosonic fields [7] and $m_{\rm DM}\gtrsim 100$ eV [8]. Both
of which can be derived from the extra forces that appear in dwarf spheroidals
for the models at the limit222For the bosonic case, the occupation number is
so large that a wave description is efficient, and backreacts when the
distribution is compressed to distances below the de Broglie wavelength. More
on this later. For the fermionic case, it is the degenerate pressure that
backreacts if DM particles are compressed.. Regarding the upper bound, DM
particles can be as heavy as allowed by the our effective description of
Nature. It is customary to also mention the possibility of dark matter being
made of macroscopic bodies of astrophysically significant mass [1].

##  2 Galactic rotation curves to test ultra-light bosonic DM

When shall we describe the light from a source as a classical field or as a
collection of photons? The rough answer to this question is that the classical
field emerges when the light configuration is made of several phase-space
states occupied with large occupation numbers and uncorrelated phases [9]. For
dark matter in a galaxy (e.g. the Milky way) one can ask the same question,
assuming it is bosonic. An order of magnitude estimate can be given as
follows: the Milky Way halo is supposed to have a size of $\sim$100 kpc, and a
maximum speed of $v\sim 2\times 10^{-3}$ for virialized matter (i.e. matter
that remains bounded to the galaxy). These numbers give a total number of
states which can be filled in the relevant phase of the Milky Way of order
$N_{s}\approx 10^{75}\left(m_{\rm DM}/\rm eV\right)^{3}$. The predicted mass
of its DM halo being $M_{\rm DM}\approx 10^{12}\,M_{\odot}$ [10], and assuming
that DM is distributed universally over all the allowed states, on finds that
this mass requires the occupation number

$N_{p}=M_{\rm DM}/(N_{s}m_{\rm DM})\approx 10^{3}\left({\rm eV}/m_{\rm
DM}\right)^{4}.$ (1)

As a result, for $m_{\rm DM}\lesssim$ eV, DM is better described using
classical equations. In the simplest example of a scalar field $\phi$, the
latter will correspond to the Klein-Gordon equation coupled to gravity. To
understand the final configuration quantitatively, one can consider the
typical wavelength of these waves $\lambda$. For scales larger than $\lambda$,
the wavelike nature won’t play any role, and these DM configurations can be
described as a collection of free waves (or wavepackets) virializing through
gravity and with typical velocity $\bar{v}$, corresponding to a typical
wavelength $\bar{\lambda}\sim 1/(m\bar{v})$. As a result, for distances
$r\gg\bar{\lambda}$ (the ‘halo’ part of the galaxy), the expected field is (up
to normalization)

$\phi_{h}\sim\int_{0}^{v_{\rm
max}}\mathrm{d}^{3}ve^{-v^{2}/\sigma_{0}^{2}}e^{i(\omega
t-m\vec{v}\cdot\vec{x})+\psi_{v}}+\mathrm{h.c.},\vspace{-.2cm}$ (2)

where I included random phases $\psi_{v}$, a Maxwellian form for the
distribution with $\sigma_{0}$ ($10^{-3}$ for the Milky Way) and
$\omega=m\sqrt{1+v^{2}}$. From the previous formula, it is clear that most
modes (labelled by $v$) in this non-relativistic distribution, oscillate with
very similar frequencies, which can be considered coherent for times
$t\lesssim 1/(m\sigma_{0}^{2})$. Even if this is not part of this talk, I want
to stress that these quasi-coherent oscillators generate _unique_ phenomena in
the galactic halo, such as an extra source of gravitational heating [11, 12,
13]. Closer to the galactic center, the backreaction from trying to compress
waves below their de Broglie wavelength $\lambda_{dB}=/(mv)$ (sometimes called
‘quantum pressure’) becomes efficient. This effect can be summarized by an
extra term appearing in the Euler equation for the ‘fluid’ degrees of freedom
of the field configuration $\phi$, $\rho$ and $\vec{v}$ (see [7]),

$\dot{\vec{v}}+\left(\vec{v}\cdot\frac{\nabla}{a}\right)\vec{v}=-\frac{\nabla}{a}\left(V-\frac{1}{2m^{2}a^{2}}\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}\right),$
(3)

where $V$ is the gravitational potential, and the last term is the one
reacting when the field is contracted below $\lambda_{dB}\approx(10^{-22}{\rm
eV}/m_{\rm DM})(10^{-3}/v)\,$kpc. Below these scales, numerical
simulations333These simulations were also very relevant to show that the rest
of the DM halo has a standard Navarro-Frenk-White (NFW) profile. have shown
the existence of a coherent structure at the center of galaxies, that have
received the named of ‘soliton’ [14, 15]. Quite remarkably, one can understand
the main features of these structures by considering a field configuration of
the form $\phi\approx\frac{1}{\sqrt{2}m}e^{-im_{\rm DM}t}e^{-i\gamma
t}\chi(r)+\mathrm{h.c.},$, where we have assumed a spherically symmetric
configuration and $\gamma\ll m$ is a parameter tat controls the non-
relativistic limit. The final system of equations satisfied by $\chi(r)$
correspond to a self-gravitating configuration, and are not fixed by simply
assuming that the field vanishes at large distances, see e.g. [5]. In fact,
the scaling $\chi_{\lambda}(r)=\lambda^{2}\chi_{1}(\lambda r)$ and
$\gamma_{\lambda}=\lambda^{2}\gamma$ generates a new solution. By fitting the
configurations of the solitonic part of the halos in [14, 15], one can see
that they are fitted by the $\chi(r)$ functions found in the previous limit
for a particular value of $\gamma$ for a given DM halo. The numerical
simulations yield a very intriguing/interesting relation that allows one to
predict the size of the soliton for any DM halo up to some uncertainties. The
first of these relations was noticed in [16], and in [5] we expressed it as
(this is valid up to factors of a few)

$\left.\frac{E}{M}\right|_{\text{soliton}}\approx\left.\frac{E}{M}\right|_{\text{halo
}},$ (4)

where $E/M$ is the total energy (including gravitational energy) divided by
the mass of either the soliton or the halo. This relation suggests some kind
of equilibrium between the degrees of freedom of the halo and those of the
soliton. The simulations performed in [15] found that a whole variety of
relations is possible, and it is is still not clear if this corresponds to the
fact that the configuration is not relaxed. In any case, what seems to be
universal is the existence of the soliton in the center of galactic DM halos
of considerable size and that changes the profile of DM at scales smaller than
$\lambda_{dB}$, generating a core where otherwise one would expect a a cuspy
NFW profile [1].

If the DM profile is modified, its corresponding gravitational potential will
also change. As a consequence, the properties of the visible matter connected
to the average value of the DM potential (e.g. the rotation curves) may be
modified. To model this, one can take a simplistic NFW profile at the
outskirts of a galaxy, characterized by a concentration and a radius
parameter, and use Eq. (4) to find the corresponding soliton of this
particular DM halo. The final gravitational potential will be a sum of
$\Phi_{\rm halo}$ and $\Phi_{\rm soliton}$ and the expected orbital velocities
will be modified as $V_{\text{circ }}^{2}=r\partial_{r}(\Phi_{\rm
halo}+\Phi_{\rm soliton})$. A particularly interesting aspect is that the
soliton profile generates a second peak in the velocity of rotation at small
distances of the same size as the one generated by the halo part at large
distances from the galactic center (see left panel of Fig. 1). With this
observation, one can now try to find the effect of the potential from the
soliton into the rotation curves of galaxies. For this, in [5], we used the
SPARC sample [17], which provides high resolution rotation curves of late-type
disc galaxy data. Furthermore, the baryonic effects can also be estimated from
the provided photometric data. As a result of not finding the inner peak
associated to the soliton (see left panel of Fig. 1), or the corresponding
changes in rotation curves produced by it (when the peak distance to the
center is too small to be resolved), we concluded that it ULDM of masses below

$m_{\rm DM}\lesssim 10^{-21}\,\rm eV\vspace{-.3cm}$ (5)

are in tension with data. Since the size of the soliton scales as $1/m_{\rm
DM}$, smaller masses are hard to probe, though, as also shown in [5], our own
galaxy may probe relevant in this case, since the gravitational potential of
the Milky Way has been probed until $10^{-4}$pc. A couple of points are
relevant: first, using a sample of hundreds of galaxies is very relevant to
make sure that we are not looking at outlayers when searching for the soliton.
Our method in [5] used 175 galaxies for this purpose. We also considered the
possible influence of stars and supermassive black holes in the centers of the
galaxies, and concluded that our bound in (5) is robust against it (see also
[18, 19, 20]).

##  3 Dynamical friction to test light fermionic DM

As I discussed above, the study of the velocity dispersion of dwarf
spheroidals and their size, can be translated into an argument similar to the
one behind (1) to claim that, because the occupation number of fermionic
states cannot surpass 1, the bound $m\gtrsim 100$ eV (known as Tremaine-Gunn
bound) emerges [8, 6]. A very important aspect that is not always emphasized
is that, as the mass is reduced, and the phase space fills, it is not clear
what the distribution should be. This is because it should correspond to a
collection of (almost) degenerate fermions, out of equilibrium and interacting
through gravity. These two points make numerical simulations quite
challenging, as one should resort to quantum Vlasov equations, which I have
not seen solved in the context of dark matter444A similar problem has been
addressed in condensed matter. As I mentioned in my talk, this is a problem
that interests me, if the reader wants to reach out to discuss it further..
However, the relevance of this calculation seems limited, since the effects of
degenerate DM as far as the distribution is concerned will be likely described
by an effective temperature, related to the virial velocity, whose effect can
be easily studied.

The question we posed ourselved in [21] was what happens when you move
_across_ this degenerate medium of fermions, and interact with it
gravitationally. Would aspects of dynamical friction, tidal striping, or other
dynamical processes where the DM halo interacts gravitationally (see [2, 3])
be modified to observable levels? A similar question was asked in [22] for
ULDM, and several interesting directions were identified. In this contribution
I only discuss dynamical friction (DF). The reason why we expect it to be
modified in ULDM and light fermionic DM is the following. DF arises from the
wake (a field overdensity) left behind as a perturber (e.g. a collection of
stars, as a globular cluster) travels through a medium (e.g. the DM) with
which it interacts gravitationally. For ULDM, if this overdensity is smaller
than $\lambda_{dB}$, the final term in (3) responds to the contraction, and
reduces the size of the wake, modifying DF. It turns out that the final answer
is more complex, and still being investigated, se e.g [23, 13, 24, 25, 26].
For the degenerate fermions, one expects that degenerate pressure counteracts
when one tries to compress it below the scale where fermions should occupy the
same state (recall that for a galaxy, we assume a virial distribution, which
as a related Fermi momentum). To realize this intuition, in [21] (see also
[27]) we followed a Fokker-Planck approach where a single perturber with
momentum $p$ and distribution $f_{1}$ interacts with a medium with
distribution $f_{2}$ through gravity. This is encapsulated in a collision term

$\displaystyle C\left[f_{1}\right]=\frac{(2\pi)^{4}}{2E_{p}}\int
d\Pi_{k,p^{\prime},k^{\prime}}\delta^{(4)}\left(p+k-p^{\prime}-k^{\prime}\right)|\overline{\mathcal{M}}|^{2}\left[f_{1}\left(p^{\prime}\right)f_{2}\left(k^{\prime}\right)\right.$
$\displaystyle\left.~{}~{}~{}\left(1\pm f_{1}(p)\right)\left(1\pm
f_{2}(k)\right)-f_{1}(p)f_{2}(k)\left(1\pm
f_{1}\left(p^{\prime}\right)\right)\left(1\pm
f_{2}\left(k^{\prime}\right)\right)\right].$ (6)

In this formula, $E_{p}$ is the energy of the inital state,
$|\overline{\mathcal{M}}|$ represents the gravitational scattering cross-
section from states with momentum $p$ (perturber) and $k$ (medium) to
$p^{\prime}$ and $k^{\prime}$, and the fermionic/bosonic nature of the medium
is represented by the relative signs in front of $f_{2}$. The final ‘friction’
has the following effect in the velocity $\mathbf{V}$ of the perturber of mass
$M$:

$\displaystyle\frac{d\mathbf{V}}{dt}=-\frac{4\pi G^{2}\rho
M\ln\Lambda}{V^{3}}F(V/\sigma)\mathbf{V},$ (7)

where $\rho$ is the density of the medium, $M$ its mass, $G$ is the
gravitational Newton constant, $\ln\Lambda$ represents a cut-off scale
generating an $O(1)$ contribution and $F$ captures the details of the medium.
For instance, for a Maxwellian distribution of particles with velocity
dispersion $\sigma$ one gets $F\sim 1$ for $V\gg\sigma$ and
$F\sim\frac{\sqrt{2}}{3\sqrt{\pi}}\frac{V^{3}}{\sigma^{3}}$ for $V\ll\sigma$.
For a degenerate dark matter (DDM) gas characterized by a Fermi velocity
$v_{F}$, one finds that

$F_{\mathrm{DDM}}\rightarrow 1,\ {\rm for}\ V\gg
v_{F};~{}~{}~{}F_{\mathrm{DDM}}\rightarrow\frac{V^{3}}{v_{F}^{3}},\ {\rm for}\
V\ll v_{F}.\vspace{-.2cm}$ (8)

This effect generates a reduction of dynamical friction for slow enough
perturbers. Furthermore, the degenerate gas of fermions will generate a core
structure that may reduce $\sigma$ and increase $\rho$. This core extends to
distances $r_{c}\sim
681\left(\frac{\rho_{0}}{10^{7}\mathrm{M}_{\odot}/\mathrm{kpc}^{3}}\right)^{-\frac{1}{6}}\left(\frac{gm^{4}}{2\times(120\mathrm{eV})^{4}}\right)^{-\frac{1}{3}}\mathrm{pc}$
(see also footnote 2).

This last phenomena (modification of $\rho$ and/or $\sigma$) may also happen
in ULDM, models with self-interacting dark matter (SIDM), or where baryonic
feedback is modified. In [21] we considered different models (SIDM, degenerate
dark matter with $m_{\mathrm{DM}}\approx 135\,\mathrm{eV}$, different baryonic
feedback processes) that can fit the data connected to kinematics, and study
if they can be differentiated with dynamical friction. For this, we focused on
data from Fornax. This very luminous nearby ($\sim 147$ kpc away) dwarf
satellite has $\sim 4\times 10^{7}M_{\odot}$ stars in a $\sim$ kpc size. It is
a dark matter dominated galaxy, with 5-6 globular clusters (GCs) with masses
$\sim 10^{5}M_{\odot}$, two of which at locations that seem very improbable,
if dynamical friction has the size predicted by (7) with $F\sim 1$. To
understand this, one estimate the time scale for a perturber to plunging to
the center of a galaxy (we assume $F\sim 1$)

$\tau\equiv\frac{|\mathbf{V}|}{|d\mathbf{V}/dt|}\sim
1.8\left(\frac{V}{12\mathrm{~{}km}/\mathrm{s}}\right)^{3}\frac{2\times
10^{7}\frac{M_{\odot}}{\mathrm{kpc}^{3}}}{\rho}\mathrm{Gyr},$ (9)

where we have included the values typical of Fornax. This time scale is short
enough to be suspicious of GCs at locations corresponding to time scales
shorter than $\sim 10$ Gpc.

By studying the data on the line-of-sight velocity dispersion of $\sim$ 2500
stars from Fornax, we are able to fit them with the aforementioned DM models
(see [21, 27]). More interestingly, these very same models produce
_significantly_ different dynamical friction, since the medium properties are
different. Still, all SIDM, more compact DM profiles and degenerate dark
matter allows to increase $\tau$ to levels of $O(10)$ Gpc (the final answer
being model dependent), and also agree with the kinematic properties of Fornax
(see right panel of Fig. 1). As a result, one could invoke them to solve this
‘GC’ Fornax timing problem555In [22] a similar reasoning was used to suggest a
mass of DM of $m_{\rm DM}\sim 10^{-22}$eV., though it is hard to distinguish
between them.

Before considering this possibility, one has to face a question about initial
conditions of the GCs. Indeed, the probability to end up with GCs at the
positions we measure depends on the initial distribution. Unfortunately, the
data from a single galaxy (Fornax) is not enough to make a statement, since
the initial conditions are only known statistically. Even assuming that Fornax
is a good representive of the initial distribution of GCs, we showed that the
standard DM effect does still produce GCs in the positions observed in Fornax
with a 25% probability. Hence, it is not clear to which extend Fornax poses a
problem to standard DM models, though the idea of using DF to differentiate
among DM models is still sound. In particular, the most natural way forward is
to more data from systems where GCs or similar trackers are resolved in DM
halos.

##  4 Summary and outlook

This contribution is devoted to new ideas to use galactic dynamics to
differentiate between DM models. The final goal would be to implement the
dynamics of different DM models into all the phenomenology related to the DM
halo. This should include rotation curves, but also many other dynamical
phenomena, already known for more standard DM models [1, 2, 3].

Figure 1: Left panel: rotation curves for NGC1560 vs the peak predicted by a
soliton of $m_{\rm DM}\sim 10^{-22}$ eV, from [5]. Right panel: typical time
scale associated to dynamical friction for different models of DM, and GCs
located at different radii $x$ of Fornax. One would like to have all these
times above few Gyrs.

I have focused on two particular examples. In Sec. 2, I discussed some aspects
in which the distribution of ULDM differs from that of standard DM. In
particular, the existence of ‘quantum pressure’ implies that the field reacts
when trying to squeeze it below $\lambda_{dB}$. As a result, a coherent
structure (soliton) appears in the center of DM halos. This soliton modifies
the rotation curves to levels that seem to be in tension with data for $m_{\rm
DM}\lesssim 10^{-21}$eV. Higher masses requires measuring the properties of
the galactic halo at distances below 0.1 kpc, and one could think about using
the dynamical data from the Milky Way to explore them. This seems a good point
to mention that other properties of ULDM may modify this conclusion, e.g.
self-interaction or the presence of supermassive black holes in the cases of
higher masses. Also, it is important to recall that [15] found a condition
between the solitons and the halo differing from (4). Relevant for ULDM, the
halo itself has patches with wave behaviour, which modifies gravitational
heating or dynamical friction.

In the second part of this contribution I focused on degenerate DM fermions.
For these, an optimistic bound of $m_{\rm DM}\gtrsim 100$ eV arises from
galactic dynamics. In the limiting case, small galaxies develop a core from
the degenerate pressure. Furthermore, since the DM distribution is degenerate,
one expects that it is harder to exchange momentum with it (since the ‘out-
going’ momentum state may be occupied). As a result, the effective coupling of
a perturber with the DM halo should decrease, and enhance dynamical friction
time scales. We showed in [21] that this intuition is correct, and that the
time scales for perturbers that we see today in the DM halo to fall to the
galactic center can be increased to levels of order $\sim$10 Gpc. However, the
same conclusion can be achieved with other DM models. Furthermore, the mass
considered $m_{\mathrm{DM}}\approx 135\mathrm{eV}$, is almost excluded by
galactic dynamics and hard to reconcile this small mass with other
cosmological bounds [21, 28]. As a possible future direction, I find it
interesting to understand which distribution is generated by a collection of
degenerate fermions in ‘virial equilibrium’ generated by gravity.

Acknowledgements. I thank the organizers of the 1st General Meeting and 1st
Training School of the COST Action COSMIC WSIPers for their kind invitation.
It is also a pleasure to thank all my great collaborators on the topics I
discussed, N. Bar, K. Blum, H. Kim and S. Sibiryakov. I am also grateful to A.
Caputo for comments on a previous draft. The research leading to these results
has received funding from the Spanish Ministry of Science and Innovation
(PID2020–115845GB-I00/AEI/10.13039/501100011033). IFAE is partially funded by
the CERCA program of the Generalitat de Catalunya. I also acknowledge the
support from the Departament de Recerca i Universitats de la Generalitat de
Catalunya al Grup de Recerca i Universitats from Generalitat de Catalunya to
the Grup de Recerca 00649 (Codi: 2021 SGR 00649).

## References

  * [1] G. Bertone, J. Silk, B. Moore, J. Diemand, J. Bullock, M. Kaplinghat, L. Strigari, Y. Mellier, D. Merritt and J. Bekenstein, et al. “Particle Dark Matter: Observations, Models and Searches,” Cambridge Univ. Press, 2010.
  * [2] J. Binney and S. Tremaine, “Galactic Dynamics: Second Edition,” Princeton Series in Astrophysics, 2008.
  * Mo et al. [2010] Mo, H., van den Bosch, F. C., & White, S. 2010, “Galaxy Formation and Evolution”, Cambridge, UK: Cambridge University Press, 2010.
  * [4] G. G. Raffelt, Ann. Rev. Nucl. Part. Sci. 49 (1999), 163-216 [arXiv:hep-ph/9903472 [hep-ph]].
  * [5] N. Bar, D. Blas, K. Blum and S. Sibiryakov, Phys. Rev. D 98 (2018) no.8, 083027 [arXiv:1805.00122 [astro-ph.CO]].
  * [6] J. Alvey, N. Sabti, V. Tiki, D. Blas, K. Bondarenko, A. Boyarsky, M. Escudero, M. Fairbairn, M. Orkney and J. I. Read, Mon. Not. Roy. Astron. Soc. 501 (2021) no.1, 1188-1201 [arXiv:2010.03572 [hep-ph]].
  * [7] E. G. M. Ferreira, Astron. Astrophys. Rev. 29 (2021) no.1, 7 [arXiv:2005.03254 [astro-ph.CO]].
  * [8] S. Tremaine and J. E. Gunn, Phys. Rev. Lett. 42 (1979), 407-410
  * [9] A. Duncan, “The Conceptual Framework of Quantum Field Theory,” Oxford University Press, 2012.
  * [10] L. Posti, A. Helmi A&A 621, A56 (2019) [arXiv:1805.01408 [astro-ph.GA]].
  * [11] D. J. E. Marsh and J. C. Niemeyer, Phys. Rev. Lett. 123 (2019) no.5, 051103 [arXiv:1810.08543 [astro-ph.CO]].
  * [12] N. Dalal and A. Kravtsov, Phys. Rev. D 106 (2022) no.6, 063517 [arXiv:2203.05750 [astro-ph.CO]].
  * [13] B. Bar-Or, J. B. Fouvry and S. Tremaine, Astrophys. J. 871 (2019) no.1, 28 [arXiv:1809.07673 [astro-ph.GA]].
  * [14] H. Y. Schive, T. Chiueh and T. Broadhurst, Nature Phys. 10 (2014), 496-499 [arXiv:1406.6586 [astro-ph.GA]].
  * [15] H. Y. J. Chan, E. G. M. Ferreira, S. May, K. Hayashi and M. Chiba, Mon. Not. Roy. Astron. Soc. 511 (2022) no.1, 943-952 [arXiv:2110.11882 [astro-ph.CO]].
  * [16] H. Y. Schive, M. H. Liao, T. P. Woo, S. K. Wong, T. Chiueh, T. Broadhurst and W. Y. P. Hwang, Phys. Rev. Lett. 113 (2014) no.26, 261302 [arXiv:1407.7762 [astro-ph.GA]].
  * [17] F. Lelli, S. S. McGaugh and J. M. Schombert, Astron. J. 152 (2016), 157 [arXiv:1606.09251 [astro-ph.GA]].
  * Chan et al. [2018] Chan, J. H. H., Schive, H.-Y., Woo, T.-P., et al. 2018, MNRAS, 478, 2686 [1712.01947[astro-ph.GA]].
  * [19] N. Bar, K. Blum, J. Eby and R. Sato, Phys. Rev. D 99 (2019) no.10, 103020 [arXiv:1903.03402 [astro-ph.CO]].
  * [20] N. Bar, K. Blum, T. Lacroix and P. Panci, JCAP 07 (2019), 045 [arXiv:1905.11745 [astro-ph.CO]].
  * [21] N. Bar, D. Blas, K. Blum and H. Kim, Phys. Rev. D 104 (2021) no.4, 043021 [arXiv:2102.11522 [astro-ph.GA]].
  * [22] L. Hui, J. P. Ostriker, S. Tremaine and E. Witten, Phys. Rev. D 95 (2017) no.4, 043541 [arXiv:1610.08297 [astro-ph.CO]].
  * [23] L. Lancaster, C. Giovanetti, P. Mocz, Y. Kahn, M. Lisanti and D. N. Spergel, JCAP 01 (2020), 001 [arXiv:1909.06381 [astro-ph.CO]].
  * [24] R. Buehler and V. Desjacques, Phys. Rev. D 107 (2023) no.2, 023516 [arXiv:2207.13740 [astro-ph.CO]].
  * [25] Y. Wang and R. Easther, Phys. Rev. D 105 (2022) no.6, 063523 [arXiv:2110.03428 [gr-qc]].
  * [26] R. Vicente and V. Cardoso, Phys. Rev. D 105 (2022) no.8, 083008 [arXiv:2201.08854 [gr-qc]].
  * [27] D. Blas, [arXiv:2205.00289 [hep-ph]].
  * [28] M. Carena, N. M. Coyle, Y. Y. Li, S. D. McDermott and Y. Tsai, Phys. Rev. D 106 (2022) no.8, 083016 doi:10.1103/PhysRevD.106.083016 [arXiv:2108.02785 [hep-ph]].


$\frac{C_{b}}{C_{c}}=\frac{1/C_{f}^{\prime}-1}{K_{R}}$, we have
$C_{e}=C_{e}^{\prime}$ and $C_{f}=C_{f}^{\prime}$ and
$1/K_{R}=C_{f}^{2}+\frac{C_{f}^{4}}{C_{e}(2-C_{e})}$. Then it follows from
$\overset{b}{\geq}$ in (58) that

$\displaystyle\frac{1}{K_{A}n_{A}}\sum_{k=1}^{K_{A}}\sum_{i=1}^{n_{A}}\mathcal{L}(\mathbf{W}\bm{h}_{k,i},\mathbf{y}_{k})\geq-
C_{g}(1+1/K_{R})\sqrt{E_{H}E_{W}}+C_{d}=L_{2}.$ (63)

On the other hand, consider the solution $(\mathbf{H},\mathbf{W})$ that
satisfies

$\displaystyle\begin{aligned}
&\mathbf{w}_{k}=g_{N}(x_{0})\mathbf{P}_{A}\left[\frac{x_{0}}{\sqrt{(K_{A}-1)K_{A}}}(K_{A}\mathbf{y}_{k}-\mathbf{1}_{K_{A}})+\frac{1-x_{0}}{\sqrt{K_{A}}}\mathbf{1}_{K_{A}}\right],\quad
k\in[K_{A}],\\\
&\mathbf{w}_{k}=-\frac{C_{e}(2-C_{e})}{C_{f}^{2}K_{A}}\mathbf{P}_{A}\sum_{k=1}^{K_{A}}\mathbf{w}_{k},\quad
k\in[K_{A}+1:K],\\\
&\bm{h}_{k,i}=\frac{\sqrt{(1+1/K_{R})E_{H}}}{\|\mathbf{w}_{i}+\frac{C_{e}}{C_{f}K_{A}}\sum_{k=1}^{K_{A}}\mathbf{w}_{k}\|}\mathbf{P}_{A}\left[\mathbf{w}_{i}+\frac{C_{e}}{C_{f}K_{A}}\sum_{k=1}^{K_{A}}\mathbf{w}_{k}\right],\quad
k\in[K_{A}],~{}i\in[n_{A}],\\\ &\bm{h}_{k,i}=\mathbf{0}_{p},\quad
k\in[K_{A}+1:K],~{}i\in[n_{B}],\end{aligned}$

where $\mathbf{y}_{k}\in\mathbb{R}^{K}$ is the vector containing one in the
$k$-th entry and zero elsewhere and $\mathbf{P}_{A}\in\mathbb{R}^{p\times
K_{A}}$ is a partial orthogonal matrix such that
$\mathbf{P}^{\top}_{A}\mathbf{P}_{A}=\mathbf{I}_{K_{A}}$. We have
$\exp\left(\bm{h}_{k,i}^{\top}\mathbf{w}_{k}\right)=g_{a}(x_{0})$ for
$i\in[n_{A}]$ and $k\in[K_{A}]$,
$\exp\left(\bm{h}_{k,i}^{\top}\mathbf{w}_{k^{\prime}}\right)=g_{b}(x_{0})$ for
$i\in[n_{A}]$ and $k,{k^{\prime}}\in[K_{A}]$ such that $k\neq{k^{\prime}}$,
and $\exp\left(\bm{h}_{k,i}^{\top}\mathbf{w}_{k^{\prime}}\right)=g_{c}(x_{0})$
for $i\in[n_{A}]$, $k\in[K_{A}]$, and ${k^{\prime}}\in[K_{B}]$. Moreover,
$(\mathbf{H},\mathbf{W})$ can achieve the equality in (63). Finally, following
a same argument as Case (ii), we have that (1) $L_{2}$ is the global minimum
of (52); (2) any minimizer satisfies that $\mathbf{h}_{k,i}=\mathbf{0}_{p}$
for all $k\in[K_{A}+1:K]$ and $i\in[n_{B}]$; (3) for any feasible solution
$\left(\bm{H}^{\prime},\mathbf{W}^{\prime}\right)$, if there exist
$k,{k^{\prime}}\in[K_{A}+1:K]$ such that
$\left\|\mathbf{w}_{k}-\mathbf{w}_{k^{\prime}}\right\|=\varepsilon>0$, then
(62) holds.

Combining the three cases, we obtain Lemma 4, completing the proof. ∎

######  Lemma 5.

For any constants $C_{a}>0$, $C_{b}>0$, and $C_{c}>0$, define
$C_{a}^{\prime}:=\frac{C_{a}}{C_{a}+(K_{A}-1)C_{b}+K_{B}C_{c}}\in(0,1)$,
$C_{b}^{\prime}:=\frac{C_{b}}{C_{a}+(K_{A}-1)C_{b}+K_{B}C_{c}}\in(0,1)$, and
$C_{c}^{\prime}:=\frac{C_{c}}{C_{a}+(K_{A}-1)C_{b}+K_{B}C_{c}}\in(0,1)$,
$C_{d}:=-C_{a}^{\prime}\log(C_{a}^{\prime})-C_{b}^{\prime}(K_{A}-1)\log(C_{b}^{\prime})-K_{B}C_{c}^{\prime}\log(C_{c}^{\prime})$,
$C_{e}:=\frac{K_{A}C_{b}}{K_{A}C_{b}+K_{B}C_{c}}\in(0,1)$,
$C_{f}:=\frac{K_{B}C_{c}}{K_{A}C_{b}+K_{B}C_{c}}\in(0,1)$, and
$C_{g}:=\frac{K_{A}C_{b}+K_{B}C_{c}}{C_{a}+(K_{A}-1)C_{b}+K_{B}C_{c}}>0$. For
any feasible solution $(\bm{H},\mathbf{W})$ of (52), the objective value of
(52) can be bounded from below by:

$\displaystyle\frac{1}{K_{A}n_{A}}\sum_{k=1}^{K_{A}}\sum_{i=1}^{n_{A}}\mathcal{L}(\mathbf{W}\bm{h}_{k,i},\mathbf{y}_{k})$
(64) $\displaystyle\overset{a}{\geq}$
$\displaystyle-\frac{C_{g}}{K_{A}}\sqrt{KE_{H}}\sqrt{\sum_{k=1}^{K_{A}}\left\|C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k}\right\|^{2}}+C_{d}$
$\displaystyle\overset{b}{\geq}$
$\displaystyle-\frac{C_{g}}{K_{A}}\sqrt{KE_{H}}\sqrt{KE_{W}\\!-K_{A}\left(1/K_{R}-C_{f}^{2}-\frac{C_{f}^{4}}{C_{e}(2-C_{e})}\right)\|\mathbf{w}_{B}\|^{2}-\\!\sum_{k=K_{A}+1}^{K}\left\|\mathbf{w}_{k}-\mathbf{w}_{B}\right\|^{2}}+C_{d},$

where $\mathbf{w}_{A}:=\frac{1}{K_{A}}\sum_{k=1}^{K_{A}}\mathbf{w}_{k}$,
$\mathbf{w}_{B}:=\frac{1}{K_{B}}\sum_{k=K_{A}+1}^{K}\mathbf{w}_{k}$, and
$K_{R}:=\frac{K_{A}}{K_{B}}$. Moreover, the equality in $\overset{a}{\geq}$
hold only if $\mathbf{h}_{k,i}=\mathbf{0}_{p}$ for all $k\in[K_{A}+1:K]$.

###### Proof of Lemma 5.

For $k\in[K_{A}]$ and $i\in[n_{k}]$, we introduce
$\mathbf{z}_{k,i}=\mathbf{W}\bm{h}_{k,i}$. Because that
$C_{a}^{\prime}+(K_{A}-1)C_{b}^{\prime}+K_{B}C_{c}^{\prime}=1$,
$C_{a}^{\prime}>0$, $C_{b}^{\prime}>0$, and $C_{c}^{\prime}>0$, by the
concavity of $\log(\cdot)$, we have

$\displaystyle\begin{aligned}
&-\log\left(\frac{\exp(\mathbf{z}_{k,i}(i))}{\sum_{{k^{\prime}}=1}^{K}\exp(\mathbf{z}_{{k^{\prime}},i}(k))}\right)\\\
=&-\mathbf{z}_{k,i}(k)+\log\left(C_{a}^{\prime}\left(\frac{\exp(z_{k,i}(k))}{C_{a}^{\prime}}\right)+\sum_{{k^{\prime}}=1,~{}{k^{\prime}}\neq
k}^{K_{A}}C_{b}^{\prime}\left(\frac{\exp(z_{k,i}({k^{\prime}}))}{C_{b}^{\prime}}\right)+\sum_{{k^{\prime}}=K_{A}+1}^{K}C_{c}^{\prime}\left(\frac{\exp(z_{k,i}({k^{\prime}}))}{C_{c}^{\prime}}\right)\right)\\\
\geq&-\mathbf{z}_{k,i}(k)+C_{a}^{\prime}\mathbf{z}_{k,i}(k)+C_{b}^{\prime}\sum_{{k^{\prime}}=1,~{}{k^{\prime}}\neq
k}^{K_{A}}\mathbf{z}_{k,i}({k^{\prime}})+C_{C}^{\prime}\sum_{{k^{\prime}}=K_{A}+1}^{K}\mathbf{z}_{i,j}(k)+C_{d}\\\
=&C_{g}C_{e}\left(\frac{1}{K_{A}}\sum_{{k^{\prime}}=1}^{K_{A}}\mathbf{z}_{k,i}({k^{\prime}})-\mathbf{z}_{k,i}(k)\right)+C_{g}C_{f}\left(\frac{1}{K_{B}}\sum_{{k^{\prime}}=K_{A}+1}^{K}\mathbf{z}_{k,i}({k^{\prime}})-\mathbf{z}_{k,i}(k)\right)+C_{d}.\end{aligned}$

Therefore, integrating (A.2.2) with $k\in[K_{A}]$ and $i\in[n_{A}]$, recalling
that $\mathbf{w}_{A}=\frac{1}{K_{A}}\sum_{k=1}^{K_{A}}\mathbf{w}_{k}$ and
$\mathbf{w}_{B}=\frac{1}{K_{B}}\sum_{k=K_{A}+1}^{K}\mathbf{w}_{k}$, we have

$\displaystyle\frac{1}{K_{A}n_{A}}\sum_{k=1}^{K_{A}}\sum_{i=1}^{n_{A}}\mathcal{L}(\mathbf{W}\bm{h}_{k,i},\mathbf{y}_{k})$
(65) $\displaystyle\geq$
$\displaystyle\frac{1}{K_{A}n_{A}}\sum_{k=1}^{K_{A}}\sum_{i=1}^{n_{A}}C_{g}\left[C_{e}(\bm{h}_{k,i}\mathbf{w}_{A}-\bm{h}_{k,i}\mathbf{w}_{k})+C_{f}(\bm{h}_{k,i}\mathbf{w}_{B}-\bm{h}_{k,i}\mathbf{w}_{k})\right]+C_{d}$
$\displaystyle\overset{a}{=}$
$\displaystyle\frac{C_{g}}{K_{A}}\sum_{k=1}^{K_{A}}\bm{h}_{k}^{\top}(C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k})+C_{d},$

where in $\overset{a}{=}$, we introduce
$\bm{h}_{k}:=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\bm{h}_{k,i}$ for $k\in[K]$, and
use $C_{e}+C_{f}=1$. Then it is sufficient to bound
$\sum_{k=1}^{K_{A}}\bm{h}_{k}^{\top}(C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k})$.
By the Cauchy–Schwarz inequality, we have

$\displaystyle\sum_{k=1}^{K_{A}}\bm{h}_{k}^{\top}(C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k})\geq$
$\displaystyle-\sqrt{\sum_{k=1}^{K_{A}}\|\bm{h}_{k}\|^{2}}\sqrt{\sum_{k=1}^{K_{A}}\left\|C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k}\right\|^{2}}$
$\displaystyle\overset{a}{\geq}$
$\displaystyle-\sqrt{\sum_{k=1}^{K_{A}}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\|\bm{h}_{k,i}\|^{2}}\sqrt{\sum_{k=1}^{K_{A}}\left\|C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k}\right\|^{2}}$
$\displaystyle\overset{b}{\geq}$
$\displaystyle-\sqrt{KE_{H}}\sqrt{\sum_{k=1}^{K_{A}}\left\|C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k}\right\|^{2}},$
(66)

where $\overset{a}{\geq}$ follows from Jensens’s inequality
$\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\|\bm{h}_{k,i}\|^{2}\geq\bm{h}_{k}$ for
$k\in[K_{A}]$ and $\overset{b}{\geq}$ uses the constraint that
$\frac{1}{K}\sum_{k=1}^{K}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\left\|\bm{h}_{k,i}\right\|^{2}\leq
E_{H}$. Moreover, we have
$\sum_{k=1}^{K_{A}}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\left\|\bm{h}_{k,i}\right\|^{2}=E_{H}$
only if $\mathbf{h}_{k,i}=\mathbf{0}_{p}$ for all $k\in[K_{A}+1:K]$. Plugging
(66) into (65), we obtain $\overset{a}{\geq}$ in (64).

We then bound
$\sum_{k=1}^{K_{A}}\left\|C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k}\right\|^{2}$.
We have

$\displaystyle\frac{1}{K_{A}}\sum_{k=1}^{K_{A}}\left\|C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}-\mathbf{w}_{k}\right\|^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{K_{A}}\sum_{k=1}^{K_{A}}\|\mathbf{w}_{k}\|^{2}-2\frac{1}{K_{A}}\sum_{k=1}^{K_{A}}\mathbf{w}_{k}\cdot(C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B})+\|C_{e}\mathbf{w}_{A}+C_{f}\mathbf{w}_{B}\|^{2}$
$\displaystyle\overset{a}{=}$
$\displaystyle\frac{1}{K_{A}}\sum_{k=1}^{K_{A}}\|\mathbf{w}_{k}\|^{2}-2C_{f}^{2}\mathbf{w}_{A}^{\top}\mathbf{w}_{B}-C_{e}(2-C_{e})\|\mathbf{w}_{A}\|^{2}+C_{f}^{2}\|\mathbf{w}_{B}\|^{2}.$
(67)

where $\overset{a}{=}$ uses
$\sum_{k=1}^{K_{A}}\mathbf{w}_{k}=K_{A}\mathbf{w}_{A}$. Then using the
constraint that $\sum_{k=1}^{K}\|\mathbf{w}_{k}\|\leq KE_{W}$ yields that

$\displaystyle\frac{1}{K_{A}}\sum_{k=1}^{K_{A}}\|\mathbf{w}_{k}\|^{2}-2C_{f}^{2}\mathbf{w}_{A}^{\top}\mathbf{w}_{B}-C_{e}(2-C_{e})\|\mathbf{w}_{A}\|^{2}+C_{f}^{2}\|\mathbf{w}_{B}\|^{2}$
(68) $\displaystyle\leq$
$\displaystyle\frac{K}{K_{A}}E_{W}^{2}-\frac{1}{K_{A}}\sum_{k=K_{A}+1}^{K}\\!\|\mathbf{w}_{k}\|^{2}-C_{e}(2-C_{f})\left\|\mathbf{w}_{A}+\frac{C_{f}^{2}}{C_{e}(2-C_{e})}\mathbf{w}_{B}\right\|^{2}\\!\\!+\\!\left(C_{f}^{2}+\frac{C_{f}^{4}}{C_{e}(2-C_{e})}\right)\|\mathbf{w}_{B}\|^{2}$
$\displaystyle\overset{a}{=}$
$\displaystyle\frac{K}{K_{A}}E_{W}^{2}-\left(1/K_{R}-C_{f}^{2}-\frac{C_{f}^{4}}{C_{e}(2-C_{e})}\right)\|\mathbf{w}_{B}\|^{2}-\frac{1}{K_{A}}\sum_{k=K_{A}+1}^{K}\\!\left\|\mathbf{w}_{k}-\mathbf{w}_{B}\right\|^{2},$

where $\overset{a}{\geq}$ applies
$\sum_{k=K_{A}+1}^{K}\|\mathbf{w}_{k}\|^{2}=K_{B}\|\mathbf{w}_{B}\|^{2}+\sum_{k=K_{A}+1}^{K}\left\|\mathbf{w}_{k}-\mathbf{w}_{B}\right\|^{2}$.
Plugging (A.2.2) and (68) into $\overset{a}{\geq}$ in (64), we obtain
$\overset{b}{\geq}$ in (64), completing the proof. ∎

##  Appendix B Additional Results

##### Comparison of Oversampling and Weighted Adjusting.

Oversampling and weight adjusting are two commonly-used tricks in deep
learning [JK19]. Both of them actually consider the same objective as (16),
but applies different optimization algorithms to minimize the objective. It
was observed that oversampling is more stable than weight adjusting in
optimization. As a by product of this work, we compare the two algorithms
below and shows that the variance of updates for oversampling will be
potentially much smaller than that of weight adjusting. It was well-known in
stochastic optimization field that the variance of the updates decides the
convergence of an optimization algorithm (see e.g, [BCN18, FLLZ18, FLZ19]).
Thus we offer a reasonable justification for the stability of the oversampling
technique. We simply consider sampling the training data without replacement.
It slightly differs from the deep learning training methods in practice.
Besides, we only consider sampling a single data in each update. The analysis
can be directly extended to the mini-batch setting.

We first introduce the two methods. The weight adjusting algorithm in each
update randomly samples a training data, and updates the parameters
$\bm{W}_{\textnormal{full}}$ by the Stochastic Gradient Descent algorithm as

$\displaystyle\bm{W}_{\textnormal{full}}^{t+1}=\bm{W}_{\textnormal{full}}^{t}-\eta_{w}\mathbf{v}_{w}^{t},\quad
t=0,1,2,\dots,$ (69)

where $\bm{W}_{\textnormal{full}}^{t}$ denotes the parameters at iteration
step $t$, $\eta_{w}$ is a positive step size, and the stochastic gradient
$\mathbf{v}_{w}^{t}$ satisfies that

$\mathbf{v}_{w}^{t}=\begin{cases}\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k}),&k\in[K_{A}],i\in[n_{A}],\text{~{}with
probability~{}}\frac{1}{K_{A}n_{A}+K_{B}n_{B}},\\\
w_{r}\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k}),&k\in[K_{A}+1:K_{B}],i\in[n_{B}],\text{~{}with
probability~{}}\frac{1}{K_{A}n_{A}+K_{B}n_{B}}.\end{cases}$

We have

$\displaystyle\mathbb{E}\left[\mathbf{v}_{w}^{t}\mid\bm{W}_{\textnormal{full}}^{t}\right]$
(70) $\displaystyle=$
$\displaystyle\frac{1}{n_{A}K_{A}+n_{B}K_{B}}\left[\sum_{k=1}^{K_{A}}\sum_{i=1}^{n_{A}}\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k})+w_{r}\\!\\!\sum_{k=K_{A}+1}^{K}\\!\sum_{i=1}^{n_{B}}\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k})\right],$

and

$\displaystyle\mathbb{E}\left[\|\mathbf{v}_{w}^{t}\|^{2}\mid\bm{W}_{\textnormal{full}}^{t}\right]=$
$\displaystyle\frac{1}{n_{A}K_{A}+n_{B}K_{B}}\sum_{k=1}^{K_{A}}\sum_{i=1}^{n_{A}}\left\|\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k})\right\|^{2}$
$\displaystyle+\frac{w_{r}^{2}}{n_{A}K_{A}+n_{B}K_{B}}\sum_{k=K_{A}+1}^{K}\sum_{i=1}^{n_{B}}\left\|\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k})\right\|^{2}.$
(71)

For the oversampling method, the algorithm in effect duplicates the data by
$w_{r}$ times and runs Stochastic Gradient Descent on the “whole” data.
Therefore, the update goes as

$\displaystyle\bm{W}_{\textnormal{full}}^{t+1}=\bm{W}_{\textnormal{full}}^{t}-\eta_{s}\mathbf{v}_{s}^{t},\quad
t=0,1,2,\dots,$ (72)

where $\mathbf{v}_{s}^{t}$ satisfies that

$\mathbf{v}_{s}^{t}=\begin{cases}\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k}),&k\in[K_{A}],i\in[n_{A}],\text{~{}with
probability~{}}\frac{1}{K_{A}n_{A}+K_{B}w_{r}n_{B}},\\\
\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k}),&k\in[K_{A}+1:K_{B}],i\in[n_{B}],\text{~{}with
probability~{}}\frac{w_{r}}{K_{A}n_{A}+K_{B}w_{r}n_{B}}.\end{cases}$

We obtain

$\displaystyle\mathbb{E}\left[\mathbf{v}_{s}^{t}\mid\bm{W}_{\textnormal{full}}^{t}\right]=$
$\displaystyle\frac{1}{n_{A}K_{A}+w_{r}n_{B}K_{B}}\sum_{k=1}^{K_{A}}\sum_{i=1}^{n_{A}}\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k})$
$\displaystyle+\frac{w_{r}}{n_{A}K_{A}+w_{r}n_{B}K_{B}}\sum_{k=K_{A}+1}^{K}\sum_{i=1}^{n_{B}}\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k}),$

and

$\displaystyle\mathbb{E}\left[\|\mathbf{v}_{s}^{t}\|^{2}\mid\bm{W}_{\textnormal{full}}^{t}\right]=$
$\displaystyle\frac{1}{n_{A}K_{A}+w_{r}n_{B}K_{B}}\sum_{k=1}^{K_{A}}\sum_{i=1}^{n_{A}}\left\|\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k})\right\|^{2}$
$\displaystyle+\frac{w_{r}}{n_{A}K_{A}+w_{r}n_{B}K_{B}}\sum_{k=K_{A}+1}^{K}\sum_{i=1}^{n_{B}}\left\|\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k})\right\|^{2}.$
(73)

We suppose the two updates in expectation are in a same scale. That means we
assume
$\eta_{w}=\frac{n_{A}K_{A}+w_{r}n_{B}K_{B}}{n_{A}K_{A}+n_{B}K_{B}}\eta_{s}$.
Then
$\eta_{w}\mathbb{E}\left[\mathbf{v}_{w}^{t}\mid\bm{W}_{\textnormal{full}}^{t}\right]=\eta_{s}\mathbb{E}\left[\mathbf{v}_{s}^{t}\mid\bm{W}_{\textnormal{full}}^{t}\right]$.
In fact, if $K_{A}\asymp 1$, $K_{B}\asymp 1$, $n_{A}\gg n_{B}$, and $1\ll
w_{r}\lesssim\left(n_{A}/n_{B}\right)$, we have
$\frac{n_{A}K_{A}+w_{r}n_{B}K_{B}}{n_{A}K_{A}+n_{B}K_{B}}\asymp 1$ and so
$\eta_{w}\asymp\eta_{s}$. Now by comparing (B) with (B), we obtain that the
second moment of $\eta_{w}\mathbf{v}_{w}^{t}$ is much smaller than that of
$\eta_{s}\mathbf{v}_{s}^{t}$ since the order of $w_{r}$ for the latter is
larger by $1$. For example, let us assume that all the norms of the gradients
are in a same order, i.e.,
$\left\|\nabla_{\bm{W}_{\textnormal{full}}}\mathcal{L}(f(\mathbf{x}_{k,i};\bm{W}_{\textnormal{full}}^{t}),\mathbf{y}_{k})\right\|\asymp
a$ for all $k$ and $i$, where $a>0$. Then (B) implies that
$\mathbb{E}\left[\|\mathbf{v}_{s}^{t}\|^{2}\mid\bm{W}_{\textnormal{full}}^{t}\right]\asymp\eta_{s}^{2}a^{2}$.
However, (B) reads that
$\mathbb{E}\left[\|\mathbf{v}_{w}^{t}\|^{2}\mid\bm{W}_{\textnormal{full}}^{t}\right]\asymp\eta_{s}^{2}\frac{n_{A}K_{A}+w_{r}^{2}n_{B}K_{B}}{n_{A}K_{A}+w_{r}n_{B}K_{B}}a^{2}$.
Furthermore, if we set $w_{r}\asymp n_{A}/n_{B}$, then
$\mathbb{E}\left[\|\mathbf{v}_{w}^{t}\|^{2}\mid\bm{W}_{\textnormal{full}}^{t}\right]\asymp\eta_{s}^{2}w_{r}a^{2}$.
Thus the second moment for $\eta_{w}\mathbf{v}_{w}^{t}$ is around $w_{r}$
times of that for $\eta_{s}\mathbf{v}_{s}^{t}$. And this fact also holds for
the variance because
$\left\|\eta_{s}\mathbb{E}\left[\mathbf{v}_{s}^{t}\mid\bm{W}_{\textnormal{full}}^{t}\right]\right\|\asymp\eta_{s}a$
and the property that
$\mathbb{E}\|\mathbf{x}-\mathbb{E}[\mathbf{x}]\|^{2}=\mathbb{E}\|\mathbf{x}\|^{2}-\|\mathbb{E}[\mathbf{x}]\|^{2}$
for any random variable $\mathbf{x}$. Therefore, we can conclude that the
variance of updates for oversampling is potentially much smaller than that of
weight adjusting.

##### More Discussions on Convex Relaxation and Cross-Entropy Loss.

We show Program (7) can also be relaxed as a nuclear norm-constrained convex
optimization. The result heavily relies on the progress of matrix
decomposition, e.g. [BMP08, HV19]. We will the use the equality (see e.g.,
[BMP08, Section 2]) that for any matrix $\mathbf{Z}$ and $a>0$,

$\|\mathbf{Z}\|_{*}=\inf_{r\in\mathbb{N}_{+}}\inf_{\mathbf{U},\mathbf{V}:\mathbf{U}\mathbf{V}^{\top}=\mathbf{Z}}\frac{a}{2}\|\mathbf{U}\|^{2}+\frac{1}{2a}\|\mathbf{V}\|^{2},$
(74)

where $r$ is the number of columns for $\mathbf{U}$ and $\|\cdot\|_{*}$
denotes the nuclear norm.

For any feasible solution $\left(\bm{H},\mathbf{W}\right)$ for the original
program (7), we define

$\mathbf{h}_{k}=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\mathbf{h}_{k,i},~{}k\in[K],\quad\tilde{\bm{H}}=[\mathbf{h}_{1},\mathbf{h}_{2},\dots,\mathbf{h}_{K}]\in\mathbb{R}^{p\times
K},~{}~{}\text{and}~{}~{}\mathbf{Z}=\mathbf{W}\tilde{\bm{H}}\in\mathbb{R}^{K\times
K}.$ (75)

We consider the convex program:

$\displaystyle\min_{\mathbf{Z}\in\mathbb{R}^{K\times K}}$
$\displaystyle\sum_{k=1}^{K}\frac{n_{k}}{N}\mathcal{L}(\mathbf{Z}_{k},\mathbf{y}_{k})$
(76) $\displaystyle\mathrm{s.t.}$ $\displaystyle\|\mathbf{Z}\|_{*}\leq
K\sqrt{E_{H}E_{W}}.$

where $\mathbf{Z}_{k}$ denotes the $k$-th column of $\mathbf{Z}$ for
$k\in[K]$.

######  Lemma 6.

Assume $p\geq K$ and the loss function $\mathcal{L}$ is convex on the first
argument. Let $\mathbf{Z}^{\star}$ be a minimizer of the convex program (76).
Let $r$ be the rank of $\mathbf{Z}^{\star}$ and consider thin Singular Value
Decomposition (SVD) of $\mathbf{Z}^{\star}$ as
$\mathbf{Z}^{\star}=\mathbf{U}^{\star}\mathbf{\Sigma}^{\star}\mathbf{V}^{\star}$.
Introduce two diagonal matrices $\mathbf{\Sigma}_{1}^{\star}$ and
$\mathbf{\Sigma}_{2}^{\star}$ with the entries defined as
$\mathbf{\Sigma}_{1}^{\star}(i,i)=\sqrt{\frac{E_{W}}{E_{H}}}\sqrt{|\mathbf{\Sigma}^{\star}(i,i)|}$
and
$\mathbf{\Sigma}_{2}^{\star}(i,i)=\sqrt{\frac{E_{H}}{E_{W}}}\mathbf{\Sigma}^{\star}(i,i)/\sqrt{|\mathbf{\Sigma}^{\star}(i,i)|}$
for $i\in[r]$, respectively. Let
$\left(\mathbf{H}^{\star},\mathbf{W}^{\star}\right)$ be

$\displaystyle\mathbf{W}=\mathbf{U}^{\star}\mathbf{\Sigma}_{1}^{\star}\mathbf{P}^{\top},\quad\left[\bm{h}_{1}^{\star},\bm{h}_{2}^{\star},\dots,\bm{h}_{K}^{\star}\right]=\mathbf{P}\mathbf{\Sigma}_{2}^{\star}\mathbf{V}^{\star},$
(77) $\displaystyle\bm{h}_{k,i}^{\star}=\bm{h}_{k}^{\star},\quad
k\in[K],~{}i\in[n_{k}],$

where $\mathbf{P}\in\mathbb{R}^{p\times r}$ is any partial orthogonal matrix
such that $\mathbf{P}^{\top}\mathbf{P}=\mathbf{I}_{r}$. Then
$(\mathbf{H}^{\star},\mathbf{W}^{\star})$ is a minimizer of (7).

###### Proof of Lemma 6.

For any feasible solution $\left(\bm{H},\mathbf{W}\right)$ for the original
program (7), define $\mathbf{h}_{k}$ for $k\in[K]$, $\tilde{\bm{H}}$, and
$\mathbf{Z}$ by (75). We show $\mathbf{Z}$ is a feasible solution for the
convex program (76). In fact, by (74) with $r=K$ and $a=\sqrt{E_{H}/E_{W}}$,
we have

$\displaystyle\left\|\mathbf{Z}\right\|_{*}$
$\displaystyle\leq\frac{\sqrt{E_{H}/E_{W}}}{2}\left\|\mathbf{W}\right\|^{2}+\frac{\sqrt{E_{W}/E_{H}}}{2}\left\|\tilde{\bm{H}}\right\|^{2}$
$\displaystyle\overset{a}{\leq}\frac{\sqrt{E_{H}/E_{W}}}{2}\sum_{k=1}^{K}\|\mathbf{w}_{k}\|^{2}+\frac{\sqrt{E_{W}/E_{H}}}{2}\sum_{k=1}^{K}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\left\|\bm{h}_{k,i}\right\|^{2}$
$\displaystyle\leq K\sqrt{E_{H}E_{W}},$ (78)

where $\overset{a}{\leq}$ applies Jensen’s inequality as:

$\left\|\tilde{\bm{H}}\right\|^{2}=\sum_{k=1}^{K}\|\mathbf{h}_{k}\|^{2}\leq\sum_{k=1}^{K}\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\left\|\bm{h}_{k,i}\right\|^{2}.$

Let $L_{0}$ be the global minimum of the convex problem (76). Since
$\mathcal{L}$ is convex on the first argument, by the same argument as
(A.2.1), we obtain, for any feasible solution
$\left(\bm{H},\mathbf{W}\right)$,

$\displaystyle\frac{1}{N}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\mathcal{L}(\mathbf{W}\bm{h}_{k,i},\mathbf{y}_{k})$
$\displaystyle=\sum_{k=1}^{K}\frac{n_{k}}{N}\left[\frac{1}{n_{k}}\sum_{k=1}^{n_{k}}\mathcal{L}(\mathbf{W}\bm{h}_{k,i},\mathbf{y}_{k})\right]$
$\displaystyle\geq\sum_{k=1}^{K}\frac{n_{k}}{N}\mathcal{L}(\mathbf{W}\bm{h}_{k},\mathbf{y}_{k})=\sum_{k=1}^{K}\frac{n_{k}}{N}\mathcal{L}(\mathbf{Z}_{k},\mathbf{y}_{k})\geq
L_{0}.$ (79)

On the other hand, for the solution
$\left(\bm{H}^{\star},\mathbf{W}^{\star}\right)$ defined in (77) with
$\mathbf{Z}^{\star}$, we can verify that
$\left(\bm{H}^{\star},\mathbf{W}^{\star}\right)$ is a feasible solution for
(7) and

$\frac{1}{N}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\mathcal{L}(\mathbf{W}^{\star}\bm{h}_{k,i}^{\star},\mathbf{y}_{k})=\sum_{k=1}^{K}\frac{n_{k}}{N}\mathcal{L}(\mathbf{Z}_{k}^{\star},\mathbf{y}_{k})=L_{0}.$
(80)

Combining (B) and (80), we have that $L_{0}$ is the global minimum of (7) and
$(\mathbf{H}^{\star},\mathbf{W}^{\star})$ is a minimizer. ∎

######  Property 1.

For the cross-entropy loss, we have the following properties.

  1. (A)

Any minimizer $\mathbf{Z}^{\star}$ of (76) satisfies that
$\|\mathbf{Z}\|_{*}=\sqrt{E_{H}E_{W}}$.

  2. (B)

Any minimizer $(\bm{H}^{\star},\mathbf{W}^{\star})$ of (7) satisfies

$\frac{1}{K}\sum_{k=1}^{K}\frac{1}{n}\sum_{i=1}^{n}\left\|\bm{h}_{k,i}^{\star}\right\|^{2}=E_{H},\quad\text{and}\quad\quad\frac{1}{K}\sum_{k=1}^{K}\left\|\mathbf{w}_{k}^{\star}\right\|^{2}=E_{W}.$

  3. (C)

Any minimizer $\mathbf{X}^{\star}$ of (14) satisfies that

$\frac{1}{K}\sum_{k=1}^{K}\mathbf{X}^{\star}(k,k)=E_{H},\quad\text{and}\quad\quad\frac{1}{K}\sum_{k=K+1}^{2K}\mathbf{X}^{\star}(k,k)=E_{W}.$

###### Proof of Property 1.

We first prove (A). Let $\mathbf{Z}^{\star}$ be any minimier of (76). Then by
the Karush–Kuhn–Tucker conditions, there is a pair $(\lambda,\mathbf{\xi})$
with $\lambda\geq 0$ and $\mathbf{\xi}\in\partial\|\mathbf{Z}^{\star}\|_{*}$
such that

$\nabla_{\mathbf{Z}}\left[\sum_{k=1}^{K}\frac{n_{k}}{N}\mathcal{L}(\mathbf{Z}_{k}^{\star},\mathbf{y}_{k})\right]+\lambda\mathbf{\xi}=\mathbf{0}^{K\times
K},$

where $\partial\|\mathbf{Z}\|_{*}$ denotes the set of sub-gradient of
$\|\mathbf{Z}\|_{*}$. For the cross-entropy loss, one can verify that
$\nabla_{\mathbf{Z}}\left[\sum_{k=1}^{K}\frac{n_{k}}{N}\mathcal{L}(\mathbf{Z}_{k},\mathbf{y}_{k})\right]\neq\mathbf{0}^{K\times
K}$ for all $\mathbf{Z}$. So $\lambda\neq 0$. By the complementary slackness
condition, we have $\mathbf{Z}$ will reach the boundary of the constraint,
achieving (A).

For (B), suppose there is a minimizer $(\bm{H}^{\star},\mathbf{W}^{\star})$ of
(7) such that
$\frac{1}{K}\sum_{k=1}^{K}\frac{1}{n}\sum_{i=1}^{n}\left\|\bm{h}_{k,i}^{\star}\right\|^{2}<E_{H}$
or $\frac{1}{K}\sum_{k=1}^{K}\left\|\mathbf{w}_{k}^{\star}\right\|^{2}<E_{W}$.
Letting $\mathbf{Z}^{\star}$ defined as (75), it follows from (B) that
$\mathbf{Z}^{\star}$ is a minimizer of (76). However, by (B), we have
$\|\mathbf{Z}^{\star}\|_{*}<\sqrt{E_{H}E_{W}}$, which is contradictory to (A).
We obtain (B).

For (C), suppose there is a minimizer $\mathbf{X}^{\star}$ of (14) such that
$\frac{1}{K}\sum_{k=1}^{K}\mathbf{X}^{\star}(k,k)<E_{H}$ or
$\frac{1}{K}\sum_{k=K+1}^{2K}\mathbf{X}^{\star}(k,k)<E_{W}$. Then letting
$(\bm{H}^{\star},\mathbf{W}^{\star})$ defined in (15),
$(\bm{H}^{\star},\mathbf{W}^{\star})$ is a minimizer of (7) from Theorem 1.
However, we have
$\frac{1}{K}\sum_{k=1}^{K}\frac{1}{n}\sum_{i=1}^{n}\left\|\bm{h}_{k,i}^{\star}\right\|^{2}<E_{H}$
or $\frac{1}{K}\sum_{k=1}^{K}\left\|\mathbf{w}_{k}^{\star}\right\|^{2}<E_{W}$,
which is contradictory to (B). We complete the proof.

∎


# BSRA: Block-based Super Resolution Accelerator with Hardware Efficient Pixel
Attention

Dun-Hao Yang, and Tian-Sheuan Chang, Senior Member, IEEE Institute of
Electronics, National Yang Ming Chiao Tung University,
Hsinchu, Taiwan

###### Abstract

Increasingly, convolution neural network (CNN) based super resolution models
have been proposed for better reconstruction results, but their large model
size and complicated structure inhibit their real-time hardware
implementation. Current hardware designs are limited to a plain network and
suffer from lower quality and high memory bandwidth requirements. This paper
proposes a super resolution hardware accelerator with hardware efficient pixel
attention that just needs 25.9K parameters and simple structure but achieves
0.38dB better reconstruction images than the widely used FSRCNN. The
accelerator adopts full model block wise convolution for full model layer
fusion to reduce external memory access to model input and output only. In
addition, CNN and pixel attention are well supported by PE arrays with
distributed weights. The final implementation can support full HD image
reconstruction at 30 frames per second with TSMC 40nm CMOS process.

###### Index Terms:

Convolution neural network, deep learning accelerators, pixel attention
mechanism, super resolution

##  I Introduction

Convolutional neural network (CNN) based super resolution (SR) [1] is getting
popular in recent years because of better reconstructed high resolution images
over traditional interpolation methods. These SR models are getting deeper and
wider and use complicated structure to get better performance[2, 3, 4, 5] as
shown in Fig. 1. However, SR models suffer from high computational complexity
and memory bandwidth because their feature sizes are not getting smaller over
layers and they need to process large input size. Thus, hardware acceleration
is demanded for real-time applications.

Various hardware accelerators have been proposed[6, 7, 8]. However, current
designs only use plain networks like the widely used FSRCNN[9] instead of
recent large size and complicated models due to the high model complexity,
which limits their reconstructed image quality. In addition, most of the
designs use layer-by-layer processing that will need to store intermediate
data to DRAM and load it back for each layer or large buffer size, which is
extremely bad for SR applications due to their large feature size.

To address the above problems, we propose a full model block convolution based
super resolution accelerator (BSRA) with hardware efficient pixel attention
(HPAN) model. The HPAN model uses pixel attention for better image quality
than FSRCNN but keeps the structure simple for small model size and low
complexity. The accelerator uses block convolution for the whole model instead
of certain layers as in [10] to enable the whole model fusion to reduce the
external memory bandwidth to model input and output only and needs a small on-
chip buffer size. The required convolution and pixel attention are well
supported by the proposed distributed weight PE array. The final
implementation can achieve real-time full HD image throughput with TSMC 40nm
CMOS process.

##  II Proposed SR Model with Hardware Efficient Pixel Attention

Figure 1: Image quality vs model size for different SR models Figure 2: The
proposed HPAN model

###  II-A Related Work

Since the first CNN based SR model, SRCNN [1], was proposed, many deep
learning skills are added to the super resolution model to get better image
quality. CARN [3] cascades the residual in residual blocks to get multilevel
representations. RCAN [2], which is over 400 layers, uses residual blocks and
channel attention to learn more channel-wise information. SAN [4] proposes
non-locally enhanced residual group to gain long distance and structure
message. These models are getting deeper and wider for better images.

Fig. 1 the scatter chart of the SR model size and its image quality. This
figure shows that the quality of the reconstructed image is increasing with
model size, but the increment of parameters is not efficient at million scale.
Take FSRCNN [9] and RCAN [2] for example, the difference in performance on
Set5 between them is only 1.3 dB, but the model size of the latter one is
about 1300 larger than the former one. The model size increment is very
inefficient so that most of hardware designs [6, 7, 8] use plain and shallow
networks just like SRCNN [1] and FSRCNN [9]. For example, in [8], they uses a
classification network to judge whether image tiles are with more high
frequency signal and send them to appropriate plain networks. This mechanism
can reduce the computation overhead of the SR part, but the classification
part is much bigger than the SR part. Therefore, this paper tries to find a
better architecture which has a trade-off between model size and reconstructed
image quality for single image super resolution.

###  II-B Proposed HPAN Model

Fig. 2 shows the proposed HPAN model based on convolution pixel attention
mechanism and the number of parameters is just about 26K. The model design is
based on the following observations. First, for the overall model structure,
skip or dense connection is popular in most of the models to learn the
residual between high resolution (HR) images and SR images, but this will need
extra memory access for hardware implementation. Thus, instead of the residual
architecture, we take the attention mechanism into consideration. In Fig. 1,
RCAN [2] and SAN [4] are the best two models in performance with different
attention methods. However, the trend line shows that pixel attention (PAN)
[5] is better. Thus, we use pixel attention to construct our model. However,
[5] has a local skip connection in the basic block and a global interpolated
image addition. Besides, in the up-sampling stage, the final image is not
generated immediately because the interpolated features still need to be tuned
by some convolution layers and thus needs large computation. These are not
hardware friendly.

Based on the above observations, we propose the hardware-friendly pixel
attention model with 26K parameters as shown in Fig. 2. The model consists of
three stages: feature extraction, mapping, and reconstruction. The feature
extraction stage uses one convolution layer. The mapping stage is based on the
Clamping Pixel Attention Block (CPAB) in (4) to extract features. After
cascading some CPABs, we use transpose convolution to reconstruct the image
directly to prevent the problem mentioned before.

The representation of CPAB in Fig. 2 is shown as follows.

$Clamp(f_{in})=\left\\{\begin{array}[]{lr}+255,&f_{in}\geq+255\\\
-255,&f_{in}\leq-255\\\ f_{in},&\text{otherwise}\end{array}\right.$ (1)

$D(f_{in})=sigmoid(f_{in}/256)$ (2)

$PA(f_{in})=M_{ew}(f_{in},Clamp(D(CL_{1}(f_{in}))))$ (3)

$CPAB(f_{in})=CL_{3}(PA(CL_{1}(f_{in})))$ (4)

where $f_{in}$ is input feature, and $Clamp()$, $D()$, $M_{ew}(,)$, $PA()$,
and $CL_{i}()$ denote clamping function, divisor, element-wise multiplication,
pixel attention and convolution with kernel size $i$, respectively. In this
model, with (2), the weight and activation distribution will be more
reasonable after training because the range of pixel intensity is from 0 to 1
or 0 to 255. This is applied to every convolution to ensure the proper range.
Besides, since the value of the attention map is smaller than 1, the absolute
value of multiplication results will not exceed 255. The divisor in (2)
chooses 256 for hardware-friendly design as a simple shifting operation.

##  III Proposed Architecture

###  III-A Overview

Fig. 3 shows the proposed system architecture. This design gets weights and
input images by accessing external memory and stores them into weight and
feature SRAM buffers. Both have 32 banks to store kernels of different layers
and intermediate results, respectively, for the computing core. This core
consists of 32 processing element (PE) arrays to process 32 channels of input.
Then the outputs of different channels are summed in a local accumulator and
stored in the partial sum buffer. These partial sums will be further
accumulated in the selective adder for convolution. The results will be
further processed by sigmoid and clamping function according to the proposed
model. The results will be sent to the multiplexer for attention mechanism or
feature memory for the next layer.

Figure 3: System architecture Figure 4: PE arrays

###  III-B PE Array with Distributed Weight

Fig. 4 shows the proposed one PE array, which consists of six PE lines (PEL0
to PEL5) with six PEs in each PELx. Each PE is a multiplier with registers. In
the figure, each PE array is for one channel of convolution in a layer. In
this way, the input feature maps and weights share the same spatial location
in the same PE array. Therefore, the multipliers in each PE get different
features according to their location in the spatial space. In other words, 36
input pixels are allocated to 36 PEs. As for the multiplicands in one PE
array, there are two choices, one is convolutional kernel weights and the
other is pixel-wise attention mask. If the multiplicands are kernel weights,
the weights will broadcast vertically. If the multiplicands are an attention
mask, the mask values will be distributed in PE array. With this distributed
weight scheme, this design can support convolution as well as pixel attention.
After completing the multiplications, the 36 results are transferred to the
next stage to be accumulated.

Figure 5: The three-stage accumulator.

###  III-C Three-Stage Accumulator

Fig. 5 shows the accumulator to accumulate the partial sum, which is
partitioned into three-stage pipelines for minimizing the critical path and
attention mechanism. In which, for our 32 channel computation, eight partial
sums at the same spatial location but different channels are summed up at the
first stage, and their results are further accumulated to complete a 32
channel summation at the second stage. More channel summations are done at the
third stage for the final convolutional results. The selective adder at the
third stage will select the proper input based on convolution or attention
operation for the correct output.

Figure 6: An example of the proposed data flow. In which, for data allocation,
the same colored geometric shapes share the same weights. If kernel size is
$k$, after $k$ clock cycle, the cached feature should be updated. Figure 7:
Comparison of reconstructed images from the Set14 and Urban100 datasets. The
images from left to right are original ground truth, enlarged ground truth,
bicubic interpolation, our floating-point result and quantization result with
block convolution, respectively.

###  III-D Data Flow of Convolution

Fig. 6 shows an example data flow with one PE array with 4 PEs in each PEL for
simplicity. In the proposed data flow, the kernel weights are broadcast
vertically along the PE array row and the feature inputs are transferred to
PEs depending on their spatial locations. With this allocation, the results of
PEs, $P^{sum}_{i,j}$, are belong to the output pixel $O_{i,j}$ at the first
cycle. At the next cycle, the weights in the PE array are circularly shifted
right. The results $P^{sum}_{i,j}$ are belong to $O_{i-1,j}$. At the third
cycle, the weights are shifted again and the results are belonging to
$O_{i-2,j}$. When the first row of weights are fully used, the second row of
weights are loaded into the PE array. By repeating this processing flow, all
partial sums can be generated. These partial sums will be accumulated in the
accumulator and stored in the partial sum buffer. After nine clock cycles, the
feature input will be shifted down by four pixels and will be loaded for
another round of computation until the bottom of this tile. After that, the
feature input will be shifted right by six pixels for computation. With this,
the feature map can be reused to reduce the memory bandwidth due to the large
input in the SR model.

###  III-E Block Convolution for Whole Model Fusion Execution

Because the feature size of intermediate layers is as large as the input of
the SR model, the required on-chip buffer size will grow up quickly as the
input size is increased. This problem cannot be solved well by commonly used
tile-based processing [8] due to large boundary data storage or computation
between tiles. Thus, this paper adopts the block convolution [10] by splitting
the input into nonoverlapped tiles and processing them with suitable padding
without boundary information. However, the previous approach in [10] limits
this to only certain layers and needs a large tile size due to information
loss at the boundary. In contrast, the proposed network with pixel level
attention could reduce this impact to minimum. We apply this to all layers
with small tile size ($40,48$) , and enable the whole model fused execution by
processing one tile from input to output directly. All intermediate data can
be stored in small on-chip buffers without external DRAM access as other SR
accelerators. With this, the external DRAM access can be reduced to the model
input and output only.

##  IV Experimental Results

###  IV-A Model Simulation Result

Fig.7 shows the simulation results and comparisons. The results with bicubic
interpolation will lose a lot of detail and become blur, especially at the
stripe boundary. In contrast, the proposed network can reconstruct images
correctly and the generated results are more clear than the traditional ones.
The results are almost kept the same even when we quantize the weights and
activations and include the full model block convolution without handling
boundary information.

Table I shows performance comparisons with other works. In which, VDSR and PAN
have better performance than ours, but they also need 10 and 25 times of model
parameters than ours. The performance difference between VDSR and ours is
marginal for several datasets. In contrast, our model size is close to the
lightweight models but has much better performance. Our model has achieved a
good balance between model size and performance, which has also been
illustrated in Fig. 1.

TABLE I: Comparison of several software approaches on commonly used datasets (Urb.100 is Urban100 and M.109 is Manga109). Algorithm | Param. | Set5 | Set14 | B100 | Urb.100 | M.109
---|---|---|---|---|---|---
Bicubic | - | 33.66 | 30.25 | 29.57 | 26.89 | 30.86
SRCNN-Ex | 57,184 | 36.66 | 32.45 | 31.36 | 29.50 | 35.60
FSRCNN | 12,464 | 37.00 | 32.63 | 31.50 | 29.88 | 36.67
This work | 25,920 | 37.38 | 32.91 | 31.69 | 30.29 | 37.00
VDSR | 665,000 | 37.53 | 33.05 | 31.90 | 30.77 | 37.22
PAN | 261,000 | 38.00 | 33.59 | 32.18 | 32.01 | 38.70
TABLE II: Comparison with other hardware works

Work | [7] | [8] | This work
---|---|---|---
Method | Improved IDN | S/L FSRCNN | Pixel Attention
Resolution | FHD | FHD | FHD
Process | 32 nm | 65 nm | 40nm
Frequency | 200 MHz | 200 MHz | 471 MHz
Parameters | 12,726 | 16,66011footnotemark: 1 | 25,920
PSNR (dB) | 36.96 | 37.06 | 37.18
Precision(W/A) | \- / 12 | 8 / 16 | 11 / 18
Memory Size | - | 572 KB | 232 KB
Gate Count | 2614K | - | 2837K
External BW | I/O + | I/O + partial | I/O only
| intermediate | intermediate |
FPS | 60 | 31.8 | 30

11footnotemark: 1

SR model only without classification model

###  IV-B Hardware Implementation

Table II shows hardware implementation results of TSMC 40nm CMOS process and
comparison. This design can achieve real-time full HD image reconstruction
with 30 frames per second when running at 471 MHz. Compared to other designs,
this design can save significant intermediate feature map I/O to DRAM with
smaller on-chip buffer size and better reconstruction image quality.

##  V Conclusion

This paper proposes an SR accelerator with hardware efficient pixel attention
model to reconstruct images with 37.38 dB PSNR on Set5 using 25.9 K
parameters, which is the first hardware work without just plain network.
Besides, we also reduce the external memory access to model I/O only with full
model block convolution and layer fusion. Such convolution and pixel attention
computation have been well supported by PE arrays with distributed weights.
The final implementation can achieve real-time full HD image reconstruction
with smaller buffer size and bandwidth but better image quality than other SR
accelerators.

## Acknowledgment

This work was supported by the Ministry of Science and Technology, Taiwan,
under Grant 109-2634-F-009-022, 109-2639-E-009-001, 110-2221-E-A49-148-MY3 and
110-2622-8-009-018-SB, and TSMC.

## References

  * [1] C. Dong, C. C. Loy, K. He, and X. Tang, “Image super-resolution using deep convolutional networks,” 2015.
  * [2] Y. Zhang, K. Li, K. Li, L. Wang, B. Zhong, and Y. Fu, “Image super-resolution using very deep residual channel attention networks,” 2018.
  * [3] N. Ahn, B. Kang, and K.-A. Sohn, “Fast, accurate, and lightweight super-resolution with cascading residual network,” 2018.
  * [4] T. Dai, J. Cai, Y. Zhang, S.-T. Xia, and L. Zhang, “Second-order attention network for single image super-resolution,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 11 065–11 074.
  * [5] H. Zhao, X. Kong, J. He, Y. Qiao, and C. Dong, “Efficient image super-resolution using pixel attention,” 2020.
  * [6] S. Lee, S. Joo, H. K. Ahn, and S.-O. Jung, “CNN acceleration with hardware-efficient dataflow for super-resolution,” _IEEE Access_ , vol. 8, pp. 187 754–187 765, 2020.
  * [7] P.-W. Yen, Y.-S. Lin, C.-Y. Chang, and S.-Y. Chien, “Real-time super resolution CNN accelerator with constant kernel size winograd convolution,” _2020 2nd IEEE International Conference on Artificial Intelligence Circuits and Systems (AICAS)_ , pp. 193–197, 2020.
  * [8] J. Lee, J. Lee, and H.-J. Yoo, “SRNPU: An energy-efficient CNN-based super-resolution processor with tile-based selective super-resolution in mobile devices,” _IEEE Journal on Emerging and Selected Topics in Circuits and Systems_ , vol. 10, no. 3, pp. 320–334, 2020.
  * [9] C. Dong, C. C. Loy, and X. Tang, “Accelerating the super-resolution convolutional neural network,” 2016.
  * [10] G. Li, Z. Liu, F. Li, and J. Cheng, “Block convolution: Towards memory-efficient inference of large-scale CNNs on FPGA,” 2021.