\section{Introduction}
 One of the major areas of interest in hadron physics and QCD in recent days is to understand the mechanical properties like mass, angular momentum and pressure distribution inside the nucleon in terms of quarks and gluons \cite{Polyakov:2018zvc,Lorce:2018egm,PhysRevD.99.094026}. The structure of the nucleon can be studied by bombarding it with a high energy probe like an electron. The mechanical properties of the nucleon are encoded in the gravitational form factors (GFFs), which are the form factors of the energy-momentum tensor.  The GFFs are functions of the square of the momentum transfer in the process ($Q^2$). They give information on how matter couples to gravity. These GFFs are related to the generalized parton distributions (GPDs)
that can be accessed in exclusive electron-proton scattering process, for example deeply virtual Compton scattering (DVCS). The extensive experimental programs at HERA \cite{H1:1999pji,H1:2001nez,ZEUS:2003pwh}, HERMES \cite{HERMES:2001bob}, COMPASS \cite{ATLAS:2018zdn} and JLAB \cite{CLAS:2001wjj,CLAS:2020yqf} in the recent past have given valuable data as inputs for the extraction of the GPDs. The upcoming electron-ion collider (EIC) will also probe the GPDs as well as give inputs for  understanding of the mechanical properties of the nucleon.  There are  four GFFs of the proton, $A(Q^2), B(Q^2), C(Q^2)$ and $\overline{C}(Q^2)$. The GFFs $A(Q^2)$ and $B(Q^2)$ are related to the mass and angular momentum distributions of the proton. Conservation of the energy-momentum tensor constraints the GFFs  $A(Q^2), \  B(Q^2)$, and $\overline{C}(Q^2)$; however, $C(Q^2)$, also known as the D-term, is not related to any Poincare generator and is unconstrained by such conservation laws. This form factor, also known as the D-term, contributes to the DVCS process  when the skewness $\xi$ is non-zero, or when there is a longitudinal momentum transfer from the initial state proton to the final state proton. The pressure distribution inside the proton is related to the D-term. The pressure distribution is usually defined in the Breit frame, and it is subject to  relativistic corrections.  In Ref.\cite{Lorce:2018egm} the pressure distribution is defined in the infinite momentum frame, or equivalently in the light-front formalism \cite{Freese:2021czn}. As the transverse boosts in the light-front formalism are Galilean, one can describe the two-dimension pressure distribution of a relativistic system using overlaps of light-front wave functions (LFWFs) for the GFFs. The D-term is zero for a free fermion and $-1$ for a boson. For an interacting system, in order to ensure stability, the D-term has to be negative.  Recently the D-term has been extracted from the JLab data; the pressure distribution inside the nucleon obtained from fits of JLab data is found to be repulsive near the center and attractive towards the periphery of the nucleus \cite{Burkert:2018bqq}. At the core, the pressure is even higher than that inside a neutron star, which is the most dense object in the universe. The JLab result generated a lot of interest in this field and in the recent past, quite a lot of theoretical studies have been performed on the shear and pressure distributions in the nucleon in different models for example in the Bag model \cite{Neubelt:2019sou}, chiral quark soliton model \cite{Goeke:2007fp,Wakamatsu:2007uc,Schweitzer:2002nm,Goeke:2007fq,Wakamatsu:2006dy}, AdS/QCD motivated diquark model \cite{Chakrabarti:2015lba} as well as a simple multipole model \cite{Lorce:2018egm}. Lattice results are also available \cite{Hagler:2003jd,Gockeler:2003jfa,LHPC:2010jcs,LHPC:2007blg,QCDSF-UKQCD:2007gdl,Deka:2013zha}.

Most of these models are phenomenological models for the proton and they do not incorporate any gluonic  degree of freedom. However, the D-term in particular, which is related to the pressure distribution inside the nucleon; involves the light-cone transverse  component of the energy-momentum tensor, that  depends on the quark-gluon interactions and hence the gluons are expected to play a major role here.  While  it is non-trivial to incorporate gluons in models for the nucleon,  which is a bound state, it is often interesting to replace the hadron state by  a simpler relativistic spin-$1/2$ state like a quark dressed with a gluon at one loop in QCD. This may be thought of as a perturbative model having a gluonic degree of freedom. The dressed quark state can be expanded in Fock space in terms of multi-parton LFWFs. These LFWFs can be calculated analytically using the light-front Hamiltonian. Upto one loop, this model incorporates the full quark-gluon interaction.
The GPDs can be expressed in terms of overlaps of these LFWFs. Previously the GPDs as well as the DVCS amplitudes  have been calculated in a similar QED model \cite{PhysRevD.71.014038}. The  Wigner functions and spin distributions have been investigated in \cite{More:2017zqq,Mukherjee:2014nya}. In this work, we investigate the GFFs and the pressure distributions for a quark dressed with a gluon at one loop. The twist two GPDs as well as the GFFs $A(Q^2), \  B(Q^2)$ require a calculation of the `good' light-cone component of the energy-momentum tensor, whereas the remaining two GFFs, $C(Q^2)$ and $\overline{C} (Q^2)$ requires  a calculation of some of the `bad' components, which include the interaction terms. In order to calculate the matrix elements of these terms, we use the two-component formulation of light-front QCD developed in \cite{Zhang:1993dd} in light-cone gauge, $A^+=0$, where $A^\mu$ is the gluon field. The `bad' component of the fermion field in this case is constrained  and can be eliminated from the operator expressions, which would then be written only in terms of the dynamical quark and gluon fields.  Matrix elements of these interaction terms in the operator structures would involve particle number changing off-diagonal contributions, when expressed in terms of the multi-parton LFWFs. There would be diagonal contributions as well, that does not change the particle number.  In this work, we keep all terms upto $O(g^2)$, where $g$ denotes the quark-gluon coupling. A similar perturbative model was used in \cite{Metz:2021lqv} in QED, where the D-term  and the pressure and shear distribution were calculated for an electron dressed with a photon at one loop in QED using the Feynman diagram approach and in the Breit frame.

The plan of the paper is as follows.  In section \ref{model} we introduce the gravitational form factors in a light-front dressed quark model, in section  \ref{Dterm} we discuss the $D$- term and the pressure distributions. Finally in section \ref{con} we present the summary and conclusion. Details of the calculation and useful formulas are given in the appendices.

\section{Gravitational form factors in a Light-front dressed quark model}
\label{model}
\subsection{Light-front wave function}
In this section, we outline the formalism used to calculate the gravitational form factors (GFFs) for a dressed quark state, which we call dressed quark model (DQM).  Instead of a proton state, we take  a quark dressed with a gluon.  This is a composite spin-$1/2$ state, which in the light-front (LF) formalism is treated fully relativistically. To list a few merits of this formalism :
\begin{itemize}
    \item Due to the presence of gluon dressing, the model employs a  gluonic degree of freedom \cite{Harindranath:1996sj,Harindranath:1998pd}.
    \item The dressed quark state can be expanded in terms of light-front wavefunctions (LFWFs).  Although the LFWF of a bound state, like a proton, cannot be calculated analytically, the LFWF for a dressed quark can be calculated analytically in perturbation theory \cite{PhysRevD.63.045006}.
    \item  One can write the LFWFs in terms of relative momenta that are frame independent \cite{Brodsky:1997de}. Thus LFWFs are boost invariant.
    \end{itemize}

In LF Hamiltonian formalism, the dressed quark state can be expanded in Fock space, where the contributions come from a  single quark state, a quark and a gluon state, quark and two gluons state, and so on. However, we truncate the Fock space expansion upto the two-particle sector. Such truncation, in LF framework,  is boost invariant. The dressed quark state is written as   \cite{Harindranath:1998pd, Harindranath:1998fm}
\begin{eqnarray} \label{state}
|P,\lambda \rangle
&=& \psi_1 (P, \lambda) b^{\dagger}_{\lambda}(P)|0 \rangle + \sum_{\lambda_1, \lambda_2}\int[k_1][k_2]   \sqrt{2(2\pi)^3P^+} \delta^3 (P-k_1-k_2)\ \psi_2 (P,\lambda|k_1,\lambda_1;k_2,\lambda_2) b^{\dagger}_{\lambda_1}(k_1) a^{\dagger}_{\lambda_2}(k_2)|0\rangle,  \nonumber\\
\text{where}~~ [k]&=& \frac{dk^+ d^2{\boldsymbol k}^{\perp}} {\sqrt{2(2\pi)^3{k^{+}}}}.
\end{eqnarray}
In Eq. (\ref{state}), $\psi_1(P, \lambda)$ in the first term, corresponds to a single particle state with momentum (helicity) $P$ ($\lambda$) and is  also responsible for the wavefunction normalization. The two-particle LFWF, $\psi_2(P,\lambda|k_1,\lambda_1;k_2,\lambda_2)$ is related to the probability amplitude of finding two particles namely a quark and a gluon with helicity $\lambda_1$ and $\lambda_2$ inside the dressed quark state.
$b^\dagger$ and $a^\dagger$ correspond to the creation operator of quark and gluon respectively.

In LF Hamiltonian framework \cite{Zhang:1993dd}, one uses the constraint equations in the light-cone gauge to eliminate the redundant degree of freedom and express the fields in terms of physically independent degrees of freedom. The quark fields are decomposed as \cite{Zhang:1993dd}
\begin{equation}
\psi=\psi_+\ +\ \psi_-,
\end{equation}
where $\psi_{\pm}=\Lambda_{\pm} \psi$ and $\Lambda_{\pm}$ are the projection operators.

We use the  two-component framework developed in \cite{Zhang:1993dd}, where using a suitable  representation of the gamma matrices one can write:
\begin{align}
 	\psi_+= \begin{bmatrix}
 		\xi\\0
 	\end{bmatrix}, ~~~~\psi_-=\begin{bmatrix}
 	0\\ \eta
 \end{bmatrix},
 \end{align}
