# QCD axion dark matter and the cosmic dipole problem

Chengcheng Han<EMAIL_ADDRESS>School of Physics, Sun Yat-Sen
University, Guangzhou 510275, China  Asia Pacific Center for Theoretical
Physics, Pohang 37673, Korea

###### Abstract

There is increasing evidence suggesting a discrepancy between the cosmic
dipole observed in the number count of distant galaxies and the one derived
from the cosmic microwave background (CMB). In this study, we investigate the
possibility that the cosmic dipole problem can be addressed by considering the
QCD axion, a hypothetical particle that arises from the spontaneous breaking
of the Peccei-Quinn symmetry and is postulated to constitute the dark matter
in our universe.

##  I Introduction

Dark matter remains one of the enigmas of our universe, constituting the
dominant component of observed matter and playing a pivotal role in the
formation of large-scale structures. Despite its fundamental importance, the
true nature of dark matter remains elusive, except for its gravitational
interactions with ordinary matter. Numerous experiments have been dedicated to
unraveling the particle nature of dark matter; however, thus far, no
compelling evidence has emerged. Given that astronomical observations have
been instrumental in providing evidence for the existence of dark matter, it
is reasonable to expect that further insights from astronomy may shed light on
its properties.

Interestingly, recent findings have revealed a significant discrepancy between
the cosmic dipole measured from galaxy number-counts and the dipole derived
from the cosmic microwave background (CMB) Wright  _et al._ (2010). The CMB
dipole arises from the motion of the solar system relative to the CMB
background and indicates a velocity of the solar system of approximately 370
km/s relative to the CMB Aghanim  _et al._ (2020). Consequently, this velocity
should induce a dipole in the measurements of galaxy number counts.

Comparing the cosmic dipole obtained from distant galaxy number counts with
the dipole inferred from the CMB has long been employed as a means to test the
cosmological principle Ellis and Baldwin (1984); Singal (2011, 2019); Tiwari
and Nusser (2016); Rubart and Schwarz (2013); Bengaly  _et al._ (2018);
Siewert  _et al._ (2021); Gibelyou and Huterer (2012). Studies conducted by
Ellis and Baldwin Ellis and Baldwin (1984) have demonstrated that for radio
sources with identical flux density spectra described by
$S\propto\nu^{-\alpha}$ (where $\nu$ denotes the frequency and $\alpha$
represents the spectral index) and an integral source count above a flux
density threshold $S_{*}$ characterized by $dN(>S_{*})/d\Omega\propto
S_{*}^{-x}$, a dipole in the galaxy number counts arises as follows:

$\frac{\Delta N}{N}=\mathbf{d_{k}}\cdot\hat{n}~{},$ (1)

where $\mathbf{d_{k}}=[2+x(1+\alpha)]\mathbf{\beta}$, and
$\mathbf{\beta}=\mathbf{v}/v_{c}$ represents the relative velocity of the
solar system with respect to the “matter rest frame”. By employing the
velocity $\mathbf{v}$ inferred from CMB measurements in Eq. (1) to estimate
the dipole in galaxy number counts, one can compare it with the dipole
obtained from actual measurements based on galaxies, thereby assessing the
consistency of the cosmological principle. It is worth noting that for a
comprehensive overview of previous dipole measurements, references such as
Peebles (2022); Aluri  _et al._ (2022) offer valuable resources.

The analysis of the CatWISE catalogue, which examines over a million quasars,
has recently unveiled a discrepancy between the amplitude of the number-count
dipole and the dipole derived from the CMB. This deviation has been reported
to be approximately 4.9$\sigma$ in references Secrest  _et al._ (2021, 2022).
Such findings raise doubts regarding the validity of the cosmological
principle, a fundamental assumption stating that the universe should exhibit
homogeneity and isotropy on sufficiently large scales. Consequently,
alternative principles have been contemplated as potential replacements for
this assumption Krishnan  _et al._ (2022). However, before entirely discarding
this fundamental principle, it is crucial to investigate whether the cosmic
dipole problem can be addressed within the framework of cosmic perturbation
theory. Interestingly, recent discoveries indicate that the inconsistency
between the two dipoles could be alleviated if a significant dark matter
isocurvature perturbation exists at super-horizon scales Domènech  _et al._
(2022).

In the case of isocurvature perturbations originating from the density
fluctuations of dark matter, they cannot be attributed to WIMP (Weakly
Interacting Massive Particle) dark matter. This is due to the fact that WIMPs
would have thermally equilibrated with photons in the early universe, leading
to the eradication of any initial isocurvature perturbations Weinberg (2004,
2003). However, axion dark matter presents a viable candidate in this context.
Axions are hypothetical particles arising from the spontaneous breaking of the
Peccei-Quinn symmetry, proposed as a solution to the strong CP problem Peccei
and Quinn (1977a, b). It has been long known that axion dark matter can
generate isocurvature perturbations with an amplitude on the order of
$\frac{H}{\pi f_{a}\theta_{0}}$, where $H$ represents the Hubble parameter
during inflation, $f_{a}$ is the axion decay constant, and $\theta_{0}$
denotes the initial misalignment angle of the axion field Weinberg (1978);
Wilczek (1978); Dine and Fischler (1983); Abbott and Sikivie (1983); Preskill
_et al._ (1983); Kim and Carosi (2010); Wantz and Shellard (2010); Marsh
(2016); Di Luzio  _et al._ (2020); Sakharov and Khlopov (1994); Sakharov  _et
al._ (1996); Khlopov  _et al._ (1999). However, it is important to note that
the non-observation of isocurvature perturbations in the CMB places stringent
constraints on the ratio of the Hubble parameter and the axion decay constant,
yielding $H/f_{a}\lesssim 10^{-5}$ for $\theta_{0}\sim\mathcal{O}(1)$, as
reported by the Planck data Akrami  _et al._ (2020). Therefore, if axions
constitute the dark matter and are responsible for the isocurvature
perturbations, their properties must satisfy this constraint.

On the other hand, to address the cosmic dipole problem, a significant
isocurvature perturbation is required, which contradicts the constraints
imposed by CMB observations. However, if we can generate a large isocurvature
perturbation at super-horizon scales that becomes sufficiently small at the
scale of recombination $k_{\rm dec}$, it may be possible to reconcile this
contradiction. In the following, we will present an axion model that satisfies
this requirement.

This paper is organized as follows. In Section II, we will provide a brief
overview of the cosmic dipole problem and present the necessary conditions to
address it. In Section III, we will introduce an axion model that has the
potential to explain the cosmic dipole problem. The numerical calculations
will be presented in Section IV. Finally, in Section V, we will draw our
conclusions based on the findings.

##  II The cosmic dipole problem

Firstly, let us present some statistics on the calculation of the power
spectrum of the cosmic microwave background (CMB). The temperature
anisotropies in the CMB can be decomposed into spherical harmonics as follows:

$\frac{\Delta
T}{T}(\hat{\mathbf{n}})=\sum_{l,m}a_{lm}Y_{lm}(\hat{\mathbf{n}})~{},$ (2)

where $a_{lm}$ are the spherical harmonic coefficients given by the integral:

$a_{lm}=\int d\Omega\frac{\Delta
T}{T}(\hat{\mathbf{n}})Y_{lm}(\hat{\mathbf{n}})~{}.$ (3)

The two-point correlation function for the temperature fluctuations is
expressed as:

$\left\langle\frac{\Delta T}{T}(\hat{\mathbf{n}})\frac{\Delta
T}{T}(\hat{\mathbf{n^{\prime}}})\right\rangle=\frac{1}{4\pi}\sum_{l=1}^{\infty}(2l+1)C_{l}P_{l}(\hat{\mathbf{n}}\cdot\hat{\mathbf{n^{\prime}}})~{},$
(4)

where we ignore the monopole term ($l=0$). In a homogeneous and isotropic
universe, we have,

$\left\langle
a_{lm}^{*}a_{lm^{\prime}}\right\rangle=\delta_{ll^{\prime}}\delta_{mm^{\prime}}C_{l}~{}.$
(5)

Typically, the dipole corresponds to $l=1$ and the quadrupole to $l=2$. The
measured CMB dipole can be written as:

$\frac{\Delta
T}{T}(\hat{\mathbf{n}})=d^{\text{CMB}}\hat{\mathbf{n}}\cdot\hat{\mathbf{n}}_{0}~{},$
(6)

where $d^{\text{CMB}}=(1.23357\pm 0.00036)\times 10^{-3}$ is the magnitude of
the dipole and $\hat{\mathbf{n}}_{0}$ specifies its direction at
$l=264^{\circ}$ and $b=48^{\circ}$. It is commonly assumed that this dipole
arises from the Doppler effect due to the motion of the solar system with
respect to the rest frame of the CMB. The relative velocity is given by:

$v_{o}=d^{\text{CMB}}\times v_{c}=369.82\pm 0.11,\text{km/s}~{},$ (7)

where $v_{c}$ is the speed of light. For the CMB dipole, we have the relation:

$\left\langle\frac{\Delta T}{T}(\hat{\mathbf{n}})\frac{\Delta
T}{T}(\hat{\mathbf{n^{\prime}}})\right\rangle=\frac{3}{4\pi}C_{1}\hat{\mathbf{n}}\cdot\hat{\mathbf{n^{\prime}}}.$
(8)

Using

$\left\langle(\hat{\mathbf{n}}\cdot\hat{\mathbf{n}}_{0})(\hat{\mathbf{n^{\prime}}}\cdot\hat{\mathbf{n}}_{0})\right\rangle=\frac{1}{3}\hat{\mathbf{n}}\cdot\hat{\mathbf{n^{\prime}}}~{},$
(9)

we get

$d^{\text{CMB}}=\sqrt{\frac{9C_{1}}{4\pi}}~{}.$ (10)

These statistics provide the necessary background for understanding the CMB
dipole and its relation to the power spectrum.

In recent studies, the measurement of kinematic dipoles from distant radio
galaxies or quasars has revealed that these dipoles are 2-3 times larger than
the expected value inferred from the velocity $v_{o}$ obtained from the CMB
dipole Secrest  _et al._ (2021, 2022). This result suggests the existence of a
preferred velocity given by:

$v_{o}^{\prime}=797\pm 87,\text{km/s}~{},$ (11)

which points in a similar direction. Notably, this velocity deviates from the
velocity $v_{o}$ inferred from the CMB dipole by approximately $4.9\sigma$. At
first glance, one might be inclined to conclude that the cosmological
principle should be discarded. However, before abandoning this fundamental
assumption, it is crucial to investigate whether this inconsistency can be
explained within the framework of the perturbed Friedmann-Lemaître-Robertson-
Walker(FLRW) model. Several studies have demonstrated that if a large
isocurvature perturbation exists at super-horizon scales, a portion of the CMB
dipole could be intrinsic Turner (1991); Langlois and Piran (1996); Erickcek
_et al._ (2008). Therefore, a possible solution to reconcile the dipole
discrepancy is to assume that our velocity is $v_{o}^{\prime}$, and an
intrinsic dipole cancels out a part of the kinematic dipole from the CMB,
i.e.,

$d^{\text{CMB}}=d_{\text{kin}}^{\text{CMB}}+D_{1}^{\text{CMB}}=1.23357\times
10^{-3}~{},$ (12)

where $d_{\text{kin}}^{\text{CMB}}=v_{o}^{\prime}/c$ is the kinematic dipole
and the $z$-direction is chosen as the dipole direction. $D_{1}^{\text{CMB}}$
represents the intrinsic dipole of the CMB arising from isocurvature
perturbations in the opposite direction. We can approximate it as:

$D_{1}^{\text{CMB}}\approx-1.4\times
10^{-3}-({v_{o}^{\prime}-797\text{km/s}})/{v_{c}}~{}.$ (13)

To account for the discrepancy within $2\sigma$, it would require
$|D_{1}^{\text{CMB}}|>8\times 10^{-4}$. This scenario of cancellation
involving an intrinsic dipole provides a potential explanation for the
observed discrepancy between the CMB and kinematic dipoles.

