# Supplementary Material: Collective transition quenching in the presence of
multiple competing decay channels

Wai-Keong Mok  Institute for Quantum Information and Matter, California
Institute of Technology, Pasadena, CA 91125, USA  Stuart J. Masson  Department
of Physics, Columbia University, New York, New York 10027, USA  Dan M.
Stamper-Kurn  Department of Physics, University of California, Berkeley,
California 94720  Challenge Institute for Quantum Computation, University of
California, Berkeley, California 94720  Materials Sciences Division, Lawrence
Berkeley National Laboratory, Berkeley, California 94720  Tanya Zelevinsky
Department of Physics, Columbia University, New York, New York 10027, USA  Ana
Asenjo-Garcia  Department of Physics, Columbia University, New York, New York
10027, USA

###### Contents

  1. I Collective jump operators
  2. II Derivation of the rate equation
  3. III Effect of non-collective decay channels
  4. IV Degenerate case of identical decay rates
  5. V Fluctuations of the ground state population density
  6. VI Estimate of the convergence time to steady state
  7. VII Integrability and conservation laws of the rate equation
    1. VII.1 Dominant transition ratio
    2. VII.2 Continuum approximation
  8. VIII Molecular parameters for photoassociation

##  I Collective jump operators

Here, we provide a simple derivation for the collective jump operators
$\hat{A}_{\mu}$ corresponding to the transition $\ket{e}\to\ket{g_{\mu}}$. In
the permutation-symmetric subspace, we label the basis states as
$\ket{n_{1},n_{2},\ldots,n_{d}}$, where $n_{\mu}$ denotes the number of
emitters in the ground state $\ket{g_{\mu}}$. In terms of the product basis
kets,

$\ket{n_{1},n_{2},\ldots,n_{d}}=\begin{pmatrix}N\\\ n_{1}\,n_{2}\ldots
n_{d}\,n_{e}\end{pmatrix}^{-1/2}\left(\ket{\underbrace{g_{1}\ldots
g_{1}}_{n_{1}}\ldots\underbrace{g_{d}\ldots g_{d}}_{n_{d}}\underbrace{e\ldots
e}_{n_{e}}}+\text{permutations}\right),$ (1)

where $n_{e}=N-\sum_{\mu=1}^{d}n_{\mu}$ and $N\choose{n_{1}\ldots n_{e}}$ is
the multinomial coefficient. The action of the collective jump operators on
$\ket{n_{1},n_{2},\ldots,n_{d}}$ reads,

$\begin{split}\hat{A}_{\mu}\ket{n_{1},n_{2},\ldots,n_{\mu},\ldots
n_{d}}&=\begin{pmatrix}N\\\ n_{1}\,n_{2}\ldots
n_{d}\,n_{e}\end{pmatrix}^{-1/2}\hat{A}_{\mu}\left(\ket{\underbrace{g_{1}\ldots
g_{1}}_{n_{1}}\ldots\underbrace{g_{d}\ldots g_{d}}_{n_{d}}\underbrace{e\ldots
e}_{n_{e}}}+\text{permutations}\right)\\\
&\propto\ket{n_{1},n_{2},\ldots,n_{\mu}+1,\ldots n_{d}}.\end{split}$ (2)

The proportionality constant can be worked out via a simple counting argument.
Firstly, $\hat{A}_{\mu}$ acting on a product basis ket in Eq. (1) produces a
superposition of $n_{e}$ product basis kets. We can thus write

$\begin{split}\hat{A}_{\mu}\ket{n_{1},n_{2},\ldots,n_{\mu},\ldots
n_{d}}&=\begin{pmatrix}N\\\ n_{1}\,n_{2}\ldots
n_{d}\,n_{e}\end{pmatrix}^{-1/2}\left[\text{sum of}\quad\begin{pmatrix}N\\\
n_{1}\,n_{2}\ldots n_{d}\,n_{e}\end{pmatrix}n_{e}\quad\text{product basis
kets}\right].\end{split}$ (3)

However, this is overcounting since there are only
$N!/(n_{1}!\ldots(n_{\mu}+1)!\ldots n_{d}!(n_{e}-1)!)$ unique product basis
kets. Since the state is permutationally-symmetric, each unique product basis
ket is repeated the same number of times in the sum. Hence,

$\begin{split}\hat{A}_{\mu}\ket{n_{1},n_{2},\ldots,n_{\mu},\ldots
n_{d}}&=\begin{pmatrix}N\\\ n_{1}\,n_{2}\ldots
n_{d}\,n_{e}\end{pmatrix}^{-1/2}\times\frac{\begin{pmatrix}N\\\
n_{1}\,n_{2}\ldots n_{d}\,n_{e}\end{pmatrix}n_{e}}{\begin{pmatrix}N\\\
n_{1}\,n_{2}\ldots(n_{\mu}+1)\ldots n_{d}(n_{e}-1)\end{pmatrix}}\\\
&\times\left(\ket{\underbrace{g_{1}\ldots
g_{1}}_{n_{1}}\ldots\underbrace{g_{\mu}\ldots
g_{\mu}}_{n_{\mu}+1}\underbrace{g_{d}\ldots g_{d}}_{n_{d}}\underbrace{e\ldots
e}_{n_{e}-1}}+\text{permutations}\right).\end{split}$ (4)

Since

$\ket{n_{1},n_{2},\ldots,n_{\mu}+1,\ldots n_{d}}=\begin{pmatrix}N\\\
n_{1}\,n_{2}\ldots(n_{\mu}+1)\ldots
n_{d}(n_{e}-1)\end{pmatrix}^{-1/2}\left(\ket{\underbrace{g_{1}\ldots
g_{1}}_{n_{1}}\ldots\underbrace{g_{\mu}\ldots
g_{\mu}}_{n_{\mu}+1}\underbrace{g_{d}\ldots g_{d}}_{n_{d}}\underbrace{e\ldots
e}_{n_{e}-1}}+\text{permutations}\right),$ (5)

we have

$\begin{split}\hat{A}_{\mu}\ket{n_{1},n_{2},\ldots,n_{\mu},\ldots
n_{d}}&=\begin{pmatrix}N\\\ n_{1}\,n_{2}\ldots
n_{d}\,n_{e}\end{pmatrix}^{-1/2}\times\frac{\begin{pmatrix}N\\\
n_{1}\,n_{2}\ldots n_{d}\,n_{e}\end{pmatrix}n_{e}}{\begin{pmatrix}N\\\
n_{1}\,n_{2}\ldots(n_{\mu}+1)\ldots n_{d}\,(n_{e}-1)\end{pmatrix}}\\\
&\times\begin{pmatrix}N\\\ n_{1}\,n_{2}\ldots(n_{\mu}+1)\ldots
n_{d}\,(n_{e}-1)\end{pmatrix}^{1/2}\ket{n_{1},n_{2},\ldots,n_{\mu}+1,\ldots
n_{d}},\end{split}$ (6)

which can be simplfied as

$\hat{A}_{\mu}\ket{n_{1},n_{2},\ldots,n_{\mu},\ldots
n_{d}}=\sqrt{\left(N-\sum_{\nu=1}^{d}n_{\nu}\right)(n_{\mu}+1)}\ket{n_{1},n_{2},\ldots,n_{\mu}+1,\ldots
n_{d}}.$ (7)

Correspondingly,

$\hat{A}_{\mu}^{\dagger}\ket{n_{1},n_{2},\ldots,n_{\mu},\ldots
n_{d}}=\sqrt{\left(N-\sum_{\nu=1}^{d}n_{\nu}+1\right)n_{\mu}}\ket{n_{1},n_{2},\ldots,n_{\mu}-1,\ldots
n_{d}}.$ (8)

##  II Derivation of the rate equation

From the master equation

$\dot{\rho}=-i\left[\sum_{\mu=1}^{d}\chi_{\mu}\hat{A}_{\mu}^{\dagger}\hat{A}_{\mu},\rho\right]+\sum_{\mu=1}^{d}\Gamma_{\mu}\mathcal{D}[\hat{A}_{\mu}]\rho,$
(9)

we derive the corresponding rate equation for the populations
$P_{n_{1},\ldots,n_{d}}$. Let us expand the density matrix in the permutation
symmetric subspace as

$\rho=\sum_{\mu=1}^{d}\sum_{n_{\mu},n_{\mu}^{\prime}=0}^{N}\rho_{n_{1},\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{d}^{\prime}}\ket{n_{1},\ldots,n_{d}}\bra{n_{1}^{\prime},\ldots,n_{d}^{\prime}},$
(10)

with the basis states defined in Eq. (1). The diagonal terms correspond to the
populations

$P_{n_{1},\ldots,n_{d}}\equiv\rho_{n_{1},\ldots,n_{d}}^{n_{1},\ldots,n_{d}}$
(11)

of occupying the state $\ket{n_{1},\ldots,n_{d}}$, while the off-diagonal
terms are the coherences. Using the collective operators in Eqs. (7) and (8),
the master equation can be written as

$\begin{split}\partial_{t}\rho_{n_{1},\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{d}^{\prime}}&=-i\sum_{\mu=1}^{d}\chi_{\mu}\rho_{n_{1},\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{d}^{\prime}}\left[\braket{n_{1},\ldots,n_{d}}{\hat{A}_{\mu}^{\dagger}\hat{A}_{\mu}}{n_{1},\ldots,n_{d}}-\braket{n_{1}^{\prime},\ldots,n_{d}^{\prime}}{\hat{A}_{\mu}^{\dagger}\hat{A}_{\mu}}{n_{1}^{\prime},\ldots,n_{d}^{\prime}}\right]\\\
&+\sum_{\mu=1}^{d}\Gamma_{\mu}\Bigg{\\{}\rho_{n_{1},\ldots,n_{\mu}-1,\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{\mu}^{\prime}-1,\ldots,n_{d}^{\prime}}\braket{n_{1},\ldots,n_{d}}{\hat{A}_{\mu}}{n_{1},\ldots,n_{\mu}-1,\ldots,n_{d}}\braket{n_{1}^{\prime},\ldots,n_{\mu}^{\prime}-1,\ldots,n_{d}^{\prime}}{\hat{A}_{\mu}^{\dagger}}{n_{1}^{\prime},\ldots,n_{d}^{\prime}}\\\
&-\frac{1}{2}\rho_{n_{1},\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{d}^{\prime}}\left[\braket{n_{1},\ldots,n_{d}}{\hat{A}_{\mu}^{\dagger}\hat{A}_{\mu}}{n_{1},\ldots,n_{d}}+\braket{n_{1}^{\prime},\ldots,n_{d}^{\prime}}{\hat{A}_{\mu}^{\dagger}\hat{A}_{\mu}}{n_{1}^{\prime},\ldots,n_{d}^{\prime}}\right]\Bigg{\\}}.\end{split}$
(12)

From the definitions (7) and (8) of the collective jump operators,
$\braket{n_{1},\ldots,n_{d}}{\hat{A}_{\mu}^{\dagger}\hat{A}_{\mu}}{n_{1},\ldots,n_{d}}=(N-\sum_{\nu}n_{\nu})(n_{\mu}+1)$
and
$\braket{n_{1},\ldots,n_{d}}{\hat{A}_{\mu}}{n_{1},\ldots,n_{\mu}-1,\ldots,n_{d}}=\sqrt{(N-\sum_{\nu}n_{\nu}+1)\,n_{\mu}}$.
The master equation is thus