where, the two-component quark fields are given by
\begin{eqnarray}
\xi(y) &=& \sum_{\lambda}\chi_{\lambda}\int \frac{[k]}{\sqrt{2(2\pi)^3}}[b_{\lambda}(k)e^{-ik\cdot y}+d^{\dagger}_{-\lambda}(k)e^{ik\cdot y}],
\\
\eta(y) &=& \left(\frac{1}{i\partial^+}\right)\left[\sigma^{\perp}\cdot\left(i\partial^{\perp}+g A^{\perp}(y)\right)+im\right]\xi(y),
\end{eqnarray}
and the dynamical components of the gluon field are given by
\begin{align}
	A^{\perp}(y)= \sum_{\lambda} \int \frac{[k]}{\sqrt{2(2\pi)^3k^+}}[{\bf\epsilon}^{\perp}_{\lambda}a_{\lambda}(k)e^{-i k \cdot y}+ {\bf \epsilon}^{\perp*}_{\lambda}a^{\dagger}_{\lambda}(k)e^{i k \cdot y}].
\end{align}
Here we have suppressed the color indices.

The LFWFs are written in terms of relative momenta \cite{Harindranath:1998pd, Harindranath:1998pc} and hence are independent of the momentum of the bound state. The Jacobi momenta $x_i$, ${\boldsymbol\kappa}^\perp_i$ are defined  such that they satisfy the relation $x_1+x_2=1$ and ${\boldsymbol\kappa}^\perp_1+{\boldsymbol\kappa}^\perp_2=0$.
\begin{eqnarray}
k_i^+=x_iP^+,~~~~ {\boldsymbol k}^{\perp}_i={\boldsymbol\kappa}^\perp_i+x_i \bsP,
\end{eqnarray}
where $x_i$ is the longitudinal momentum fraction for the quark or gluon , inside the two-particle LFWF.
The boost invariant two-particle LFWF can be written as,
\begin{eqnarray}\label{BILFWF}
\phi^{\lambda a}_{\lambda_1,\lambda_2}(x_i,{\boldsymbol\kappa}^\perp_i)&=& \bigg[\frac{x(1-x)}{{\boldsymbol\kappa}^{\perp2}+m^2(1-x)^2}\bigg]\frac{g}{\sqrt{2(2\pi)^3}}\frac{T^a}{\sqrt{1-x}} \chi_{\lambda_1}^{\dagger}\nonumber\\
&\times&\bigg[\frac{-2({\boldsymbol\kappa}^\perp\cdot \epsilon_{\lambda_2}^{\perp*})}{1-x}-\frac{1}{x}(\tilde{\sigma}^{\perp}\cdot{\boldsymbol\kappa}^\perp)(\tilde{\sigma}^{\perp}\cdot \epsilon_{\lambda_2}^{\perp*})+im(\tilde{\sigma}^{\perp}\cdot \epsilon_{\lambda_2}^{\perp*})\frac{1-x}{x}\bigg]\chi_{\lambda} \psi_1^{\lambda}
\end{eqnarray}
where, $\phi^{\lambda a}_{\lambda_1,\lambda_2}(x_i,{\boldsymbol\kappa}^\perp_i )= \sqrt{P^+} \psi_2 (P,\lambda|k_1,\lambda_1;k_2,\lambda_2)$,  $g$ is the quark-gluon  coupling. $T^a$ and ${\boldsymbol \epsilon}_{\lambda_2}^\perp$ are colour SU(3) matrices and polarization vector of gluon. The quark mass and the two-component spinor for the quark are denoted by $m$ and $\chi_\lambda$ respectively, $\lambda=1,2$ correspond to helicity up/down. We have used the notation  $\tilde{\sigma}_1=\sigma_2$ and $\tilde{\sigma}_2=-\sigma_1$ \cite{Harindranath:2001rc}.
It is customary to define the four momenta in light-front as
\begin{eqnarray}
P^\mu &=&(P^+,\bsP,P^-).
\end{eqnarray}
We choose Drell-Yan frame (DYF), so the longitudinal momentum transfer $q^+=0$.  Thus, the initial and the final state four momenta will be
\begin{eqnarray}\label{initialmom}
P^{\mu} &=&\bigg(P^+, {\bf0}^{\perp}, \ \frac{M^2}{P^+}\bigg),\\
\label{finalmom}
P^{\prime\mu}&=&\bigg(P^+,\ {\boldsymbol q}^{\perp},\ \frac{{\boldsymbol q}^{\perp 2}  + M^2}{P^+}\bigg),
\end{eqnarray}
and the invariant momentum transfer
\begin{eqnarray}\label{momtranfer}
q^\mu&=&(P^{\prime}-P)^\mu=\bigg(0, \ {\boldsymbol q}^{\perp},   \frac{{\boldsymbol q}^{\perp 2}}{P^+}\bigg).
\end{eqnarray}
In our frame,  $\bsP=0$ and $q^2=-{\boldsymbol q}^{\perp 2}$.

\subsection{Energy-momentum tensor  and extraction of gravitational form factors}
We start by writing QCD Lagrangian
\begin{eqnarray}
\mathcal{L}_{QCD} &=&  \frac{1}{2}\overline{\psi}\left(i \gamma_{\mu}D ^{\mu}-m\right)\psi-  \frac{1}{4} F^{\mu \nu}_a F^a_{\mu \nu},  \end{eqnarray}
where
the covariant derivative $
iD^{\mu} = i\overleftrightarrow{\partial}^\mu+gA^{\mu} $ and
$\alpha (i\overleftrightarrow{\partial}^\mu)\beta
= \frac{i}{2}\alpha\left(\partial^{\mu}\beta\right) -\frac{i}{2}\left(\partial^{\mu}\alpha\right) \beta$. The field strength tensor for non-Abelian gauge theory is
\begin{eqnarray}
F^{\mu \nu}_a = \partial^{\mu} A^{\nu}_a - \partial^{\nu} A^{\mu}_a +  g \ f^{abc} A^{\mu}_b A^\nu_c.
\end{eqnarray}
$\psi$ and $A^{\mu}$ is the fermion and boson field respectively.
The above Lagrangian can be used to obtain the gauge invariant symmetric  energy-momentum tensor \cite{Harindranath:1997kk}
\begin{eqnarray}
\theta^{\mu \nu} = \frac{1}{2}\overline{\psi}\ i\left[\gamma^{\mu}D^{\nu}+\gamma^{\nu}D^{\mu}\right]\psi - F^{\mu \lambda a}F_{\lambda a}^{\nu} + \frac{1}{4} g^{\mu \nu} \left( F_{\lambda \sigma a}\right)^2 - g^{\mu \nu} \overline{\psi} \left(
i\gamma^{\lambda}D_{\lambda} -m
\right)
\psi. \label{emtqcd}
\end{eqnarray}
The last term in Eq. \ref{emtqcd} goes to zero when one writes Lagrangian equation of motion of fermions. We are interested only in the fermionic part of the energy-momentum tensor for this work. So, Eq. \ref{emtqcd} reduces to
\begin{eqnarray}
\theta^{\mu \nu}_\textsf{\scriptsize Q} = \frac{1}{2}\overline{\psi}
i\left[\gamma^{\mu}D^{\nu}+\gamma^{\nu}D^{\mu}\right]
\psi \label{emtqcdforquark}
\end{eqnarray}
The matrix element of totally symmetric energy-momentum tensor (EMT) encodes the information of GFFs. In order to calculate the matrix elements, we use the dressed quark state defined in Eq. \ref{state} and extract the GFFs from the quark part of the symmetric energy-momentum tensor $\theta^{\mu \nu}_Q$.

For a spin half system, the standard parameterization used to obtain the four GFFs from the symmetric energy-momentum tensor up to $O(q^2)$ is \cite{Ji:2012vj,Harindranath:2013goa}
\begin{eqnarray}
	\langle P^{\prime},S^{\prime}|	\theta^{\mu\nu}_i(0)|P,S \rangle &=&\overline{U}(P^{\prime},S^{\prime})\bigg[-B_i(q^2)\frac{\overline{P}^{\mu}\ \overline{P}^{ \nu}}{M}+\left(A_i(q^2)+B_i(q^2)\right)\frac{1}{2}(\gamma^{\mu}\overline{P}^{\nu}+\gamma^{\nu}\overline{P}^{\mu})\nonumber\\
	&+& C_i(q^2)\frac{q^{\mu}q^{\nu}-q^2g^{\mu\nu}}{M}+\overline{C}_i(q^2)M\ g^{\mu\nu}\bigg]U(P,S),
	\label{e36}