However, before drawing this conclusion, it is necessary to investigate
whether super-horizon isocurvature perturbations also impact the galaxy
number-count dipole. In a recent study by Domènech  _et al._ (2022), the
authors examined this issue and found that isocurvature perturbations at the
super-horizon scale primarily affect the CMB dipole and have only a minor
influence on the galaxy number-count dipole. Consequently, the authors
interpret the cosmic dipole problem as originating from a single-mode
isocurvature perturbation at the super-horizon scale.

However, it is important to note that in realistic models, the power spectrum
of perturbations typically exhibits a continuum rather than being dominated by
a single mode. In this work, we propose an alternative explanation by
attributing the cosmic dipole problem to the continuum isocurvature
perturbations associated with axion dark matter. It is worth mentioning that
while some articles attempt to employ axion strings to explain the CMB
hemispherical anomaly, this is a distinct problem from the cosmic dipole
problem we are addressing Yang  _et al._ (2017, 2019).

According to the contributions of the isocurvature perturbation, neglecting
the baryonic matter contribution, we have Domènech  _et al._ (2022); Erickcek
_et al._ (2008):

$\displaystyle\left.\frac{\Delta
T}{T}\right|_{\mathbf{k}}=-\frac{1}{3}S^{ini}_{\mathbf{k}}~{},$ (14)

where $S^{ini}_{\mathbf{k}}$ represents the initial isocurvature perturbation
long before recombination. This isocurvature perturbation also sources the
adiabatic perturbation:

$\displaystyle\Phi_{\mathbf{k}}=-\frac{1}{5}S^{ini}_{\mathbf{k}}~{}.$ (15)

It is well-known that the adiabatic perturbation generates temperature
fluctuations through the Sachs-Wolfe effect:

$\displaystyle\left.\frac{\Delta
T}{T}\right|^{\text{adia}}_{\mathbf{k}}=\frac{1}{3}\Phi_{\mathbf{k}}.$ (16)

Summing these contributions, we obtain,

$\displaystyle\left.\frac{\Delta
T}{T}\right|_{\mathbf{k}}=-\frac{2}{5}S^{ini}_{\mathbf{k}}~{}.$ (17)

However, the super-horizon adiabatic perturbation also induces a velocity of
us, which contributes to the dipole. Remarkably, its effect precisely cancels
the dipole contribution from the adiabatic perturbation Turner (1991);
Langlois and Piran (1996); Erickcek  _et al._ (2008). This explains why only
the super-horizon adiabatic perturbation alone cannot generate a CMB dipole.

The power spectrum of the isocurvature perturbation, denoted as
$\mathcal{P}_{S}(k)$, can be defined as follows:

$\displaystyle\langle
S_{\mathbf{k}}S_{\mathbf{k}^{\prime}}\rangle=2\pi^{2}\frac{\mathcal{P}_{S}(k)}{k^{3}}\delta(\mathbf{k-k^{\prime}})\Theta(k-k_{\text{min}})~{},$
(18)

where $k_{\text{min}}$ is the minimum comoving wave number at which the
isocurvature perturbation becomes relevant. Using this power spectrum, we can
express the angular power spectrum $C_{l}$ as,

$\displaystyle
C_{l}=g_{l}\int_{k_{\text{min}}}^{\infty}\frac{dk}{k}\mathcal{P}_{S}(k)j_{l}^{2}(kr{\text{dec}})~{},$
(19)

where $g_{l}=\frac{4\pi}{9}$ for $l=1$ and $g_{l}=\frac{16\pi}{25}$ for $l>1$.
Here, $k$ is the comoving wave number and $r_{\text{dec}}\approx 14.1$ Gpc is
the comoving distance to the photon last scattering surface. To explain the
cosmic dipole problem, we require

$\displaystyle\sqrt{\frac{9C_{1}}{4\pi}}\gtrsim 8\times 10^{-4}\Rightarrow
C_{1}\gtrsim 9\times 10^{-7}~{}.$ (20)

Constraints from multipole observations, particularly the quadrupole, provide
additional constraints on the power spectrum of the isocurvature perturbation.
According to the Planck 2018 results Aghanim  _et al._ (2020), the temperature
variance for the quadrupole is measured to be $(200^{+540}_{-120})\mu K^{2}$.
Taking the 3$\sigma$ upper limit, we have,

$\displaystyle\frac{2(2+1)C_{2}}{2\pi}\lesssim 2.5\times 10^{-10}\Rightarrow
C_{2}\lesssim 2.6\times 10^{-10}.$ (21)

Thus we have

$\displaystyle\frac{C_{2}}{C_{1}}\lesssim 3\times 10^{-4}~{}.$ (22)

It’s important to note that this condition is derived from the quadrupole
limit and the constraints from other multipoles are much weaker.

Before presenting a specific model, let’s assume that the power spectrum of
the isocurvature perturbation, $\mathcal{P}_{S}(k)$, follows a power law given
by $\mathcal{P}_{S}(k)=A(k/k{\rm min})^{n-1}$, where $n<1$. For $n>-1$, we can
approximate the expression for the multipole moments $C_{l}$ as:

$\displaystyle C_{l}\approx g_{l}A(k_{\rm min}r_{\rm dec})^{1-n}c(n,l)~{},$
(23)

where $g_{l}$ and $c(n,l)$ are defined as before. However, this expression is
only valid when $n-2+2l>-1$. If $n-2+2l<-1$, the main contribution to the
integration comes from the region around $k_{\rm min}$. By using the
properties of $j_{l}(x)\sim\frac{x^{l}}{(2l+1)!!}$ for $x\ll 1$, it can be
shown that:

$\displaystyle C_{l}\approx\frac{g_{l}A}{|n-1+2l|((2l+1)!!)^{2}}(k_{\rm
min}r_{\rm dec})^{2l}~{}.$ (24)

This expression is valid when $n-2+2l<-1$.

Taking the example where $n=-2$ and $k_{\rm min}r_{\rm dec}=0.01$, we find
that

$\displaystyle\frac{C_{2}}{C_{1}}\approx 1.3\times 10^{-3}~{}.$ (25)

This value is too large to satisfy the limit on the quadrupole. To address
this, we can either decrease $k_{\rm min}r_{\rm dec}$ or make $n$ smaller. For
instance, if we take $k_{\rm min}r_{\rm dec}=0.0023$, we obtain

$\displaystyle\frac{C_{2}}{C_{1}}=3\times 10^{-4}~{}.$ (26)

Assuming $A=1$, we have $C_{1}=8.4\times 10^{-7}$, which marginally explains
the cosmic dipole discrepancy.

##  III Large isocurvature from the axion

In this section, we will demonstrate that a significant isocurvature
perturbation can be generated by the axion dark matter. The energy density of
the axion can be approximated as

$\displaystyle\rho_{a}=\frac{1}{2}m_{a}^{2}f^{2}\theta_{0}^{2}~{},$ (27)

where $f$ is the decay constant of the axion and $\theta_{0}$ is the initial
displacement angle. In the case of the QCD axion, $f$ can typically range from
$10^{11}$ GeV to $10^{15}$ GeV for a $\theta_{0}$ varying from
$\mathcal{O}(10^{-2})$ to $\mathcal{O}(1)$. The axion mass can be estimated
using the following formula,

$\displaystyle m_{a}=5.7\left(\frac{10^{9}~{}\rm GeV}{f}\right)\rm meV~{}.$
(28)

For an axion-like particle (ALP), the relation between the mass and $f$ can be
more relaxed.

The isocurvature perturbation of the axion can be defined as

$\displaystyle
S\equiv\frac{\delta\rho_{a}}{\rho_{a}}-\frac{3}{4}\frac{\delta\rho_{r}}{\rho_{r}}\simeq\frac{\delta\rho_{a}}{\rho_{a}}=\frac{2\delta\theta}{\theta_{0}},$
(29)

where $\delta\rho_{a}$ and $\delta\rho_{r}$ are the perturbations in the axion
and radiation energy densities, respectively. The amplitude of the axion
perturbation $\delta\theta$ can be related to the Hubble parameter $H$ and the
axion field value at a given comoving scale as

$\displaystyle\mathcal{P}^{1/2}_{\delta\theta}=\frac{H}{2\pi\varphi(k)}~{}.$
(30)

Using this relation, we can express the amplitude of the isocurvature
perturbation as

$\displaystyle\mathcal{P}^{1/2}_{S}=\frac{H}{\pi\varphi(k)\theta_{0}}~{}.$
(31)

It is important to note that if the perturbation $\delta\theta$ is comparable
to the background value $\theta_{0}$, we should use the effective angle
$\theta_{\rm eff}=\sqrt{\delta\theta^{2}+\theta_{0}^{2}}$ to calculate the
power spectrum. As a result, we always have $\mathcal{P}^{1/2}_{S}\lesssim 1$.
From the equation above, it is evident that a large isocurvature requires
$\varphi\sim H$ at the super horizon scale, while $\varphi$ should be
sufficiently large at $k_{\rm dec}$ to avoid the isocurvature limit imposed by
CMB observations. If the value of $\varphi$ rapidly increases to a very large
value during the early stage of inflation, the tension with the isocurvature
limit from CMB observations can be alleviated, as illustrated in Figure 1. In
the following, we will present a model that realizes this mechanism.

The model consists of a complex scalar field $\Phi$ with a non-vanishing
Peccei-Quinn charge and a potential

$\displaystyle V(\Phi)=\lambda(\Phi\Phi^{\dagger}-f^{2}/2)^{2}~{}.$ (32)

We define the axion field as the angular mode of $\Phi$,

$\displaystyle\Phi=\frac{1}{\sqrt{2}}\varphi\exp(i\frac{a}{f})~{}.$ (33)

Assuming that initially $\varphi$ is located close to the origin but with a
displacement around $H$. Ignoring the motion of angular direction, the
equation of motion for $\varphi$ is then given by,

$\displaystyle\ddot{\varphi}+3H\dot{\varphi}+V^{\prime}(\varphi)=0~{}.$ (34)

Figure 1: Illustration for a varying isocurvature perturbation during the
early universe. $\varphi$ evolves from a small value to a large value during
inflation.

The assumption that the field $\Phi$ is initially located around the origin of
the potential is motivated by the presence of a coupling between the inflaton
field $I$ and $\Phi$. The coupling term, $V_{I}=g{\Phi^{\dagger}\Phi I^{2}}$,
provides an effective mass for the $\Phi$ field during inflation due to the
large value of the inflaton $I$. For a positive coupling constant $g$, this
effective mass is larger than the Hubble parameter, leading to a minimum of
the potential at $\Phi=0$ and the preservation of the Peccei-Quinn symmetry.

During inflation, as the inflaton $I$ rolls from a large value to a smaller
value, the effective mass term of $\Phi$ from $V_{I}$ becomes sub-dominant
compared to the negative mass term from the potential in Eq. (32). At a
certain stage, when $I$ becomes sufficiently small, the negative mass from the
potential dominates, causing the potential to develop a non-vanishing vacuum
expectation value. Since the potential is flat at the top, the field $\Phi$
initially undergoes random walking and is then displaced from the origin with
an amplitude around $H$. Subsequently, $\Phi$ evolves according to the
equation of motion described by Eq. (34), eventually reaching the value of $f$
after a few e-folds of inflation. If this occurs at the super horizon scale
$k_{\rm min}$, the isocurvature perturbation is large at the super horizon
scale and small at $k_{\rm dec}$.

##  IV Numerical calculation

For the numerical calculation, since

$\displaystyle\mathcal{P}^{1/2}_{S}=\frac{H}{\pi\varphi(k)\theta_{0}}~{},$
(35)

thus

$\displaystyle\mathcal{P}_{S}(k)\approx\mathcal{P}_{S}(k_{\rm
min})\left(\frac{\varphi(k)}{\varphi(k_{\rm min})}\right)^{-2}~{},$ (36)

where $\mathcal{P}_{S}(k_{\rm min})\simeq 1$. The evolution of $\varphi(k)$
follows the equation of motion of Eq. (34). Given that $k=k_{\rm
min}\exp(Ht)$, we have