$\begin{split}\partial_{t}\rho_{n_{1},\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{d}^{\prime}}&=-i\sum_{\mu=1}^{d}\chi_{\mu}\left[\left(N-\sum_{\nu}n_{\nu}\right)(n_{\mu}+1)-\left(N-\sum_{\nu}n_{\nu}^{\prime}\right)(n_{\mu}^{\prime}+1)\right]\rho_{n_{1},\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{d}^{\prime}}\\\
&+\sum_{\mu=1}^{d}\Gamma_{\mu}\Bigg{\\{}\sqrt{\left(N-\sum_{\nu}n_{\nu}+1\right)n_{\mu}\left(N-\sum_{\nu}n_{\nu}^{\prime}+1\right)n_{\mu}^{\prime}}\,\,\rho_{n_{1},\ldots,n_{\mu}-1,\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{\mu}^{\prime}-1,\ldots,n_{d}^{\prime}}\\\
&-\frac{1}{2}\left[\left(N-\sum_{\nu}n_{\nu}\right)(n_{\mu}+1)+\left(N-\sum_{\nu}n_{\nu}^{\prime}\right)(n_{\mu}^{\prime}+1)\right]\rho_{n_{1},\ldots,n_{d}}^{n_{1}^{\prime},\ldots,n_{d}^{\prime}}\Bigg{\\}}.\end{split}$
(13)

Note that the populations are decoupled from the coherences, which decay to
zero at long times. Thus, for the purposes of obtaining the steady state
populations, we can ignore the dynamics of the coherences. The diagonal terms
from setting $n_{\mu}^{\prime}=n_{\mu}$ read

$\dot{P}_{n_{1},\ldots,n_{d}}=\sum_{\mu=1}^{d}\Gamma_{\mu}\left[\left(N-\sum_{\nu}n_{\nu}+1\right)n_{\mu}P_{n_{1},\ldots,n_{\mu}-1,n_{d}}-\left(N-\sum_{\nu}n_{\nu}\right)(n_{\mu}+1)P_{n_{1},\ldots,n_{d}}\right],$
(14)

which leads to the rate equation for the populations.

##  III Effect of non-collective decay channels

Non-collective decay channels can be modelled as a leakage out of the subspace
spanned by $\ket{e}$ and the $d$ ground states
$\ket{g_{1}},\ldots,\ket{g_{d}}$. We model non-collective decay channels by
introducing an additional ground state $\ket{g_{0}}$ for each emitter. The
transition $\ket{e}\to\ket{g_{0}}$ occurs at a rate of $\Gamma_{\text{leak}}$
which is not collectively enhanced. Defining the lowering operator
$\hat{\sigma}_{i}^{-}=(\ket{g_{0}}\bra{e})_{i}$, the non-collective decay
contributes a dissipative term
$\Gamma_{\text{leak}}\sum_{i}\mathcal{D}[\hat{\sigma}_{i}^{-}]\rho$ to the
master equation. We label states by $\ket{n_{0};n_{1},\ldots,n_{d}}$, where
$n_{\nu}$ is the number of emitters in the ground state $\ket{g_{\nu}}$,
$\nu=0,1,\ldots,d$. The number of emitters in the excited state is thus
$n_{e}=N-\sum_{\nu=0}^{d}n_{\nu}$. However, due to the presence of independent
decay, the states $\ket{n_{0};n_{1},\ldots,n_{d}}$ cannot be interpreted as a
symmetrized state, since the local jump operators $\hat{\sigma}_{i}$ do not
preserve permutation symmetry, although the density matrix remains symmetric
if $\Gamma_{\text{leak}}$ is the same for all $i\in\\{1,\ldots,N\\}$ [1].
Instead, we modify the notation to use $\ket{n_{0};n_{1},\ldots,n_{d}}$ as a
shorthand to denote any state containing $n_{\mu}$ emitters in the ground
state $\ket{g_{\mu}}$, with $\mu=\\{0,\ldots,d\\}$.

Consider the state $\ket{n_{0};n_{1},\ldots,n_{d}}$. The probability of a
collective transition with the operator $\hat{A}_{\mu}$ is

$p_{\mu}\propto\Gamma_{\mu}\braket{n_{0};n_{1},\ldots,n_{d}}{\hat{A}_{\mu}^{\dagger}\hat{A}_{\mu}}{n_{0};n_{1},\ldots,n_{d}}=\Gamma_{\mu}\left(N-\sum_{\nu=0}^{d}n_{\nu}\right)(n_{\mu}+1),$
(15)

and the total leakage probability is

$p_{\text{leak}}\propto\Gamma_{\text{leak}}\sum_{i=1}^{N}\braket{n_{0};n_{1},\ldots,n_{d}}{\hat{\sigma}_{i}^{+}\hat{\sigma}_{i}^{-}}{n_{0};n_{1},\ldots,n_{d}}=\Gamma_{\text{leak}}\left(N-\sum_{\nu=0}^{d}n_{\nu}\right).$
(16)

The point to note here is that the collective transition probabilities and the
total leakage probability depend only on the occupation numbers
$\\{n_{\mu}\\}$ and hold even if $\ket{n_{0};n_{1},\ldots,n_{d}}$ is not a
permutationally symmetric state. Thus, for our purposes of calculating the
DTR, the exact form of $\ket{n_{0};n_{1},\ldots,n_{d}}$ is not relevant.

In the continuum approximation, this is modelled as a fluid flow in
$\mathbb{R}_{+}^{d+1}$, with the velocity vector field
$\vec{v}(\vec{x})=(v_{0},v_{1},\ldots,v_{d})^{T}$ having components

$v_{0}(\vec{x})=\Gamma_{\text{leak}}\left(N-\sum_{\nu=0}^{d}x_{\nu}\right),$
(17)

and

$v_{\mu}(\vec{x})=\Gamma_{\mu}\left(N-\sum_{\nu=0}^{d}x_{\nu}\right)(x_{\mu}+1),\quad\mu=1,\ldots,d.$
(18)

Dividing the components of $\vec{v}(\vec{x})$, we have (for $\mu=1,\ldots,d$)

$\frac{\text{d}x_{0}}{\text{d}x_{\mu}}=\frac{\Gamma_{\text{leak}}}{\Gamma_{\mu}}\frac{1}{x_{\mu}+1},$
(19)

which can be integrated to obtain

$x_{0}=\frac{\Gamma_{\text{leak}}}{\Gamma_{\mu}}\ln(x_{\mu}+1).$ (20)

In the steady state, $x_{0}+\sum_{\mu=1}^{d}x_{\mu}=N$, so

$x_{0}+\sum_{\mu=1}^{d}\left(e^{\Gamma_{\mu}x_{0}/\Gamma_{\text{leak}}}-1\right)=N.$
(21)

Solving this for $x_{0}$ exactly is difficult. Instead, for large $N$, we use
dominant balance to get

$e^{\Gamma_{1}x_{0}/\Gamma_{\text{leak}}}\approx N\implies
x_{0}\approx\frac{\Gamma_{\text{leak}}}{\Gamma_{1}}\ln N.$ (22)

It remains to check that the omitted terms in Eq. (21) are subdominant to $N$.
The largest omitted term is

$e^{\Gamma_{2}x_{0}/\Gamma_{\text{leak}}}\approx e^{r_{2}\ln N}=N^{r_{2}}\ll
N$ (23)

for a sufficiently large $N$. This is the same regime of validity as for our
approximate DTR formula. Thus, for large $N$, the various $x_{\mu}\sim
N^{r_{\mu}}$ scale as a power law with $N$, while $x_{0}\sim\ln N$ only scales
logarithmically with $N$. We thus conclude that non-collective decay does not
affect the asymptotic behavior of our approximate DTR.

In the special case of $d=1$, the system reduces to the well-studied problem
of Dicke superradiance in two-level systems, with additional leakage that is
not collectively enhanced. There is no competition between the various
collective decay channels. As a result, the mean population density of the
ground state $\ket{g_{1}}$ scales as

$\frac{\overline{n_{1}}}{N}\approx\frac{x_{1}}{N}\sim
1-\frac{\Gamma_{\text{leak}}}{\Gamma_{1}}\frac{\ln N}{N}$ (24)

which always converges to unity, even if $\Gamma_{1}<\Gamma_{\text{leak}}$.

##  IV Degenerate case of identical decay rates

In the main text, we analyze the collective transition quenching in the non-
degenerate case, i.e, $\Gamma_{1}>\Gamma_{2}>\ldots>\Gamma_{d}$. Here, we
analyze the degenerate case where some of the $\Gamma_{\mu}$ are identical. If
all $\Gamma_{\mu}$ are the same, the steady state probability distribution
$P_{n_{1},\ldots,n_{d}}$ is the uniform distribution supported on points
$(n_{1},\ldots,n_{d})$ satisfying the constraint $n_{1}+\ldots+n_{d}=N$ [2].
More generally, we consider $d$ ground states split into $L$ degenerate
blocks, with degeneracies $\\{m_{1},m_{2},\ldots,m_{L}\\}$ and
$\sum_{\mu=1}^{L}m_{\mu}=d$. The decay rates are ordered as
$\Gamma_{1}>\Gamma_{2}>\ldots>\Gamma_{L}$ with the corresponding ratios
$1>r_{2}>\ldots>r_{L}$. Note that the degeneracy here refers to identical
decay rates (e.g., $\Gamma_{1}=\Gamma_{2}=\ldots=\Gamma_{m_{1}}$) and not
identical ground state energies.

Following the continuum and single-particle approximation explained in the
main text, we have

$x_{\mu}^{\Gamma_{\nu}}\approx x_{\nu}^{\Gamma_{\mu}},$ (25)

where $x_{\mu}$ is interpreted as the population of any of $m_{\mu}$ ground
states in the $\mu$-th degenerate block. Then, $x_{\mu}\approx
x_{L}^{r_{\mu}/r_{L}}$. The steady state condition is written as

$\sum_{\mu=1}^{L}m_{\mu}x_{\mu}=N\implies\sum_{\mu=1}^{L}m_{\mu}x_{L}^{r_{\mu}/r_{L}}\approx
N.$ (26)

Using dominant balance, we get

$x_{L}\approx\left(\frac{N}{m_{1}}\right)^{r_{L}}.$ (27)

Equation (25) then gives the other ground state populations as

$x_{\mu}\approx\left(\frac{N}{m_{1}}\right)^{r_{\mu}},\quad\mu=2,\ldots,L.$
(28)

Plugging this into the steady state condition,

$m_{1}x_{1}+\sum_{\mu=2}^{L}m_{\mu}\left(\frac{N}{m_{1}}\right)^{r_{\mu}}\approx
N.$ (29)

Defining the dominant transition ratio (DTR) as $x_{1}/N$, we obtain

$\text{DTR}\approx\frac{1}{m_{1}}-\sum_{\mu=2}^{L}\frac{m_{\mu}}{m_{1}^{1+r_{\mu}}N^{1-r_{\mu}}},$
(30)

which recovers the result in the main text when $m_{1}=m_{2}=\ldots=m_{L}=1$.
The factor of $1/m_{1}$ in the first term reflects the fact that as
$N\to\infty$, the population density is shared equally between the $m_{1}$
dominant ground states with the largest decay rate $\Gamma_{1}$.

##  V Fluctuations of the ground state population density

We numerically simulate the rate equation for $d=2$ ground states. From the
steady state probability distribution of the dominant ground state $P(n_{1})$,
we find the mean and variance of the distribution, to compute the relative
fluctuation $\delta n_{1}/\overline{n_{1}}$, where $\delta n_{1}$ is the
standard deviation for the population of $n_{1}$. We then fit the relative
fluctuation to the phenomenological power law $\sim N^{-\beta}$ for $N$ up to
$10^{4}$, to extract the fluctuation power-law index $\beta$. The fit is done
for small values of $r_{2}=\Gamma_{2}/\Gamma_{1}$, since larger values of
$r_{2}$ require a larger $N$ to accurately capture the asymptotic behavior.
From our fit, we estimate $\beta=0.930-0.842r_{2}$. In the main text, we show
that the variance $\text{Var}(n_{1}/N)\leq(1-\text{DTR})\text{DTR}$. Thus,