\label{FF}
\end{eqnarray}
where $\overline{P}^{\mu}=\frac{1}{2}(P^{\prime}+P)^{\mu}$ is the average nucleon four momentum. $\overline{U}(P',S')$, $U(P,S)$  are the Dirac spinors and $M$ is the mass of the target state, $i \equiv(Q,G)$. $A_i$, $B_i$, $C_i$ and $\overline{C}_i$ are the quark or gluon form factors. In the Drell-Yan frame, as discussed above, there is no momentum transfer in the longitudinal direction and from Eq. (\ref{FF}), one can easily extract the form factors $A_\textsf{\scriptsize Q}(q^2)$ and $B_\textsf{\scriptsize Q}(q^2)$ from the diagonal component of the energy-momentum tensor $\theta^{++}_{\textsf{\scriptsize Q}}$. $A_i(q^2) \ (B_i(q^2))$ is obtained from the sum of the helicity conserving (flip) states. The outline of the calculation is briefly described below.
The matrix element that one needs to calculate can be written in the compact form as
\begin{eqnarray} \label{matrixelement}
\mathcal{M}^{\mu \nu }_{SS'} = \frac{1}{2}\left[\langle P'  ,S'|	\theta^{\mu \nu }_\textsf{\scriptsize Q}(0)|P,S \rangle \right]
\end{eqnarray}
where the Lorentz indices $(\mu, \nu)\ \equiv \{+,-,1,2\}$, $(S, S') \equiv \{ \uparrow,\downarrow \}$  is the helicity of the initial and final state. $\uparrow(\downarrow)$ positive (negative) spin projection along $z-$ axis. Using Eq. \ref{matrixelement} we obtain
\begin{eqnarray}\label{rhsA}
 \mathcal{M}^{++}_{\uparrow \uparrow} + \mathcal{M}^{++}_{\downarrow \downarrow} &=& 2\  (P^+)^2A_\textsf{\scriptsize Q}(q^2), \\
 \label{rhsB}
\mathcal{M}^{++}_{\uparrow \downarrow} + \mathcal{M}^{++}_{\downarrow \uparrow} &=& \frac{ i q^{(2)}}{M} \ (P^+)^2 B_\textsf{\scriptsize Q}(q^2) .
\end{eqnarray}
The most apt way to extract GFFs $C_\textsf{\scriptsize Q}(q^2)$ and $\overline{C}_\textsf{\scriptsize Q}(q^2)$ is using transverse components $(1, 2)$ of the energy-momentum tensor
\begin{eqnarray}\label{rhsC}
\mathcal{M}^{11}_{\uparrow \downarrow} + \mathcal{M}^{11}_{\downarrow \uparrow} - \mathcal{M}^{22}_{\uparrow \downarrow} - \mathcal{M}^{22}_{\downarrow \uparrow} &=& i\bigg[\frac{B_\textsf{\scriptsize Q}(q^2)}{4M}- \frac{C_Q(q^2)}{M}\bigg]\left((q^{(1)})^2 \ q^{(2)}-(q^{(2)})^3\right), \\
\label{rhsCbar}
\mathcal{M}^{11}_{\uparrow \downarrow} + \mathcal{M}^{11}_{\downarrow \uparrow} + \mathcal{M}^{22}_{\uparrow \downarrow} + \mathcal{M}^{22}_{\downarrow \uparrow} &=&
i\bigg[B_\textsf{\scriptsize Q}(q^2)\frac{q^2}{4M}-C_\textsf{\scriptsize Q}(q^2)\frac{3q^2}{M}+\overline{C}_\textsf{\scriptsize Q}(q^2) 2M\bigg]q^{(2)}.
\end{eqnarray}
Here $q^{(i)}$ are the components of $\bf{q}^\perp$.  By knowing the form of $B_\textsf{\scriptsize Q}(q^2)$ from Eq. \ref{rhsB} and from Eq. \ref{rhsC} we obtain the analytical expression for $C_\textsf{\scriptsize Q}(q^2)$. We can now easily calculate the form of $\overline{C}_\textsf{\scriptsize Q}(q^2)$.
The final expression for the four independent GFFs are as follows:
\begin{eqnarray}\label{gffAexp}
A_\textsf{\scriptsize Q}(q^2)&=& 1+ \frac{g^2\ C_F}{2 \pi ^2}\left[\frac{11}{10} -\frac{4}{5}\left(1+\frac{2 m^2}{q^2}\right) \frac{f_2}{f_1}-\frac{1}{3} \log \left(\frac{\Lambda ^2}{m^2}\right)\right]\\
B_\textsf{\scriptsize Q}(q^2)&=& \frac{g^2 C_F}{12\pi^2}\ \frac{m^2}{q^2}\ \frac{f_2}{f_1},\label{gffBexp}\\
D_\textsf{\scriptsize Q}(q^2)&=&  \frac{5 g^2 C_F}{6\pi^2}\  \frac{m^2}{q^2} \ \bigg(1-f_1f_2\bigg)= 4 \ C_\textsf{\scriptsize Q}(q^2),\label{gffDexp}\\
\overline{C}_\textsf{\scriptsize Q}(q^2)&=& \frac{g^2 C_F}{72\pi^2} \ \left(29-30 \ f_1 \ f_2+3~\mathrm{log}\bigg(\frac{\Lambda^2}{m^2}\bigg)\right),
\label{gffCbarexp}\end{eqnarray}
where
\begin{eqnarray}
f_1 :&=&\frac{1}{2}\sqrt{1+\frac{4 m^2}{q^2}}, \\ f_2:&=&\log\left(1+\frac{q^2\left(1+2f_1\right)}{2 m^2}\right)
\end{eqnarray}
$C_F$ is the colour factor and $\Lambda$ is the ultra-violet cut-off.
The expression for all the form factors in terms of the overlap of LFWFs are given in Appendix \ref{appA}.
The necessary forms of the integrals required for momentum integration in the transverse plane are also shown in Appendix  ~\ref{appB}. Then it is straightforward to obtain the compact form of the GFFs using Eq. (\ref{BILFWF}) and integrating over the Jacobi momenta ${\boldsymbol\kappa}^\perp$ and $x$.
 The transverse momentum integrals of the GFFs $B_\textsf{\scriptsize Q}(q^2)$ and $C_\textsf{\scriptsize Q}(q^2)$ are convergent and hence are finite and so is the GFF $D_\textsf{\scriptsize Q}(q^2)$ as it is four times $C_\textsf{\scriptsize Q}(q^2)$. However, while calculating the transverse momentum integrals of the GFFs $A_\textsf{\scriptsize Q}(q^2)$ and $\overline{C}_\textsf{\scriptsize Q}(q^2)$ we encounter ultra-violet (UV) divergences which we regulate using the UV cut-off $\Lambda$. This is introduced as the upper limit of the transverse momentum integration. The longitudinal momentum integral involved in the calculation of the GFF $A_\textsf{\scriptsize Q}(q^2)$ receives a contribution from the wave function normalization at $x=1$ \cite{Harindranath:1998pd, Harindranath:1998pc}, which contributes in the single particle sector upto the order we are interested in.  In this work, we included the contribution from the two-particle LFWF only, and we have included a cut-off at $x=1$ in the longitudinal momentum integration when calculating $A_\textsf{\scriptsize Q}(q^2)$.

\subsection{Numerical Analysis: Gravitational Form Factors}

In this section, we show the numerical results for the four GFFs viz $A(q^2)$, $ B(q^2)$, $D (q^2)$ and $\overline{C}(q^2)$  shown in Eqs.~\ref{gffAexp} -\ref{gffCbarexp}.
We now study the analytical results by fixing the parameters for all the GFFs:
the quark mass $m = 0.3 ~\mathrm{GeV}$, the coupling constant $g=1$, the color factor $C_F = 1$ and the UV cut-off $\Lambda = 10^3~\mathrm{GeV}$. We choose the cut-off such that $\Lambda>>M$, where $M$ is the mass of the dressed quark system. In this work, we take $M=m$.
\begin{figure}[h]

	\begin{minipage}{0.45\linewidth}
		\includegraphics[scale=0.3]{Images/plot_A.pdf}
	\end{minipage}

	\begin{minipage}{0.45\linewidth}
		\includegraphics[scale=0.3]{Images/plot_B.pdf}
	\end{minipage}
	\caption{Plot of the GFFs $A_\textsf{\scriptsize Q}(q^2)$ and $B_\textsf{\scriptsize Q}(q^2)$ as function of $q^2$, here $m=0.3$ GeV and $g=1$.}
	\label{figgffAnB}
\end{figure}

In Fig.~\ref{figgffAnB}, we plot the GFFs $A_Q(q^2)$ and $B_Q(q^2)$ as functions of squared momentum transfer $q^2 $. The behavior of the form factors as $q^2 \rightarrow 0$  has received substantial interest in the literature~\cite{Lowdon_2017,
Hagler:2003jd,Dorati:2007bk,Lorce:2018egm,Kim:2012ts,Jung:2014jja}. The light-cone momenta of the parton is related to $A_i(0)$ and the total angular momenta of the parton is  $J_i(0) = \frac{1}{2}\left( A_i(0) + B_i(0)\right)$.
The conservation of momentum and total angular momentum demands that $A(0) + B(0) =1$ and $J(0) = \frac{1}{2}\left( A(0) + B(0)\right)={1\over 2}$, summed over all quarks and gluons \cite{Teryaev:111207,Brodsky:2000ii,Teryaev:1999su}.
 We observe in Fig.~\ref{figgffAnB} that both $A_Q(q^2)$ and $B_Q(q^2)$ are positive over the chosen
$q^2$ range having their maxima at $q^2 = 0$. We see that $A_\textsf{\scriptsize Q}(0) = 0.7412$ for $\Lambda = 10^3~\mathrm{GeV}$.
The GFF $B_\textsf{\scriptsize Q}(0) = 0.0084$ which can be interpreted as the value of the anomalous gravitomagnetic moment of the quark in the DQM. One can expect that the contribution coming from the gluon anomalous gravitomagnetic moment would be negative such the total $B(0)$ vanishes \cite{Brodsky:2000ii}. The calculation of the gluon part of the GFFs and verifying the sum rules is part of a future publication.
$B_\textsf{\scriptsize Q}(q^2)$ is related to the second moment of the helicity flip generalized parton distribution (GPD) $E(x,q^2)$ by the spin sum rule such that $B_i(q^2) = \int dx~ x E_i(x,q^2)$ ~\cite{Ji_1997}. In the case of QED the first moment of the GPD $E_e (x,0)$ corresponds to the Schwinger value for the anomalous magnetic moment $\int dx E_e (x,0) = F_2(0) = \frac{\alpha}{2\pi}$ \cite{PhysRevD.71.014038}. Our approach leads $B_e(0) = \frac{\alpha}{3\pi}$ for the case of electron, as expected  where $\alpha$ is the fine structure constant.

  \begin{figure}[htp!]

	\begin{minipage}{0.45\linewidth}
		\includegraphics[scale=0.3]{Images/plot_D.pdf}
	\end{minipage}

	\begin{minipage}{0.45\textwidth}
		\includegraphics[scale=0.3]{Images/plot_Cbar.pdf}
		\label{gffDDQM}
	\end{minipage}
	\caption{Plot of the GFFs $D(q^2)$ and $\overline{C}_\textsf{\scriptsize Q}(q^2)$ as function of $q^2$, here $m=0.3$ GeV and $g=1$.}
		\label{figgffDnCbar}
\end{figure}

In Fig.~\ref{figgffDnCbar} shows our result for the GFFs $D_Q(q^2)$ and $\overline{C}_\textsf{\scriptsize Q}(q^2)$ as function of $q^2$. We observe that $D_\textsf{\scriptsize Q}(q^2)$ is negative which means that the DQM behaves like a bound system for the chosen set of parameters. The value of $D_\textsf{\scriptsize Q}(0) = -\frac{5 g^2 \ C_F}{72\pi^2}$ which for our model parameter has numerical value of $D_\textsf{\scriptsize Q}(0) = -0.007$. The GFF $\overline{C}_\textsf{\scriptsize Q}(q^2)$ is found to be positive for the chosen range of $q^2$ and the value of $\overline{C}_\textsf{\scriptsize Q}(0) = 0.067$.
The sum of quark and gluon contribution for $\overline{C}(0) $ is constrained to be zero. Since we are getting a positive contribution from the quark,
 so we expect the gluon contribution to be negative in accordance with the constraint. The calculation of the gluon part is part of a future publication.  The sign of both $A_\textsf{\scriptsize Q}(q^2)$ and $\overline{C}_\textsf{\scriptsize Q}(q^2)$ is $\Lambda$ dependent.\\

\section{D-term and the Pressure Distribution}
\label{Dterm}
The fundamental information like mass and spin of a particle is encrypted in the EMT which couples gravity to matter \cite{Pagels:1966zza,Polyakov:2002yz,Polyakov:2018zvc,Polyakov:2018exb,belitsky2005unraveling,diehl2003generalized}. The EMT, also offer internal mechanical properties via  \emph{D}-term which is called as the last unknown global property \cite{Polyakov:2002yz,Polyakov:2018zvc,Polyakov:2018exb}. As mentioned in the introduction, the \emph{D}-term is not related to any Poincare generators and is unconstrained.  The recent experimental extraction of the quark \emph{D}-term  resulted in the first experimental estimate of the pressure distribution inside the proton \cite{Burkert:2021ith}.

  \begin{figure}[htp!]

	\begin{minipage}{0.45\linewidth}
		\includegraphics[scale=0.35]{Images/plot_D_compare.pdf}
	\end{minipage}
	\caption{Plot of electron GFF $D_e(q^2)$ as function of $q^2$, here we set $m=0.511$ MeV, $\alpha=\frac{1}{137}$. }
		\label{Dcompare}
\end{figure}
The electron gravitational form factor $D_e(q^2)$ was studied in \cite{Metz:2021lqv} using Feynman diagram approach. It was found in \cite{Metz:2021lqv} that the total D-term for a dressed electron is divergent at $q^2=0$, the divergent comes from the photon part of the QED energy-momentum tensor. A non-zero photon mass acts as a regulator for this divergence. In order to compare our result for $D_e(q^2)$ with Metz {\it et al.} (MPR) we use eq. [\ref{gffDexp}] and set the model parameters to QED domain such that, $m=0.511$ MeV and $g^2=4\pi \alpha$, where $\alpha$ is the fine structure constant. We observe that our result obtained in LFWF approach is in excellent agreement with MPR as seen in the plot in Fig. \ref{Dcompare}.

\subsection{Pressure and shear forces for a quark dressed with a gluon}

We obtain the $D$-term from the transverse components of EMT that in turn is related to the pressure and shear distributions  \cite{Polyakov:2002yz,Polyakov:2018zvc,Polyakov:2018exb}.
These quantities can be well studied in the impact parameter space  by performing a Fourier transformation.
We work in light-front framework, where the light-front time is $x^+$. At  $x^+=0$, the longitudinal LF coordinate $x^-$ can be integrated out from the EMT and hence we are left with two transverse coordinates \cite{Lorce:2018egm}. As the transverse boosts in the light-front framework are Galilean, one can define the pressure and shear distributions in 2D, which are free from relativistic corrections. The light-front definitions of 2D pressure and shear distributions are related to the 3D definitions in the Breit frame by Abel transformation \cite{Panteleeva:2021iip}. In the Drell-Yan frame, a two-dimensional Fourier transform of the $D$-term from the momentum space ${\boldsymbol q}^{\perp}$ to impact parameter space $\bsb$ \cite{Burkardt:2000za} gives the mechanical properties like the shear and pressure distributions in the transverse impact parameter space .
The expressions for pressure and shear distributions in two-dimensions \cite{Freese:2021czn} are

\begin{eqnarray}
\label{Prefun}
p(\bsb)&=&\frac{1}{2M \bsb}\frac{d}{d\bsb} \left[\bsb\ \frac{d}{d\bsb} D_\textsf{\scriptsize Q}(\bsb)\right]-M \ \overline{C}_\textsf{\scriptsize Q}(\bsb),\\
s(\bsb)&=&-{\bsb\over M} \frac{d}{d\bsb}\left[ \frac{1}{\bsb} \frac{d}{d\bsb} D_\textsf{\scriptsize Q}(\bsb)\right],
\end{eqnarray}

where
\begin{eqnarray}
\label{fgff}
F(\bsb)
&=&  \frac{ 1}{(2\pi)^2}~\int d^2 {\boldsymbol q}^{\perp} \ e^{-i{\boldsymbol q}^{\perp} \bsb} \mathcal{F}(q^2) \nonumber \\
&=&\frac{1}{2\pi}\int_0^{\infty} d  {\boldsymbol q}^{\perp 2} J_0\left( {\boldsymbol q}^{\perp} \bsb\right)\mathcal{F}(q^2),
\end{eqnarray}
where $\mathcal{F}=(A, B, C, \overline{C})$. $J_0$ is the Bessel function of the zeroth order and $\bsb$ represents the impact parameter. $M$ is the mass of the dressed quark state.

In order to study the spatial distribution,  we take  wave packet states  in position space centered at origin instead of plane waves. The most prevalent forms are {\it Gaussian wave packets} \cite{Diehl:2002he, Chakrabarti:2005zm}.  The dressed quark state which is confined in transverse momentum space and has definite longitudinal momentum can be written as
\begin{eqnarray}
\frac{1}{16\pi^3}\int \frac{d^2 \bsp dp^+}{p^+}\phi\left( p\right) \mid p^+,\bsp,\lambda \rangle
\end{eqnarray}
with $\phi(p)=p^+\ \delta(p^+-p_0^+)\ \phi\left(\bsp\right)$.
We choose a Gaussian shape for $\phi\left( \bsp\right)$ in transverse momentum :
\begin{eqnarray}\label{gaussian}
\phi\left(\bsp\right)
= e^{-\frac{\bspp}{2\Delta^2}}
\end{eqnarray}
where $\Delta$ is the width of Gaussian. Since we choose to work in LF, the longitudinal and transverse momentum integration can be done independently.

\subsection{Numerical analysis: pressure and shear forces}

In Fig.~\ref{figPresuure} we have shown the plot for the pressure $2\pi \bsb \ p(\bsb)$ and shear $2\pi \bsb \ s(\bsb)$ distributions in the impact parameter space. We have chosen wave packets over the plane waves as it not only yields the Fourier transformed pressure in the impact
parameter space \cite{Diehl:2002he} but also gives smooth plots for the distribution. The spread of the distribution depends on the width of the Gaussian $\Delta$. So, we study the dependence of $\Delta$ in Fig.~\ref{figPresuure} by plotting the distributions for two different values of $\Delta$.

 In our model, pressure distribution is in accordance
with the von Laue condition \cite{Lane,Freese:2021czn} which is a stability condition such that \begin{eqnarray}\int_0^\infty d^2 \bsb ~p(\bsb) = 0.
\end{eqnarray}
The stability condition requires the presence of at-least one node in the pressure distribution which is seen in Fig.~\ref{figPresuure} (left panel).
The pressure profile shows a central positive core which becomes negative towards the tail region. The repulsive central pressure in contrast with the confining pressure in the outer region maintains the stability of the system \cite{Burkert:2018bqq}. We observe that the shear distribution is
non-negative in the b-space. This behavior of the shear force distributions is also observed in stable hydrostatic systems \cite{Polyakov:2018zvc}.
We also observe that as the width of the Gaussian wave packet, $\Delta$ increases the peak value of the distribution decreases and shifts away from the center in the impact parameter space.

The combination of pressure and shear defines the normal and the tangential forces experienced by a ring of radius $\bsb$ :
\begin{eqnarray}
\label{Fn-Ft}
F_n(\bsb)\es2\pi \bsb\left( p(\bsb) + \frac{1}{2} s(\bsb)\right),
\\
F_t(\bsb)&=& 2\pi \bsb \left( p(\bsb) - \frac{1}{2} s(\bsb)\right).
\end{eqnarray}
In Fig.~\ref{figForce} we show the plot for the normal force  $F_n$, and the tangential forces $F_t$ in the impact parameter space.
We see that the normal force  is positive. The positive nature of the normal component of the force ensures stability against collapse \cite{Polyakov:2018zvc}.
We observe that the tangential force is positive for smaller $\bsb$ and has a negative peak at around $\bsb=0.03$~fm   ($\bsb=0.1$~fm) for $\Delta=0.2$ ($\Delta=0.5$), while for larger $\bsb$ it is zero. The presence of both positive and negative regions maintains stability in the tangential direction.

 The pressure and shear distributions for a dressed electron in QED in  Breit frame was calculated in \cite{Lorce:2018egm}. The pressure and shear distributions coming from the fermionic part of the QED EMT is negative and approaches zero near the periphery. The 2D pressure and shear distribution for the proton in the Drell-Yan frame has been calculated in \cite{Chakrabarti:2020kdc} in a light-front diquark model with inputs from AdS/QCD. The behavior is similar as seen here in the DQM. The qualitative behavior is also similar as calculated in  the light cone sum rule approach \cite{Anikin:2019kwi}  and that extracted from JLab data for the nucleon \cite{Burkert:2018bqq}. Qualitative behavior of the normal and tangential force is  similar to that observed in phenomenological models for the nucleon \cite{Lorce:2018egm,Chakrabarti:2020kdc}.

\begin{figure}[h]

	\begin{minipage}{0.45\linewidth}
		\includegraphics[scale=0.3]{Images/PressurePlots_new_v2.pdf}
	\end{minipage}

	\begin{minipage}{0.45\textwidth}
		\includegraphics[scale=0.3]{Images/ShearPlots_new_v2.pdf}
					\end{minipage}
	\caption{Plots of (left panel) the pressure distribution $2\pi \bsb p(\bsb)$ and (right panel) the shear force distribution $2\pi \bsb s(\bsb)$ as a function of $\bsb$.}
	\label{figPresuure}
\end{figure}
\begin{figure}[htp!]

	\begin{minipage}{0.45\linewidth}
		\includegraphics[scale=0.3]{Images/NormalForcePlots_new_v2.pdf}
	\end{minipage}

	\begin{minipage}{0.45\textwidth}
		\includegraphics[scale=0.3]{Images/TangentialForcePlots_new_v2.pdf}

	\end{minipage}
	\caption{Plots of (left panel) the normal forces  $F_n(\bsb)$, and (right panel) the tangential forces $F_t(\bsb)$ as a function of $\bsb$.}
	\label{figForce}
\end{figure}
\subsection{The energy density and pressure distributions in DQM}
We also calculate the Galilean energy density, radial pressure, tangential pressure, isotropic pressure, and pressure anisotropy, which are defined in Ref ~\cite{Lorce:2018egm}
\begin{eqnarray} \label{Genergy}
\mu_i(\bsb) &=& M \left [\frac{1}{2} A_i(\bsb)+ \overline{C}_i(\bsb) + \frac{1}{4M^2}\frac{1}{\bsb}\frac{d}{d\bsb} \left( \bsb\frac{d}{d\bsb}\left[ \frac{1}{2}B_i(\bsb) - 4C_i(\bsb) \right]\right) \right], \\
\label{radialP}
\sigma^r_i(\bsb) &=& M  \left[
-\overline{C}_i(\bsb) + \frac{1}{M^2} \frac{1}{\bsb} \frac{d C_i(\bsb)}{d\bsb}
\right],  \\
\label{tangentialP}
\sigma^t_i(\bsb)
&=& M  \left[
-\overline{C}_i(\bsb) + \frac{1}{M^2} \frac{d^2C_i(\bsb)}{d \bsbb}
\right] \label{mp}, \\
\label{totalP}
\sigma_{i}(\bsb) &=& M  \left[
-\overline{C}_i(\bsb) + \frac{1}{2}\frac{1}{M^2} \frac{1}{\bsb} \frac{d}{d\bsb}\left(\bsb \frac{d \ C_i(\bsb)}{d \bsb} \right)
\right] , \\
\label{shearlike}
\Pi_i(\bsb)&=&  M  \left[
-\frac{1}{M^2} \bsb \frac{d}{d\bsb}\left(\frac{1}{\bsb} \frac{dC_i(\bsb)}{d \bsb} \right)\right].
\end{eqnarray}

In this section, we study five 2-dimensional distributions in impact parameter space $\bsb$. We first take the Fourier transform of all the GFFs from momentum space to impact parameter space as defined in Eq. \ref{fgff}. For all the two-dimensional distributions we choose a suitable value for the Gaussian width $\Delta=0.2$.

The radial pressure, tangential pressure and isotropic pressure are related to the form factor $C$ and $\overline{C}$ only. The relation between radial pressure, tangential pressure and isotropic pressure is evident from Eqs. \ref{radialP} to \ref{totalP}
\begin{eqnarray}\label{pressurerelation}
\sigma_\textsf{\scriptsize Q}&=& \frac{(\sigma^r_\textsf{\scriptsize Q}+\sigma^t_\textsf{\scriptsize Q})}{2},
\end{eqnarray}
We can also infer from Eq. \ref{radialP}, \ref{tangentialP} and \ref{shearlike} the relation between pressure anisotropy, radial pressure and tangential pressure as
\begin{eqnarray}
\label{shearrelation}
\Pi_\textsf{\scriptsize Q}&=& \sigma^r_\textsf{\scriptsize Q}-\sigma^t_\textsf{\scriptsize Q}.
\end{eqnarray}
Therefore, from Eq. \ref{pressurerelation} we can corroborate that three pressure must show similar behavior.
One can clearly anticipate that from Eq. \ref{shearrelation} that the nature of pressure anisotropy will have a bell shape.

\subsection{Numerical Analysis: The energy density and pressure distributions}
In this section we discuss, the distributions listed in Eqs. \ref{Genergy}-\ref{shearlike}.
\begin{figure}[h]

		\includegraphics[scale=0.3]{Images/EnergyDensityPlots1.pdf}
	\caption{Plot the 2D Galilean energy density $(2\pi \bsb \mu_\textsf{\scriptsize Q} (\bsb))$ as a function of $\bsb$.
We choose the Gaussian width $\Delta=0.2$.}
\label{energydensity}
\end{figure}
The Galilean energy density is nothing but the combination of diagonal components (from the way we extract the GFFs) of the EMT which encapsulates all the GFFs that can be seen from Eqs. \ref{rhsA} - \ref{rhsCbar} and Eq. \ref{Genergy}. In Fig \ref{energydensity}, we observe that the Galilean energy density increases with a peak value around $\bsb = 0.04~ \mathrm{fm}$ and then decreases at the periphery.

\begin{figure}[h]

	\includegraphics[scale=0.3]{Images/pressurePlotsMerged1.pdf}
	\caption{Plot of the 2D radial pressure $ (2\pi \bsb \sigma^r_\textsf{\scriptsize Q}$), tangential pressure  $(2\pi \bsb \sigma^t_\textsf{\scriptsize Q})$,  isotropic pressure  $(2\pi \bsb \sigma_\textsf{\scriptsize Q})$, as a function of $\bsb$. We choose the Gaussian width $\Delta=0.2$.}
	\label{allpressure}
\end{figure}
\begin{figure}[h]

		\includegraphics[scale=0.3]{Images/PresssureAnistropyPlots1.pdf}

	\caption{Plot of the 2D pressure anisotropy $(2\pi \bsb \Pi_{\textsf{\scriptsize Q}})$ as a function of $\bsb$. We choose the Gaussian width $\Delta=0.2$.}\label{anisotropy}
\end{figure}
The energy density shows unphysical behavior (negative) for  $\bsb \le 0.02\ \mathrm{fm}$, which is different from phenomenological models for the nucleon \cite{Lorce:2018egm, Chakrabarti:2020kdc}.
The tangential pressure flips the sign slightly away from the center. The total tangential force near the central region up to about $0.02$ fm is repulsive but becomes attractive after that till the outer region. This gives the signature that {\it quark contribution} dominates in the central region while the {\it gluon contribution} can be manifested in the small negative region away from the central region towards the periphery \cite{Lorce:2018egm}.

The radial pressure is always positive and is governed by the quark contribution. The isotropic pressure shows a slight negative region.
In Fig. \ref{anisotropy}, we show the plot of pressure anisotropy. Typically, it is expected that the total pressure anisotropy should vanish at the center that is one of the reasons for the radial pressure to be greater than tangential pressure. However, in our model, for the quark case, the pressure anisotropy shows a slight negative contribution around the central region. To understand this, it will be interesting to study the effect of the gluonic part of the EMT. Then, it becomes zero slightly away from the center and further, we observe that anisotropy increases monotonically till about $0.03$ fm and then again decreases. Overall, for a broad range of $\bsb$,  the 2D pressure anisotropy behaves in a similar way as in a diquark model and in a multipole model for the nucleon \cite{Lorce:2018egm}.

\section{Conclusion}
\label{con}

In this work we have investigated the GFFs and the mechanical properties like the pressure, shear and the energy distributions for a composite spin $1/2$ system, namely a quark dressed with a gluon at one loop in QCD. The GFFs and the mechanical properties of the proton in terms of its quark and gluon content have attracted a lot of interest  in recent days. There are a lot of theoretical studies in phenomenological models of the nucleon, there are also lattice calculations. Most of these theoretical models for the nucleon do not include gluons. The C form factor, or the D-term, that gives the pressure distribution inside the nucleon, in particular, is related to the non-light-cone (plus)  component of the energy-momentum tensor, and  therefore depends on the quark-gluon interaction.  In this work we replace the nucleon state by a composite spin-$1/2$ state, namely a quark dressed with a gluon at one loop in QCD, in order to calculate the GFFs. This may be thought of as a perturbative model that incorporates quark-gluon interaction. Our approach based on light-front Hamiltonian perturbation theory treats such a state fully relativistically. The LFWFs of the state   can be calculated analytically in perturbation theory, and the GFFs are expressed in terms of overlaps of LFWFs. We have used the LF two-component formalism to eliminate the constrained fields in light cone gauge.     In order to verify the correctness of the approach, we have compared with the existing results for the D-term for a dressed electron at one loop, in QED.  The pressure and shear distributions as well as the 2D energy density are calculated in the Drell-Yan frame when there is no momentum transfer in the longitudinal direction. We have compared our results with the results for the nucleon calculated in phenomenological models. Overall we have observed a qualitatively similar behavior, however the magnitude of the results in DQM  cannot be compared with that of a nucleon, which is expected. In this work we presented the GFF and the mechanical properties due to the fermionic part of the EMT. An interesting study would be to calculate the contributions coming from the gluonic part, this is part of ongoing work.


\section{Introduction}\label{sec:introduction}
K-nearest neighbor graphs (K-NNG), which contain the identities of the closest datapoints for each datum, are fundamental to many data science and machine learning techniques. The rapid growth of attainable data has therefore increased the need for practical K-NNG algorithms and efficient implementations thereof.\\
This paper presents an optimized implementation of the NN-Descent algorithm by Wei Dong et al. \cite{nndescent}. NN-Descent is an iterative, randomized heuristic which improves a random guess of the K-NNG.\\
A particularly popular implementation of NN-Descent is PyNNDescent \cite{pynndescent} which is implemented in Python and uses the Numba JIT copiler \cite{noauthor_numba_nodate} to achieve sufficient runtime performance. Particular strengths of PyNNDescent are its ease of use and support of custom distance metrics. PyNNDescent is used in the sci-kit learn compatible implementation of UMAP \cite{mcinnes_lmcinnesumap_2020} \cite{umap}. UMAP is currently gaining popularity as an alternative dimensionality reduction technique to t-SNE due to favorable embedding properties and because it is faster than current implementations of t-SNE \cite{becht_dimensionality_2019}.\\
The use of NN-Descent in practice makes performance optimized implementations especially desirable.\\
A major challenge in optimizing the performance of NN-Descent implementations is the irregular memory access pattern stemming from unordered input data together with the neighbor-of-neighbor heuristic explained in Section \ref{sec:nn_descent_algo}.\\
This paper presents a single-core implementation which is limited to the l2-distance metric at the benefit of a vast runtime reduction when compared to PyNNDescent.\\
Optimization techniques are discussed which performance for both low or high-dimensional input respectively. In particular a novel heuristic is introduced which increases locality by improving the otherwise irregular memory access pattern.\\
The overall performance advantage over PyNNDescent motivates the use of specialized implementations of NN-Descent for commonly used distance metrics and use cases.\\

Algorithms for K-NNG computation are still an active area of research due to their enormous importance. Alternative approaches include "Goldfinger" \cite{Ruas} which is based on bit-comparisons of hashes as well as methods building upon the same ideas as NN-Descent \cite{baron}.

\section{K-Nearest Neighbor Graph and NN-Descent}\label{sec:nn_descent_algo}
In this section we present an overview of the K-Nearest Neighbor Graph (K-NNG) problem and the NN-Descent algorithm optimized in our implementation.\\

A K-NNG is a directed simple graph $G = (V,E)$ such that we have a directed edge $(u,v) \in E$ iff $v$ is one of the k-nearest neighbors of $u$ according to a given distance metric. Let $V$ be the set of datapoints in $\mathbb{R}^{d}$ space where $|V| = n$. The dimensionality of all datapoints is denoted by $d$. The $k$ nodes adjacent to some node $u$ in $G$, denoted as $\text{adj}_{G}(u)$, are thus the $k$ nearest neighbors according to the given metric.

One trivial way of computing the KNN-Graph is to compute all $\frac{n(n-1)}{2}$ mutual distances and selecting the best $k$ points. This is not computationally viable for most datasets of practical interest, especially in data science applications due to the $O(n^{2})$ asymptotic number of distance evaluations.\\
NN-Descent by Wei Dong et al.\cite{nndescent} is a randomized, iterative heuristic which computes an approximation of the K-NNG. The NN-Descent algorithm requires far fewer distance evaluations (empirical cost $O(n^{1.14})$ \cite{nndescent}) at the expense of the quality of the resulting KNN-Graph.\\

The key insight used is that "a neighbor of a neighbor is also likely to be a neighbor"\cite{nndescent} (Figure \ref{fig:neighbor}).

\begin{figure}\centering
  \includegraphics[width=0.3\textwidth]{images/nn.pdf}
  \caption{A neighbor of a neighbor is likely to be also a neighbor \cite{nndescent}. In the current estimate of the K-NNG, N is considered a neighbor of NN and also of C. NN-Descent evaluates the distance between NN and C updates the current approximation of the K-NNG accordingly by replacing edge $($C,N$)$ with $($C,NN$)$. \label{fig:neighbor}}
\end{figure}

Using this assumption allows NN-Descent to iteratively improve the KNN-Graph estimate by evaluating distances primarily for promising pairs of vertices.\\
NN-Descent begins with a random initialization of $G$, so that for every vertex the set of $k$ nearest neighbors are uniformly sampled from $V$. Multiple iterations of the NN-Descent algorithm are then executed. Every iteration improves the K-NNG estimate by performing the following two steps for every vertex in $V$:
\begin{enumerate}
    \item Select candidates: find a suitable set $S$ with vertices that are neighbors-of-neighbors
    \item Calculate and update: calculate all pairwise distances in $S$ and update $G$ if closer vertices were found.
\end{enumerate}
We will refer to the first step as \textit{Selection step} and the second one as \textit{Calculation step}. Iterations are performed until the number of changes to $G$ falls below a specified threshold. In this case the estimated K-NNG is a reasonable approximation of the true K-NNG.\\
Section \ref{subsec:selection_step} contains more details on the selection of candidates. For more details on the algorithm itself refer to \cite{nndescent}.\\
As the NN-Descent algorithm is a heuristic it is necessary to validate both the quality of the returned K-NNG as well as the computational cost. Recall is used to measure how close the K-NNG approximation is to the true K-NNG. Our implementation achieved a recall of over $99\%$ on all examined datasets. Multiple parameters could if desired be altered to change the runtime-quality trade-off.\\
To calculate the number of floating point operations, the number of distance evaluations is counted. This enables performance comparisons between different code versions which have small but varying amounts of additional floating point comparisons. Together with the dimensionality of the datapoints the number of operations can be computed. For each l2-distance evaluation $d$ subtractions, $d$ multiplications, and $d-1$ additions are performed.

\section{Optimization of NN-Descent implementation}\label{sec:yourmethod}
Numerous possible optimizations were implemented and evaluated to yield a fast single-core implementation of NN-Descent. Some approaches aim at making the selection step highly efficient while others aim to improve the locality or optimize the mutual distance calculations. The relative impact of these optimizations differs depending on parameters of the problem, such as the dimensionality $d$ of the datapoints.\\

We started out with a C implementation\cite{Eskomski19} which adheres closely to the pseudo code found in \cite{nndescent} (NNDescentFull). We then studied the most popular implementation of NN-Descent,  PyNNDescent \cite{pynndescent}, which was introduced in Section \ref{sec:introduction}.
Although asymptotically the two implementations behave the same, the Python implementation was considerably faster than the straightforward C implementation. We identified the main differences and adopted the improvements

\subsection{Selection Step}\label{subsec:selection_step}
An important improvement found in PyNNDescent is the selection step, which is performed in each iteration for every vertex (see Section \ref{sec:nn_descent_algo}). These improvements were then expanded upon for an additional speed-up.

In the selection step of every iteration we need to find the neighborhood $N(u)$ for every node $u$ in the current KNN-graph approximation. The neighborhood $N(u)$ contains every node $v \in V$ for which $(u,v) \in E$, i.e. node $v$ is one of the k-nearest neighbors according to the current approximation. More interestingly, it also contains every node $w$ for which $(w,u) \in E$, so every node $w$ which has node $u$ as one of its k-nearest neighbors.
In the pseudo code as presented \cite{nndescent} such nodes $w$ are found by first inverting all the directed edges in the current KNN-graph $G$, resulting in a graph $G' = (V, E')$ for which $E' = \{ (u,v) | (v,u) \in E\}$. After the \textit{reverse} step, we then proceed by setting $N(u) = \text{adj}_{G}(u) \cup \text{adj}_{G'}(u)$ for all $u$ which we refer to as \textit{union} step. Note that while $\text{adj}_{G}(u)$ is bounded in size by k, $\text{adj}_{G'}(u)$ can contain up to $n$ elements, which requires the usage of a dynamically growing data structure.

Finally we \textit{sample} the neighborhood $N(u)$ to contain $\rho \cdot k$ elements. This reduces the number of the pairwise distance calculations of the neighborhood $N(u)$ which is quadratic in $|N(u)|$. The whole process of finding a suitably sized neighborhood for every node can be regarded as a composition of three functions: \textit{reverse, union} and \textit{sample}. In the basic implementation the intermediate results are stored in memory and three full passes over the entire K-NNG are required.

PyNNDescent improves over this by only doing one pass over the data but introducing a heap in the sampling process. Instead of building each neighborhood $N(u)$ and then sampling from it, PyNNDescent does both simultaneously in one pass over the KNN-graph. We adopted this change. For each edge $r = (u,v)$ a weight $r_e$ is drawn uniformly at random (u.a.r.) in $[0,1]$. Both $N(u)$ and $N(v)$ are implemented as heaps and we insert the node $v$ in $N(u)$ with weight $r_e$, we do the symmetric thing for $N(v)$. This corresponds to both the reverse and union step of the naïve selection implementation. Because the heaps are bounded in size, every neighborhood $N(v)$ ends up containing at most $\rho \cdot k$ elements. Selecting a subset of size $\rho \cdot k$ is equivalent to assigning a random weight u.a.r. to each element and selecting the $\rho \cdot k$ elements with the smallest weights. This gives a considerable speedup; on our synthetic dataset we observed a 16x runtime speedup (see Section \ref{subsec:exp_selection}).

Because the heaps incurred a lot of cache misses we further optimized the fused selection function described above in order to get rid of the heaps.
Upon every update of the KNN-graph we keep track of how large the neighborhood of every node $v$ is. Since when doing these updates we access the relevant data structures anyway, we do not incur any additional cache misses by these modifications. Knowing how large each neighborhood is allows us to simplify the sampling process: for every edge $e = (u,v)$ we insert $v$ into $N(u)$ with probability $\frac{\rho \cdot k}{|N(u)|}$. In expectation this is equivalent to the previous sampling procedure, but it works without heaps. This gives a small speedup of around 1.12x.

The optimized one-pass sampling step reduced the number of memory accesses by eliminating the intermediate results of the sequentially applied functions (\textit{reverse, union} and \textit{sample}).
Furthermore, we avoided the difficulty of controlling the (previously unbounded) size of the reverse graph $G'$.
This is important since for many relevant input sizes the KNN-graph does not fit into the caches found in commodity hardware. Multiple passes over the KNN-graph then lead to many cache misses.

\subsection{Greedy Reordering Heuristic}\label{subsec:greedy}
Datapoints which are close in the dataspace are frequently accessed together but the underlying data is not usually located closely together. This leads to a difficulty in exploiting spatial locality. Our roofline model analysis (cf. Section \ref{subsec:exp_roofline}) indicates that our implementation is memory-bound for low-dimensional inputs. The major difficulty of exploiting locality in spite of the non-uniform memory access pattern is the primary issue to be solved. Introducing an assumption on the data space distribution of the input data allows the development of a heuristic approach to tackle this problem.

Without any assumptions about the input distribution, our access pattern is irregular and we cannot improve locality. This is due to the tight relationship of temporal locality of two nodes and their distance in data space. We proceeded by assuming our input is clustered, meaning for every node all its $k$ nearest neighbors are within the same cluster (\textit{clustered assumption}). This assumption will allow us to partly recover those clusters from an early approximation of the K-NNG. After reordering memory such that the clusters are close together, we proceed with the remaining iterations of NN-Descent. Experiments on a synthetic dataset and on real world data set are promising (consult section \ref{sec:exp} for more details).

Recall that during the selection step we iterate over all edges $e \in E$ of the current KNN-graph approximation $G = (V,E)$. For one edge $e = (u,v)$ we will access both $\text{adj}_{G}(u)$ and $\text{adj}_{G}(v)$. Those two lists are likely to be in completely different locations in memory. Since the edge $e$ is part of our current KNN-graph, $u$ and $v$ are likely close in data space according to the given metric. This is why, especially after the initial iteration when our KNN-graph approximation becomes more accurate, closeness in data-space and temporal locality in the access pattern are highly correlated. For the remainder of this section we will consider clustered inputs (\textit{clustered assumption}).

Intuitively, after the first iteration a nodes nearest neighbor is likely to be in the same cluster. Recall that we start with a randomly initialized approximation, meaning every node has $k$ u.a.r. chosen nearest neighbors. The probability that within those $k$ nodes we do not have any node from the same cluster can be bounded from above by (for $c$ clusters, and $k$ neighbors):
$$\Pr[\text{all $k$ nodes not in cluster $i$}] \leq \Big (\frac{c-1}{c} \Big )^k$$

For a wide range of practical $k$ and $c$ this probability is sufficiently small. Armed with the intuition that the nearest neighbor of every node in our approximation is likely to be within the same cluster, we may now try to exploit that to reorder our memory. As a first step we want an algorithm that:
\begin{itemize}
    \item may only use the existence of those clusters, the input is \textit{not} ordered in any way revealing information about the structure of those clusters
    \item recovers most of the clusters. Moreover, it should output a permutation $\sigma: [n] \rightarrow [n]$ which we may use in the end to permute our data in memory all at once to bring the clusters together.
    \item makes at most one pass over the KNN-graph
\end{itemize}

The above requirements inform the design of our greedy clustering heuristic. In the pseudo code (Algorithm \ref{alg:greedy}), the permutations $\sigma$ and $\sigma^{-1}$ are modeled as $n$-dimensional arrays. We initialize them with the identity function: for each $i \in [n]$ we have $\sigma(i) = i$. We proceed by looking at node $i=0$. In each iteration of the outermost loop we would like to find a good candidate for the spot $i+1$, meaning whichever node permutation $\sigma$ maps onto $i+1$, it should be close in data space to node $i$. To achieve that we now sort the adjacency list of $i$ by distance, so $a_i[j]$ contains the identifier of the $j$'th closest node. Now we check whether the spot assigned to by $a_i[j]$ by permutation $\sigma$ is smaller than our current position, if so, we assume that $a_i[j]$ already has a good spot where it is close to its (data-space) neighbors in memory space. If not we check whether $a_i[j]$ already occupies the spot we would like to have it at ($\sigma(a_i[j]) = i+1$) - then we conclude the search and break out of the inner loop.

Otherwise if $\sigma(a_i[j]) > i+1$, we would like to set $\sigma$ such that $\sigma(a_i[j]) = i+1$, such that the node $a_i[j]$ will occupy spot $i+1$. Note that this is the desired outcome, the nodes on spots $i$ and $i+1$ are close together in data space. This specific sequence of two swaps in the permutation $\sigma$ and its inverse turn out to give the desired result. By creating and updating both the permutation and its inverse at the same time, we save ourselves a costly inversions of the permutation at several steps. This way we can satisfy the second requirement of only doing one pass through the KNN-graph.

The algorithm then returns a permutation $\sigma$. We proceed by permuting all of our data in memory using that permutation. Afterwards we continue with the remaining iterations of NN-Descent using the permuted memory layout. The copying itself is done all at once using $\sigma$.

When the clustered assumption is given, it is intuitive that our heuristic will succeed in clustering most of the data. Consult section \ref{sec:exp} for an experimental evaluation on a synthetic data set. More surprisingly we have even seen a small speedup on real world datasets where the clustered assumption does not hold.

\begin{algorithm}
\SetAlgoLined
\KwResult{Permutation $\sigma$ }
 $\sigma \gets id$\;
 $\sigma^{-1} \gets id$\;
    \For{$i\gets0$ \KwTo $n-1$ }{
        $a_i \gets \text{sorted}(\text{adj}_{G}(i))$\;
        \For{$j\gets0$ \KwTo $k-1$ }{
            \uIf{$\sigma(a_i[j]) < i+1$}{
                continue\;
            }
            \uElseIf{$\sigma(a_i[j]) = i+1$}{
                break\;
            }
            \ElseIf{$\sigma(a_i[j]) > i+1$}{
                swap in $\sigma$ entries $a_i[j]$ and $\sigma^{-1}(i+1)$ \;
                swap in $\sigma^{-1}$ entries $\sigma(a_i[j])$ and $i+1$ \;
                break\;
            }
        }
    }

 \caption{Greedy Clustering Heuristic}\label{alg:greedy}
\end{algorithm}

\subsection{Compute Step}\label{subsec:highdimension}

For higher dimensional datasets each l2-distance evaluation becomes more costly while the overhead of sampling and updating data structures remains constant. In such cases optimizing the distance evaluations becomes significantly more important than the optimizations described in previous sections.\\
A single l2 distance evaluation for two vectors is computed by summing the component-wise distances and taking the square root. As the actual value of the l2-distance is unimportant, the square root is omitted and the implementation uses the squared l2-distance. This improvement is not significantly impactful as for high dimensional vectors the cost of computing and summing the differences is dominant.\\
We decided to limit our implementation to vector dimensions which are divisible by 8 in order to simplify the use of AVX2 SIMD intrinsics. As each AVX2 register can hold 8 single-precision floating point numbers this alleviates the need for auxiliary code to handle the last (fewer than 8) components. The real-world datasets described in Section \ref{sec:exp} happen to fulfill this requirement without modification.\\
We use a AVX2 vector of accumulators for each distance evaluation and process 8 components at a time by subtracting to compute the difference and using an \emph{fmadd} instruction to square the difference and add it to the accumulator vector (cf. tag \textit{l2intrinsics} in Section \ref{sec:exp}).\\
We noticed that the restriction of the dimensionality to divisibles by 8 allowed for an easy modification to the way the datapoints are stored in memory. This modification allocates the data neatly aligned to 256 bits at the cost of at most additional 192 bits. This change significantly improves the performance since the \emph{loadu} instructions become faster. We did not observe a speedup by replacing \emph{loadu} intrinsic instructions by the equivalent \emph{load} intrinsic (cf. tag \textit{mem-align}).
The compute step for a single node can be further improved by blocking in order to compute the mutual distances between multiple vectors simultaneously. Modifications to the basic NN-Descent algorithms sampling step lead to a neighbourhood never exceeding a fixed size (in our implementation 50 nodes). For these, all the mutual distances have to be computed which can effectively be blocked.
We use a blocksize of 5 by 5 vectors (not scalars) in our implementation (Figure \ref{fig:matrix}).

\begin{figure}\centering
  \includegraphics[width=0.5\textwidth]{images/matrix.pdf}
  \caption{Illustration of blocked distance evaluations. For the neighborhood containing $v_{1}$ through $v_{10}$ the red cells represent distance evaluation. All distance evaluations within a 5 by 5 block are computed simultaneously. \label{fig:matrix}}
\end{figure}

For each such block we allocate an 256 bit AVX2 accumulator for every distance evaluation (either 10 or 25 in total). We then proceed with the computation of the squared l2-distance 8 features at a time for all combinations. In a block where 25 mutual distances are computed simultaneously only 10 AVX2 vectors of data are loaded per 8 dimensions. Without blocking each of theses would be loaded once for every distance evaluation. (1 vs. 25 loads per component). For high dimensional vectors especially this drastic reduction of data loads has a major impact on the performance (cf. tag \textit{blocked}).\\
Choosing a blocksize of 5x5 allows for each of the 25 accumulators to be allocated to a register which can be seen by inspecting the assembly. If the number of vectors to be compared to one-another is not divisible by 5 a more flexible but slower function is used for the remaining pairs.\\

\section{Experimental Results}\label{sec:exp}
This section discusses the experimental performance of the optimizations discussed in the previous section.
The main questions we aim to answer are the following:
\begin{enumerate}
    \item Does our clustering heuristic successfully cluster our data in memory?
    \item How does performance scale in terms of dataset size $n$ and dataset dimension $d$?

    \item Which optimizations are beneficial specifically in low or high dimensional settings?
    \item How does our final implementation perform to the numba based Python implementation PyNNDescent on real world data?
\end{enumerate}
 For all experiments in this section we used the squared euclidean distance and $k=20$ on single precision floats.\\

\mypar{Experimental setup} All the experiments were conducted on a computer running Ubuntu 20.04 LTS. The CPU is a Intel Core i7-9700K CPU @3.60GHz (turbo boost disabled) with the cache sizes L1: 256 KiB, L2: 2 MiB and L3: 12 MiB. The GCC version 9.3.0 was used in all experiments with the following flags: \emph{03}, \emph{ffast-math}, \emph{march=native}, \emph{flto}. The flags \emph{p} and \emph{pg} were used to add run time instrumentation to the code for specific profiling tasks.
\\

Four datasets were used to evaluate the performance and recall of our modified implementation. To investigate scaling behaviours, synthetically generated datasets were used. Additionally real-world datasets including the MNIST handwritten digits \cite{mnist} were used to verify the generalization of the performance from the synthetic datasets.
\begin{description}
\item[Synthetic Gaussian Dataset] The input parameter for the generation of this data set is the number of dimensions and the number of points. For the \textit{Single Gaussian Dataset} all points are drawn from one gaussian distribution centered at the origin. In the non-single variant, for each dimension a gaussian is created and centered around the canonical basis vector. For all evaluations the covariance is $2\cdot I_d$.
\item[Synthetic Clustered Dataset] A dataset designed to fulfill our clustered assumption. For every cluster we draw its points from a multivariate Gaussian. Mean and covariance are chosen such that the clustered assumption holds with high probability.
\item[MNIST Dataset \cite{mnist}] The MNIST database contains 70'000 images of handwritten digits given as 784 dimensional vectors of pixel intensity values.
\item[Audio Dataset] The audio dataset used in the NN-Descent publication by Wei Dong et al. \cite{nndescent}. Each of the 54'387 points consists 192 features which were extracted from recordings of English sentences.
\end{description}

\subsection{Selection Step}\label{subsec:exp_selection}
Every iteration of NN-Descent consists of a selection and computation step (cf. Section \ref{sec:nn_descent_algo}). Here we evaluate the improvements detailed in Section \ref{subsec:selection_step}. As previously mentioned the PyNNDescent inspired sampling is up to 16 times faster than the naïve sampling implementation in C. Our improved sampling, which we call \textit{turbosampling}, gives a speedup on top of that of up to 1.12x. The speedup was measured in terms of runtime since the flop count varies across those three implementations and a performance plot would be misleading. We used the \textit{Synthetic Gaussian Dataset} with parameters $n=16'384$, $d=8$ for said measurements. Due to space constraints, we omit further evaluations of improvements inspired by PyNNDescent and focus on our own contributions.

\subsection{Roofline model}\label{subsec:exp_roofline}
\begin{figure}\centering
  \includegraphics[width=0.4\textwidth]{images/roofline.pdf}
  \caption{Roofline plot with peak performance $\pi = 24$ [flops/cycle] and bandwidth $\beta = 4.77$ [bytes/cycle]. The \textit{Synthetic Gaussian Dataset} was used with $n=131'072$ and the dimensions $8$ and $256$ respectively.}
  \label{fig:roofline}
\end{figure}
Our optimizations in Section \ref{subsec:greedy} and \ref{subsec:highdimension} are motivated by the fact that the bottleneck of our algorithm changes as we change the dimension. For low dimensional data the computation is memory bound whereas for higher dimension we actually become compute bound as shown in figure \ref{fig:roofline}. In the following the derivation of the roofline model parameters is explained.

As stated in section \ref{sec:nn_descent_algo} the work $W(n)$ is calculated from the number of distance evaluations and the dimensionality of the datapoints.
For the sake of the roofline model analysis, we additionally need to reason about the data movement $Q(n)$, the bandwidth $\beta$ as well as the peak performance $\pi$.

\mypar{Bandwidth} Using the stream benchmark tool \cite{streambenchmark}, we measured the bandwidth $\beta = 4.77$[bytes/cycle] which seems reasonable as the maximal bandwidth from the manual is rarely attainable ($41.6 GB/s$ \cite{intelcpumanual} results in $12.41$[bytes/cycle]).

\mypar{Peak performance}
Utilizing 8-way AVX2 SIMD instructions for single precision floating point numbers the peak performance is $32$ [flops/cycle] on the Coffee Lake processor based on the Skylake microarchitecture.
This bound does not account for the instruction mix which is about $50\%$ 8-way subtractions and $50\%$ 8-way FMA's leading to a bound of $24$ [flops/cycle] which is used as our peak performance $\pi$.

\mypar{Data movement} The number of bytes transferred from memory to the cache far exceeds the size of the cache, which necessitates repeated loading and rewriting of data. To reason about the bytes transferred, cachegrind (Valgrind extension)\cite{cachegrind} was used. Cachegrind enables examination of the code's cache behaviour by simulating the memory behaviour in terms of first and last level (LL) caches. Using this tool we measured the LL cache data read and write misses. Table \ref{table:cachegrind} summarizes the cachegrind output on our \textit{Synthetic Clustered Dataset} with $n=131072$ and $16$ clusters.

\begin{table}\centering
\begin{tabular}{@{}ccc@{}}\toprule
& LL read misses & LL write misses  \\ \midrule
no-heuristic ($d=8$) & 122'150'286 & 14'777'070 \\
greedyheuristic ($d=8$) & 69'653'838 & 12'328'994 \\
no-heuristic ($d=256$) & 450'209'609 & 20'438'131\\
\bottomrule
\end{tabular}
\caption{Cachegrind results for versions of our implementation with (\textit{greedyheuristic}) or without (\textit{no-heuristic}) memory reordering on the \textit{Synthetic Clustered Dataset} ($n=131'072$ and $16$ clusters) for the specified dimension. Increasing $d$ by a factor of 32 increases the last-level read misses by a smaller factor.}
\label{table:cachegrind}
\end{table}

A major improvement of locality when using the greedy clustering heuristic can be seen in terms of the cache misses. Our heuristic nearly halves the number of LL cache data read misses. The increase in operational intensity by the reduction in cache misses moves the computation to the right in the roofline plot (\emph{no-heuristic, dim8} to \emph{greedyheuristic, dim8} in figure \ref{fig:roofline}).

If we increase the dimension we do more work $W(n)$ and achieve a higher operational intensity because the last-level read misses increase by a smaller factor. Since features for a single datapoint are continuous in memory, we attribute this to spatial locality.

\subsection{Greedy Clustering Heuristic Evaluation}
\begin{figure}\centering
  \includegraphics[width=0.4\textwidth]{images/clusters.pdf}
  \caption{Each line represents the fraction of datapoints belonging to single cluster in a 2000 sized window (y-axis) starting at the position given in the x-axis. Plot represents the cluster distribution \textit{greedyclustering} reordering retrieved on our \textit{Synthetic Clustered Dataset} with parameters $n=16'384$, $d=8$ and 8 clusters.\label{fig:clustering}}
\end{figure}

\mypar{Quality of Clustering with Greedy Heuristic}
This section investigates both how well the greedy clustering algorithm recovers the clusters of the underlying dataset and how this causes the speedup.

To test this, we ran our greedy clustering algorithm on the \textit{Synthetic Clustered Dataset}. We then reordered our data according to the permutation $\sigma$ found by the greedy clustering algorithm.
Figure \ref{fig:clustering} shows that the algorithm successfully recovers many of the clusters in the beginning of the dataset, where the heuristic can chose from many unassigned nodes.

Towards the end of the dataset all of the relative frequencies of the eight clusters are around $\frac{1}{8}$ - our heuristic stops working. This is expected since our algorithm is restricted by design to a single pass through the data. If a certain cluster was already handled and most of the points belonging to it were moved to some location, then the missing few will end up at the end of our simplified memory layout.

\begin{figure}[h!]
  \includegraphics[width=0.4\textwidth]{images/barplot.pdf}
  \caption{Time spent on each iteration on {Synthetic Clustered Dataset} with 16'384 points, 16 clusters and 8 dimensions. While the first iteration is slower due to the overhead of the heuristic, we profit in later iterations.\label{plot:iterations}}
\end{figure}

\mypar{Empirical Evaluation on Synthetic Dataset}
To test whether our heuristic leads to a speedup when the clustered assumption is satisfied we created such a dataset (\textit{Synthetic Clustered Dataset}).
We ran our algorithm twice on this dataset, once with greedy clustering and subsequent reordering of the data in memory and once without those changes. As reasoned in the previous section, we expect a shorter running time when doing the reordering. One can see in Figure \ref{plot:iterations} that the first iteration takes slightly longer to finish in the greedy clustering version of our algorithm due to overhead, but we can profit from the reordering in all subsequent iterations where we are consistently faster than the non reordered version. In total this amounts to a speedup of 18.46\% over all iterations.

\begin{table}\centering
\begin{tabular}{@{}ccc@{}}\toprule
runtime & MNIST & Audio \\ \midrule
no-heuristic & 12.12s & 4.78s\\
greedyheuristic & 11.45s & 4.53s\\
PyNNDescent & 24.41s & 14.47s\\
\bottomrule
\end{tabular}
\caption{Runtimes on the real-world \textit{MNIST} and \textit{Audio} datasets. Note how even though the clustered assumption might not hold, we still get a speedup on both datasets using our greedy clustering and subsequent reordering of data in memory. Our final implementation \textit{greedyclustering} is significantly faster on both datasets than PyNNDescent.\label{tbl:clusteredrealworld}}
\end{table}

\mypar{Real World Data Evaluation}
As previously hinted at we observed a speedup on real world data when applying our clustering heuristic. Consult Table \ref{tbl:clusteredrealworld} to see the runtimes on the \textit{MNIST} and \textit{Audio} datasets. According to our observations, the clustering heuristic does not manage to cluster MNIST semantically. Still, the reordering does improve locality which makes sense considering we move nodes together which are close in data space.

\subsection{Performance}
To show how the performance scales with number of input points $n$, we used the \textit{Synthetic Gaussian Dataset} with dimension $256$. Every line in Figure \ref{plot:performanceplot256} corresponds to a specific version of our code, building on top of the improvements of the lines below.
\\
\\
\textbf{turbosampling} Selection step improvements (Section \ref{subsec:selection_step})\\
\textbf{l2intrinsics} Optimized L2 distance calculation using intrinsics (Section \ref{subsec:highdimension})\\
\textbf{mem-align} Data aligned to 256bits (Section \ref{subsec:highdimension})\\
\textbf{blocked} Blocking of the L2 distance calculations (Section \ref{subsec:highdimension})\\
\textbf{greedyheuristic} Reordering of memory using heuristic (Section \ref{subsec:greedy})\\

Even for small $n$ the performance rapidly degrades as cache sizes are exceeded.
Note how all of the improvements lead to speedups, resulting in a performance gain of 1.5x. Considering the difficult setting of a randomized approximation algorithm, this is a remarkable improvement from our \textit{turbosampling} baseline.

\begin{figure}\centering
  \includegraphics[width=0.4\textwidth]{images/perfpd256.pdf}
  \caption{Performance of the discussed improvements on the \textit{Synthetic Gaussian Dataset}. The dimensions is fixed at $256$, while the $n$ increases along the x-axis. Lines correspond to specific tags of our code where every versions contains the improvements of previous versions to showcase the effect of specific changes.\label{plot:performanceplot256}}
\end{figure}

\subsection{Impact of Dimensionality}
We previously argued that increasing the dimension of our input data makes our implementation compute bound. In low dimensions, our algorithm spends a large portion of its runtime in the selection step. For higher dimensions, the focus changes to the calculation step. Some improvements can therefore only unfold their true potential when we increase the dataset dimensionality. In Figure \ref{plot:dimensionality} we kept the number of datapoints $n$ fixed at $16'384$ and varied the dimension from $8$ up to $3144$ in increments of $64$. We used our \textit{Synthetic Single Gaussian Dataset} for this evaluation. Note how \textit{turbosampling} only sees a 3.52x performance gain when increasing the dimension while our optimizations targeted at high-dimensional data profit much more from the higher dimension - \textit{blocked} sees a 8.90x speedup going from $d=8$ to $d=3144$. Those results indicate that the impact of our optimizations varies with the dimension as expected.
\begin{figure}\centering
  \includegraphics[width=0.4\textwidth]{images/dimensionality.pdf}
  \caption{Performance for the various versions of our code on the \textit{Synthetic Single Gaussian Dataset} with $n=16'384$ and increasing dimension along the x-axis.\label{plot:dimensionality}}
\end{figure}

\section{Conclusion}
NN-Descent is an efficient algorithm for computing K-NNGs.
The presented fast single-core C implementation incorporates numerous orthogonal optimizations which improve the runtime significantly in low-dimensional, as well as high dimensional usecases.\\
The presented optimizations to the selection step which reduce the number of memory passes and simplify the datastructures are particularly primarily for low dimensional data.\\
For high dimensional input data the implementation becomes compute bound and beneficial memory alignment and the use of blocked distance evaluations become paramount. Blocking is made possible by the restriction to the l2-distance while restricing the dimensionality makes memory alignment cheaper.\\
An interesting topic for future work would be to further explore heuristics for reordering the data. The evaluation shows that such heuristics can improve performance even on real world data. Such reordering heuristics are made possible by the iterative nature of the NN-Descent algorithm.\\
The significantly lower runtime on real world datasets compared to the popular PyNNDescent implementation illustrates the value of specialized implementations for common usecases.

\printbibliography

\end{document}