$\displaystyle C_{l}=g_{l}\int_{k_{\rm
min}}^{\infty}\frac{dk}{k}\mathcal{P}_{S}(k)j_{l}^{2}(kr_{\rm dec})~{}.$ (37)

In practically we can define $\tau=Ht$, now above equation becomes,

$\displaystyle C_{l}=g_{l}\int_{0}^{\infty}{d\tau}\mathcal{P}_{S}(k_{\rm
min}e^{\tau})j_{l}^{2}(k_{\rm min}e^{\tau}r_{\rm dec})~{}.$ (38)

Figure 2: The motion for the radial mode during inflation. The blue, red and
green curves correspond to $\lambda=10^{-9},2\times 10^{-9},4\times 10^{-9}$.

In Figure 2, we present the evolution of $\varphi$ for different values of the
parameter $\lambda$. The blue, red, and green curves correspond to
$\lambda=10^{-9}$, $2\times 10^{-9}$, and $4\times 10^{-9}$, respectively. The
initial conditions are chosen as $\dot{\varphi}=0$, $\varphi(k_{\rm
min})=H/\pi$, and $f=10^{5}H$. As observed in the figure, for smaller values
of $\lambda$, the field $\varphi$ remains at a small value for a longer
duration. However, after a few e-folds, it eventually reaches the value of
$f$. It is important to note that if $\lambda$ is too small, the field may
remain at the small value for an extended period, resulting in a large
isocurvature perturbation at recombination. Conversely, if $\lambda$ is too
large and combined with a small $k_{\rm min}$, the field quickly reaches the
value of $f$, and it becomes challenging to explain the cosmic dipole problem.

In order to address the cosmic dipole problem, we require $C_{1}>9\times
10^{-7}$. However, the Planck data imposes a limit on the dark matter
isocurvature perturbation, $\mathcal{P}_{S}\lesssim 0.8\times 10^{-10}$ at
$k_{\rm low}=0.002~{}{\rm Mpc}^{-1}$. To evade this limit, it is necessary to
have $\varphi(k_{\rm low})\gtrsim f/3$ for $f=10^{5}H$.

In Figure 3, we illustrate the parameter space required to explain the cosmic
dipole anomaly in the $1/(k_{\rm min}r_{\rm dec})$ versus $\lambda$ plane. The
two solid curves represent $C_{1}=9\times 10^{-7}$ and $C_{1}=5.6\times
10^{-6}$, which correspond to the lower and upper limits of the intrinsic CMB
dipole, respectively. The gray region is excluded due to the presence of a
too-large quadrupole in the CMB, while the red region is excluded due to a
significant isocurvature contribution at $k_{\rm low}=0.002~{}{\rm Mpc}^{-1}$.

Figure 3: Parameter space to explain the comic dipole problem. The two real
lines correspond $C_{1}=9\times 10^{-7}$ and $C_{1}=5.6\times 10^{-6}$. The
gray region is excluded because of the too large quadrupole at CMB. The red
region is excluded due to too large isocurvature.

From the plot, we observe that for $\lambda<7\times 10^{-10}$, regardless of
the value of $1/(k_{\rm min}r_{\rm dec})$, the predicted quadrupole in the CMB
would be too large. On the other hand, for $\lambda>4\times 10^{-9}$, values
of $1/(k_{\rm min}r_{\rm dec})<250$ are necessary to generate a sufficient
contribution to the cosmic dipole. It is worth noting that as $\lambda$
increases, the curve $C_{1}=5.6\times 10^{-6}$ eventually reaches the gray
region ($C_{2}\lesssim 2.6\times 10^{-10}$), which sets an upper limit on
$\lambda$. Numerically, we find this occurs at $\lambda=6.6\times 10^{-5}$ and
$1/(k_{\rm min}r_{\rm dec})=14.6$. As a benchmark point in the model to
explain the dipole problem, we consider $\lambda=2\times 10^{-9}$, $1/(k_{\rm
min}r_{\rm dec})=200$, and $\mathcal{P}_{S}(k_{\rm min})\approx 1$. This
yields $C_{1}=1.6\times 10^{-6}$ and $C_{2}=6.2\times 10^{-12}$.

Including the baryon matter contribution to the matter component introduces an
additional factor
$\left(\frac{\Omega_{cd}}{\Omega_{cd}+\Omega_{b}}\right)^{2}$ in the
calculation of $C_{l}$. While the inclusion of this factor can affect the
precise numerical values of the predicted CMB multipoles, the essential result
and the overall conclusion of our analysis should remain unchanged.

##  V Conclusion

In this work, we propose an axion model as a potential solution to the cosmic
dipole problem. Our model predicts a large isocurvature perturbation at the
super horizon scale while maintaining a small perturbation at recombination to
comply with the isocurvature limit derived from CMB observations. This is
achieved by allowing the radial mode of the axion field to evolve from a small
value to a significantly larger value during the inflationary period. Through
our analysis, we find that by choosing appropriate parameters such as
$f=10^{5}H$ and $7\times 10^{-10}<\lambda<6.6\times 10^{-5}$, along with a
suitable value of $k_{\rm min}$, it is possible to explain the cosmic dipole
problem within the framework of our axion model.

It is worth noting that the cosmic dipole problem, if successfully addressed
by our model, could serve as the first evidence supporting the existence of
axion dark matter. This provides additional motivation for further exploration
and investigation of axion physics in the context of cosmology and particle
physics. Overall, our proposed axion model offers a promising avenue to
resolve the cosmic dipole problem and potentially sheds light on the nature of
dark matter through the presence of axions.

Acknowledgment.—CH acknowledges Tsutomu T. Yanagida, Misao Sasaki, Kazunori
Kohri and Shi Pi for helpful discussions. CH acknowledges support from the Sun
Yat-Sen University Science Foundation and the Fundamental Research Funds for
the Central Universities, Sun Yat-sen University under Grant No. 23qnpy58.
This research is supported by the Associate Fellow Program from the APCTP.

##  Appendix A Galaxy number-count dipole from super-horizon perturbation

In this appendix we give a rough proof that the super-horizon perturbation
does not give sizable galaxy number-count dipole. In the estimation of the
galaxy number-count dipole, we adopt the results from previous studies Bonvin
and Durrer (2011); Di Dio  _et al._ (2013); Nadolny  _et al._ (2021); Maartens
_et al._ (2018), which are also consistent with the findings of the authors in
Domènech  _et al._ (2022). The expression for the dipole is given by

$\displaystyle\Delta({\bf
n},z,m_{*})=-\left(2+\frac{\mathcal{H}^{\prime}}{\mathcal{H}^{2}}+\frac{2-5s}{r_{S}\mathcal{H}}-f_{\rm
evo}\right){\bf V_{0}\cdot n}$ $\displaystyle+b\delta+(f_{\rm
evo}-3)\mathcal{H}V/k-3\Phi$
$\displaystyle+(1+5s)\Phi+\Psi+\frac{1}{\mathcal{H}}\left[\Phi^{\prime}+\partial_{r}({\bf
V\cdot n})\right]$
$\displaystyle+\left(\frac{\mathcal{H}^{\prime}}{\mathcal{H}^{2}}+\frac{2-5s}{r_{S}\mathcal{H}}+5s-f_{\rm
evo}\right)\left(\Psi+{\bf V\cdot n}\right.$
$\displaystyle\left.+\int_{0}^{r_{S}}dr(\Phi^{\prime}+\Psi^{\prime})\right)$
$\displaystyle+\frac{2-5s}{2r_{S}}\int_{0}^{r_{S}}dr\left[2-\frac{r_{S}-r}{r}\Delta_{\Omega}\right](\Phi+\Psi)~{}.$
(39)

In the equation above, ${\bf V_{0}\cdot n}$ term represents the velocity
induced by super-horizon fluctuations. This term only contributes to the
dipole component and is typically neglected in studies such as Bonvin and
Durrer (2011); Nadolny  _et al._ (2021). The quantities $f_{\rm evo}$ and $s$
are defined as follows:

$\displaystyle f_{\rm evo}$ $\displaystyle=$
$\displaystyle\frac{\partial\ln\left(a^{3}N(z,L>L_{*})\right)}{\mathcal{H}\partial\tau_{S}}~{},$
(40) $\displaystyle s$ $\displaystyle=$
$\displaystyle-\frac{2}{5}\frac{\partial\ln N(z,L_{S}>L_{*})}{\partial
L_{S}}~{},$ (41)

where $N(z,L>L_{*})$ represents the number of galaxies with luminosity greater
than $L_{*}$ at redshift $z$, and $L_{S}$ is the luminosity associated with
the comoving scale $r_{S}$. In addition, $\delta$ is the matter density
contrast in the co-moving frame, and $b$ is the linear bias parameter that
connects the galaxy number-count fluctuations to the dark matter fluctuation.
In this case, we assume $b=1$ for simplicity. $V,\Phi$ and $\Psi$ are the
gauge invariant velocity and Bardeen potentials.

The gauge-invariant velocity $V$, $\Phi$ and $\Psi$ and $\delta$ are related
by the following expressions:

$\displaystyle\Psi$ $\displaystyle=$
$\displaystyle-\frac{3}{4}\left(\frac{k_{\rm
eq}}{k}\right)^{2}\frac{1+y}{y^{2}}\delta~{},$ (42) $\displaystyle V$
$\displaystyle=$
$\displaystyle-\frac{3}{k}\mathcal{H}\frac{1+y}{4+3y}\left[y\frac{d\delta}{dy}-\frac{1}{1+y}\delta\right]~{},$
(43)

where $y=a/a_{\rm eq}$ is a parameter related to the scale factor and $k_{\rm
eq}$ represents the scale of equality between matter and radiation. It is
worth noting that for super-horizon perturbations ($k\ll H_{0}$), the
anisotropy of the stress can be ignored, leading to the relation $\Phi=\Psi$.

For adiabatic perturbation, assuming initial primordial power spectrum
$\mathcal{P}(k)$ defined as

$\displaystyle\langle\Psi^{ini}_{\bf k}\Psi^{ini}_{\bf
k^{\prime}}\rangle=2\pi^{2}\frac{\mathcal{P}(k)}{k^{3}}\delta({\bf
k-k^{\prime}})~{}.$ (44)

Then the power spectrum $C^{N}_{l}$ can be calculated by:

$\displaystyle C^{N}_{l}(z)=4\pi\int\frac{dk}{k}\mathcal{P}(k)F^{2}_{l}~{},$
(45)

where

$\displaystyle F_{l}$ $\displaystyle=$ $\displaystyle
j_{l}(kr_{S})\left[bT_{\delta}+\left(\frac{\mathcal{H}^{\prime}}{\mathcal{H}^{2}}+\frac{2-5s}{r_{S}\mathcal{H}}+5s-f_{\rm
evo}+1\right)T_{\Psi}\right.$ (46) $\displaystyle+$
$\displaystyle\left.(-2+5s)T_{\Phi}+\mathcal{H}^{-1}T^{\prime}_{\Phi}\right]$
$\displaystyle+$
$\displaystyle\left[j^{\prime}_{l}(kr_{S})\left(\frac{\mathcal{H}^{\prime}}{\mathcal{H}^{2}}+\frac{2-5s}{r_{S}\mathcal{H}}+5s-f_{\rm
evo}\right)+j^{\prime\prime}_{l}(kr_{S})\frac{k}{\mathcal{H}}\right.$
$\displaystyle+$ $\displaystyle\left.(f_{\rm
evo}-3)j_{l}(kr_{S})\frac{\mathcal{H}}{k}\right]T_{V}(r_{S})$ $\displaystyle-$
$\displaystyle
j^{\prime}_{l}(0)\left(2+\frac{\mathcal{H}^{\prime}}{\mathcal{H}^{2}}+\frac{2-5s}{r_{S}\mathcal{H}}-f_{\rm
evo}\right)T_{V}(0)$ $\displaystyle+$
$\displaystyle\frac{2-5s}{2r_{S}}\int_{0}^{r_{S}}\left[(T_{\Phi}+T_{\Psi})\left(l(l+1)\frac{r_{S}-r}{r}+2\right)\right]$
$\displaystyle+$
$\displaystyle\int_{0}^{r_{S}}(T^{\prime}_{\Phi}+T^{\prime}_{\Psi})\left(\frac{\mathcal{H}^{\prime}}{\mathcal{H}^{2}}+\frac{2-5s}{r_{S}\mathcal{H}}+5s-f_{\rm
evo}\right)~{},$