$\frac{\delta
n_{1}}{\overline{n_{1}}}=\frac{\sqrt{\text{Var}(n_{1}/N)}}{\text{DTR}}\leq\sqrt{\frac{1-\text{DTR}}{\text{DTR}}}\sim
N^{(-1+r_{2})/2}$ (31)

for $d=2$. Comparing this with the asymptotic behavior $\delta
n_{1}/\overline{n_{1}}\sim N^{-\beta}$, this gives the lower bound

$\beta\geq\frac{1}{2}(1-r_{2}).$ (32)

Our estimated $\beta$ from the fitting satisfies this bound for all
$r_{2}\in(0,1)$.

(a)  Figure 1: Power law index $\beta$ from fitting the tail of the relative
fluctuation $\delta n_{1}/\overline{n_{1}}$ of the dominant ground state
population to $\sim N^{-\beta}$. The dashed line indicates the
phenomenological fit $\beta=0.930-0.842r_{2}$. The inset shows the population
distribution of the dominant ground state for $r_{2}=0.2$, plotted between
$0.994\leq n_{1}/N\leq 1$, with a width of order $N^{-\beta}$.

##  VI Estimate of the convergence time to steady state

From the continuum approximation, we have the velocity field
$\vec{v}(\vec{x})$ with components

$v_{\mu}\equiv\frac{\text{d}x_{\mu}}{\text{d}t}=\left(N-\sum_{\mu=1}^{d}x_{\mu}\right)(x_{\mu}+1).$
(33)

In the single-particle approximation, the trajectory $\mathcal{C}$ in
$\mathbb{R}_{+}^{d}$ is defined by the set of equations

$(x_{\mu}+1)^{\Gamma_{\nu}}=(x_{\nu}+1)^{\Gamma_{\mu}}$ (34)

for all pairs $\mu,\nu\in\\{1,\ldots,d\\}$, in the region
$\sum_{\mu=1}^{d}x_{\mu}\leq N$. The convergence time can be written as a line
integral

$T=\int_{\mathcal{C}}\frac{\text{d}x_{1}}{v_{1}}=\frac{1}{\Gamma_{1}}\int_{\mathcal{C}}\frac{\text{d}x_{1}}{(N-x_{1}-\sum_{\mu>1}x_{\mu})(x_{1}+1)}$
(35)

from $x_{1}=0$ to $x_{1}=N-\sum_{\mu>1}N^{r_{\mu}}$. As in the main text, we
consider the non-degenerate case of $r_{\mu}$ being all different for
simplicity. Performing the integral directly is not trivial. Instead, we can
obtain an approximation by splitting the trajectory $\mathcal{C}$ into two
parts [3]: $\mathcal{C}_{1}$, where the particle travels along $\mathcal{C}$
from $x_{1}=0$ to $x_{1}=N^{a}$ for some $0<a<1$, and $\mathcal{C}_{2}$, where
the particle travels along $\mathcal{C}$ from $x_{1}=N^{a}$ to
$x_{1}=N-\sum_{\mu>1}N^{r_{\mu}}$. There is no loss of generality in setting
the intermediate point to be at $x_{1}=N^{a}$, since we can always choose a
point on $\mathcal{C}$ with $1\ll x_{1}\ll N$ and define $a$ appropriately.

For the trajectory $\mathcal{C}_{1}$, $N\gg\\{x_{1},x_{\mu>1}\\}$, so

$T_{1}\equiv\int_{\mathcal{C}_{1}}\frac{\text{d}x_{1}}{v_{1}}\approx\frac{1}{N\Gamma_{1}}\int_{0}^{N^{a}}\frac{dx_{1}}{x_{1}+1}\approx\frac{a}{\Gamma_{1}}\frac{\ln
N}{N}.$ (36)

For the trajectory $\mathcal{C}_{2}$, $x_{1}\gg x_{\mu>1}$, so

$\begin{split}T_{2}\equiv\int_{\mathcal{C}_{2}}\frac{\text{d}x_{1}}{v_{1}}&\approx\frac{1}{\Gamma_{1}}\int_{N^{a}}^{N-\sum_{\mu>1}N^{r_{\mu}}}\frac{\text{d}x_{1}}{(N-x_{1})(x_{1}+1)}\approx\frac{1}{N\Gamma_{1}}\ln\left(\frac{x_{1}+1}{N-x_{1}}\right)\bigg{|}_{N^{a}}^{N-\sum_{\mu>1}N^{r_{\mu}}}\\\
&\approx\frac{1}{N\Gamma_{1}}\ln\left(\frac{N^{2}}{N^{r_{2}+a}}\right)=\frac{2-r_{2}-a}{\Gamma_{1}}\frac{\ln
N}{N}.\end{split}$ (37)

This gives the convergence time $T$ as

$T=T_{1}+T_{2}\approx\frac{2-r_{2}}{\Gamma_{1}}\frac{\ln N}{N}.$ (38)

The final result for $T$ does not depend on $a$, which is consistent since the
intermediate point is chosen arbitrarily. Thus, adding competing collective
decay channels does not affect the convergence timescale of $\sim\ln N/N$, but
only changes the prefactor.

##  VII Integrability and conservation laws of the rate equation

Since the steady state is non-unique and depends on the initial conditions,
the rate equation encodes conservation laws that we now derive. Our goal is to
find a function $h(n_{1},\ldots,n_{d})$ that satisfies the conservation law
$\partial_{t}\overline{h(n_{1},\ldots,n_{d})}\equiv\sum_{n_{1},\ldots,n_{d}}h(n_{1},\ldots,n_{d})\dot{P}_{n_{1},\ldots,n_{d}}=0$,
where $\dot{P}_{n_{1},\ldots,n_{2}}$ is defined by the rate equation.
Explicitly,

$\begin{split}0&=\sum_{n_{1},\ldots,n_{d}}\sum_{\mu=1}^{d}\Gamma_{\mu}\left[-\left(N-\sum_{\nu}n_{\nu}\right)\left(n_{\mu}+1\right)P_{n_{1},\ldots,n_{\mu},\ldots,n_{d}}+\left(N-\sum_{\nu}n_{\nu}+1\right)n_{\mu}P_{n_{1},\ldots,n_{\mu}-1,\ldots,n_{d}}\right]h(n_{1},\ldots,n_{d})\\\
&=\sum_{n_{1},\ldots,n_{d}}\sum_{\mu=1}^{d}\Gamma_{\mu}\left(N-\sum_{\nu}n_{\nu}\right)(n_{\mu}+1)\left[h(n_{1},\ldots,n_{\mu}+1,\ldots,n_{d})-h(n_{1},\ldots,n_{\mu},\ldots,n_{d})\right]P_{n_{1},\ldots,n_{d}}.\end{split}$
(39)

Subsituting the separable ansatz

$h(n_{1},\ldots,n_{d})=\prod_{\mu=1}^{d}h_{\mu}(n_{\mu})$ (40)

into the conservation law, and dividing through by $h(n_{1},\ldots,n_{d})$, we
have

$\sum_{\mu}\Gamma_{\mu}(n_{\mu}+1)\left(\frac{h_{\mu}(n_{\mu}+1)}{h_{\mu}(n_{\mu})}-1\right)=0,$
(41)

using the fact that the probabilities $P_{n_{1},\ldots,n_{d}}$ are arbitrary.
Since each term in the sum above only depends on $n_{\mu}$ that are arbitrary,
we must have

$\Gamma_{\mu}(n_{\mu}+1)\left(\frac{h_{\mu}(n_{\mu}+1)}{h_{\mu}(n_{\mu})}-1\right)=c_{\mu}$
(42)

where $c_{\mu}$ are constants satisfying the constraint

$\sum_{\mu=1}^{d}c_{\mu}=0.$ (43)

Equation (42) is a recurrence relation for $h_{\mu}(n_{\mu})$, which has the
general solution

$h_{\mu}(n_{\mu})\propto\frac{\tilde{\Gamma}(1+n_{\mu}+c_{\mu}/\Gamma_{\mu})}{n_{\mu}!},$
(44)

where $\tilde{\Gamma}(z)$ is the gamma function. Multiplying the $h_{\mu}$’s
together and choosing the normalization factor such that $\overline{h}=1$ for
the fully excited state (described by the probability distribution
$P_{n_{1},\ldots,n_{d}}=\prod_{\mu}\delta_{n_{\mu},0}$), we obtain the
generating function

$h(n_{1},\ldots,n_{d})=\prod_{\mu=1}^{d}\frac{\tilde{\Gamma}(1+n_{\mu}+c_{\mu}/\Gamma_{\mu})}{n_{\mu}!\tilde{\Gamma}(1+c_{\mu}/\Gamma_{\mu})}.$
(45)

By construction, the steady state distribution satisfies $\overline{h}=1$ with
$\sum_{\mu}n_{\mu}=N$. The generating function $h(n_{1},\ldots,n_{d})$ is a
polynomial of degree $n_{\mu}$ in the variable $c_{\mu}$, which can be seen by
expanding out Eq. (45) using the definition
$\tilde{\Gamma}(1+z)=z\tilde{\Gamma}(z)$ of the gamma function iteratively,

$h(n_{1},\ldots,n_{d})=\prod_{\mu=1}^{d}\frac{1}{n_{\mu}!}\prod_{j=1}^{n_{\mu}}\left(j+\frac{c_{\mu}}{\Gamma_{\mu}}\right).$
(46)

This expression holds for $n_{\mu}>0$. If $n_{\mu}=0$, the product
$\prod_{j}(j+c_{\mu}/\Gamma_{\mu})$ is replaced by $1$. Thus,
$h(n_{1},\ldots,n_{d})$ encodes $N+1$ independent conserved quantities, which
are the coefficients of the polynomial. A trivial conserved quantity can be
obtained by setting all $c_{\mu}=0$,

$1=\overline{h}=\sum_{n_{1},\ldots,n_{d}}P_{n_{1},\ldots,n_{d}},$ (47)

which is simply the conservation of probability. By choosing suitable values
of $c_{\mu}$, we seek non-trivial conserved quantities.

###  VII.1 Dominant transition ratio

Let us choose $c_{\mu}=\Gamma_{\mu}$ for some particular $\mu>1$, and all
other $c_{\nu>1}=0$. From the constraint $\sum_{k=1}^{d}c_{k}=0$, this fixes
$c_{1}=-\Gamma_{\mu}$. Substituting this into the conservation law
$\overline{h}=1$, we obtain

$\bigg{\langle}\frac{\tilde{\Gamma}(1+n_{1}-r_{\mu})\tilde{\Gamma}(2+n_{\mu})}{n_{1}!n_{\mu}!}\bigg{\rangle}=\tilde{\Gamma}(1-r_{\mu}).$
(48)

The average on the LHS is taken with respect to the steady state probability
distribution. We denote $r_{\mu}=\Gamma_{\mu}/\Gamma_{1}<1$. Assuming that the
main contribution from $n_{\mu}$ comes from $1\ll n_{\mu}\ll N$, we use
Stirling’s approximation to get

$\overline{n_{1}^{-r_{\mu}}n_{\mu}}\approx\tilde{\Gamma}(1-r_{\mu}).$ (49)

Since $r_{\mu}<1$, the gamma function is always finite and positive.
Substituting the steady state constraint $n_{1}=N-\sum_{\nu>1}n_{\nu}$, and
Taylor expanding in powers of $1/N$, we have

$\overline{n_{\mu}}+\frac{r_{\mu}}{N}\sum_{\nu>1}\overline{n_{\mu}n_{\nu}}\approx\tilde{\Gamma}(1-r_{\mu})N^{r_{\mu}}.$
(50)

To proceed further, we assume $\sum_{\nu>1}\overline{n_{\mu}n_{\nu}}\ll
N\overline{n_{\mu}}$ to finally arrive at

$\overline{n_{\mu}}\approx\tilde{\Gamma}(1-r_{\mu})N^{r_{\mu}}$ (51)

which gives the DTR as

$\text{DTR}=\frac{\overline{n_{1}}}{N}=1-\frac{1}{N}\sum_{\mu>1}\overline{n_{\mu}}=1-\sum_{\mu>1}\frac{\tilde{\Gamma}(1-r_{\mu})}{N^{1-r_{\mu}}}.$
(52)

The assumption $\sum_{\nu>1}\overline{n_{\mu}n_{\nu}}\ll N\overline{n_{\mu}}$
can be thought of as a mean-field approximation
$\overline{n_{\mu}n_{\nu}}\approx\overline{n_{\mu}}\,\,\overline{n_{\nu}}$.
Then, using Eq. (51), we have
$\sum_{\nu>1}\overline{n_{\mu}n_{\nu}}\approx\overline{n_{\mu}}\sum_{\nu>1}\overline{n_{\nu}}\approx\overline{n_{\mu}}\sum_{\nu>1}N^{r_{\nu}}\ll
N\overline{n_{\mu}}$ which is self-consistent. Since the population density
$n_{1}/N\in[0,1]$ is bounded, we can use the Bhatia-Davis inequality [4] to
bound the variance by

$\text{Var}(n_{1}/N)\leq(1-\text{DTR})\text{DTR}\sim\sum_{\mu>1}\frac{\tilde{\Gamma}(1-r_{\mu})}{N^{1-r_{\mu}}}.$
(53)

Although Eq. (51) is finite for any $r_{\mu}<1$, one has to be careful about
the divergence of the gamma function as $r_{\mu}\to 1^{-}$. Since
$\tilde{\Gamma}(z)\sim 1/z$ near $z=0$, we get

$\overline{n_{\mu}}\approx\frac{N^{r_{\mu}}}{1-r_{\mu}},\quad r_{\mu}\to
1^{-}.$ (54)

To observe collective transition quenching, we demand $\overline{n_{\mu}}\ll
N$ which gives the condition

$N\gg y^{y},\quad y=\frac{1}{1-r_{\mu}}.$ (55)

Physically, this means that for $r_{\mu}$ close to $1$, we require a larger
system size to suppress the sub-dominant collective transitions. As an
estimate, if $r_{\mu}=0.9$, then our theoretical predictions hold for $N\gg
10^{10}$.

###  VII.2 Continuum approximation

We can also derive the conservation laws in the continuum approximation.
Applying the continuum approximation to the difference equation (39), we have

$\sum_{\mu=1}^{d}\Gamma_{\mu}(x_{\mu}+1)\frac{\partial}{\partial
x_{\mu}}h(\vec{x})=0.$ (56)

Using the separable ansatz

$h(\vec{x})=\prod_{\mu=1}^{d}h_{\mu}(x_{\mu}),$ (57)

we obtain

$\sum_{\mu=1}^{d}\Gamma_{\mu}(x_{\mu}+1)\frac{h_{\mu}^{\prime}(x_{\mu})}{h_{\mu}(x_{\mu})}=0.$
(58)

Each term in the sum depends on an independent coordinate $x_{\mu}$, so the
terms in the sum must all be equal to some constant $c_{\mu}$

$\Gamma_{\mu}(x_{\mu}+1)\frac{h_{\mu}^{\prime}(x_{\mu})}{h_{\mu}(x_{\mu})}=c_{\mu}$
(59)

with the constraint $\sum_{\mu=1}^{d}c_{\mu}=0$. Solving these equations, we
get the generating function

$h(\vec{x})=\prod_{\mu=1}^{d}(x_{\mu}+1)^{c_{\mu}/\Gamma_{\mu}}.$ (60)

The fully excited initial state is represented by the Dirac delta distribution
$P(\vec{x})=\delta^{d}(\vec{x})$, which satisfies $\overline{h(\vec{x})}=1$.
Equation (60) is the continuous analog of Eq. (45). Similar to the discrete
case, we can choose $c_{\mu}=\Gamma_{\mu}$ for some $\mu>1$ and all other
$c_{\nu>1}=0$, which constraints $c_{1}=-\Gamma_{\mu}$. Assuming that the main
contribution comes from $1\ll x_{\mu}\ll N$ and using
$x_{1}=N-\sum_{\mu>1}x_{\mu}$, we have

$\overline{x_{\mu}}\approx N^{r_{\mu}},$ (61)

which gives the approximate DTR as

$\text{DTR}\approx 1-\sum_{\mu>1}\frac{1}{N^{1-r_{\mu}}}.$ (62)

This agrees with the approximate formula derived in the main text using the
single-particle approximation.

##  VIII Molecular parameters for photoassociation

Vibrational quantum number $\nu$ | Rotational quantum number $J$ | Branching factor | Binding energy, MHz$/h$
---|---|---|---
-1 | 0 | 0.54 | 137
-1 | 2 | 0.27 | 67
-2 | 0 | 0.13 | 1400
-2 | 2 | 0.06 | 1227
Table 1: Predominant decay paths from the excited
$(1)1_{u}(\nu^{\prime}=-1,J^{\prime}=1)$ state of a strontium dimer [5]. There
are four states with significant branching ratios, those with rotational
quantum number $J=\\{0,2\\}$ and vibrational quantum number $\nu=\\{-1,-2\\}$
(where negative vibrational numbers count down from the dissociation
threshold) [6]. All other states have branching ratios of $<0.01$ and are
ignored in the calculations. The minimum separation between these states is
for the two $\nu=-1$ states, where the separation between the $J=0$ and $J=2$
states is $~{}70$ MHz [7], significantly larger than the linewidth of the
excited state which is $\sim 15$ kHz [8].

## References

  * Shammah  _et al._ [2018] N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori, Open quantum systems with local and collective incoherent processes: Efficient numerical simulations using permutational invariance, Phys. Rev. A 98, 063815 (2018).
  * Piñeiro Orioli  _et al._ [2022] A. Piñeiro Orioli, J. K. Thompson, and A. M. Rey, Emergent dark states from superradiant dynamics in multilevel atoms in a cavity, Phys. Rev. X 12, 011054 (2022).
  * Malz  _et al._ [2022] D. Malz, R. Trivedi, and J. I. Cirac, Large-$n$ limit of Dicke superradiance, Phys. Rev. A 106, 013716 (2022).
  * Bhatia and Davis [2000] R. Bhatia and C. Davis, A better bound on the variance, Am. Math. Mon. 107, 353 (2000).
  * [5] M. Borkowski, private communication, .
  * [6] I. Majewska, PhD Thesis, 2021, University of Warsaw. Theoretical description of ultracold strontium molecules in an optical lattice: Control of photodissociation and interpretation of molecular clock experiments, .
  * [7] K. H. Leung, PhD Thesis, 2023, Columbia University. The strontium molecular lattice clock: Vibrational spectroscopy with hertz-level accuracy, .
  * Zelevinsky  _et al._ [2006] T. Zelevinsky, M. M. Boyd, A. D. Ludlow, T. Ido, J. Ye, R. Ciuryło, P. Naidon, and P. S. Julienne, Narrow line photoassociation in an optical lattice, Phys. Rev. Lett. 96, 203201 (2006).


# Generalized Neighbor Search using Commodity Hardware Acceleration

Durga Mandarapu  nnnn-nnnn-nnnn-nnnn Department of Computer SciencePurdue
UniversityStreet1 Address1West LafayetteIN47906USA<EMAIL_ADDRESS>,  Vani
Nagarajan  School of Electrical and Computer EngineeringPurdue UniversityWest
LafayetteINUSA<EMAIL_ADDRESS>and  Milind Kulkarni  School of Electrical
and Computer EngineeringPurdue UniversityWest LafayetteINUSA<EMAIL_ADDRESS>
###### Abstract.

Tree-based Nearest Neighbor Search (NNS) is hard to parallelize on GPUs.
However, newer Nvidia GPUs are equipped with Ray Tracing (RT) cores that can
build a spatial tree called Bounding Volume Hierarchy (BVH) to accelerate
graphics rendering. Recent work proposed using RT cores to implement NNS, but
they all have a hardware-imposed constraint on the type of distance metric,
which is the Euclidean distance. We propose and implement two approaches for
generalized distance computations: filter-refine, and monotone transformation,
each of which allows non-euclidean nearest neighbor queries to be performed in
terms of Euclidean distances. We find that our reductions improve the time
taken to perform distance computations during the search, thereby improving
the overall performance of the NNS.

##  1\. Introduction

Nearest Neighbor Search (NNS) is the problem of finding points similar to a
query point in terms of a desired distance metric. The most popular metric is
the Euclidean distance ($L^{2}$ norm), which is assumed by many high-
performance realizations of the NNS problem [1; 2].

NNS with non-Eucledean metrics has various applications.

  * •

Celestial objects, such as stars, are positioned using 3D spherical or
Cartesian coordinates. The metric used to measure the proximity of such
neighborhood stars to Earth is the Angular distance, determined by the stellar
parallax method (parallax, ).

  * •

In a spatial data application like street maps, finding the nearest points of
interest involves ordering them according to their Manhattan distance from the
query point location (manhattan, ) since the data points in a city usually
adhere to taxicab geometry.

An incomplete list of popular non-Euclidean metrics include Manhattan distance
($L^{1}$ norm), Chebyshev distance ($L^{\infty}$ norm), Minkowski distance
($L^{p}$ norm), angular or cosine distance, Jacquard distance, and, Hamming
distance.

Besides the choice of the metric, another important detail for a high-
performance NNS problem implementation is dimensionality. NNS in lower
dimensions has several applications, such as point cloud registration [5],
spatial data query processing, facial recognition, low dimensional embedding
[6], etc. Low-dimensional embeddings are used in various machine learning
applications. One way of producing them is through dimension reduction
followed by NNS, which efficiently extracts the prominent features (embedding,
).

Depending on dimensionality, efficient implementations of NNS employ
specialized indexing techniques. At lower dimensions, traditional approaches
to fast indexing utilize tree-based data structures. Better performance than
trees at higher dimensions [7] can be achieved by other approaches to fast
indexing, like hashing, quantization, or graphs.

Non-tree based indexing can be efficiently implemented on GPUs [8; 9; 7],
while tree-based approaches do not enjoy the same property. In particular,
GPU-accelerated tree-based implementations suffer from memory coalescing and
thread divergence [10]. Prior work for making tree traversals faster on GPUs
[10; 11] is not satisfactory, because it makes programmer-end optimization
attempts and does not bring tree-based approaches on par with the state-of-
the-art in the non-tree based implementations.

A promising way to accelerate tree traversals on GPU employs ray-tracing cores
(RT cores), as found in recent Nvidia GPUs marked with RTX. RT cores offer
hardware support for specific visualization tasks but can be repurposed for
more generic tree problems. With the right choice of scene and rays, RT cores
can be utilized to perform non-ray tracing applications, such as locating a
point in an unstructured tetrahedral mesh [12], drawing a force-directed graph
[13], and performing nearest neighbor searches [1; 2].

The key idea in all generic applications of RT cores is reducing the target
application to ray tracing. To achieve such reduction, one produces a scene
with objects and rays that encode the input to the given generic application.
The main operation in RT cores, ray-bounding-box intersection, determines
whether a ray of a given length with a given origin and direction can
intersect with any of the objects in a scene. A successful reduction ensures
that hardware-accelerated ray-bounding-box intersection built by the RT cores
provides a partial or complete solution to the generic problem of interest.