where $j_{l}(x)$ is the spherical Bessel function, $T_{\delta}$, $T_{\Psi}$,
$T_{\Phi}$, and $T_{V}$ are transfer functions, and $r_{S}$ represents the
comoving distance of the galaxy sources. The equation above provides a way to
compute the power spectrum $C^{N}_{l}$ at a given redshift $z$ by integrating
over the primordial power spectrum $\mathcal{P}(k)$ and evaluating the
function $F_{l}$ using the transfer functions $T_{\delta}$, $T_{\Psi}$,
$T_{\Phi}$, and $T_{V}$.

The evaluation of the transfer functions needs numerical calculation. For
super-horizon modes with $k<H_{0}$ and $y\gg 1$, the transfer functions can be
analytically approximated as follows (from reference Hu and Sugiyama (1995)):

$\displaystyle T_{\Psi}$ $\displaystyle=$ $\displaystyle\frac{9}{10}~{},$
$\displaystyle T_{\Phi}$ $\displaystyle=$ $\displaystyle T_{\Psi}~{},$
$\displaystyle T_{V}$ $\displaystyle=$ $\displaystyle\frac{1}{3}k\eta
T_{\Psi}~{},$ $\displaystyle T_{\delta}$ $\displaystyle=$
$\displaystyle-\frac{4}{3}y\left(\frac{k}{k_{\rm eq}}\right)^{2}T_{\Psi}~{}.$
(47)

These expressions provide analytical approximations for the transfer functions
$T_{i}$ ($i=\delta,V,\Phi,\Psi$) during the matter-dominated stage for super-
horizon modes.

To check the result from the paper Domènech  _et al._ (2022), we take $z=1$
(corresponding to $r_{S}=$ 3.4 Gpc) as an example, since the average redshift
of quasars used in Wright  _et al._ (2010) is around 1. We fix $k$ to satisfy
$kr_{\rm dec}=0.1$ and take $l=1$. We find that

$\displaystyle F_{l}(l=1)=-10^{-5}-3\times 10^{-7}f_{\rm evo}-10^{-5}s~{}.$
(48)

There must be large cancellations between different terms, as many terms are
orders of magnitude larger than $F_{l}$. For example, the value of the
$j^{\prime}_{l}(kr_{S})$ term in the third line of Eq. (46) is,

$\displaystyle(j^{\prime}_{l}(kr_{S})~{}{\rm term})(l=1)=0.02-0.0076f_{\rm
evo}-0.02s~{}.$

As demonstrated by the authors Domènech  _et al._ (2022), for a single super-
horizon mode, the cancellation occurs at the order of $kr_{S}$, and the main
contribution is approximately at the order of $(kr_{S})^{3}$.

In the case of a super-horizon isocurvature perturbation, as described in Hu
and Sugiyama (1995), we have the following relations,

$\displaystyle\Phi_{\bf k}$ $\displaystyle=$ $\displaystyle\Psi_{\bf
k}=-\frac{1}{5}S_{\bf k}~{},$ $\displaystyle V_{\bf k}$ $\displaystyle=$
$\displaystyle-\frac{1}{15}k\eta S_{\bf k}~{},$ $\displaystyle\delta_{\bf k}$
$\displaystyle=$ $\displaystyle\frac{4}{15}y\left(\frac{k}{k_{\rm
eq}}\right)^{2}S_{\bf k}~{}.$ (50)

To calculate $C^{N}_{l}$, we need to replace $\mathcal{P}$ with
$\mathcal{P}_{S}$, and the transfer functions become,

$\displaystyle T_{\Psi}$ $\displaystyle=$ $\displaystyle-\frac{1}{5}~{},$
$\displaystyle T_{\Phi}$ $\displaystyle=$ $\displaystyle T_{\Psi}~{},$
$\displaystyle T_{V}$ $\displaystyle=$ $\displaystyle\frac{1}{3}k\eta
T_{\Psi}~{},$ $\displaystyle T_{\delta}$ $\displaystyle=$
$\displaystyle-\frac{4}{3}y\left(\frac{k}{k_{\rm eq}}\right)^{2}T_{\Psi}~{}.$
(51)

Since the relation of $\Phi,V,\delta$ with $\Psi$ does not change, the
cancellation in $F_{l}$ still exists. As shown in Domènech  _et al._ (2022),
such cancellation also holds for a single mode even if the dark energy-
dominated stage is considered. For the continuum power spectrum, a full
numerical calculation of the galaxy number-count dipole including the dark
energy-dominated stage is needed, but it is beyond the scope of this paper and
we leave it for future study. In the end, we provide the predicted galaxy
number-count dipole for our benchmark point, assuming the universe is only
matter-dominated. By considering only the super-horizon mode ($k<H_{0}$),
assuming $f_{\rm evo}=1,s=1$ we find $C_{1}^{N}=6\times 10^{-14}$,
corresponding to a dipole around $6\times 10^{-8}$. Note that the sub-horizon
mode also contributes to the dipole, but since
$\mathcal{P}_{S}<\mathcal{P}_{R}$ (which is approximately $2.5\times 10^{-9}$)
and it has been shown that the adiabatic perturbation $\mathcal{P}_{R}$
contributes only around $0.001$ to the galaxy number-count dipole at $z=1$
Nadolny  _et al._ (2021), we believe the contribution of the galaxy number-
count dipole from the sub-horizon isocurvature mode is negligible in our
model.

## References

  * Wright  _et al._ (2010) E. L. Wright _et al._ , Astron. J. 140, 1868 (2010), arXiv:1008.0031 [astro-ph.IM] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A1 (2020), arXiv:1807.06205 [astro-ph.CO] .
  * Ellis and Baldwin (1984) G. F. R. Ellis and J. E. Baldwin, Monthly Notices of the Royal Astronomical Society 206, 377 (1984), https://academic.oup.com/mnras/article-pdf/206/2/377/18187025/mnras206-0377.pdf .
  * Singal (2011) A. K. Singal, Astrophys. J. Lett. 742, L23 (2011), arXiv:1110.6260 [astro-ph.CO] .
  * Singal (2019) A. K. Singal, Phys. Rev. D 100, 063501 (2019), arXiv:1904.11362 [physics.gen-ph] .
  * Tiwari and Nusser (2016) P. Tiwari and A. Nusser, JCAP 03, 062 (2016), arXiv:1509.02532 [astro-ph.CO] .
  * Rubart and Schwarz (2013) M. Rubart and D. J. Schwarz, Astron. Astrophys. 555, A117 (2013), arXiv:1301.5559 [astro-ph.CO] .
  * Bengaly  _et al._ (2018) C. A. P. Bengaly, R. Maartens, and M. G. Santos, JCAP 04, 031 (2018), arXiv:1710.08804 [astro-ph.CO] .
  * Siewert  _et al._ (2021) T. M. Siewert, M. Schmidt-Rubart, and D. J. Schwarz, Astron. Astrophys. 653, A9 (2021), arXiv:2010.08366 [astro-ph.CO] .
  * Gibelyou and Huterer (2012) C. Gibelyou and D. Huterer, Mon. Not. Roy. Astron. Soc. 427, 1994 (2012), arXiv:1205.6476 [astro-ph.CO] .
  * Peebles (2022) P. J. E. Peebles, Annals Phys. 447, 169159 (2022), arXiv:2208.05018 [astro-ph.CO] .
  * Aluri  _et al._ (2022) P. K. Aluri _et al._ , (2022), arXiv:2207.05765 [astro-ph.CO] .
  * Secrest  _et al._ (2021) N. J. Secrest, S. von Hausegger, M. Rameez, R. Mohayaee, S. Sarkar, and J. Colin, Astrophys. J. Lett. 908, L51 (2021), arXiv:2009.14826 [astro-ph.CO] .
  * Secrest  _et al._ (2022) N. J. Secrest, S. von Hausegger, M. Rameez, R. Mohayaee, and S. Sarkar, Astrophys. J. Lett. 937, L31 (2022), arXiv:2206.05624 [astro-ph.CO] .
  * Krishnan  _et al._ (2022) C. Krishnan, R. Mondol, and M. M. Sheikh-Jabbari, (2022), arXiv:2209.14918 [astro-ph.CO] .
  * Domènech  _et al._ (2022) G. Domènech, R. Mohayaee, S. P. Patil, and S. Sarkar, JCAP 10, 019 (2022), arXiv:2207.01569 [astro-ph.CO] .
  * Weinberg (2004) S. Weinberg, Phys. Rev. D 70, 083522 (2004), arXiv:astro-ph/0405397 .
  * Weinberg (2003) S. Weinberg, Phys. Rev. D 67, 123504 (2003), arXiv:astro-ph/0302326 .
  * Peccei and Quinn (1977a) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn (1977b) R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977b).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Kim and Carosi (2010) J. E. Kim and G. Carosi, Rev. Mod. Phys. 82, 557 (2010), [Erratum: Rev.Mod.Phys. 91, 049902 (2019)], arXiv:0807.3125 [hep-ph] .
  * Wantz and Shellard (2010) O. Wantz and E. P. S. Shellard, Phys. Rev. D 82, 123508 (2010), arXiv:0910.1066 [astro-ph.CO] .
  * Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Di Luzio  _et al._ (2020) L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Sakharov and Khlopov (1994) A. S. Sakharov and M. Y. Khlopov, Phys. Atom. Nucl. 57, 485 (1994).
  * Sakharov  _et al._ (1996) A. S. Sakharov, D. D. Sokoloff, and M. Y. Khlopov, Phys. Atom. Nucl. 59, 1005 (1996).
  * Khlopov  _et al._ (1999) M. Y. Khlopov, A. S. Sakharov, and D. D. Sokoloff, Nucl. Phys. B Proc. Suppl. 72, 105 (1999).
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Turner (1991) M. S. Turner, Phys. Rev. D 44, 3737 (1991).
  * Langlois and Piran (1996) D. Langlois and T. Piran, Phys. Rev. D 53, 2908 (1996), arXiv:astro-ph/9507094 .
  * Erickcek  _et al._ (2008) A. L. Erickcek, S. M. Carroll, and M. Kamionkowski, Phys. Rev. D 78, 083012 (2008), arXiv:0808.1570 [astro-ph] .
  * Yang  _et al._ (2017) Q. Yang, Y. Liu, and H. Di, Phys. Rev. D 96, 083516 (2017), arXiv:1612.03708 [astro-ph.CO] .
  * Yang  _et al._ (2019) Q. Yang, H. Yu, and H. Di, Phys. Dark Univ. 26, 100407 (2019), arXiv:1712.09917 [astro-ph.CO] .
  * Bonvin and Durrer (2011) C. Bonvin and R. Durrer, Phys. Rev. D 84, 063505 (2011).
  * Di Dio  _et al._ (2013) E. Di Dio, F. Montanari, J. Lesgourgues, and R. Durrer, JCAP 11, 044 (2013), arXiv:1307.1459 [astro-ph.CO] .
  * Nadolny  _et al._ (2021) T. Nadolny, R. Durrer, M. Kunz, and H. Padmanabhan, JCAP 11, 009 (2021), arXiv:2106.05284 [astro-ph.CO] .
  * Maartens  _et al._ (2018) R. Maartens, C. Clarkson, and S. Chen, JCAP 01, 013 (2018), arXiv:1709.04165 [astro-ph.CO] .
  * Hu and Sugiyama (1995) W. Hu and N. Sugiyama, Phys. Rev. D 51, 2599 (1995), arXiv:astro-ph/9411008 .


# Per Clip Lagrangian Multiplier Optimisation for HEVC

Daniel J Ringis, François Pitié and Anil Kokaram; Sigmedia Group, Electronic
and Electrical Engineering Dept., Trinity College Dublin; Dublin, Ireland