Crucially, the reductions of NNS to RT cores devised so far [1; 2] are
inherently limited to $L^{2}$ distance. This is due to RT cores themselves,
which arrange the objects in a scene following the Euclidean distance. As a
result, one cannot merely use $L^{2}$ NNS implementation to work with other
distance metrics. In particular, an object being at a particular Euclidean
distance from the point of interest says nothing about a non-$L^{2}$ distance
metric like the angular distance between them (Figure 1).

In this paper, we show how to extend reductions of nearest neighbour search
computation to RT cores to work in the case of non-Ecludean metrics and make
the following

### Contributions:

  1. (1)

We define two approaches for generalized distance computations: filter and
refine, and monotone transformations, each allowing non-Euclidean nearest
neighbor queries to be performed in terms of $L^{2}$ distances.

  2. (2)

We show how to implement these generalized distance computations on RT
hardware, including specifically optimized realizations for several different
metrics: $L^{1}$ norm, $L^{\infty}$ norm, general $L^{p}$ norms, and cosine
similarity.

  3. (3)

We evaluate our reductions by implementing them as stand-alone applications
and find that our optimized approaches are more efficient and scalable than
various baselines.

A notable limitation of our approach to efficient NNS implementation comes
from the tree-oriented nature of the computations that RT cores provide. This
makes our approach worthwhile in lower dimensions only.

Figure 1. Euclidean and Angular distances: $a$ and $b$ are data points, $q$ is
a query point, and $O$ is the point of reference.
$L^{2}(q,a)<L^{2}(q,b)\centernot\implies\beta<\alpha$

##  2\. Overview

#### Generalized nearest neighbor search

finds the $k$ nearest points from a training set $P$ to a query point $q$
according to some distance measure $D$. Designing a generalized tree-based NNS
solution involves supporting a broader set of parameters along two axes: the
dimension of the data and the distance metric. As a first step towards a truly
general hardware-accelerated NNS, we provide and accelerate generalized tree-
based NNS algorithm 2 for popular distances besides Euclidean distance,
despite the RT cores’ inherent constraint of using Euclidean distances.

Our key observation is that the $L^{2}$ reduction from prior work [1; 2] does
not directly solve the nearest neighbor problem. Instead, the reduction to RT
cores accelerates an $r$-bounded distance query: find all points within
distance $r$ of query point $q$, according to their Euclidean distance. So the
question is not how to solve NNS problems using other distance measures but
instead how to reduce the nearest neighbor searches for other distance
measures to the $r$-bounded Euclidean distance query. This paper introduces
two more-general reductions for this task, thus enabling acceleration of NNS
for other distance measures via RT cores: Filter and refine, and Monotone
transformations.

#### Filter-Refine

is a traditional two-step search process.

We identify FR as a general strategy for NNS for which the reductions from
prior papers are a special case. By performing the $r$-bounded distance query
in $L^{2}$, all the points within a radius $r$ of the query point are
collected and are then filtered and sorted to find the closest $k$ points. We
identify that many distance measures have an inclusion property with respect
to $L^{2}$: there exists some function $f:\mathbb{R}\rightarrow\mathbb{R}$
such that all the points that are within $r$ of the query point according to
the target distance measure are guaranteed to be within $f(r)$ of the query
point in $L^{2}$. We can thus adapt the filter-refine process to filter
according to $f(r)$ and then refine following the desired target distance
measure.

We further identify several norms whose geometric properties allow us to
create ray tracing scenes that implement more accurate filters, allowing for a
more precise initial filter of the points that are within $r$ of the query
point. Section 4.1 discusses this general approach in more detail and defines
specific reductions for $L^{1}$, $L^{\infty}$, and generalized $L^{p}$ norms.

#### Monotone Transformations

Not all distance metrics satisfy the inclusion property mentioned above. For
certain distance metrics, such as cosine similarity or angular distance, we
show that we can construct a monotonic spatial transformation
$f_{m}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ such that the distances
between the input points and the query point $q$ according to the desired
metric are in the same order as the distances between the output points and
$q$ according to $L^{2}$ norm. This order-preserving transformation of NNS in
terms of desired metric to NNS in $L^{2}$, allows us to use the filter-refine
approach.

Section 4.6 discusses this general approach in more detail and defines
specific reductions for cosine similarity and 2D distance. We demonstrate the
power of all of our reductions by implementing them as stand-alone
applications using RT cores of the NVIDIA GeForce RTX 2060 GPU.

##  3\. Background

###  3.1. $k$-Nearest Neighbour Search

There are several flavors of NNS, so we start by defining the NNS problem Def
3.1. The naive way of performing $k$-NNS is computing the distance between
every pair of points and ordering them, which takes $O(kn)$ or $O(nlog(k))$
per query point, which means it scales up linearly with the number of queries.
It’s only natural that fast and scalable $k$-NNS algorithms focus on finding
approximate nearest neighbors due to the time-consuming computation. This
inexact search of approximate nearest neighbors produces a set of neighbors
$T^{\prime}$ which contains a correctly found subset of true nearest neighbors
and a small subset of false nearest neighbors. Recall 1 is a performance
metric that defines the search quality and correctness of an approximate
$k$-NNS solution. It is defined as the ratio of the number of correctly found
nearest neighbors by the inexact search to the number of true nearest
neighbors from the exact search.

###### Definition 3.0.

Generalized nearest neighbor search. Given a query point $q\in\mathbb{R}^{d}$
(a point in $d$ dimensions), a training set of points,
$A\subseteq\mathcal{A}(\mathbb{R}^{d})$, $n=|A|$, a value $k\in\mathbb{N}$,
and a distance measure $D$, a function
$\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ that computes a
distance between two $d$-dimensional points, the generalized nearest neighbor
problem finds a result set of points, $T\subseteq A$, that contains the
closest $k$ points to $q$ according to $D$.

(1) $T^{\prime}=T_{1}\cup T_{2},T_{1}\subset T,T_{2}\not\subset T.\
Recall=\frac{|T^{\prime}\cap T|}{|T|}=\frac{|T_{1}|}{|T|}$

####  $L^{p}$ norm

A set of distance functions that are defined for points in 3-dimensional space
as

$L^{p}(a,b)={(|a_{x}-b_{x}|^{p}+|a_{y}-b_{y}|^{p}+|a_{z}-b_{z}|^{p})}^{\frac{1}{p}},p\in\mathbb{R}\geq
1$ $L^{\infty}(a,b)=max(|a_{x}-b_{x}|,|a_{y}-b_{y}|,|a_{z}-b_{z}|)$

where $|.|$ represents the absolute value. When $p=1$, the norm is called
Taxicab or Manhattan distance which corresponds to the distance travelled
between points along the lines in a rectangular grid. When $p=2$, it is the
well-known euclidean distance. And $p\to\infty$, the norm value reaches the
maximum norm which is the maximum of the coordinates. In general, $L^{p}$
norms have several applications in statistical computations, penalized
regressions, information theory, etc. For example, in path finding algorithms
where movement is restricted to certain directions, Manhattan distance helps
useful.

#### Cosine distance

Angular distance, cosine distance, inner product, or dot product convey
information about distances in spherical geometry. Angular distance is the
shorter angle between two vectors, while cosine distance or cosine similarity
is the cosine function applied to this angle. The inner product or dot product
is the same for vectors, and these are in turn same as cosine distance when
the vectors are of unit length. Because cosine distance measures how similar
two vectors are, it is highly useful in recommendation systems.

###  3.2. Related work on NNS

Tree-based NNS algorithms are only efficient at lower dimensions, so they are
mostly specialized for certain applications. An optimized KD-tree [5] is
proposed to aid in point cloud registration, however, this optimization is
particular to the registration pipeline. For a generic tree traversal,
Treelogy(autoropes, ) proposes several optimizations to reduce the memory and
thread divergence caused by the threads in GPU that are traversing different
parts of tree. They simply port these optimizations to perform NNS using
$L^{2}$ norm. The state-of-the-art exact NNS library [7], FAISS[8] uses two-
level quantization to index. They further use Tensorflow to provide several
optimizations and as a framework for their GPU implementation. The state-of-
the-art approximate NNS library [7], SONG[9] uses HNSW[14] graph indexing and
accelerating the distance computations by efficiently using GPU memory
allocations. cite several trees pointnet sweetknn tigris

###  3.3. Ray Tracing Cores

Ray tracing is a graphics rendering algorithm where several rays are modeled
from the data point as a source and followed till they hit/miss various
portions of the object in the scene to determine the way the object should be
presented. The object’s volume is bounded by the axis-aligned boxes to form a
hierarchical tree of such bounding volumes, with leaf nodes being boxes that
bound a single object and internal nodes being boxes that engulf more than one
bounding box. Since the goal of the RT cores is to accelerate ray-triangle
intersection tests of the ray-tracing pipeline, they need the object to be
expressed in terms of 3D triangles. For example, a tetrahedron is a volume
engulfed by 4 triangular faces. RT cores also support user geometries by
building axis-aligned bounding boxes (AABB) around them, with a downside being
that boxes may engulf more volume than just the user geometry. Though the RT
cores do not accelerate the intersection test algorithm, they help reduce the
number of tests performed. When a ray traverses the BVH, it tests for
intersection with bounding volumes that enclose multiple objects rather than
individual objects. bounding box choices

Figure 2. RTX Pipeline for user geometry (top) and triangles (bottom)

###  3.4. Programming and Execution Model

Optix provides an interface with both the traditional Streaming Multi-
processors (SMs) and RT cores, allowing the user to write traditional shader
programs that are executed on the SMs and leverage the RT hardware for BVH
build, traversal, and, if applicable, intersection testing. It supports
triangles, curves, and user-defined geometries. Figure 2 shows the stages of
the pipeline that we use in our implementations. The more general set of
programs for different stages of the ray tracing pipeline that Optix allows
are as follows:

RayGen:

The user defines ray parameters such as origin, direction and length of the
ray. RayGen creates the ray and proceeds to call for BVH traversal and
intersection testing.

Intersection:

If the BVH is constructed using user-defined geometries, the user is required
to provide a custom intersection test for ray-primitive intersections.

AnyHit:

The AnyHit program allows the user to record information about the
intersection and choose whether to continue or terminate the BVH traversal.

ClosestHit:

The user can specify a ClosestHit program to determine the primitive closest
to the ray origin (i.e) it finds the closest intersection between a primitive
and the ray.

Miss:

After the ray completes BVH traversal, if there are no intersections between
the ray and primitives in the scene, the user can specify a miss program to
handle shading in such scenarios.

###  3.5. Related Work on RT-Accelerated Applications

Wald et al.[12] were the first to use RT cores for non-ray tracing
applications. They modeled a query point as a ray to locate the tetrahedron
the point is in amongst all the tetrahedrons in the mesh. Zellman et al.[13]
showed how to use RT cores to perform graph drawing, where they draw solid
spheres around all points and then shoot an infinitesimal ray to let it hit
the spheres to identify the nearest neighbors. They used the force exerted by
the nearest neighbors to direct their graph drawing algorithm. Evangelou et
al. [1] used RT cores to perform photon mapping by finding the set points in a
fixed-radius neighborhood of a query point. Zhu [2] proposed optimizations
such as point reordering and query partitioning to improve ray coherency for
RT-accelerated nearest neighbor applications.

###  3.6. NNS on RT cores

Although RT architecture is intended to aid in visualization tasks, recent
research (section 3.5) has shown that with the right choice of the scene and
rays, RT cores can be utilized to perform non-ray tracing applications. The
key to all of these approaches is reducing the target application to ray-
tracing: producing a scene with objects and rays such that the hardware-
accelerated ray-bounding-box intersection provided by the RT cores provides a
partial or complete solution to the problem.

##  4\. Design

###  4.1. Filter-Refine

The reduction of the Nearest Neighbor Search problem to a Ray Tracing problem
that was introduced by prior work can be generalized to a Filter-Refine
approach, which is a two-step selection strategy for search problems. First,
we filter a subset of possible candidates from data points and then refine
this subset to produce a result that exactly matches the result of the
baseline search algorithm. We call this approach of filter-refining NNS on RT
cores the $L^{2}$ reduction. To find the nearest $k$ points within a radius
$r$ of a query point $q$,

  1. (1)

We first find all the points that are within an Euclidean distance of $r$, and
this boils down to solving a $r$-bound query problem. We model all the points
that fall inside this bound as the result of ray-object intersection, where
the ray source represents the query point, and the object represents a volume
of radius $r$ around the data point.

  2. (2)

Once the candidates are filtered, we compute the $L^{2}$ distances between
query point $q$ and each of the candidates. They are then sorted to find the
$k$ nearest points. Figure 3 shows this pictorially.

Formulating $L^{2}$ reduction as two phases is the first step towards
proposing a generalized distance computation for $k-$NNS. Solving the
$r$-bound query problem is where the RT cores come in. The hardware
accelerates the search process of candidates by posing them as a part of the
ray tracing scene. The second phase of refining the candidates is no different
from using GPU threads to perform the same work. These two phases can be
interleaved — find a point that is within $r$ distance, dynamically update the
list of $k$ nearest neighbors, and go back to search for candidates again.

###  4.2. $L^{2}$ Reduction

By defining solid spheres of radius $r$ centered around all data points and
casting rays of infinitesimal length from the query points as ray origins, we
solve the $r$-bound query problem. We construct spheres because the locus of
points that are within the $L^{2}$ distance of $r$ from a point $q$ is a solid
sphere of radius $r$ centered at $q$. Since a sphere is a user-geometry, the
RT cores encompass this sphere with an AABB and index the AABB, using BVH, as
a proxy for the sphere. When a ray is launched in RayGen and hits an indexed
bounding box that is encompassing a geometrical object, this intersection is
reported to the user. If a ray of finite length $r$ is cast from an origin
point $q$ into a scene containing several translucent objects with a
refractive index $0$, all the AABBs that are within a Euclidean distance of
$r$ along the ray will be returned as intersections in the AnyHit primitive.
Any AABB that contains $q$ will be reported in the intersection primitive,
irrespective of the direction of the ray. Using a finite-length ray requires
us to pick a direction to cast the ray. However, a sphere is uniform in all
directions, and picking particular directions leads to not choosing to find
spheres in the directions that are not picked. Observe that if an
infinitesimal ray cast from the point source $q$ intersects a solid sphere
with a radius $r$ and center $a$, it implies that the point $a$ is within the
Euclidean distance of $r$ to $q$. The rays of infinitesimal length cast from
the query points as ray origins report every indexed bounding box encompassing
a sphere no further than $r$ distance as an intersection.

When the RT cores report a ray-object intersection since the intersection is
with the AABB that contains the sphere, there is a possibility that the ray
hit dead space — the volume outside the sphere but inside the AABB. We observe
that we do not have to process the points that lead to such intersections. We
only have to consider the data points that fall inside the sphere centered
around the query point. Solving the $r$-bound query problem before computing
the distances helps prune the false positives that led to such dead space
intersections. Figure 3 shows this pictorially.

Figure 3. filter-refine: 1) AABB is applied by RT cores, 2) a sphere is
applied by user filter of $L^{2}$ reduction, 3) candidates are refined to
select top-$k$ neighbors.

###  4.3. $L^{p}$ Reduction

RT cores build BVH by hierarchically combining the closest euclidean distant
AABBs. By using the AABBs that are spatially ordered according to
$L^{2}$-norm, can the objects be ordered to fit a different distance metric?
Yes, as long as the metric holds the following inclusion property.

###### Definition 4.0 (Inclusion Property).

For a given finite radius $r>0$ in distance metric $D$, $\exists$ a function
$f:\mathbb{R}\rightarrow\mathbb{R}$ in $L^{2}$ such that all points
$T^{\prime}$ within $r$ of query point $q$ according to the target distance
metric $D$ are within $f(r)$ in $L^{2}$.

Our filter-refine generalization works for many different distance metrics.
For example, all points within $r=1$ according to $L^{\infty}$ in two
dimensions are within $r=\sqrt{2}$ in $L^{2}$. Thus, we can set $r^{\prime}$
to ensure that all of the nearest points according to the desired norm $D$ are
found. We present $L^{p}$ reduction Algorithm 1 to generalize all $L^{p}$
distance computations.

Algorithm 1 Filter-Refine Algorithm for $L^{p}$ Reduction

1:$AABB.bounds\leftarrow center\pm r^{\prime}$ $\triangleright$ AABB of side
$2r^{\prime}$

2:if an intersection is reported then $\triangleright$ Hardware Filter

3: $X\leftarrow object.center$ $\triangleright$ Data point

4: $Q\leftarrow ray.origin$ $\triangleright$ Query point

5: if $(X.x-Q.x)^{2}+(X.y-Q.y)^{2}+(X.z-Q.z)^{2}\leq(r^{\prime})^{2}$ then

6: $w\leftarrow abs(X.x-Q.x)^{p}+abs(X.y-Q.y)^{p}\newline \hskip
39.83368pt+abs(X.z-Q.z)^{p}$

7: if $w\leq r^{p}$ then $\triangleright$ User Filter

8: if $w<neighbors[k].distance$ then $\triangleright$ Refine

9: Topk.insert($X,w$)

In line 1 of algorithm 1, an AABB is defined to have at least a side of length
$2r^{\prime}$ to fit the sphere of radius$r^{\prime}$. It is also up to the
user to decide how big of a bounding box is needed to render the desired
geometry. Lines 2-9 are executed for every AABB-ray intersection that RT cores
detect. In line 5, we perform preliminary filtering by making use of the
inclusion property to remove the data points whose representative spheres do
not contain the query point. In line 6, we compute the $L^{p}$ norm distance
between the data point which is the center of the locus, and the query point
which is the ray origin. In line 7, we discard the data points whose $L^{p}$
norm volumes do not contain the query point, thus removing false positives by
solving the $r$ bound query problem. Lines 8-9 refine the selected candidates
and store the top k closest points.

###  4.4. Correctnes of Our Reductions

We start by proving $L^{2}$ reduction, which we later on use to prove $L^{p}$
reduction. Proving $L^{2}$ reduction means proving that the following filter
and refine steps correctly return K nearest neighbors within radius $r$ of a
query point.

Hardware Filter: AABB-ray Intersection:

All points that could be at a euclidean distance $\leq r$ are filtered by RT
cores.

User Filter: Sphere-ray Intersection:

Only the points that are at a euclidean distance $\leq r$ are filtered by user
filter.

Refine: Ordering.:

Correctly picks only the K-nearest neighbors of all the candidates screened
out by the user filter.

######  $L^{2}$ reduction Proof-Sketch.

Hardware and User filters solve the $r$-bound query problem. Proving the
correctness of filters boils down to proving that the filter step doesn’t
incorrectly remove candidates for K-nearest neighbors within radius $r$. An
intersection at this stage implies that the distance between the ray
origin(query point) and the center of the sphere(data point) is $\sqrt{3}r$ to
$0$ in 3D space. The set of points that lead to these intersections is a
superset of the points that are $\leq r$. The sphere-ray intersection test
further reduces the upper bound on distances to $r$. If a point falls outside
the sphere in the AABB, it’s farther than any point that is in/on the sphere.
So any point that doesn’t fall inside the sphere need not be considered for
ordering in the refining step. In the refining process, we compare the
distance between the query and data points resulting from every new
intersection with that of the distance between the farthest of top $k$ known
neighbors to the query point. ∎

Since we already proved the correctness of the filter and refine steps, all
that is left to prove the correctness of $L^{p}$ reduction is to prove that
the inclusion property is not removing the candidates that could be a
$k$-nearest neighbor. If distance metric $D$ satisfies the inclusion property,
we can find $k$-nearest neighbors according to $D$ by performing $L^{p}$
reduction with $r^{\prime}$.

######  $L^{p}$ reduction Proof-Sketch.

For a finite $r$, there exists a surface $S_{1}$ such that $D(q,t)\leq
r,\forall t\in T$ around the points following the distance metric $D$. A
sphere $S_{2}$ of radius $r^{\prime}=max(\\{L^{2}(q,t))|\forall t\in T\\})$
encompasses this surface. The candidate set filtered by sphere $S_{2}$ is a
superset of those filtered by $S_{1}$. Since $S_{2}$ has points that are
outside $S_{1}$, points that are farther than $r$ could be in the candidate
set. If the number of points inside $S_{1}$ is $\geq k$, then these points
that are outside $S_{1}$ are refined in line8. On the other hand, if the
number of points inside $S_{1}$ is $<k$, then they are discarded by the user
filter in line7. ∎

###  4.5. Custom geometry Enhanced Reductions

When the RT cores report a ray-object intersection, since the intersection is
with the AABB that contains the sphere, there is a possibility that the ray
hit dead space — the volume outside the sphere but inside the AABB. We observe
that we do not have to process the points that lead to such intersections by
building custom geometries. For a given search radius ($r$) and distance
metric ($D:=L^{p}$), we can build a customized $L^{p}$ reduction (Algorithm 2)
by building locus geometries instead of spheres and guarantee that the
candidate set is a superset of the points that metric $D$ needs to consider to
solve $r$-bound query problem.

Algorithm 2 Enhanced $L^{p}$ reduction

1:$AABB.bounds\leftarrow center\pm r^{\prime}$ $\triangleright$ AABB of side
$2r^{\prime}$

2:if an intersection is reported then $\triangleright$ Hardware Filter

3: $X\leftarrow object.center$

4: $Q\leftarrow ray.origin$

5: $w\leftarrow abs(X.x-Q.x)^{p}+abs(X.y-Q.y)^{p}\newline \hskip
28.45274pt+abs(X.z-Q.z)^{p}$

6: if $w\leq r^{p}$ then $\triangleright$ User Filter with Custom Geometry

7: if $w<neighbors[k].distance$ then $\triangleright$ Refine

8: Topk.insert($X,w$)

#### Case of $1\leq p<2$.

All $L^{p},p<2$ norm volumes are already inside the $L^{2}$ norm sphere. By
the inclusion property, all spheres with radius $r^{\prime}\geq r$ allow us to
use $L^{2}$ reduction. Sphere of radius $r^{\prime}=r$ tightly circumscribes
all the $L^{p},p<2$ norms of radius $r$. In particular, $L^{1}$ norm is a
square bi-pyramid, which is two pyramids on either face of a square base, with
the six vertices on the three axes, on either side of the origin at a uniform
distance $r$. Where $(x,y)$ and $(r,\theta)$ are the Cartesian and spherical
coordinates of a point,

$L^{1}(x,y)=|x|+|y|=rcos\theta+rsin\theta\geq r=L^{2}(x,y)$

All the norms $L^{p},1<p<2$ are curves that lie between $L^{1}$ norm square
bi-pyramid and $L^{2}$ norm sphere. Since the $L^{2}$ norm sphere already
encompasses the norm surfaces, the bounds of AABB remain the same. $L^{p}$
reduction allows us to use BVH that is constructed using the $L^{2}$ norm
distances and find $k$-nearest neighbors according to $L^{1}$ norm. Instead of
constructing a sphere, we construct the exact surface bounding the data
points. In this reduction, observe that we only have to change the user filter
step of the Filter-Refine model.

#### Case of $p>2$.