###### Abstract

The majority of internet traffic is video content. This drives the demand for
video compression in order to deliver high quality video at low target
bitrates. This paper investigates the impact of adjusting the rate distortion
equation on compression performance. A constant of proportionality, k, is used
to modify the Lagrange multiplier used in H.265 (HEVC). Direct optimisation
methods are deployed to maximise BD-Rate improvement for a particular clip.
This leads to up to 21% BD-Rate improvement for an individual clip.
Furthermore we use a more realistic corpus of material provided by YouTube.
The results show that direct optimisation using BD-rate as the objective
function can lead to further gains in bitrate savings that are not available
with previous approaches.

## Introduction

Video content is predicted to rise to over 80% of internet traffic by the year
2021 [1]. This continues to drive video compression technology. A successful
video codec would provide a bitstream which is capable of being decoded at the
highest possible quality for a target bitrate, or the lowest possible bitrate
for a target quality or distortion. One of the key challenges to be solved in
the design of a practical codec is the tradeoff between rate and distortion. A
Lagrange multiplier approach was advocated by Sullivan et al [2] and that has
been the adopted view since 1998. In this approach the codec makes a number of
decisions in order to minimise a cost J as follows.

$J=D+\lambda R$ (1)

As can be seen $J$ combines both a distortion D (for a frame or macroblock)
and a rate R (the number of coded bits for that unit) through the action of
the Lagrangian multiplier $\lambda$. Different choices for $\lambda$ result in
different R/D tradeoffs. The idea is used across a wide range of codec
operations e.g. block type (Skipped/Intra/Inter), motion vectors and bit
allocation at the frame and clip levels [3].

The choice of $\lambda$ was determined experimentally in the early days of
codec research [4]. This parameter can be considered as a kind of
hyperparameter in a hybrid codec as it affects so many other decisions. This
was originally implemented for the H.263 codec and then modified for H.264 and
H.265 (HEVC) [5]. The approach was to select a value based on a cohort of
examples such that the performance was optimised on average across the set.
That set was quite small, less than 5 examples. More recently, Katsavounidis
et al [6][7] recognised that by choosing codec hyperparameters per clip
significant gains could be achieved. This is because the statistics of video
clips varies greatly over any corpus. While that work considered only the
design of a bitrate ladder for DASH streaming, it provided an impetus for per
clip optimisation throughout the codec system. Covell et al[8] also built on
that idea by optimising other codec hyperparameters (CRF) in that case for per
clip transcoding.

This paper explores the idea that a better $\lambda$ exists for an individual
video clip. Unlike previous work, $\lambda$ was adjusted using direct
optimisation techniques using BD-Rate as the objective function. In addition,
a more comprehensive and modern corpus of test data was used, provided
recently by YouTube [9]. For a given clip it was found that up to a 21%
improvement can be had by adjusting the Lagrangian multiplier in the rate
distortion equation. An example of the gains reported in this paper are shown
in Figure 1. The proposed algorithms all result in RD performance that is well
above the default behaviour in the X265 codec (in this case). This result is
discussed in more detail later, but it contrasts well with previous work that
reports up to 11% improvement. The next section presents more detailed
background and a review of recent work.

Figure 1: Rate Distortion curves for sequence: Gaming_360P-1e31. The RD curves
generated using any one of the proposed algorithms outperforms the curve for
X265 using the default $\lambda$ (black).

## Lagrangian Multiplier in Compression

There has been a tremendous amount of work done to establish what the
Lagrangian multiplier should be in the different iterations of the MPEG codecs
(H.263, H.264 and H.265). Being able to find the balance between the quality
of decoded images and channel capacity is the fundamental problem which the
rate distortion algorithm attempts to solve [10].

As alluded to in the introduction, the seminal work of Sullivan and Wiegand[2]
laid the foundation for an experimental approach to choosing an appropriate
$\lambda$. In that work, they established the relationship between
quantisation step size $Q$ (used for DCT coefficient quantisation) and the
distortion $D$ in a frame i.e. $D(Q)=aQ^{2}$, where $a$ is a constant relating
to the content. Minimising $J$ wrt $D$ then yields a relationship
$\lambda=D/b$ where $b$ is another content related constant. Combining these
two relationships then leads to a relationship between $\lambda$ and $Q$ as
follows

$\lambda=0.85\times Q^{2}$ (2)

where the constant was determined experimentally over a small corpus of 3
video clips (CIF resolution). This was explicitly presented for H.263 in 2001
[3].

Updates to these experiments then yielded similar relationships for H.264 and
H.265 (HEVC). Because of the introduction of bi-directional (B) frames the
constants are all different. In fact three different relationships were
established for each of the Intra (I), Predicted (P) and B frames as follows.

$\lambda_{I}=0.57\times 2^{(Q-12)/3}$ (3) $\lambda_{P}=0.85\times
2^{(Q-12)/3}$ (4) $\lambda_{B}=0.68\times\max(2,\min(4,(Q-12)/6))\times
2^{(Q-12)/3}$ (5)

Clearly the size of the corpus used for these experiments was quite small, so
there is room for exploration with a larger corpus. In addition, average
performance is known to be worse than individual adaptation as discussed
previously. Hence the following section outlines different methods which have
been used to determine a better Lagrangian multiplier for use in video
encoding.

### Adaptive Lagrange Multipliers

There have been a few attempts to adjust the rate distortion optimisation
within a video codec. In all cases, the authors have discovered that a better
$\lambda$ does exist per clip. The techniques are best considered as two
classes. In one class (M-Algorithms), $\lambda$ is adjusted away from the
default stated previously, using a constant of proportionality $k$ as follows

$\lambda_{\textrm{new}}=k\times\lambda_{\textrm{orig}}$ (6)

where $\lambda_{\textrm{orig}}$ is the default Lagrangian multiplier estimated
in the video codec, and $\lambda_{\textrm{new}}$ is the updated Lagrangian.
Algorithms in this class, learn a relationship between $k$ and features
extracted from the video stream. Zhang and Bull [11] used a single feature
${D}_{P}/{D}_{B}$, the ratio between the MSE of P and B frames. This feature
gives some notion of temporal complexity. Experiments based on the DynTex
database yielded a choice for $k$ as follows.

$k=a({D}_{P}/{D}_{B}+d)^{b}+c$ (7)

where $a,b,c\textrm{ and }d$ were experimentally calculated and different for
each codec tested (H.264 and H.265). In their work, $k$ was updated every GOP
and they report up to 11% improvement in BD-Rate. This approach of fitting a
simple curve was also used by Yang et al, except that they used a different
feature, a perceptual content measurement $PC$ as follows

$k=a\textrm{PC}-b$ (8)

where again $a,b$ were determined experimentally using a corpus of the Derfs
dataset. Unlike Zhang and Bull, this was a straight line fit, the loss in
complexity of the fit is compensated for by the increase in complexity of the
feature.

Ma et al [12] used four features with a SVM to determine $k$. The SVM was
determined using the DynTex dataset. They reported up to 2dB improvement in
PSNR and 0.05 improvement in SSIM at equal bitrates. The four video features
measured spatial and temporal complexity directly from texture and motion
measurements on the pixel level data itself. Hamza et al [13] also take a
Machine learning based approach but here using a semantic scene model. Their
work used the Derfs dataset and reported up to 6% BD-Rate improvement.

Since $\lambda$ is linked to $Q$ as explained previously, the second class of
algorithms (Q-algorithms) adjusted $\lambda$ implicitly by adjusting the
quantiser parameter $Q$. Papadopoulos et al[14] applied an offset to Q, in
HEVC, based on the ratio of the distortion in the P and B frames. Each QP was
updated from the previous Group of Pictures (GOP) using the following
equation:

$\textrm{QP}=a\times\textrm{MSE}_{\textrm{ratio}}-b$ (9)

where $a,b$ are constants determined experimentally. This lead to an average
BD-Rate improvement of 1.07% on the DynTex dataset. With up to 3% BD-Rate
improvement achieved for a single sequence. Taking a local, exhaustive
approach, Im and Chan [15] proposed encoding a frame multiple times in a
video. Each frame was encoded using a Quantiser Parameter from the set
($QP\in({QP,QP\pm 1,QP\pm 2,QP\pm 3,QP\pm 4})$). This led to an increased
complexity in coding in HEVC, however for a single sequence achieved a 14.15%
BD-Rate improvement.

However, both M/Q-Algorithms have overlooked the possibility that a direct
search using BD-Rate as the objective function, could lead to even more
improvements. In some way all these previous works have reported data showing
$k$ versus $BD-Rate$, but they use various summary models to simplify this
relationship across a corpus. By using instead Direct optimisation, we are
able to explore just how much more BD-Rate gains there are per-clip. This idea
is presented next.