The $L^{2}$ norm sphere inscribes all $L^{p},p>2$ norm spheres. To build a
circumscribing sphere following the inclusion property, the radius
$r^{\prime}>r$. The spheres of all $L^{p},p>2$ lie between $L^{2}$ norm sphere
and $L^{\infty}$ norm sphere that is a cube. $L^{\infty}$ norm is an axis-
aligned square and cube in 2D and 3D spaces, respectively . A cube is also an
AABB, which means we are not indexing any unnecessary volume when we render it
as a user-defined geometry.

To encompass a cube of radius $r$, we need to construct a sphere of radius
$r^{\prime}=r\sqrt{d}$, where $d$ is the number of dimensions. The AABB would
be of the size of side $2r\sqrt{d}$ without $L^{p}$ reduction because AABB
encompasses the sphere, which encompasses the cube of side $2r$. But with our
enhanced $L^{p}$ reduction, we do not need the cube to be encompassed by a
sphere, so an AABB with a side of $2r$ suffices. The volume difference between
the inscribing and circumscribing spheres is why we can expect custom geometry
enhanced $L^{p}$ reduction to work faster than $L^{p}$ reduction.

Again, We only change the user filter step to filter all the points inside the
$L^{p},p>2$ volume instead of just the points inside the sphere. In
particular, for the $L^{\infty}$ norm, the user filter step is non-existent
since the user filter has to return all the points filtered out by the
hardware filter as the cube is an AABB.

###  4.6. Monotone Transforms

While $L^{p}$ norms are a large class of helpful distance metrics, not all
metrics satisfy the inclusion property. For example, consider cosine distances
between data points in Figure 1b; one can not draw a surface without including
all the points. If the candidate set is the same as the set of original
points, the reduction is not useful. No $L^{2}$ radius will help us take
advantage of RT cores. Inclusion property is not enough to generalize all the
distances of our interests. We can exploit $L^{p}$ reduction if only we can
transform the points such that their distances according to the $L^{2}$ norm
have the exact ordering as their distances according to the target metric.

Consider performing $k$-NNS on 2D points. It is impossible to come up with a
radius $r^{\prime}$ in $L^{2}$ such that the inclusion property is satisfied.
If we have two $3D$ distances $dist_{3}$, $dist_{3}^{\prime}$ and their
corresponding projected $2D$ distances $dist_{2}$, $dist_{2}^{\prime}$ then
$dist_{2}>dist_{2}^{\prime}\centernot\implies dist_{3}>dist_{3}^{\prime}$. But
consider this particular monotone transformation
$f_{m}:(a_{x},a_{y})\to(a_{x},a_{y},0)$ that zeroes out the $z$ coordinate of
the $3D$ points. Now the $3D$ distance of the transformed points is the same
as their $2D$ distance. The relative ordering of $dist_{2}(q,a)$ in $2D$ space
is the same as their relative ordering of $dist_{3}(f_{m}(q),f_{m}(a))$ in
$3D$. The monotone property of this transformation allows us to apply the
earlier $L^{p}$ reduction using RT cores on the transformed points to get
$k$-NN of original points from domain space. As long as we can find a
monotonic transformation on a distance metric $D$ to $L^{2}$ norm, the order
of the points in the original domain space will be preserved in the target
range space.

###### Definition 4.0.

Monotone Property: Given a query point $q$, and two distance metrics, $D(q,a)$
and $L^{2}(q,a)$, a monotone transformation is a function
$f_{m}:\mathbb{R}^{m}\to\mathbb{R}^{3},m\in\mathbb{N}$ that has the property:
for two points $a_{1}$ and $a_{2}$, and $D(q,a_{1})<D(q,a_{2})$, then
$L^{2}(f(q),f(a_{1}))<L^{2}(f_{m}(q),f_{m}(a_{2}))$.

#### Correctness of Monotone Reduction

A distance metric $D$ that satisfies monotone property, can use $L^{p}$
reduction to correctly compute $k$-nearest neighbors.

######  Proof-Sketch.

The result of $k$-NNS using $L^{p}$ reduction is an ordered set of $k$ points
that are closest to query point $q$. By the monotone property, we know that
this ordering is the same as the ordering of these points according to the
distance metric $D$. ∎

Figure 4. Cosine Reduction: 1) Building sphere user-geometry around points, 2)
Normalizing points for cosine reduction.

###  4.7. Cosine Reduction

As depicted in Figure 1, the cosine distance between arbitrary vectors from
the origin does not correlate with the Euclidean distance between their
endpoints. However, when they are normalized, they fall on the circumference
of a circle with a radius equal to the normalized length, say $R$. The angular
distance between two vectors is the same as the angular distance between their
normalized versions that fall on this circle which is the angle they make at
the center. Let $\alpha$ be the angle between the query point $q$ and data
point $b$. Then we would end up with the following relation between Euclidean
distance and cosine/angular distance.

$r(q,b)=\sqrt{||q||^{2}+||b||^{2}-2||q||||b||cos(\alpha)}$
$=\sqrt{R^{2}+R^{2}-2R^{2}cos(\alpha)}=R\sqrt{2-2cos(\alpha)}$

In the normalized space, the shorter the Euclidean distance between the
points, the closer the points on the circle, the smaller the angle they make
at the center, and the longer the cosine distance between them. Therefore, the
Cosine distance decrease as the Euclidean distance between two data points
increases. Angular distance ordering is positively preserved by Euclidean
distances while Cosine distance ordering is negatively preserved (reversed) by
Euclidean distances. Normalizing the input space is the monotonic
transformation
$f_{m}:(a_{x},a_{y},a_{z})\to(\frac{a_{x}}{\omega},\frac{a_{y}}{\omega},\frac{a_{z}}{\omega},\omega=\sqrt{a_{x}^{2},a_{y}^{2},a_{z}^{2}})$
that let us map the cosine/angular distance ordering to $L^{2}$ distance.
Figure 4 shows working of cosine reduction.

###  4.8. Composability

Putting both the filter-refine and the monotonic transformation stand-alone
reductions together, we can compose multiple monotone transformations together
to chain from one metric to another, and then pass it into any filter/refine
reduction. For example, consider performing $k$-NNS on $2D$ points according
to $L^{1}$ norm. We can first apply a monotonic transformation to points in
$2D$ and map them to $3D$, and then apply $L^{p},p=1$ reduction.

###  4.9. Other Distances

Jaccard Distance:

Jaccard similarity($JS$) between two sets $A,B$ is the $A\cap B/A\cup B$ and
Jaccard distance is $1-JS$. The reduction depends on the type of data and the
application for which the data is being ranked. Let’s say we represent set $A$
as a bit vector$Av$ where $i^{th}$ bit indicates the presence of $i^{th}$ in
the set. $JS=\frac{A\cap B}{A\cup B}=\frac{A\cap B}{{(A-A\cap B)+(B-A\cap
B)}+A\cap B}=\frac{1}{\frac{L^{1}(Av,Bv)}{dot(Av,Bv)}+1}$. It’s hard to bind
the ratio of Manhattan distance to the dot product. Preserving the order
according to Jaccard distance and mapping the distance computation is not
feasible either using Filter-Refine or Monotonic Transformation.

Hamming Distance:

Given two binary data strings, Hamming distance is the count of bit positions
in which the respective bits of the strings are different. This is simple as
computing the Manhattan distance on the binary strings. In 3D space, all the
possible binary data strings represent vertices of a unit cube, and Hamming
distance between these strings, therefore, is the number of edges that need to
be walked from one vertex to the other.

Mahalanobis Distance:

Mahalanobis distance is the distance between a point and distribution. It
allows the detection of an outlier by comparing the standard deviation of the
point to the mean of the distribution. After a particular whitening spatial
transformation, where axes are scaled to unit variance, Mahalanobis distance
is Euclidean distance (mahalanobis, ).

##  5\. Implementation

RT algorithm breaks down every object that needs rendering into triangles and
RT cores are designed to provide hardware support to efficiently bound
triangles and perform ray-triangle intersections. When an object is expressed
in terms of triangles to the ray-tracing architecture, we don’t have to define
the AABBs around the triangles. But when the same object is expressed as a
user geometry type, we have to define how the AABBs bind the geometry. And
there will always be a volume inside the AABB that is not the object volume
unless the user geometry itself is a box. So when the RT cores report an
intersection, it will also include false positives and increases the time
spent filtering.

Out of all the $L^{p}$ geometries, only $L^{1}$ norm and $L^{\infty}$ have
polygonal surfaces that can be easily broken down into triangular faces.
Because most of the optical effects RT cores optimize for involve reflections
where it is enough to report just the closest intersection, RT cores do not
provide as much hardware support for AnyHit even if the surface is
triangulated. AnyHit takes almost the same time irrespective of whether the
object is expressed as hardware-supported triangles or user geometries. That
leaves us with an interesting implementation trade-off between triangulating
the surfaces vs tightly bounding them with AABBs.

##  6\. Evaluation

Experimental Setup:

We used NVIDIA GeForce RTX 2060 GPU with 6GB memory and Intel Xeon Gold 5215
CPU with 2 NUMA nodes, 10 cores, and 512 GB of memory. To interface with the
RT cores on the GPU, we used Optix Wrapper Library[16]. In every run, for
chosen parameters, we build the index over the entire dataset and search for
neighbors once for all the query points. We perform 5 such runs to collect and
average metrics such as build time, search time, and recall (Equation 1). We
perform experiments for $k$-NNS using $L^{1}$ norm and $L^{\infty}$ norm and
expect the performance of other norms to range in between because $L^{1}$ norm
and $L^{\infty}$ norm are the smallest and largest norm volumes, respectively.
We explain further in section 6.2.

RTX implementations:

We refer to our RT-accelerated cosine reduction implementations as "RTX" and
RT-accelerated $L^{p}$ norm implementations as $L^{p}$ or Enhanced $L^{p}$.
$L^{p}$ implementations follow the algorithm 1, while Enhanced $L^{p}$
implementations follow algorithm 2 that build custom specific enhanced
geometries. Enhanced $L^{p}$ implementations have tightly fit locus shapes,
allowing them to be faster in certain cases. We leave the choice of radius to
the user as a tuning parameter, like the other libraries - FAISS and SCANN.

Baseline:

For the cosine reduction, we use the state-of-the-art exact NNS GPU library
FAISS [8] as the baseline and provide a comparison with the state-of-the-art
approximate-NNS library SONG [9]. FAISS uses tensorflow-gpu to interface with
CUDA cores, so we also compare with Treelogy, a Tree-based GPU implementation
(autoropes, ). Treelogy builds a KD-Tree using $L^{2}$ norm and performs
optimizations to reduce thread divergence during the tree traversals. We
modify Treelogy to compute cosine distance and other $L^{p}$ norms. For
$k$-NNS with $L^{p}$ norms besides $L^{2}$, FAISS implements these norms only
for CPU. So, we compare only with Treelogy as a baseline. We also note that
SCANN, FAISS, and SONG are targeted toward higher dimensional data.

Datasets:

We use the point cloud datasets from KITTI[17] raw 3D data footage, which we
sample to pick the first $n$ points as dataset followed by $q$ query points.
We also use Glove-angular[18] dataset with $\approx$ 1.1M points and $1e^{4}$
query points, reduced to 3 dimensions to perform dot product experiments. We
obtain the ground truth from a naive CPU program that visits all the data
points.

(a) Dot product, #queries=$1e^{4}$, k=10

(b) $L^{\infty}$ norm, #queries=$1e^{4}$, k=10

(c) $L^{\infty}$ norm, KITTI 3M, k=10

Figure 5. Comparison of our reductions with baselines

(a) k=10, #queries=$1e^{4}$

(b) k=10, #queries=$1e^{4}$

(c) radius=$1.5e^{6}$, #queries=$1e^{4}$

Figure 6. $L^{p}$ reduction performance analysis on KITTI 3M

(a) 3D-Glove 1.1M, #queries=$5e^{4}$, k=10

(b) 3D-Glove 1.1M, #queries=$1e^{4}$, k=10

(c) KITTI 3M, #queries=$5e^{4}$, k=10

Figure 7. Cosine reduction performance comparison

###  6.1. Comparison to Baselines

Figure 5(a) shows hardware vs algorithmic optimization advances in dot product
based $k$-NNS solutions for 3 different datasets. Except for RTX, all the
others use CPU to build index, so we only compare their search times. Reported
times are for an iso-recall of 0.997 for SCANN, SONG, and RTX implementations
while CUDA is an exact implementation. For our RTX implementation, We report a
speedup of 100x, 43x and 49x over FAISS for Glove, KITTI 3M, and KITTI 6M
datasets respectively.

The general trend we notice is that GPU accelerated and optimized libraries
outperform their counterparts, except for Treelogy. This is to be expected
from the naive GPU implementation of tree traversals. The larger size of the
dataset makes the tree complex, which further exacerbates the irregularity of
the memory accesses by the threads. Since the other GPU libraries FAISS and
SONG do not use trees and the optimizations they further employed make them
faster than SCANN and Treelogy. Finally, our RTX implementation which is only
accelerated by RT cores without any software optimizations outperforms all of
our baselines. The RT cores schedule threads so that they access closer parts
of tree and memory.

Figures 5(b) and 5(c) show the performance of our RTX implementation against
Treelogy, the CUDA baseline for $L^{\infty}$ norm reduction as the number of
queries increases from $1e^{4}$ to $1e^{5}$. We show only the performance of
$L^{\infty}$ norm reduction against Treelogy because as we will see in later
plots LABEL:fig:lp_radius,fig:lp_k that performance gap between $L^{1}$ norm
reduction and its enhanced version is very narrow. Figure 5(b) shows that the
build times of our RTX implementations against Treelogy for three different
datasets. As expected, we see that build time increases for both the
implementations as the dataset size increases. Further we note, that RTX
implementations are much faster compared to Treelogy. RT cores not only
accelerate BVH traversal but also the tree construction itself. Whereas
Treelogy builds the tree on the CPU and copies it to GPU. RT cores’ interface
also employs a parallel memory copy to the GPU from the CPU.

We next investigate the search time performance of RTX $L^{\infty}$
implementations and Treelogy as shown in figure 5(c). Treelogy implements
optimizations to tame the thread and memory divergence during the tree
traversal, however we see that its still prominent when compared to a hardware
accelerator that addresses the same irregular tree traversals.

###  6.2. Impact of parametric radius $r$

Figure 6(a) shows the search time trend for various of our $L^{p}$ norm
implementations as the radius $r$ changes. We see that search time increases
with an increase in the radius of the norm volume, which is also linearly
proportional to the width of AABB. In general, as $r$ increases, the time
taken by the search kernel increases in the order of $O(r^{3})$. This increase
in search time could be due to two reasons: (i) the number of possible
intersections increases as the AABB grows to include more data points, and
(ii) the overlap between the AABBs increases, making the BVH less useful. We
further observe that Enhanced $L^{\infty}$ norm reduction is faster compared
to rest of the implementations. For $p>2$, Enhanced reductions (sec 4.5) are
more efficient as the custom-defined geometries reduce the size of AABB
compared to their counterpart $L^{p}$ reduction. Whereas in the case of $p<2$,
they remain the same. We do not notice much difference in build times even
though bigger bounding boxes are being constructed.

Figure 6(b) shows how the increase in radius impact the recall that the
$L^{p}$ norm implementation achieves. We observe that recall increases with
radius exponentially, which is to say that the recall is proportional to the
volume of AABB more than the radius itself. As the size of the AABB increases,
the number of ray-AABB intersection tests performed during BVH traversal also
increases, which produces results with higher recall at the cost of higher
search times. For $p<2$, even though the AABB itself is not bigger, the
circumscribed sphere allows for more intersection tests compared to the custom
geometries.

Increase in radius of the AABB increases the number of points screened, so we
expect the recall to keep increasing with radius until it reaches 1. On the
contrary, the figure shows that recall reaches a saturation even though the
radius keeps increasing. This could be due to the presence of outliers that
might need very large radius. The presence and position of saturation point
depends on the dataset; however, this allows us to choose a lower radius to
achieve near-perfect recall with lower search times. After the saturation
point, increasing the radius from $r$ to $r+\delta$ increases the search time
proportional to $\delta^{3}$ but recall barely increases if at all. A wiser
choice would be to choose a radius at which we observe saturation point.

###  6.3. Sensitivity to $k$

We next investigate the performance of our reductions as the workload per
query increases. We increase the number of nearest neighbors that are searched
for each query point, keeping all the parameters same. We notice in figure
LABEL:fig:6c that the search time increases as $k$ increases. We also observe
that enhanced reductions are always faster. Moreover, the performance gap
between $L^{1}$ norm and its enhanced version, and $L^{\infty}$ norm and its
enhanced version is not the same. This again signifies the impact of tightly
fit custom geometries. For our RTX implementations, the performance depends on
how many points the Filter phase forwards to the refining process, as we do
not optimize for performance in the Refine phase. Since the radius remains the
same, the points that are forwarded to refine phase remain the same as $k$ is
increased, but the probability that a candidate can be a top-$k$ nearest
neighbor increases. With the custom-defined geometries, even fewer points are
passed to the refine process further increasing the candidate probability, and
improving performance. The performance gain is proportional to the ratio of
the volumes of custom geometry to that of the sphere. So the benefit of using
tightly-fit AABB of user geometry is maximal for $L^{\infty}$, where the
custom geometry itself is the AABB.

###  6.4. Comparison of Cosine Reduction and FAISS

Figure 7(a) shows the performance of our cosine reduction against FAISS as $k$
increases. We observe that our RTX reduction is outperformed by FAISS when
searching for larger number of nearest neighbors. Even though FAISS search is
faster than RTX implementation for a given $k$, the time FAISS would take to
achieve the same recall as RTX code is higher, which is evident from Figure
7(b). While RTX takes much less time to build than search, FAISS spends
significantly more time in indexing than RT cores. So the overall run
time(build time + search time) of both the implementations is comparable.

We next investigate the time that both the implementations take to finish when
achieving a certain recall and present the results as figure 7(b) . Before
discussing the trends, RTX and FAISS implementations a tunable parameter. The
recall of RTX implementation increases with an increase in radius $r$ till it
reaches a saturation point for cosine reduction. FAISS has a control parameter
nlist which is inversely proportional to recall and positively proportional to
search time. We do not show the performance trends of the respective control
variable (radius for RTX, nlist for FAISS) versus search time. However, in
Figure 7(b), we compare the index and search times of both implementations as
they achieve different recalls when the control parameters are varied. We note
that RTX implementation is faster than FAISS both in terms of index and search
times, in achieving a certain recall.

In figure 7(c), we plot the search times that RTX implementation, FAISS, and
Treelogy take as the number of queries are varied. The reported search times
are for an iso-recall of 0.98. In general, it is well known that search time
increases as the number of queries increases for any NNS algorithm which is
also evident in the plot. We further note that our cosine reduction is faster
than FAISS and Treelogy at any given number of queries. We believe our RTX
reduction is more scalable than FAISS, as FAISS does not run beyond $70000$
queries and it is not straightforward to make it run for a larger number of
queries the way RTX implementation does. We note that our RTX implementations
can run till 25M query size and probably beyond depending on available GPU
memory size.

##  7\. Conclusion

Irregular problems like tree traversals are ubiquitous, especially queries
like nearest neighbor search that has applications in domains such as point
cloud registration in computer vision, data compression, similarity scoring,
DNA sequencing, etc. Tree-based nearest neighbor search is naturally difficult
to scale up using purely software approaches on massively parallel commodity
hardware. Even though ray tracing cores are specialized hardware to cater to
graphics applications, by providing a class of reductions to the ray tracing
scene, we show that this specialized hardware can be generalized to accelerate
tree operations in other domains.

Without our reductions, distance metric computations such as $L^{p}$ norm and
cosine distance take significantly longer to complete or cannot be run on RT
cores. We want to emphasize that our goal is to show the wide applicability of
our reductions rather than scheduling or optimizing the applications. While RT
cores accelerate tree traversals through BVH construction, this tree structure
is not accessible to the user and is limited to 3D space. Availability and
programmability of the spatial tree itself would be more helpful in using RT
cores for general applications.

## References

  * [1] I. Evangelou, G. Papaioannou, K. Vardis, and A. A. Vasilakis.  Fast radius search exploiting ray tracing frameworks.  Journal of Computer Graphics Techniques (JCGT), 10(1):25–48, February 2021.
  * [2] Yuhao Zhu.  Rtnn: Accelerating neighbor search using hardware ray tracing.  In Proceedings of the 27th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, PPoPP ’22, page 76–89, New York, NY, USA, 2022. Association for Computing Machinery.
  * [3] Mark J. Reid and Karl M. Menten.  The first stellar parallaxes revisited.  Astronomische Nachrichten, 341(9):860–869, nov 2020.
  * [4] Paul E. Black.  Manhattan distance, in dictionary of algorithms and data structures ed., 11 February 2019.
  * [5] Tiancheng Xu, Boyuan Tian, and Yuhao Zhu.  Tigris: Architecture and algorithms for 3d perception in point clouds.  In Proceedings of the 52nd Annual IEEE/ACM International Symposium on Microarchitecture, MICRO ’52, page 629–642, New York, NY, USA, 2019. Association for Computing Machinery.
  * [6] Neil D. Lawrence.  A Unifying Probabilistic Perspective for Spectral Dimensionality Reduction: Insights and New Models.  arXiv e-prints, page arXiv:1010.4830, October 2010.
  * [7] Martin Aumüller, Erik Bernhardsson, and Alexander Faithfull.  Ann-benchmarks: A benchmarking tool for approximate nearest neighbor algorithms.  Information Systems, 87:101374, 2020.
  * [8] J. Johnson, M. Douze, and H. Jegou.  Billion-scale similarity search with gpus.  IEEE Transactions on Big Data, 7(03):535–547, Jul 2021.
  * [9] Weijie Zhao, Shulong Tan, and Ping Li.  Song: Approximate nearest neighbor search on gpu.  In 2020 IEEE 36th International Conference on Data Engineering (ICDE), pages 1033–1044, 2020.
  * [10] Michael Goldfarb, Youngjoon Jo, and Milind Kulkarni.  General transformations for gpu execution of tree traversals.  In Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, SC ’13, New York, NY, USA, 2013. Association for Computing Machinery.
  * [11] Marius Muja and David G. Lowe.  Scalable nearest neighbor algorithms for high dimensional data.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 36(11):2227–2240, 2014.
  * [12] Ingo Wald, Will Usher, Nathan Morrical, Laura Lediaev, and Valerio Pascucci.  RTX Beyond Ray Tracing: Exploring the Use of Hardware Ray Tracing Cores for Tet-Mesh Point Location.  In Markus Steinberger and Tim Foley, editors, High-Performance Graphics - Short Papers. The Eurographics Association, 2019.
  * [13] Stefan Zellmann, Martin Weier, and Ingo Wald.  Accelerating force-directed graph drawing with rt cores.  In 2020 IEEE Visualization Conference (VIS), pages 96–100, 2020\.
  * [14] Yu A. Malkov and D. A. Yashunin.  Efficient and robust approximate nearest neighbor search using hierarchical navigable small world graphs.  IEEE Trans. Pattern Anal. Mach. Intell., 42(4):824–836, apr 2020\.
  * [15] Prasanta Chandra Mahalanobis.  On the generalised distance in statistics., 1936.
  * [16] Owl: A node graph "wrapper" library for optix 7.
  * [17] The kitti vision benchmark suite.
  * [18] Glove: Global vectors for word representation.