Figure 2: BD Rate(%) vs $k$ using two clips ((NewsClip_720P-7745 and
CoverSong_720P-3dca). Top : Multi-Res search, where X’s are the initial search
results, the solid line is the fit spline interpolation and +’s are the
subsequent searches. Bottom : Brents Method $\Diamond$ and Golden Section
$\Box$.

## Direct Optimisation

The approach adopted here falls under the M-algorithm class discussed above.
Hence an estimate of $k$ is required such that
$\lambda_{\textrm{new}}=k\times\lambda_{\textrm{orig}}$. In this work however,
we directly minimise BD-Rate with respect to $k$ using a range of optimisation
strategies. Furthermore, we deploy the optimisation strategies using a modern
dataset of User Generated content containing realistic material provided by
YouTube [9]. The new objective function is now directly the BD-Rate [16] as a
function of $k$, $B_{r}(k)$ as follows.

$B_{r}(k)=\int_{D_{1}}^{D_{2}}(R_{1}(D)-R_{k}(D))dD$ (10)

where the integral is evaluated over the quality range $D_{1}..D_{2}$.
$R_{1}(D),R_{k}(D)$ are the RD curves corresponding to the default ($k=1$) and
the evaluated multiplier $k$ respectively. Each RD operating point is
generated at a constant bitrate within a range which matches typical streaming
media use cases. The overall optimisation process is then as follows.

  1. 1.

Generate an RD curve using $\lambda_{orig}$ and $R_{t}=256k:7M$. We use 11
operating points in this range.

  2. 2.

Define our BD-Rate objective function as specified above in equation 10. We
use the same polynomial-log fit for evaluating the integral as recommended in
[16].

  3. 3.

Starting from $k=1.0$, minimise $B_{r}(k)$, BD-Rate, wrt $k$ using any typical
optimisation routine.

Note that for every evaluation of $B_{r}(k)$ , eleven (11) encodes are
required as well as the subsequent BD-Rate calculation itself.

### Optimisation Algorithms

Three simple direct optimisation algorithms were tested. The first used a
multi-resolution search (multi-res). This is a simple brute force strategy in
which successively refined grid positions in $k$ were searched. The initial
grid was $k=0.2:0.4:3.0$. The value of $k$ which yielded the minimum
$B_{r}(k)$, $k_{\textrm{opt}}^{1}$ was calculated based on a spline curve fit
to the sparse points, evaluated at steps of $0.01$. Successive refinements
searching positions of $k_{\textrm{opt}}^{n}\pm\Delta$ were then performed
with $\Delta=0.2,0.1,0.05$ in succession. Figure 2 shows this process in
action using two clips. X’s indicate the initial evaluated positions, and the
solid line is the initial curve fit. The subsequent refinement grids are
indicated with +. This process therefore always required 15 evaluations of the
objective function.

Golden Section Search and Brent’s method [17] were used as the other two
optimisation schemes. These are smarter direct search methods using local
curve fits. As the multi-res method required fifteen (15) evaluations of the
objective function, these optimisation routines were terminated after fifteen
iterations and a tolerance of 0.02% in the objective function. This allowed a
fair comparison.

An example of the search results ($k$ vs BD-Rate) for the methods can be seen
in Figure 2 for two particular clips. Implementation details are in the next
section, but one key observation here is the very different appearance of each
curve. This shows that a summary model (e.g. used by Zhang and Bull) may be
unable to exploit all the gains available. Furthermore, we note that the
points discovered by the direct search methods deviate from the initial smooth
curve fit near the minimum. This shows the complexity of the surface in the
region of interest. Finally, note that the initial iterations of Brent’s and
Golden searches cover a very similar curve as the Multi-res approach. This
means that the complexity of the curves is actually limited to the region near
the minimum.

### Implementation details

For this work VideoLAN’s [18] implementation of H.265 was used ,
x265111Version: 3.0+28-gbc05b8a91. The codebase was modified to take as an
argument, the multiplier $k$ and hence alter the Lagrangian constraint 6 used
in RD control throughout the codec.

In order to generate Rate Distortion curves (RD-curve), each video was encoded
using the default input parameters for x265 at a constant bitrate. The target
bitrate was set within the range of 256k to 7M. Both SSIM and PSNR were used
as the distortion metric. The command invocation used for each video encode
was as follows:

x265 --input SEQ.y4m --bitrate <XXXX> \--tune-<YYYY> \--<YYYY> \--csv-log-
level 2 --csv DATA.csv --output OUT.mp4

Where SEQ, OUT and DATA are the filenames for the raw input file, output
encoded video and the log file for an encode respectively. XXXX is the target
bitrate and YYYY is the distortion metric used, either PSNR or SSIM.

### Dataset

Previous work in this space [11][12][14] has used generally low resolution
video with content that does not reflect current use. A sub-sample of videos
from the recently released YouTube-UGC Dataset222media.withyoutube.com [9] was
therefore selected for use in this paper333The specific location of the chunks
in the clips provided is available: http://www.mee.tcd.ie/ sigmedia/EI2020.
This should provide a better representation of widely used and viewed video
content. Each video clip used was 150 frames (5 seconds) in length to match
the chunks typically used by Dynamic Adaptive Streaming over HTTP (DASH)
systems. The selected clips contained Forty-four (44) at 720p resolution and
thirty-three (33) clips at 360p resolution, all at 30 fps. These clips
represented eleven (11) classes of video as specified by the YouTube team. A
sample of each class of the clips can be seen in Figure 3.

Figure 3: Example frames from each class in the YouTube dataset.

## Results

Figure 4 shows the average BD-Rate improvement for each sequence class for
each of the three direct optimisation methods. The bar chart shows that in
almost all cases, Brent’s and Golden methods are better than multi-res. Per
class gains in BD-Rate are within 0.5%-5%. For three classes (Live, Sports and
Vlog) multi-res is best. This is probably due to the restrictions on maximum
allowed iterations in the competing methods. The best BD-Rate improvement for
each class is listed in Table 1. Up to 21% is observered in an animation clip.
The corresponding BD-Distortion is less than 0.01dB. This highlights that the
gains are all in bitrate without loss in quality, and that there are
significant gains to be had per clip.

Figure 4: Average BD-Rate improvement (%) for each class of video. Both Brent’s Method and Golden Section Search require less iterations than Multi-Res search and often provide better BD-Rate improvement. Table 1: Table 1: Best BD-Rate (%) improvement for each class of sequence using the three direct optimisation methods Class | Brent | Golden | Multi-Res
---|---|---|---
Animation | 21.667 | 21.667 | 13.109
CoverSong | 1.839 | 2.279 | 2.377
Gaming | 1.858 | 1.762 | 1.665
HowTo | 5.023 | 5.239 | 5.037
LiveMusic | 1.570 | 1.810 | 1.776
LyricVideo | 15.879 | 15.879 | 10.383
MusicVideo | 6.417 | 6.760 | 0.516
NewsClip | 2.481 | 2.371 | 2.251
Sports | 8.391 | 7.583 | 8.378
TelevisionClip | 10.533 | 10.238 | 9.486
Vlog | 0.351 | 0.496 | 0.590

On average golden section and Brent’s search took 12.6 and 9.7 iterations
respectively. This was with a 0.02% BD-Rate improvement tolerance used as the
stopping citerion. Figure 1 shows a comparison over the RD curves used for a
single clip as an example. The gains are larger at lower bitrates as expected.
In this example, the BR-Rate gain was 2.5% and Brent’s method took the lowest
number of iterations : 8.

Videos which were encoded without the stopping criterion on the direct
optimiser had an additional improvement in BD-Rate. While a further increase
in compression performance is welcome, it comes with a heavy computational
cost. Figure 5 shows an example of convergence in the direct optimisers in
this case without the stringent stopping criterion. As can be seen the BD-Rate
gains over 8% after 35 iterations. But the lower plot shows the amount of
gains achieved after each iteration as a fraction of the final converged rate.
Roughly 80% of the potential maximum BD-Rate improvement can be achieved
within the first fifteen(15) iterations, therefore the stopping criterion are
not that detrimental to performance especially in the light of computational
load.

Figure 5: These two plots illustrate the convergence of the optimisation
methods used for a particular clip (Sports_720P-0b9e). Top : BD-Rate. Bottom :
Normalised BD-Rate improvement. The normalised plot shows the BD-Rate
improvement as a % of the maximum BD-Rate possible for this clip. Hence it
shows that 80% of the possible BD-Rate improvement is attained within the
first 15 iterations.

In practice, we are interested in the rate gains at a particular quality. For
example a 40dB quality target is a good rule of thumb for a streaming media
service. BD-Rate measures instead the average gain across the RD performance
curve. We can show our results in this light. Figure 6 is a cumulative
histogram of the percent bitrate improvement at 40 dB for the 77 clips encoded
in this paper. For 80% of the clips we can have up to 5% improvement at this
quality. This is a useful metric and further validates the proposed methods.

Figure 6: CDF vs bitrate improvement (%) at 40dB. Top to bottom: Multi-res,
Golden Section Search, Brents Method. At a target quality of 40 dB, 80% of the
clips had gains up to 5% improvement for Golden Section and Brents Method

Table 2 compares the results from this work with the results from Zhang and
Bull [11] applied to the dataset in this paper. The comparison is not strictly
fair because we do not adapt the $k$ every GOP and instead use a single GOP
for all the clips. However the results further reiterate the difference in the
shapes of the curves per clip. As can be seen our proposed algorithm estimates
values of $k$ that are quite different, and subsequently the BD-Rate
improvement is much better than that using $k$ as equation 7.

Table 2: Table 2: Comparison of the BD-Rate improvement using an estimated $k$ from [11] vs $k$ determined by our direct search (Brent’s method). As can be seen our estimated values are quite different and lead to much improved performance. However our adaptation scheme is different from that in Zhang and Bull. Negative BD-R implies worse performance than the unmodified codec Clip | $\textrm{k}_{\textrm{ZB}}$ | BD-R | $\textrm{k}_{\textrm{Ours}}$ | BD-R
---|---|---|---|---
Animation4268 | 3.173 | -5.81 | 0.777 | 0.591
CoverSong7360 | 3.444 | -4.23 | 1.071 | 0.178
Gaming6403 | 3.132 | -6.18 | 1.438 | 0.683
HowTo21c6 | 2.938 | -3.15 | 0.438 | 5.239
LiveMusic6452 | 3.632 | -0.85 | 1.145 | 0.234
LyricVideo739a | 2.491 | 0.89 | 1.699 | 1.412
MusicVideo3285 | 4.298 | -10.21 | 1.120 | 0.068
NewsClip7745 | 2.571 | -2.67 | 0.547 | 2.482
Sports0b9e | 3.108 | 1.89 | 1.784 | 8.391
TVClip02b8 | 2.112 | -3.12 | 0.587 | 2.887
Vlog11c5 | 3.579 | -2.12 | 1.096 | 0.166

## Conclusion

This paper has presented a direct optimisation technique for Lagrange
multiplier estimation in a codec. The approach in some sense represents an
upper bound for what can be achieved by adjusting the Lagrangian multiplier.
We report gains up to $21$% on a particular clip and $0.5-5$% per class
without affecting the quality. At a target quality of 40 dB, 80% of the clips
had gains up to 5%. Furthermore, the work here has used a much more realistic
database of content provided by YouTube and for a larger number of clips. The
results also show that the relationship between $k$ and BD-Rate varies
substantially between clips, and that the region near the minimum is not a
simple surface. These observations point to further work in per clip
optimisation using more sophisticated machine learning techniques that exploit
deeper features in the content.

## References

  * [1] S. Cass, “The age of the zettabyte cisco: the future of internet traffic is video [dataflow],” IEEE Spectrum, vol. 51, no. 3, pp. 68–68, 2014.
  * [2] G. J. Sullivan and T. Wiegand, “Rate-distortion optimization for video compression,” IEEE signal processing magazine, vol. 15, no. 6, pp. 74–90, 1998.
  * [3] T. Wiegand and B. Girod, “Lagrange multiplier selection in hybrid video coder control,” in Proceedings 2001 International Conference on Image Processing (Cat. No. 01CH37205), vol. 3, pp. 542–545, IEEE, 2001.
  * [4] A. Ortega and K. Ramchandran, “Rate-distortion methods for image and video compression,” IEEE Signal processing magazine, vol. 15, no. 6, pp. 23–50, 1998.
  * [5] G. J. Sullivan, J.-R. Ohm, W.-J. Han, and T. Wiegand, “Overview of the high efficiency video coding (hevc) standard,” IEEE Transactions on circuits and systems for video technology, vol. 22, no. 12, pp. 1649–1668, 2012.
  * [6] A. Aaron, Z. Li, M. Manohara, J. De Cock, and D. Ronca, “Netflix Technology Blog - Per-Title Encode Optimization.”  https://medium.com/netflix-techblog/per-title-encode-optimization-7e99442b62a2.
  * [7] I. Katsavounidis and L. Guo, “Video codec comparison using the dynamic optimizer framework,” in Applications of Digital Image Processing XLI, vol. 10752, p. 107520Q, International Society for Optics and Photonics, 2018.
  * [8] M. Covell, M. Arjovsky, Y.-c. Lin, and A. Kokaram, “Optimizing transcoder quality targets using a neural network with an embedded bitrate model,” Electronic Imaging, vol. 2016, no. 2, pp. 1–7, 2016.
  * [9] Y. Wang, S. Inguva, and B. Adsumilli, “Youtube ugc dataset for video compression research,” arXiv preprint arXiv:1904.06457, 2019.
  * [10] T. Wiegand, M. Lightstone, D. Mukherjee, T. G. Campbell, and S. K. Mitra, “Rate-distortion optimized mode selection for very low bit rate video coding and the emerging h. 263 standard,” IEEE Transactions on Circuits and Systems for Video Technology, vol. 6, no. 2, pp. 182–190, 1996.
  * [11] F. Zhang and D. R. Bull, “Rate-distortion optimization using adaptive lagrange multipliers,” IEEE Transactions on Circuits and Systems for Video Technology, pp. 1–1, 2018.
  * [12] C. Ma, K. Naser, V. Ricordel, P. Le Callet, and C. Qing, “An adaptive lagrange multiplier determination method for dynamic texture in hevc,” in 2016 IEEE International Conference on Consumer Electronics-China (ICCE-China), pp. 1–4, IEEE, 2016.
  * [13] A. M. Hamza, A. Abdelazim, and D. Ait-Boudaoud, “Parameter optimization in h. 265 rate-distortion by single frame semantic scene analysis,” Electronic Imaging, vol. 2019, no. 11, pp. 262–1, 2019.
  * [14] M. A. Papadopoulos, F. Zhang, D. Agrafiotis, and D. Bull, “An adaptive qp offset determination method for hevc,” in 2016 IEEE International Conference on Image Processing (ICIP), pp. 4220–4224, Sep. 2016.
  * [15] S. Im and K. Chan, “Multi-lambda search for improved rate-distortion optimization of h.265/hevc,” in 2015 10th International Conference on Information, Communications and Signal Processing (ICICS), pp. 1–5, Dec 2015\.
  * [16] G. Bjontegaard, “Calculation of average psnr differences between rd curves; vceg-m33,” tech. rep., 2001.
  * [17] B. P. Flannery, W. H. Press, S. A. Teukolsky, and W. Vetterling, “Numerical recipes in c,” Press Syndicate of the University of Cambridge, New York, vol. 24, p. 78, 1992.
  * [18] V. Organization, “x265, the free H.265/HEVC encoder.”  https://www.videolan.org/developers/x265.html.

## Author Biography

Daniel J. Ringis is a PhD researcher at Trinity College Dublin. He is from
Trinidad and Tobago, and previously studied Electrical and Computer
Engineering at the University of the West Indies, St. Augustine. His current
research is focused on video compression.

François Pitié is an Ussher Assistant Professor in Media Signal Processing in
the School of Engineering and part of the ADAPT Research Centre. His research
is concerned with media signal processing applications related to film
postproduction and film restoration. The aim of his research is to develop
tools that can enhance the cinematic experience of immersion.

Anil Kokaram received the Ph.D. degree from Cambridge University, U.K., in
1993. Then, he was a Research Fellow with the Cambridge University Engineering
Department and Churchill College Cambridge. In 1998, he established the
Sigmedia Group (www.sigmedia.tv) at Trinity College Dublin. In 2011 a company
he founded, Green Parrot Pictures, operating in the consumer image processing
space, was acquired by Google. He is currently working as a Tech. Lead in the
Chrome Media Group at Google, Mountain View CA. He continues to hold his
Professorship with the Department of Electronic Engineering, Trinity College
Dublin and is a Fellow of the College. He is the author of over 100 refereed
papers in various conferences and journals over the years. His research
interests include the general area of digital video processing and algorithms
for stereoscopic quality improvement and video processing in the life
sciences. He is best known for his work in automated motion picture
restoration.


#  Quantixar: High-Performance Vector Data Management System

Gulshan Yadav  Department of Information Technology
A. P. Shah Institute of Technology
Thane, Maharashtra
<EMAIL_ADDRESS>RahulKumar Yadav  Department of Information
Technology
A. P. Shah Institute of Technology
Thane, Maharashtra
<EMAIL_ADDRESS>Mansi Viramgama  Department of Information Technology
A. P. Shah Institute of Technology
Thane, Maharashtra
<EMAIL_ADDRESS>Mayank Viramgama  Department of Information Technology
A. P. Shah Institute of Technology
Thane, Maharashtra
<EMAIL_ADDRESS>Apeksha Mohite  Department of Information Technology
A. P. Shah Institute of Technology
Thane, Maharashtra
<EMAIL_ADDRESS>
###### Abstract

Traditional database management systems need help efficiently represent and
querying the complex, high-dimensional data prevalent in modern applications.
Vector databases offer a solution by storing data as numerical vectors within
a multi-dimensional space. This enables similarity-based search and analysis,
such as image retrieval, recommendation engine generation, and natural
language processing. This paper introduces Quantixar, a vector database
project designed for efficiency in high-dimensional settings. Quantixar
tackles the challenge of managing high-dimensional data by strategically
combining advanced indexing and quantization techniques. It employs HNSW
indexing for accelerated ANN search. Additionally, Quantixar incorporates
binary and product quantization to compress high-dimensional vectors, reducing
storage requirements and computational costs during search. The paper delves
into Quantixar’s architecture, specific implementation, and experimental
methodology.

###### Index Terms:

Vector database, retrieval, storage, large language models, approximate
nearest neighbor search.

##  I Introduction

The explosion of data produced by modern applications, from image repositories
to text corpora, presents challenges for traditional database management
systems. These systems, designed primarily for structured data, need help with
the complexity and high dimensionality inherent in representations commonly
used in machine learning and artificial intelligence tasks. Vector databases
provide a compelling solution. Vector databases unlock similarity-based search
by representing data points as numerical vectors within a multi-dimensional
space. This capability revolutionizes tasks such as image retrieval,
recommendation engine generation, and natural language processing. This paper
introduces Quantixar, a vector database project purpose-built for efficiency
and scalability. Quantixar tackles the critical challenge of managing high-
dimensional data, focusing on advanced indexing and quantization techniques.
At its core lies the utilization of HNSW (Hierarchical Navigable Small World)
indexing, a robust method for approximate nearest neighbor (ANN) search,
crucial for accelerating similarity queries [1, 6]. To further optimize
performance, Quantixar incorporates binary and product quantization. These
techniques compress high-dimensional vectors, reducing storage requirements
and computational costs during search processes. The optimizations made by
Quantixar are optional and can be configured by a user based on his
requirements. The effectiveness of vector databases hinges on their ability to
handle the “curse of dimensionality.” As the number of dimensions in a dataset
increases, the efficiency of similarity searches degrades[2, 3]. Quantixar’s
design addresses this challenge by employing cosine similarity for distance
calculations as a default method. This technique has proven to be more
resilient to the curse of dimensionality, where the effectiveness of Euclidean
distance degrades as the number of dimensions increases [4]. The remainder of
this paper delves deeper into Quantixar’s architecture, its specific
implementation of indexing and quantization methods, and the experimental
methodology used with the ANN Benchmark. We present a thorough analysis of the
results, highlighting the strengths and potential areas for future development
in the Quantixar system.

##  II Background

###  II-A Curse of Dimensionality

In vector databases and high-dimensional data analysis, the “curse of
dimensionality” represents a collection of phenomena that degrade the
performance and meaningfulness of traditional distance-based similarity
searches as the number of dimensions increases[4]. Understanding this curse is
crucial, as it highlights the core challenge systems like Quantixar are
designed to address.

####  II-A1 Theoretical Underpinning

Consider a scenario where data points are randomly distributed within a unit
hypercube (a cube where each side has a length of 1). As the dimensionality
“d” of the space increases, the volume of this hypercube grows exponentially.
Imagine a smaller hypercube inscribed within the larger one, with its sides
being some fraction ’r’ of the outer hypercube’s sides. Surprisingly, as
dimensionality increases, even a substantial fraction like r=0.8 will cause
the inner hypercube’s volume to shrink towards zero relative to that of the
outer cube.

The consequence is that in high-dimensional spaces, the vast majority of the
volume is concentrated within a thin ’shell’ near the boundaries of the space.
This leads to a counterintuitive result: the distances between randomly
distributed points tend to become very similar. Mathematically, the variance
of these distances decreases, making it increasingly difficult to discriminate
between what we would intuitively consider “near” and “far” neighbors.

####  II-A2 Implications for Search Efficiency

Traditional distance metrics such as Euclidean distance, which many algorithms
like KNN (k-nearest neighbors) rely on, lose their effectiveness in high
dimensions[5]. Imagine trying to find the nearest neighbors to a query point.
As dimensionality increases, the search process needs to cover an
exponentially growing volume of space while the distances between points
become less informative. This combination leads to several challenges:

  * •

Computational Cost: Nearest neighbor search becomes computationally demanding,
potentially rendering real-time similarity-based tasks infeasible.

  * •

Data Sparsity: The volume of space proliferates that even large datasets
become sparse, hindering the discovery of meaningful relationships.

###  II-B Approaches to Address High-Dimensionality Challenges

####  II-B1 Hierarchical Navigable Small World (HNSW)

HNSW is a graph-based indexing technique renowned for its efficiency in
approximate nearest neighbor (ANN) search with high dimensional spaces. It
excels in scenarios where exact nearest neighbor results can be sacrificed for
significantly improved query speeds[6]. The core principles behind HNSW are:

  * •

Graph Construction: HNSW builds a multi-layered graph. Links between graph
nodes are established based on proximity and maintaining a ”small-world”
property, with long-range connections between nodes that might be distant but
share similar neighborhoods. It avoids the need to exhaustively explore every
local neighborhood, which would be time-consuming in high dimensions[6].

  * •

Hierarchical Structure: Lower graph layers contain increasingly refined
subsets of the dataset. This hierarchy allows search algorithms to narrow down
relevant portions of the data space quickly[6].

  * •

Greedy Traversal: During a query, HNSW employs a greedy search algorithm that
traverses the graph. Starting at the top, it determines which general region
of the space is most likely to contain the query’s nearest neighbors. With
each layer down, it zeroes in on an increasingly specific subset of data
points. The hierarchy lets HNSW rapidly discard vast, irrelevant portions of
the dataset, improving efficiency[6].

####  II-B2 Quantization: Compressing Vectors for Efficiency

Quantization methods strategically reduce the precision of high-dimensional
vectors, leading to significant gains in storage efficiency and potential
search speed improvements[7]. Quantixar focuses on two primary quantization
techniques:

  * •

Product Quantization (PQ):

Input:

    * –

Original vector: $\mathbf{x}\in\mathbb{R}^{d}$

    * –

Number of sub-vectors: $m$

    * –

Codebook size for each sub-vector: $k$

Steps:

    1. 1.

Partitioning: Divide $\mathbf{x}$ into $m$ sub-vectors:
$\mathbf{x}=[\mathbf{x}^{(1)},\mathbf{x}^{(2)},...,\mathbf{x}^{(m)}]$

    2. 2.

Codebook Construction: For each sub-space, learn a codebook:
$\mathcal{C}^{(i)}=\\{\mathbf{c}_{1}^{(i)},\mathbf{c}_{2}^{(i)},...,\mathbf{c}_{k}^{(i)}\\}\subset\mathbb{R}^{d/m}$

    3. 3.

Encoding: Quantize each sub-vector by finding its nearest centroid:

$q(\mathbf{x}^{(i)})=\arg\min_{\mathbf{c}\in\mathcal{C}^{(i)}}||\mathbf{x}^{(i)}-\mathbf{c}||^{2}$

    4. 4.

Quantized Representation:
$\hat{\mathbf{x}}=[q(\mathbf{x}^{(1)}),q(\mathbf{x}^{(2)}),...,q(\mathbf{x}^{(m)})]$

  * •

Binary Quantization (BQ):

Input:

    * –

Original vector: $\mathbf{x}\in\mathbb{R}^{d}$

    * –

Number of hyperplanes: $m$

Steps:

    1. 1.

Hyperplane Learning: Learn a set of hyperplanes with normal vectors
$\mathbf{u}_{1},\mathbf{u}_{2},...\mathbf{u}_{m}\in\mathbb{R}^{d}$

    2. 2.

Encoding: Generate a binary code $\mathbf{b}=[b_{1},b_{2},...,b_{m}]$ where:

$b_{i}=\begin{cases}1,&\text{if }\mathbf{u}_{i}^{T}\mathbf{x}\geq 0\\\
0,&\text{otherwise}\end{cases}$

    3. 3.

Binary Representation: The binary code $\mathbf{b}$ is the compact
representation of $\mathbf{x}$

####  II-B3 Advanced Techniques: Harnessing Hardware with SIMD

Modern processors offer a powerful tool for accelerating the computations at
the heart of vector databases: SIMD (Single Instruction, Multiple Data)
instruction sets. Found in extensions like SSE (Streaming SIMD Extensions) and
AVX (Advanced Vector Extensions), these instructions let the processor perform
the same operation on multiple pieces of data simultaneously, offering
significant performance gains[8, 9].

  * •

SSE (Streaming SIMD Extensions)

    * –

XMM Registers: The heart of SSE lies in its introduction of eight 128-bit XMM
registers (XMM0 through XMM7)[10]. These specialized registers are designed to
hold multiple packed data elements. For example, each XMM register can store:

      * *

Four 32-bit floating-point numbers

      * *

Eight 16-bit integers

      * *

Sixteen 8-bit integers

    * –

Parallelism Within Registers: SSE instructions operate on these packed data
elements simultaneously. Let us imagine we have two XMM registers containing
four 32-bit floating-point numbers each:

      * *

XMM0: [A1, A2, A3, A4]

      * *

XMM1: [B1, B2, B3, B4]

      * *

A single SSE addition instruction can perform the following in one operation:
[A1+B1, A2+B2, A3+B3, A4+B4]

  * •

AVX2 (Advanced Vector Extensions 2)

    * –

YMM Registers: AVX2 extends the XMM registers to a wider 256-bit format,
creating sixteen YMM registers (YMM0 through YMM15). These larger registers
double the potential data processing throughput[12, 11].

    * –

Fused Multiply-Add (FMA): AVX2 introduces support for FMA instructions. These
instructions combine a multiplication and addition operation into a single
step (e.g., D = (A * B) + C) [12, 11]. This type of calculation is common in
many vector and matrix operations, providing a performance boost.

    * –

Register Relationship: It is important to note that YMM registers encompass
the XMM registers. The lower 128 bits of any YMM register correspond to the
respective XMM register (e.g., the lower half of YMM0 is the same as XMM0).
This allows for some flexibility in code that mixes SSE and AVX2 instructions.

##  III Methodology

Figure 1: Quantixar System Architecture.

This section delves into Quantixar’s architectural design, providing a
comprehensive overview of its core components and their interplay. Illustrated
in Fig. 1, Quantixar’s architecture comprises four principal modules: Query
processing module, Quantization module, Indexing, and Storage Engine module.
Each module is critical to the system’s overall functionality and performance.

###  III-A Query Processing :

In this subsection, we introduce the concept of an “entity” as it pertains to
the Quantixar Engine and discuss the types of queries supported by this
engine, the various similarity functions employed, and the application
interfaces designed for user interaction.

Concept of Entity: Within the Quantixar system, an “entity“ is the basic unit
of data representation and is crucial to understanding the database’s
operation. An entity is typically a high-dimensional vector that serves as an
abstract representation of complex data items. These data items, such as
images, texts, audio files, or other forms of multimedia, are transformed into
vectors through feature extraction methods. The vector encapsulates the key
attributes of the original data, enabling efficient execution of operations
like similarity searches or clustering within the database.

Query Types: Quantixar supports two primitive query types:

  * •

Vector Query: The vector query is a foundational query type in Quantixar,
primarily focusing on vector similarity search. In this context, each data
entity is represented as a single vector. The primary function of a vector
query is to identify and retrieve entities that are most similar to a given
query vector[13]. This process involves computing similarity based on specific
metrics or distance measures, determining how closely entities match the query
vector. The vector query is critical in scenarios where the similarity between
entities, represented as vectors, is the primary criterion for data retrieval.

  * •

Metadata-Enhanced Vector Search (MEVS): It is a query method that augments
vector-based searches by incorporating attribute-based constraints. In this
approach, each vector is associated with metadata containing key attributes
and contextual details[14]. MEVS applies these attribute constraints to the
metadata, effectively narrowing the dataset to only those entities meeting the
specified criteria. The system conducts a vector similarity search within this
selectively filtered subset. This sequential method, which starts with
metadata-based filtering and then proceeds to vector similarity analysis,
ensures more efficient and pertinent query results. By focusing the vector
search within a refined subset, MEVS achieves faster response times and
heightened accuracy in matching[15, 16, 17, 18].

###  III-B Quantization

The Data Preprocessing / Quantization module in Quantixar is the critical
first step, laying the groundwork for efficient similarity search within the
high-dimensional vector space. This section focuses on transforming dense
vector data into a sparse vector. The configuration provided by the user
controls this transformation, which is done using quantization methods,
including product and binary quantization.

###  III-C Indexing

Quantixar’s indexing architecture ensures efficient and scalable retrieval of
high-dimensional vectors. It offers a dual indexing strategy, HNSW and Flat
Index, to accommodate diverse deployment scenarios and performance
optimizations. The choice between HNSW and flat indexing hinges on carefully
considering dataset size, desired query speed, and the tolerance for
approximation in search results. Flat Indexing, while simple, offers the
guarantee of finding the actual exact nearest neighbors to a query vector.
This makes it suitable for smaller datasets where computational cost is
manageable or in use cases where the utmost precision is non-negotiable.
However, as data points grow, the linear search time required for flat
Indexing quickly becomes a significant performance bottleneck. Flat Indexing
becomes impractical if the application demands scalability to large datasets.
This is where HNSW indexing excels. By constructing a hierarchical graph-based
structure on the vector dataset, HNSW enables approximate nearest neighbor
(ANN) search with substantially improved efficiency. Search times scale
logarithmically with the dataset size, allowing HNSW to handle datasets with
orders of magnitude larger than those feasible with flat Indexing. However, it
is essential to note the ”approximate” nature of the results. HNSW introduces
a degree of trade-off between absolute precision and search speed. With
careful tuning, HNSW can achieve remarkably high accuracy while dramatically
outperforming flat Indexing in speed.

By leveraging these indexing methods, and based on query characteristics and
dataset size, the Query Processing module will determine the type of index to
utilize and passes the query to the Storage Engine for data retrieval.

###  III-D Vector Storage

In the storage architecture of Quantixar, both RocksDB and etcd play crucial
roles in catering to different user needs. RocksDB, pivotal for vector
embedding persistence, was chosen for its high-performance characteristics in
key-value storage and is specially optimized for SSDs. Its data structure
handles efficient read and write operations, supported by block-size
management, thread pool, and advanced compaction strategies[19, 20]. Moreover,
MVCC management at the cluster level makes it ideal for high-speed data
processing environments.

Etcd, on the other hand, is integrated for environments requiring distributed
data management[21, 22]. etcd’s datastore is built on BBoltDB, offering
efficient caching via a single memory-mapped file managed by the operating
system and an in-memory B-tree index for quick key-revision mapping[22]. This
structure suits distributed systems, efficiently supporting range reads and
sequential writing. etcd’s key-value store is optimal for scenarios where high
availability, consistency, and fault tolerance in a distributed environment
are paramount[21].

Quantixar offers a flexible storage solution, allowing users to opt for
RocksDB or etcd based on their specific requirements, whether high-speed data
processing in non-distributed environments or robust, distributed data
management.

##  IV Results

This experiment evaluated the performance of the Hierarchical Navigable Small
World (HNSW) indexing algorithm within the Quantixar framework. Performance
was assessed on two datasets: Fashion-MNIST and SIFT 128. Key metrics are
reported in Table I, including construction time, insertion time, search time
(at different ’ef’ values), recall rates, and more.

Table I: Performance Metrics for HNSW in Quantixar on Various Datasets Metric | Fashion-MNIST | SIFT 128
---|---|---
Construction Time (s) | 0.22887 | 26.04492
Insertion Time (s) | 38.27 | 5655.94
Search Time (s) (ef=64) | 6.05864 | 130.16474
Search Time (s) (ef=128) | N/A | 34.84278
Mean Fraction of Neighbors Returned | 1.0 | 1.0
Last Distances Ratio (ef=64) | 1.001 | 1.0003
Last Distances Ratio (ef=128) | N/A | 1.0002
Recall Rate (ef=64) | 0.978 | 0.9908
Recall Rate (ef=128) | N/A | 0.9964
Queries per Second (ef=64) | 10647.42 | 4158.96
Queries per Second (ef=128) | N/A | 2051.11

System Specifications: AWS EC2 t4g.xlarge instance (4 vCPUs, 16GB RAM)

As expected, construction and insertion time increase for the more complex
SIFT 128 dataset. Notably, the recall rate is excellent across both datasets,
even at the lower ’ef’ value, indicating HNSW’s effectiveness in approximate
nearest neighbor search.

##  V Conclusion

In conclusion, the development of Quantixar reflects our commitment to
advancing the field of vector databases by learning from and building upon the
foundations laid by established systems like Qdrant, Weaviate, and Milvus.
Each of these systems offers unique insights into handling large-scale, high-
dimensional data, which have been instrumental in shaping the design and
functionality of Quantixar. With its integrated components, such as the
Quantixar Engine, Index, and Storage Engine, Quantixar aims to enhance
operations like similarity searches and efficient data retrieval. Recognizing
this field’s dynamic and collaborative nature, our approach is not just about
creating a new solution but also about contributing to the collective
knowledge and technology base. As we move forward, our focus will be on
empirical validation and performance benchmarking, ensuring that Quantixar not
only meets but also enriches the evolving standards of vector database
technology. Through this ongoing journey, we aim to solidify Quantixar’s place
as a valuable and innovative addition to high-dimensional data processing.

## References

  * [1] Malkov, Y. & Yashunin, D. Efficient and robust approximate nearest neighbor search using Hierarchical Navigable Small World graphs. CoRR.
  * [2] Zhong, Y. Efficient Similarity Indexing and Searching in High Dimensions. (2015)
  * [3] Jin, Y., Wu, Y., Hu, W., Maggs, B., Zhang, X. & Zhuo, D. Curator: Efficient Indexing for Multi-Tenant Vector Databases. (2024)
  * [4] Altman, N. & Krzywinski, M. The curse (s) of dimensionality. Nat Methods. 15, 399-400 (2018)
  * [5] Kouiroukidis, N. & Evangelidis, G. The effects of dimensionality curse in high dimensional knn search. 2011 15th Panhellenic Conference On Informatics. pp. 41-45 (2011)
  * [6] Wang, M., Xu, X., Yue, Q. & Wang, Y. A comprehensive survey and experimental comparison of graph-based approximate nearest neighbor search. ArXiv Preprint ArXiv:2101.12631. (2021)
  * [7] Gray, R. Vector quantization. IEEE Assp Magazine. 1, 4-29 (1984)
  * [8] Hassaballah, M., Omran, S. & Mahdy, Y. A Review of SIMD Multimedia Extensions and their Usage in Scientific and Engineering Applications. The Computer Journal. 51, 630-649 (2008)
  * [9] Mustafa, D., Alkhasawneh, R., Obeidat, F. & Shatnawi, A. MIMD Programs Execution Support on SIMD Machines: A Holistic Survey. IEEE Access. 12 pp. 34354-34377 (2024)
  * [10] Russinoff, D. SSE Floating-Point Instructions. Formal Verification Of Floating-Point Hardware Design: A Mathematical Approach. pp. 241-246 (2021)
  * [11] Guide, P. Intel® 64 and ia-32 architectures software developer’s manual. Volume 3B: System Programming Guide, Part. 2, 1-64 (2011)
  * [12] Gepner, P. Using AVX2 instruction set to increase performance of high performance computing code. Computing And Informatics. 36, 1001-1018 (2017)
  * [13] Rashtchian, C., Woodruff, D. & Zhu, H. Vector-matrix-vector queries for solving linear algebra, statistics, and graph problems. ArXiv Preprint ArXiv:2006.14015. (2020)
  * [14] Jokela, S. & Others Metadata enhanced content management in media companies. (Helsinki University of Technology,2001)
  * [15] Chiu, W., Yen, G. & Juan, T. Minimum manhattan distance approach to multiple criteria decision making in multiobjective optimization problems. IEEE Transactions On Evolutionary Computation. 20, 972-985 (2016)
  * [16] Bookstein, A., Kulyukin, V. & Raita, T. Generalized hamming distance. Information Retrieval. 5 pp. 353-375 (2002)
  * [17] Leydesdorff, L. On the normalization and visualization of author co-citation data: Salton’s Cosine versus the Jaccard index. Journal Of The American Society For Information Science And Technology. 59, 77-85 (2008)
  * [18] Seal, A., Karlekar, A., Krejcar, O. & Gonzalo-Martin, C. Fuzzy c-means clustering using Jeffreys-divergence based similarity measure. Applied Soft Computing. 88 pp. 106016 (2020)
  * [19] Matsunobu, Y., Dong, S. & Lee, H. Myrocks: Lsm-tree database storage engine serving facebook’s social graph. Proceedings Of The VLDB Endowment. 13, 3217-3230 (2020)
  * [20] Kim, H., Park, J., Jung, S. & Lee, S. Optimizing RocksDB for better read throughput in blockchain systems. 2019 23rd International Computer Science And Engineering Conference (ICSEC). pp. 305-309 (2019)
  * [21] Larsson, L., Tärneberg, W., Klein, C., Elmroth, E. & Kihl, M. Impact of etcd deployment on Kubernetes, Istio, and application performance. Software: Practice And Experience. 50, 1986-2007 (2020)
  * [22] Nalawala, H., Shah, J., Agrawal, S. & Oza, P. A comprehensive study of “etcd”—an open-source distributed key-value store with relevant distributed databases. Emerging Technologies For Computing, Communication And Smart Cities: Proceedings Of ETCCS 2021. pp. 481-489 (2022)