# Target Score Matching

Valentin De Bortoli<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Google DeepMind  Michael Hutchinson
Isomorphic Labs  Peter Wirnsberger  Isomorphic Labs  Arnaud Doucet  Google
DeepMind

###### Abstract

Denoising Score Matching estimates the score of a “noised” version of a target
distribution by minimizing a regression loss and is widely used to train the
popular class of Denoising Diffusion Models. A well known limitation of
Denoising Score Matching, however, is that it yields poor estimates of the
score at low noise levels. This issue is particularly unfavourable for
problems in the physical sciences and for Monte Carlo sampling tasks for which
the score of the “clean” original target is known. Intuitively, estimating the
score of a slightly noised version of the target should be a simple task in
such cases. In this paper, we address this shortcoming and show that it is
indeed possible to leverage knowledge of the target score. We present a Target
Score Identity and corresponding Target Score Matching regression loss which
allows us to obtain score estimates admitting favourable properties at low
noise levels.

##  1 Introduction and Motivation

###  1.1 Denoising Score Identity and Denoising Score Matching

Consider a $\mathbb{R}^{d}$-valued random variable $X\sim p_{X}$ and let
$Y|(X=x)\sim p_{Y|X}(\cdot|x)$ be a “noisy” version of $X$. We denote the
joint density of $(X,Y)$ by $p_{X,Y}(x,y)=p_{X}(x)p_{Y|X}(y|x)$ and refer to
expectation and variance w.r.t. to this distribution as $\mathbb{E}_{X,Y}$ and
$\textup{Var}_{X,Y}$, respectively. We are interested in estimating the score
of the distribution of $Y$, that is $\nabla\log p_{Y}(y)$, where

$p_{Y}(y)=\int p_{X}(x)p_{Y|X}(y|x)\mathrm{d}x.$ (1)

Evaluating this score is particularly useful for denoising tasks, especially
Denoising Diffusion Models (DDM) which require estimates at different noise
levels (Sohl-Dickstein et al.,, 2015; Ho et al.,, 2020; Song et al.,, 2021). A
standard derivation shows that the identity

$\nabla\log p_{Y}(y)=\int\nabla\log p_{Y|X}(y|x)~{}p_{X|Y}(x|y)\mathrm{d}x$
(2)

holds under mild regularity assumptions, where henceforth $\nabla\log
p_{Y|X}(y|x)$ denotes the gradient with respect to $y$ and

$p_{X|Y}(x|y)=\frac{p_{X}(x)p_{Y|X}(y|x)}{p_{Y}(y)}$ (3)

is the posterior density of $X$ given $Y=y$ (see, for example, Vincent,
(2011)). We will refer to (2) as the _Denoising Score Identity_
($\mathrm{DSI}$). In scenarios where it is possible to compute $\nabla\log
p_{Y|X}(y|x)$ and $p_{X|Y}(x|y)$ is known up to a normalizing constant, the
score can then be approximated by estimating the expectation in (2) using
Markov chain Monte Carlo (MCMC); see e.g. Appendix D.3 in (Vargas et al.,
2023a, ) and (Huang et al.,, 2024).

However, in most standard generative modeling applications, we only have
access to $\nabla\log p_{Y|X}(y|x)$ and i.i.d. samples
$(X_{i},Y_{i})_{i=1}^{n}$ from $p_{X,Y}(x,y)$ and do not know $p_{X|Y}(x|y)$.
In this context, _Denoising Score Matching_ ($\mathrm{DSM}$) (Vincent,, 2011)
leverages $\mathrm{DSI}$ (2) by approximating $\nabla\log p_{Y}(y)$ using some
function $s_{Y}^{\theta}(y)$ whose parameters are obtained by minimizing the
regression loss

$\ell_{\textup{DSM}}(\theta)=\mathbb{E}_{X,Y}[||s_{Y}^{\theta}(Y)-\nabla\log
p_{Y|X}(Y|X)||^{2}],$ (4)

which is in practice approximated using samples $(X_{i},Y_{i})_{i=1}^{n}$.
$\mathrm{DSM}$ is an alternative to Implicit Score Matching (Hyvärinen,, 2005)
which only requires access to noisy samples $(Y_{i})_{i=1}^{n}$ from
$p_{Y}(y)$ but is computationally more expensive in high dimensions as
optimizing the corresponding loss requires computing the gradient of the
divergence of a $d$-dimensional vector field.

###  1.2 Limitations

Consider the case where $Y$ is obtained by adding some independent noise $W$
to $X$, i.e.

$Y=X+W,\qquad W\sim p_{W}(\cdot).$ (5)

If one can sample from $p_{X|Y}(x|y)$, $\mathrm{DSI}$ (2) suggests a Monte
Carlo estimator of the score $\nabla\log p_{Y}(y)$ obtained by averaging
$\nabla\log p_{Y|X}(y|X)$ over samples $X\sim p_{X|Y}(\cdot|y)$. In the case
of Gaussian noise, $p_{W}(w)=\mathcal{N}(w;0,\sigma^{2}I)$, we have that
$\sum_{i=1}^{d}\mathrm{Var}_{X|Y}((\nabla\log p_{Y|X}(Y|X))_{i})\sim
d\sigma^{-2}$ as $\sigma\to 0$. So the variance of such an estimator is higher
for low noise levels.

The high variance of the Monte Carlo estimator for low noise levels is an
independent issue to the high variance of the $\mathrm{DSM}$ regression loss
used to approximate this estimator. Indeed, while the original $\mathrm{DSM}$
loss also exhibits high variance at low noise levels, we can re-arrange (2) to
obtain the so-called _Tweedie identity_ (Robbins,, 1956; Miyasawa,, 1961;
Raphan and Simoncelli,, 2011)

$\mathbb{E}[X|Y=y]=y+\sigma^{2}\nabla\log p_{Y}(y).$ (6)

This identity provides us with an alternative way for computing the score by
estimating $\mathbb{E}[X|Y=y]$, using a regression loss of the form
$\mathbb{E}_{X,Y}[||x^{\theta}(Y)-X||^{2}]$. In this case, the regression
target is simply $X$ and therefore does not exhibit exploding variance as
$\sigma\to 0$. However, this approximation $x^{\theta}$ is used to compute the
score as

$s^{\theta}(y)=(x^{\theta}(y)-y)/\sigma^{2}.$ (7)

Hence the error of $x^{\theta}$ is amplified as $\sigma\to 0$. In fact, this
approach is simply equivalent to a rescaling of the $\mathrm{DSM}$ loss as
$\mathbb{E}_{X,Y}[||x^{\theta}(Y)-X||^{2}]=\sigma^{4}\ell_{\textup{DSM}}(\theta)$.
In the DDM literature, this reparameterisation is sometimes called the
$x_{0}$-prediction. Many other regression targets have been proposed in the
context of diffusion models, see Salimans and Ho, (2022) for instance. All
these parameterisations exhibit the same issues as the direct score prediction
or the $x_{0}$-prediction, since the resulting score approximation
$s_{\theta}$ exhibits exploding variance as $\sigma\to 0$.

Other techniques have also been proposed in order to derive better behaved
regression losses, see for example Wang et al., (2020) and Karras et al.,
(2022). While these works mitigate the high variance of the regression target,
we emphasize that they fail to address the fundamental variance issue of the
score estimator itself.

##  2 Target Score Identity and Target Score Matching

We will focus hereafter on scenarios where the score $\nabla\log p_{X}(x)$ of
the “clean” target $p_{X}(x)$ can be computed exactly. As mentioned earlier,
this is not the case for most generative modeling applications where
$p_{X}(x)$ is only available through samples. However, the score $\nabla\log
p_{X}(x)$ is known in many physical science applications and Monte Carlo
sampling tasks that are actively investigated using denoising techniques; see
e.g. Zhang and Chen, (2022); Arts et al., (2023); Cotler and Rezchikov,
(2023); Herron et al., (2023); Vargas et al., 2023a ; Vargas et al., 2023b ;
Wang et al., (2023); Zhang et al., (2023); Zheng et al., (2023); Akhound-
Sadegh et al., (2024); Huang et al., (2024); Phillips et al., (2024); Richter
et al., (2024).

For ease of presentation, we restrict ourselves to the additive noise model
(5) in this section and discuss a more general setup in the next section. At
low noise levels, we expect $\nabla\log p_{Y}(y)\approx\nabla\log p_{X}(y)$
yet neither the $\mathrm{DSI}$ (2) nor the $\mathrm{DSM}$ loss (4) take
advantage of this fact. On the contrary, the _Target Score Identity_
($\mathrm{TSI}$) and the _Target Score Matching_ ($\mathrm{TSI}$) loss that we
present below address this shortcoming and leverage explicit knowledge of
$\nabla\log p_{X}(x)$. Henceforth, we assume that all the regularity
conditions allowing us to differentiate the densities and interchange
differentiation and integration are satisfied.

######  Proposition 2.1.

For the additive noise model (5), the following _Target Score Identity_ holds

$\nabla\log p_{Y}(y)=\int\nabla\log p_{X}(x)~{}p_{X|Y}(x|y)\mathrm{d}x.$ (8)

######  Corollary 2.2.

By symmetry, we also have

$\nabla\log p_{Y}(y)=\int\nabla\log
p_{W}(w)~{}p_{W|Y}(w|y)\mathrm{d}w=\int\nabla\log
p_{W}(y-x)~{}p_{X|Y}(x|y)\mathrm{d}x,$ (9)

which is $\mathrm{DSI}$ for (5).

The proof of this result and all other results in the main paper are given in
Appendix A. An alternative proof of this identity for additive Gaussian noise
relying on the Fokker–Planck equation is provided in Appendix B.

This result is not new, which is to be expected given its simplicity. It is
part of the folklore in information theory and can be found, for example, in
(Blachman,, 1965)111In machine learning, it was used (and derived
independently) in Appendix C.1.3. of (De Bortoli et al.,, 2021) to establish
some theoretical properties of DDM.. However, to the recent exception of
Phillips et al., (2024), the remarkable computational implications of this
identity do not appear to have been exploited previously. As discussed
further, Akhound-Sadegh et al., (2024) also rely implicitly on this identity.
$\mathrm{TSI}$ shows that, if $p_{X|Y}(x|y)$ is known pointwise up to a
normalizing constant, then the score can be estimated by using an Importance
Sampling (IS) or MCMC approximation of $p_{X|Y}(x|y)$ to compute the
expectation (8). The integrand $\nabla\log p_{X}(x)$ will typically have much
smaller variance under $p_{X|Y}$ than the integrand $\nabla\log p_{Y|X}(y|x)$
appearing in (2) at low noise levels.

Having access to samples $(X_{i},Y_{i})_{i=1}^{n}$ from $p_{X,Y}$, we can also
can estimate the score $\nabla\log p_{Y}(y)$ by minimizing the following
_Target Score Matching_ ($\mathrm{TSM}$) regression loss.

######  Proposition 2.3.

Consider a class of regression functions
$s^{\theta}_{Y}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ for
$\theta\in\Theta$. For the additive noise model (5), we can estimate the score
$\nabla\log p_{Y}(y)$ by minimizing the following regression loss

$\ell_{\textup{TSM}}(\theta)=\mathbb{E}_{X,Y}[||s_{Y}^{\theta}(Y)-\nabla\log
p_{X}(X)||^{2}].$ (10)

Additionally, the $\mathrm{TSM}$ loss $\ell_{\textup{TSM}}(\theta)$ and the
$\mathrm{DSM}$ loss $\ell_{\textup{DSM}}(\theta)$ satisfy

$\ell_{\textup{TSM}}(\theta)=\ell_{\textup{DSM}}(\theta)+\int\|\nabla\log
p_{X}(x)\|^{2}p_{X}(x)\mathrm{d}x-\int\|\nabla\log
p_{Y|X}(y|x)\|^{2}p_{X,Y}(x,y)\mathrm{d}x\mathrm{d}y.$ (11)

Contrary to $\ell_{\textup{DSM}}(\theta)$, $\ell_{\textup{TSM}}(\theta)$ does
not require having access to $\nabla\log p_{Y|X}(y|x)=\nabla\log p_{W}(y-x)$.
This can be useful when the score of the noise distribution is not
analytically tractable; see e.g. Section 3 for applications to Riemmanian
manifolds. Finally, we note that the relationship (11) shows that
$\ell_{\textup{DSM}}(\theta)$ takes large values compared to
$\ell_{\textup{TSM}}(\theta)$ at low noise levels as both quantities are
positive and the expected conditional Fisher information,
$\mathbb{E}_{X,Y}[\|\nabla\log p_{Y|X}(Y|X)\|^{2}]$, takes large values.

######  Corollary 2.4.

In many applications, the observation model of interest is of the form

$Y=\alpha X+W$ (12)

for $\alpha\neq 0$ and $W$ independent of $X$. In this case, $\mathrm{TSI}$
becomes

$\nabla\log p_{Y}(y)=\alpha^{-1}\int\nabla\log
p_{X}(x)~{}p_{X|Y}(x|y)\mathrm{d}x,$ (13)

while $\mathrm{TSM}$ is given by

$\ell_{\textup{TSM}}(\theta)=\mathbb{E}_{X,Y}[||s_{Y}^{\theta}(Y)-\alpha^{-1}\nabla\log
p_{X}(X)||^{2}].$ (14)

For model (12), if the score is estimated through $\mathrm{TSM}$, it is thus
sensible to use a parameterization for $s^{\theta}_{Y}$ of the form
$s^{\theta}_{Y}(y)=\alpha^{-1}\nabla\log p_{X}(y)+\epsilon^{\theta}(y)$.

In practice, we can consider any convex combination of $\mathrm{DSI}$
(equation (2)) and $\mathrm{TSI}$ (equation (13)) to obtain another score
identity which we can use to derive a score matching loss222We can
alternatively consider a convex combination of the $\mathrm{DSM}$ and
$\mathrm{TSM}$ losses.. The score identity considered by Phillips et al.,
(2024) and corresponding score matching loss follow this approach. Therein one
considers a variance preserving denoising diffusion model (Song et al.,, 2021)
where $Y_{t}=\alpha_{t}X+\sqrt{1-\alpha^{2}_{t}}\epsilon$ for
$\epsilon\sim\mathcal{N}(0,I)$ and $(\alpha_{t})_{t\in[0,1]}$ a continuous
decreasing function such that $\alpha_{0}=1$ and $\alpha_{1}\approx 0$. In
this case, $\alpha^{-1}_{t}\nabla\log p_{X}(x)$ will typically be better
behaved than $\nabla\log p_{Y_{t}|X}(y|x)$ for $t$ close to $0$ and vice-versa
for $t$ close to $1$. Phillips et al., (2024) exploits this behaviour to
propose the score identity

$\nabla\log p_{Y_{t}}(y)=\int[\alpha_{t}(x+\nabla\log
p_{X}(x))-y]~{}p_{X|Y_{t}}(x|y)\mathrm{d}x,$ (15)

which is the sum of $\alpha^{2}_{t}$ times $\mathrm{TSI}$ (equation (13)) and
$1-\alpha^{2}_{t}$ times $\mathrm{DSI}$ (equation (2)). They use the integrand
in (15) to define a score matching loss. The rationale for this choice is that
the integrand will be close to the true score for $t\approx 0$ and $t\approx
1$ when $X\sim p_{X|Y_{t}}$. In practice, the “best” loss function one can
consider will be a function of the target $p_{X}$, “noise” $p_{Y|X}$ and
$\alpha$. In Appendix D, we follow the analysis Karras et al., (2022) to
derive a loss admitting desirable properties in a simplified Gaussian setting.

We finally note that the very recent score estimate proposed in (Akhound-
Sadegh et al.,, 2024) (see Eq. (8) therein) can be reinterpreted as a self-
normalized IS approximation in disguise of $\mathrm{TSI}$. It considers for
the model (5) and $W\sim\mathcal{N}(0,\sigma^{2}I)$ the IS proposal
distribution $q(x)=\mathcal{N}(x;y,\sigma^{2}I)=p_{W}(y-x)$ to approximate
$p_{X|Y}(x|y)\propto p_{X}(x)p_{W}(y-x)$. For $\sigma\ll 1$, this importance
distribution will perform well as $p_{X|Y}(x|y)$ is concentrated around $y$.
For larger $\sigma$, the variance of the resulting IS estimate could be
significant.

##  3 Extensions

Next, we present a few extensions of $\mathrm{TSI}$ and $\mathrm{TSM}$.

###  3.1 Extension to non-Additive Noise

We consider here a general noising process defined by

$p_{Y|X}(y|x)=F(\Phi(y,x)),$ (16)

where we assume that $\Phi(y,\cdot)$ is a $\mathrm{C}^{2}$-diffeomorphism for
any $y\in\mathbb{R}^{d}$ and that $\Phi$ is smooth. We denote by $\nabla_{1}$,
respectively $\nabla_{2}$, the gradient of $\Phi$ or $\Phi^{-1}$ with respect
to its first, respectively second, argument.

######  Proposition 3.1.

For the noise model (16), the following _Target Score Identity_ holds

$\displaystyle\nabla\log
p_{Y}(y)=\int[\nabla_{1}\Phi^{-1}(y,\Phi(y,x))^{\top}\nabla\log
p_{X}(x)+\nabla_{y}\log|\det(\nabla_{2}\Phi^{-1}(y,\cdot))|(\Phi^{-1}(y,x))]p_{X|Y}(x|y)\mathrm{d}x.$
(17)

For $\Phi(y,x)=y-\alpha x$ and $F(w)=p_{W}(w)$, we have
$\Phi^{-1}(y,z)=(y-z)/\alpha$ and therefore one has
$\nabla_{1}\Phi^{-1}(y,\Phi(y,x))^{\top}=\mathrm{Id}/\alpha$ and
$\nabla_{y}\log|\det(\nabla_{2}\Phi^{-1}(y,\cdot))|=0$. Hence, we recover
(13). We can thus estimate the score by minimizing the following
$\mathrm{TSM}$ loss

$\ell_{\textup{TSM}}(\theta)=\mathbb{E}_{X,Y}[||s^{\theta}_{Y}(Y)-[\nabla_{1}\Phi^{-1}(Y,\Phi(Y,X))^{\top}\nabla\log
p_{X}(X)+\nabla_{y}\log|\det(\nabla_{2}\Phi^{-1}(Y,\cdot))|(\Phi^{-1}(Y,X))]||^{2}].$
(18)

###  3.2 Extension to Lie groups

Consider a Lie group $G$ which admits a bi-invariant metric. We denote by
$\mu$ the (left) Haar measure on $G$. Let $p_{X}$ denotes the density of $X$
w.r.t. $\mu$. We assume the following additive model on the Lie group, i.e.
$Y=X+_{G}W$, where $+_{G}$ is the group addition on $G$ and $W\sim p_{W}$ with
density $p_{W}$ w.r.t. $\mu$. If $G=\mathbb{R}^{d}$, we recover the Euclidean
additive model of Section 2. For any smooth $f:\ G\to G$ and $x\in G$, we
denote $\mathrm{d}f(x):\ \mathrm{T}_{x}G\to\mathrm{T}_{f(x)}G$, the
differential operator of $f$ evaluated at $x$. Similarly, for any smooth $f:\
G\to\mathbb{R}$, we denote $\nabla f(x)$ its _Riemannian_ gradient.

######  Proposition 3.2.

For a Lie group, the following _Target Score Identity_ holds

$\nabla\log p_{Y}(y)=\int\mathrm{d}R_{x^{-1}y}(x)\nabla\log
p_{X}(x)~{}p_{X|Y}(x|y)\mathrm{d}\mu(x),$ (19)

where $R_{x}(y)=yx$ for any $x,y\in G$. In particular, if $G$ is a matrix Lie
group, we have

$\nabla\log p_{Y}(y)=\int\nabla\log
p_{X}(x)x^{-1}y~{}p_{X|Y}(x|y)\mathrm{d}\mu(x).$ (20)

$\mathrm{DSI}$ and $\mathrm{DSM}$ have been extended to Riemannian manifolds
in (De Bortoli et al.,, 2022; Huang et al.,, 2022); see e.g. Watson et al.,
(2023) for an application to protein modeling. Leveraging the Lie group
structure to obtain a tractable expression of the heat kernel defining
$p_{Y|X}$ and therefore obtain more amenable $\mathrm{DSI}$ and $\mathrm{DSM}$
was considered in (Yim et al.,, 2023; Lou et al.,, 2023; Leach et al.,, 2022).
Contrary to these works, we do not need to know the exact form of the additive
noising process $p_{W}$.

Adapting $\mathrm{DSI}$ and $\mathrm{DSM}$ to Riemmanian manifolds simply
requires replacing the Euclidean gradient by the Riemannian gradient. This is
not the case for $\mathrm{TSI}$, i.e. Proposition 3.2 is not obtained by
replacing the Euclidean gradient by the Riemannian gradient in Proposition
2.1. This is because, while $s_{Y}^{\theta}(y)\in\mathrm{T}_{y}G$, where
$\mathrm{T}_{y}G$ is the tangent space of $G$ at $y$. However, we have that
$\nabla\log p_{X}(x)\in\mathrm{T}_{x}G$ and therefore these two quantities are
not immediately comparable and we use $\mathrm{d}R_{x^{-1}y}(x)$ which
transports $\mathrm{T}_{x}G$ onto $\mathrm{T}_{y}G$. In contrast, in the case
of $\mathrm{DSI}$, both $s_{Y}^{\theta}(y)\in\mathrm{T}_{y}G$ and
$\nabla_{y}\log p_{Y|X}(y|x)\in\mathrm{T}_{y}G$. It is also possible to extend
straightforwardly Proposition 2.3 to the context of Lie groups.

######  Proposition 3.3.

Consider a class of regression functions such that for any $y\in G$,
$s_{\theta}(y)\in\mathrm{T}_{y}G$ for $\theta\in\Theta$. We can estimate the
score $\nabla\log p_{Y}(y)$ by minimizing the $\mathrm{TSM}$ regression loss

$\ell_{\textup{TSM}}(\theta)=\mathbb{E}_{X,Y}[||s_{Y}^{\theta}(Y)-\mathrm{d}R_{X^{-1}Y}(X)\nabla\log
p_{X}(X)||^{2}],$ (21)

where $R_{x}(y)=yx$ for any $x,y\in G$.

###  3.3 Extension to Bridge Matching

Let $Y$ be given by

$Y=\alpha X_{0}+(1-\alpha)X_{1}+W,$ (22)

where $X_{0}\sim p_{X_{0}}$, $X_{1}\sim p_{X_{1}}$, $W\sim p_{W}$ are
independent and $0<\alpha<1$. We are interested in evaluating the score
$\nabla\log p_{Y}(y)$. In the context of generative modeling, (22) appears
when one builds a transport map bridging $p_{X_{0}}$ to $p_{X_{1}}$; see e.g.
(Peluchetti,, 2021; Liu et al.,, 2022; Albergo et al.,, 2023; Lipman et al.,,
2023). We are again considering here a scenario where we have access to the
exact scores of $p_{X_{0}}$ and $p_{X_{1}}$.

######  Proposition 3.4.

For the model (22), the following _Target Score Identity_ holds

$\displaystyle\nabla\log p_{Y}(y)$ $\displaystyle=\alpha^{-1}\int\nabla\log
p_{X_{0}}(x_{0})~{}p_{X_{0},X_{1}|Y}(x_{0},x_{1}|y)\mathrm{d}x_{0}\mathrm{d}x_{1}$
(23) $\displaystyle=(1-\alpha)^{-1}\int\nabla\log
p_{X_{1}}(x_{1})~{}p_{X_{0},X_{1}|Y}(x_{0},x_{1}|y)\mathrm{d}x_{0}\mathrm{d}x_{1},$
(24)

where $p_{X_{0},X_{1}|Y}(x_{0},x_{1}|y)\propto
p_{X_{0}}(x_{0})p_{X_{1}}(x_{1})p_{W}(y-\alpha x_{0}-(1-\alpha x_{1}))$ is the
posterior density of $X_{0},X_{1}$ given $Y=y$.

A convex combination of (23) and (24) yields the elegant identity

$\nabla\log p_{Y}(y)=\int(\nabla\log p_{X_{0}}(x_{0})+\nabla\log
p_{X_{1}}(x_{1}))~{}p_{X_{0},X_{1}|Y}(x_{0},x_{1}|y)\mathrm{d}x_{0}\mathrm{d}x_{1}.$
(25)

We can use these score identities to compute the score using MCMC if
$p_{X_{0},X_{1}|Y}(x_{0},x_{1}|y)$ is available up to a normalizing constant.
Alternatively, given samples from the joint distribution
$p_{X_{0},X_{1},Y}(x_{0},x_{1},y)=p_{X_{0}}(x_{0})p_{X_{1}}(x_{1})p_{Y|X_{0},X_{1}}(y|x_{0},x_{1})$,
we can approximate the score by minimizing a regression loss, e.g. for (26)

$\ell_{\textup{TSM}}(\theta)=\mathbb{E}_{X_{0},X_{1},Y}[||s^{\theta}_{Y}(Y)-(\nabla\log
p_{X_{0}}(X_{0})+\nabla\log p_{X_{1}}(X_{1}))||^{2}].$ (26)

##  4 Experiments

###  4.1 Analytic estimators

We explore experimentally the benefits of these novel score estimators by
considering 1-d mixture of Gaussian targets defined by

$p_{0}(x_{0})=\sum_{i=1}^{N}\pi_{i}\mathcal{N}(x_{0};\mu_{i},\sigma_{i}^{2}).$
(27)

Motivated by DDM, the noising process is defined by a “noising” diffusion

$\mathrm{d}X_{t}=f_{t}X_{t}\mathrm{d}t+g_{t}\mathrm{d}B_{t},$ (28)

where $f_{t}=\tfrac{\mathrm{d}}{\mathrm{d}t}\log\alpha_{t}$ and
$g_{t}^{2}=\tfrac{\mathrm{d}}{\mathrm{d}t}\sigma_{t}^{2}-2f_{t}\sigma_{t}^{2}$
and $B_{t}$ is a standard Brownian motion. Initialized at $X_{0}=x_{0}$, (28)
defines the following conditional distribution of $X_{t}$ given $X_{0}=x_{0}$,
$p_{t|0}(x_{t}|x_{0})=\mathcal{N}(x_{t};\alpha_{t}x_{0},\sigma_{t}^{2})$. In
what follows, we focus on the cosine schedule where
$\alpha_{t}=\cos((\uppi/2)t)$ and $\sigma_{t}=\sin((\uppi/2)t)$. Consider
$X_{0}\sim p_{0}$ then the distribution of $X_{t}$ is given by

$p_{t}(x_{t})=\sum_{i=1}^{N}\pi_{i}\mathcal{N}(x_{t};\mu_{i,t},\sigma_{i,t}^{2}),\quad\mu_{i,t}=\alpha_{t}\mu_{i},\quad\sigma_{i,t}^{2}=\alpha_{t}^{2}\sigma_{i}^{2}+\sigma_{t}^{2}.$
(29)

The posterior distribution of $X_{0}$ given $X_{t}=x_{t}$ is given by another
mixture of Gaussians,

$p_{0|t}(x_{0}|x_{t})=\sum_{i=1}^{N}\pi_{i,t}\mathcal{N}(x_{0};\nu_{i,t},\gamma_{i,t}^{2}),$
(30)

with
$\pi_{i,t}\propto\pi_{i}\tfrac{1}{\sqrt{2\pi\sigma_{i,t}^{2}}}\exp\left(-\tfrac{(x_{t}-\mu_{i,t})^{2}}{2\sigma_{i,t}^{2}}\right),\nu_{i,t}=\mu_{i,t}+\tfrac{c_{i,t}^{2}}{\sigma_{i,t}^{2}}(x_{t}-\mu_{i,t}),\gamma_{i,t}^{2}=\sigma_{i}^{2}-\tfrac{c_{i,t}^{4}}{\sigma_{i,t}^{2}}$
and $c_{i,t}^{2}=\alpha_{t}(\sigma_{i}^{2}+\mu_{i}^{2})-\mu_{i}\mu_{i,t}$.

For convenience, we define the following quantities

$\displaystyle L_{\mathrm{DSI}}(x_{0},x_{t},t)$ $\displaystyle=\nabla\log
p_{t|0}(x_{t}|x_{0}),$ Denoising (31) $\displaystyle
L_{\mathrm{TSI}}(x_{0},x_{t},t)$ $\displaystyle=\alpha^{-1}_{t}\nabla\log
p_{0}(x_{0}),$ Target (32) $\displaystyle L_{w_{t}}(x_{0},x_{t},t)$
$\displaystyle=w_{t}L_{\mathrm{DSI}}(x_{0},x_{t},t)+(1-w_{t})L_{\mathrm{TSI}}(x_{0},x_{t},t),$
$\displaystyle\text{A }w_{t}\text{ mixture}$ (33)

where $w_{t}\in[0,1]$ for any $t\in[0,1]$ is a mixture weight. As described
before, for any of these targets denoted
$L\in\\{L_{\mathrm{DSI}},L_{\mathrm{TSI}},L_{w_{t}}\\}$, we have

$\nabla\log p_{t}(x_{t})=\int
L(x_{0},x_{t},t)p_{0|t}(x_{0}|x_{t})\mathrm{d}x_{t}.$ (34)

We can then estimate $\nabla\log p_{t}(x_{t})$ via the Monte Carlo estimate

$\nabla\log
p_{t}(x_{t})\approx\frac{1}{N}\sum_{i=1}^{N}L(X^{i}_{0},x_{t},t)\quad
X^{i}_{0}\sim p_{0|t}(\cdot|x_{t}),$ (35)

since $p_{0|t}(x_{0}|x_{t})$ given by (30) is easy to sample from. Note that
in more complicated scenarios, we would have had to use MCMC or IS. In
addition, we consider the score matching loss

$\nabla\log
p_{t}(x_{t})\approx\arg\min_{s_{\theta}}\int_{0}^{1}\lambda_{t}\mathbb{E}_{X_{0},X_{t}}\left[||L(X_{0},X_{t},t)-s^{\theta}(X_{t},t)||^{2}\right]\mathrm{d}t,$
(36)

where $\lambda_{t}$ is a weighting function over time and this approximation
holds true jointly over all $t$. Picking $L_{\text{DSI}}$ within (36) gives us
the usual DSM loss, and picking the other identities give a series of novel
score matching losses. In what follows, we explore two key quantities:

  1. 1.

The variance of the Monte Carlo estimates (35) of the score and

  2. 2.

the distribution of the estimated score matching loss functions around the
true score, estimated using Monte Carlo.

We investigate these quantities on a series of target distributions.

  1. 1.

A unit Gaussian,

  2. 2.

a gentle mixture of Gaussians,

  3. 3.

a hard mixture of Gaussians with very separated modes, with the same mode
variances,

  4. 4.

a hard mixture of Gaussians with very separated modes, with the different mode
variances.

The unit Gaussian case we use to find appropriate values for the mixture
weights $w_{t}$. Assuming that we have a target density
$\mathcal{N}(0,\sigma_{\text{data}}^{2}I)$ it is possible to compute the
$w_{t}$ that minimises the variance of the mixture estimator. For this we
recover

$\kappa_{t}:=\frac{\sigma_{t}^{2}}{\sigma_{t}^{2}+\alpha_{t}^{2}~{}\sigma_{\text{data}}^{2}},$
(37)

see Appendix C. In addition we can define another mixing weight

$\bar{\kappa}_{t}:=\frac{\sigma_{t}^{2}}{\sigma_{t}^{2}+\alpha_{t}^{2}~{}\sigma_{\text{data
mode}}^{2}}.$ (38)

The difference between these two is that one used the variance of the whole
target distribution, and the other used the mixture-weighted average of the
variance of each mode in the mixture, $\sigma_{\text{data
mode}}^{2}=\sum_{i}\pi_{i}\sigma_{i}^{2}$. While $\kappa$ is optimal for the
unimodal Gaussian, we will see that $\bar{\kappa}$ performs significantly
better for mixtures of Gaussians.

Figure 1 shows these distributions, and the mixture weightings $\kappa_{t}$
and $\bar{\kappa}_{t}$ for these distributions. In Figure 2 reports the
variance of the Monte Carlo estimators of the score (35) through time. For
each time $t$, we sample 10,000 samples from $p_{t}$. For each of these, we
estimate the score with 100 Monte Carlo samples from $p_{0|t}$. We then
compute the variance of these samples of the estimator.

Figure 1: Target distributions (top panel), and the mixture weights
$\kappa_{t}$ and $\bar{\kappa}_{t}$ through time induced by these targets
(bottom panel). Figure 2: The estimated variance of each estimator based the
score identities. Computed using 10,000 samples of the estimator. For each
estimator sample, $X_{t}$ is sampled from $p_{t}$. For each $X_{t}$, we use
100 samples from $p_{0|t}$ to estimate the score.

In this we see a few points of note. On the Unit Gaussian example, the Monte
Carlo estimates based on $\mathrm{DSI}$ and $\mathrm{TSI}$ perform exactly
oppositely, with alternatively large and small variances at $t\approx 0$ and
$t\approx 1$. The $\kappa_{t}$ and $\bar{\kappa}_{t}$ estimators perform the
same as $\sigma_{\text{data}}$ and $\sigma_{\text{data mode}}$ are the same.
They also have approximately zero variance, as for this simple case the
optimal mixture of the two estimators in fact gives exactly the score.

In the easy mixture case, we see the expected story. Again the $\mathrm{DSI}$
and $\mathrm{TSI}$ estimators perform opposite along time, with variance
blowing up at $t=0$ and $t=1$ respectively. There is however a shift in when
the one becomes better than the other towards $t=0$. The $\kappa_{t}$ and
$\bar{\kappa}_{t}$ estimators both perform well, with neither showing a
variance blowup near $t=0$ or $t=1$. The $\kappa_{t}$ mixture is not quite as
good as the $\mathrm{DSI}$ estimator for $t\geq 0.2$, but the
$\bar{\kappa}_{t}$ performs exactly as well, showing that this is the correct
optimal mixture.

In the hard mixture cases, we see again the optimal switching point between
$\mathrm{DSI}$ and $\mathrm{TSI}$ move further towards $t=0$. We again see
that the $\kappa_{t}$ mixture is not optimal, but that $\bar{\kappa}_{t}$
gives us the best of both $\mathrm{DSI}$ and $\mathrm{TSI}$. Even though the
time frame in which $\mathrm{TSI}$ is better than DSI is small, in this region
we see that the $\mathrm{DSI}$ variance blows up and would cause difficulty
estimating the score well. One thing of note is that in the hard mixture with
the same mode variance the $\bar{\kappa}_{t}$ estimator variance becomes
almost zero. In this case as the modes are so far separated that as $t\to 0$,
$p_{0|t}$ effectively becomes uni-model, bringing us back to the optimal
mixture case.

Next we investigate the variance in the score matching loss functions derived
from the score identities discussed. In this case we need to pick a weighting
scheme $\lambda_{t}$ across time for the score matching losses. We investigate
four weighting schemes here:

  * •

The weighting from Song et al., (2021),
$\lambda_{t}=\frac{1}{\sigma_{t}^{2}}$.

  * •

The weighting from Karras et al., (2022) which ensures for a Gaussian target
that the variance of the $\mathrm{DSM}$ loss at initialization is 1, given by

$\lambda_{t}=\frac{\sigma_{t}^{2}}{\sigma_{\text{data}}^{2}}(\sigma_{t}^{2}+\sigma_{\text{data}}^{2}).$

We refer to this as the DSM optimal weighting.

  * •

A new weighting derived in the spirit of Karras et al., (2022) which ensures
for a Gaussian target that the variance of the $\mathrm{TSM}$ loss is 1 at
initialization, given by

$\lambda_{t}=\frac{\alpha_{t}^{2}\sigma_{\text{data}}^{2}}{\sigma_{t}^{2}}(\sigma_{t}^{2}+\alpha_{t}^{2}\sigma_{\text{data}}^{2}).$

We refer to this as the $\mathrm{TSM}$ optimal weighting (see Appendix D).

  * •

The uniform weighting, $\lambda_{t}=1$.

In addition, we numerically normalise these schedules so that they integrate
to 1. This will not change the properties of the loss function, but will re-
scale the overall value of the loss, and is done to better compare between
losses. Figure 3 depicts the different weighting functions across time, for a
$\sigma_{\text{data}}^{2}=1$, which we ensure all the densities we employ
have. From Figure 3, we observe that the $\mathrm{DSM}$ and $\mathrm{TSM}$
optimal weightings respectively up- and down-weight times where Figure 2 shows
that the $\mathrm{DSI}$ and $\mathrm{TSI}$ estimators had larger variances.

Figure 3: The different weighting functions across time, for a
$\sigma_{\text{data}}^{2}=1$

Figure 4 is produced by looking at the $4^{3}$ combinations of the 4 targets,
4 weighting functions and 4 score matching losses. For each histogram in the
chart, we estimate the loss function at the true score setting
$s_{\theta}(x_{t},t)=\nabla\log p_{t}(x_{t})$ 10,000 times, using 100 samples
from the combined time integral (sampling $t$ uniformly) and expectations. The
distribution of these losses is plotted.

Figure 4: Comparison of the distribution of training losses for the
combinations of the 4 target densities, 4 training losses, and 4 weighting
functions.

###  4.2 Trained score models

Finally, we train our models on a $2$ dimensional Gaussian mixture model to
illustrate some of the advantages of our approach. We consider the four
different settings of $\mathrm{DSM}$, $\mathrm{TSM}$, $\kappa$ and
$\bar{\kappa}$ as described before. To parameterize $s_{\theta}$ we consider a
MLP with three layers of size 128, sinusoidal time embedding with embedding
dimension $128$. We used the ADAM optimizer with learning rate $10^{-4}$.

In Figure 5 (left), we let
$L\in\\{L_{\mathrm{DSM}},L_{\mathrm{TSM}},L_{\kappa_{t}},L_{\bar{\kappa}_{t}}\\}$
and compute the mean of the regression loss $\|s_{\theta}(t,X_{t})-L\|^{2}$.
We consider a uniform weighting strategy for all these losses. The exploding
behavior of the $\mathrm{DSM}$ loss at time $0$ and the exploding behavior of
the $\mathrm{TSM}$ loss at time $1$ is coherent with the theoretical results
of Section 1.2 and Section 2. Note that only the mixture of targets $\kappa$
and $\bar{\kappa}$ achieve a non-exploding behavior for all times. In Figure 5
(right), we show the empirical MMD distance between the empirical data
distribution and generated samples with score $s_{\theta}$ for a RBF kernel.
We emphasize the faster convergence of the mixture of targets with $\kappa$
and $\bar{\kappa}$.

Figure 5: (Left) Mean of the regression loss $\|s_{\theta}(t,X_{t})-L\|^{2}$
with $\|s_{\theta}(t,X_{t})-L\|^{2}$ with
$L\in\\{L_{\mathrm{DSM}},L_{\mathrm{TSM}},L_{\kappa_{t}},L_{\bar{\kappa}_{t}}\\}$
across training iterations. (Right) MMD distance between the empirical data
distribution and generated samples with score $s_{\theta}$ for a RBF kernel.

## References

  * Akhound-Sadegh et al., (2024) Akhound-Sadegh, T., Rector-Brooks, J., Bose, A. J., Mittal, S., Lemos, P., Liu, C.-H., Sendera, M., Ravanbakhsh, S., Gidel, G., Bengio, Y., Malkin, N., and Tong, A. (2024).  Iterated denoising energy matching for sampling from Boltzmann densities.  arXiv preprint arXiv:2402.06121.
  * Albergo et al., (2023) Albergo, M. S., Boffi, N. M., and Vanden-Eijnden, E. (2023).  Stochastic interpolants: A unifying framework for flows and diffusions.  arXiv preprint arXiv:2303.08797.
  * Arts et al., (2023) Arts, M., Garcia Satorras, V., Huang, C.-W., Zugner, D., Federici, M., Clementi, C., Noe, F., Pinsler, R., and van den Berg, R. (2023).  Two for one: Diffusion models and force fields for coarse-grained molecular dynamics.  Journal of Chemical Theory and Computation, 19(18):6151–6159.
  * Blachman, (1965) Blachman, N. (1965).  The convolution inequality for entropy powers.  IEEE Transactions on Information Theory, 11(2):267–271.
  * Cattiaux et al., (2023) Cattiaux, P., Conforti, G., Gentil, I., and Léonard, C. (2023).  Time reversal of diffusion processes under a finite entropy condition.  Annales de l’Institut Henri Poincaré (B) Probabilites et statistiques, 59(4):1844–1881.
  * Cotler and Rezchikov, (2023) Cotler, J. and Rezchikov, S. (2023).  Renormalizing diffusion models.  arXiv preprint arXiv:2308.12355.
  * De Bortoli et al., (2022) De Bortoli, V., Mathieu, E., Hutchinson, M., Thornton, J., Teh, Y. W., and Doucet, A. (2022).  Riemannian score-based generative modelling.  Advances in Neural Information Processing Systems, 35:2406–2422.
  * De Bortoli et al., (2021) De Bortoli, V., Thornton, J., Heng, J., and Doucet, A. (2021).  Diffusion Schrödinger bridge with applications to score-based generative modeling.  In Advances in Neural Information Processing Systems.
  * Herron et al., (2023) Herron, L., Mondal, K., Schneekloth, J. S., and Tiwary, P. (2023).  Inferring phase transitions and critical exponents from limited observations with thermodynamic maps.  arXiv preprint arXiv:2308.14885.
  * Ho et al., (2020) Ho, J., Jain, A., and Abbeel, P. (2020).  Denoising diffusion probabilistic models.  In Advances in Neural Information Processing Systems.
  * Huang et al., (2022) Huang, C.-W., Aghajohari, M., Bose, J., Panangaden, P., and Courville, A. C. (2022).  Riemannian diffusion models.  In Advances in Neural Information Processing Systems.
  * Huang et al., (2024) Huang, X., Dong, H., Hao, Y., Ma, Y., and Zhang, T. (2024).  Reverse diffusion Monte Carlo.  In International Conference on Learning Representations.
  * Hyvärinen, (2005) Hyvärinen, A. (2005).  Estimation of non-normalized statistical models by score matching.  Journal of Machine Learning Research, 6:695–709.
  * Karras et al., (2022) Karras, T., Aittala, M., Aila, T., and Laine, S. (2022).  Elucidating the design space of diffusion-based generative models.  In Advances in Neural Information Processing Systems.
  * Leach et al., (2022) Leach, A., Schmon, S. M., Degiacomi, M. T., and Willcocks, C. G. (2022).  Denoising diffusion probabilistic models on so(3) for rotational alignment.  In ICLR 2022 Workshop on Geometrical and Topological Representation Learning.
  * Lipman et al., (2023) Lipman, Y., Chen, R. T., Ben-Hamu, H., Nickel, M., and Le, M. (2023).  Flow matching for generative modeling.  In International Conference on Learning Representations.
  * Liu et al., (2022) Liu, X., Wu, L., Ye, M., and Liu, Q. (2022).  Let us build bridges: Understanding and extending diffusion generative models.  arXiv preprint arXiv:2208.14699.
  * Lou et al., (2023) Lou, A., Xu, M., and Ermon, S. (2023).  Scaling Riemannian diffusion models.  arXiv preprint arXiv:2310.20030.
  * Miyasawa, (1961) Miyasawa, K. (1961).  An empirical Bayes estimator of the mean of a normal population.  Bull. Inst. Internat. Statist, 38(181-188):1–2.
  * Peluchetti, (2021) Peluchetti, S. (2021).  Non-denoising forward-time diffusions.  https://openreview.net/forum?id=oVfIKuhqfC.
  * Phillips et al., (2024) Phillips, A., Dang, H. D., Hutchinson, M., De Bortoli, V., Deligiannidis, G., and Doucet, A. (2024).  Particle denoising diffusion sampler.  arXiv preprint arXiv:2402.06320.
  * Raphan and Simoncelli, (2011) Raphan, M. and Simoncelli, E. P. (2011).  Least squares estimation without priors or supervision.  Neural Computation, 23(2):374–420.
  * Richter et al., (2024) Richter, L., Berner, J., and Liu, G.-H. (2024).  Improved sampling via learned diffusions.  In International Conference on Learning Representations.
  * Robbins, (1956) Robbins, H. E. (1956).  An empirical Bayes approach to statistics.  In Proc. 3rd Berkeley Symp. Math. Statist. Probab., 1956, volume 1, pages 157–163.
  * Salimans and Ho, (2022) Salimans, T. and Ho, J. (2022).  Progressive distillation for fast sampling of diffusion models.  arXiv preprint arXiv:2202.00512.
  * Sohl-Dickstein et al., (2015) Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., and Ganguli, S. (2015).  Deep unsupervised learning using nonequilibrium thermodynamics.  In International Conference on Machine Learning.
  * Song et al., (2021) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B. (2021).  Score-based generative modeling through stochastic differential equations.  In International Conference on Learning Representations.
  * (28) Vargas, F., Grathwohl, W., and Doucet, A. (2023a).  Denoising diffusion samplers.  In International Conference on Learning Representations.
  * (29) Vargas, F., Ovsianas, A., Fernandes, D., Girolami, M., Lawrence, N. D., and Nüsken, N. (2023b).  Bayesian learning via neural Schrödinger–Föllmer flows.  Statistics and Computing, 33(1):3.
  * Vincent, (2011) Vincent, P. (2011).  A connection between score matching and denoising autoencoders.  Neural Computation, 23(7):1661–1674.
  * Wang et al., (2023) Wang, L., Aarts, G., and Zhou, K. (2023).  Generative diffusion models for lattice field theory.  arXiv preprint arXiv:2311.03578.
  * Wang et al., (2020) Wang, Z., Cheng, S., Yueru, L., Zhu, J., and Zhang, B. (2020).  A Wasserstein minimum velocity approach to learning unnormalized models.  In International Conference on Artificial Intelligence and Statistics.
  * Watson et al., (2023) Watson, J. L., Juergens, D., Bennett, N. R., Trippe, B. L., Yim, J., Eisenach, H. E., Ahern, W., Borst, A. J., Ragotte, R. J., Milles, L. F., et al. (2023).  De novo design of protein structure and function with RFdiffusion.  Nature, 620(7976):1089–1100.
  * Yim et al., (2023) Yim, J., Trippe, B. L., De Bortoli, V., Mathieu, E., Doucet, A., Barzilay, R., and Jaakkola, T. (2023).  Se(3) diffusion model with application to protein backbone generation.  In International Conference on Machine Learning.
  * Zhang et al., (2023) Zhang, D., Chen, R. T. Q., Liu, C.-H., , A., and Bengio, Y. (2023).  Diffusion generative flow samplers: Improving learning signals through partial trajectory optimization.  arXiv preprint arXiv:2310.02679.
  * Zhang and Chen, (2022) Zhang, Q. and Chen, Y. (2022).  Path integral sampler: a stochastic control approach for sampling.  In International Conference on Learning Representations.
  * Zheng et al., (2023) Zheng, S., He, J., Liu, C., Shi, Y., Lu, Z., Feng, W., Ju, F., Wang, J., Zhu, J., Min, Y., Zhang, H., Tang, S., Hao, H., Jin, P., Chen, C., Noé, F., Liu, H., and Liu, T.-Y. (2023).  Towards predicting equilibrium distributions for molecular systems with deep learning.  arXiv preprint arXiv:2306.05445.

## Appendix

In Appendix A, we present the proofs of all the Propositions presented in the
main paper. In Appendix B, we present another proof of $\mathrm{TSI}$
(equation (8)) for additive Gaussian noise based on diffusion techniques. In
Appendix C we derive the variance of Monte Carlo estimates of the score based
on $\mathrm{DSI}$ and $\mathrm{TSI}$ in the Gaussian case. In Appendix D, we
transpose the analysis of Karras et al., (2022) developed to obtain a stable
training loss for DDM to the $\mathrm{TSM}$ loss.

##  Appendix A Proofs of the Main Results

###  A.1 Proof of Proposition 2.1

For completeness, we present two proofs of this result (without any claim for
originality).

#### First proof.

We have $Y=X+W$, $X$ and $W$ being independent, so

$p_{Y}(y)=\int p_{X}(y-w)p_{W}(w)\mathrm{d}w,$ (39)

hence

$\nabla p_{Y}(y)=\int\nabla\log p_{X}(y-w)~{}p_{X}(y-w)p_{W}(w)\mathrm{d}w.$
(40)

It follows that

$\displaystyle\nabla\log p_{Y}(y)$ $\displaystyle=\int\nabla\log
p_{X}(y-w)\frac{p_{X}(y-w)p_{W}(w)}{p_{Y}(y)}\mathrm{d}w$ (41)
$\displaystyle=\int\nabla\log
p_{X}(x)\frac{p_{X}(x)p_{W}(y-x)}{p_{Y}(y)}\mathrm{d}x$ (42)
$\displaystyle=\int\nabla\log p_{X}(x)~{}p_{X|Y}(x|y)\mathrm{d}x.$ (43)

#### Second proof.

For this alternative proof, it is essential for clarity to emphasize
notationally which variable we differentiate with respect to. This proof
starts from $\mathrm{DSI}$ and shows that we can recover $\mathrm{TSI}$. We
have $p_{Y|X}(y|x)=p_{W}(y-x)$ so by the chain rule

$\nabla_{y}\log p_{Y|X}(y|x)=-\nabla_{x}\log p_{Y|X}(y|x).$ (44)

Now by Bayes rule, we have

$\log p_{Y|X}(y|x)=\log p_{X|Y}(x|y)+\log p_{Y}(y)-\log p_{X}(x).$ (45)

Hence, from (45) and (44), we obtain directly

$\nabla_{y}\log p_{Y|X}(y|x)=\nabla_{x}\log p_{X}(x)-\nabla_{x}\log
p_{X|Y}(x|y).$ (46)

Additionally, we have by the divergence theorem

$\int\nabla_{x}\log p_{X|Y}(x|y)~{}p_{X|Y}(x|y)\mathrm{d}x=0.$ (47)

The identity (8) follows directly by combining (2) with (46) and (47).

###  A.2 Proof of Proposition 2.3

We have

$\displaystyle\ell_{\textup{DSM}}(\theta)$
$\displaystyle=\int\|s_{Y}^{\theta}(y)-\nabla\log
p_{Y|X}(y|x)\|^{2}p_{X,Y}(x,y)\mathrm{d}x\mathrm{d}y$ (48)
$\displaystyle=\int\|s_{Y}^{\theta}(y)\|^{2}p_{Y}(y)\mathrm{d}y-2\int\langle
s_{Y}^{\theta}(y),\nabla\log p_{Y|X}(y|x)\rangle
p_{X,Y}(x,y)\mathrm{d}x\mathrm{d}y$ (49)
$\displaystyle\qquad\qquad\qquad+\int\|\nabla\log
p_{Y|X}(y|x)\|^{2}p_{X,Y}(x,y)\mathrm{d}x\mathrm{d}y.$ (50)

From (46) and (47), we get

$\int\langle s_{Y}^{\theta}(y),\nabla\log p_{Y|X}(y|x)\rangle
p_{X,Y}(x,y)\mathrm{d}x=\int\langle s_{Y}^{\theta}(y),\nabla\log
p_{X}(x)\rangle p_{X,Y}(x,y)\mathrm{d}x$ (51)

so it follows that

$\displaystyle\ell_{\textup{DSM}}(\theta)$
$\displaystyle=\int\|s_{Y}^{\theta}(y)\|^{2}p_{Y}(y)\mathrm{d}y-2\int\langle
s_{Y}^{\theta}(y),\nabla\log p_{X}(x)\rangle
p_{X,Y}(x,y)\mathrm{d}x\mathrm{d}y$ (52)
$\displaystyle\qquad\qquad\qquad+\int\|\nabla\log
p_{Y|X}(y|x)\|^{2}p_{X,Y}(x,y)\mathrm{d}x\mathrm{d}y$ (53)
$\displaystyle=\int\|s_{Y}^{\theta}(y)-\nabla\log
p_{X}(x)\|^{2}p_{Y}(y)\mathrm{d}y-\int\|\nabla\log
p_{X}(x)\|^{2}p_{X}(x)\mathrm{d}x$ (54)
$\displaystyle\qquad\qquad\qquad+\int\|\nabla\log
p_{Y|X}(y|x)\|^{2}p_{X,Y}(x,y)\mathrm{d}x\mathrm{d}y.$ (55)

The first term on the r.h.s. is equal to $\ell_{\mathrm{TSM}}(\theta)$ so the
result follows.

###  A.3 Proof of Proposition 3.1

Combining (1) and (16), one has

$\displaystyle p_{Y}(y)$ $\displaystyle=\int p_{X}(x)F(\Phi(y,x))\mathrm{d}x$
(56) $\displaystyle=\int
p_{X}(\Phi^{-1}(y,z))F(z)|\det(\nabla_{2}\Phi(y,\Phi^{-1}(y,z))^{-1})|\mathrm{d}z$
(57) $\displaystyle=\int
p_{X}(\Phi^{-1}(y,z))|\det(\nabla_{2}\Phi^{-1}(y,z))|F(z)\mathrm{d}z$ (58)

where we use the change of variables $z=\Phi(y,x)$. Hence, we have by the
chain rule and the change of variables $x=\Phi^{-1}(y,z)$

$\displaystyle\nabla\log p_{Y}(y)$
$\displaystyle=\int[\nabla_{1}\Phi^{-1}(y,z)^{\top}\nabla\log
p_{X}(\Phi^{-1}(y,z))+\nabla_{y}\log|\det(\nabla_{2}\Phi^{-1}(y,z))|]$ (59)
$\displaystyle\quad\quad\times
p_{X}(\Phi^{-1}(y,z))|\det(\nabla\Phi^{-1}(y,z))|F(z)/p_{Y}(y)\mathrm{d}z$
(60) $\displaystyle=\int[\nabla_{1}\Phi^{-1}(y,\Phi(y,x))^{\top}\nabla\log
p_{X}(x)$ (61)
$\displaystyle\qquad\qquad+\nabla_{y}\log|\det(\nabla_{2}\Phi^{-1}(y,\cdot))|(\Phi^{-1}(y,x))]p_{X|Y}(x|y)\mathrm{d}x.$
(62)

###  A.4 Proof of Proposition 3.2

Using that $\mu$ is left-invariant

$p_{Y}(y)=\int p_{X}(x)F(y^{-1}x)\mathrm{d}\mu(x)=\int
p_{X}(R_{x}(y))F(x)\mathrm{d}\mu(x),$ (63)

where $R_{x}(y)=yx$ for any $x,y\in G$. We have that
$\mathrm{d}R_{x}(y)\mathrm{d}R_{x^{-1}}(yx)=\mathrm{Id}$. Therefore for any
$y\in G$ and $h\in\mathrm{T}_{y}(G)$ we have

$\mathrm{d}(p_{X}\circ
R_{x})(y)(h)=\mathrm{d}p_{X}(yx)\mathrm{d}R_{x}(y)(h)=\langle\nabla
p_{X}(yx),\mathrm{d}R_{x}(y)(h)\rangle=\langle\mathrm{d}R_{x}(y)\mathrm{d}R_{x^{-1}}(yx)\nabla
p_{X}(yx),\mathrm{d}R_{x}(y)(h)\rangle.$ (64)

Finally, using that $\langle\cdot,\cdot\rangle$ is right-invariant we get for
any $y\in G$ and $h\in\mathrm{T}_{y}(G)$

$\mathrm{d}(p_{X}\circ R_{x})(y)(h)=\langle\mathrm{d}R_{x^{-1}}(yx)\nabla
p_{X}(yx),h\rangle.$ (65)

Combining this result and (63) we have

$\displaystyle\nabla p_{Y}(y)$
$\displaystyle=\int\mathrm{d}R_{x^{-1}}(yx)\nabla\log
p_{X}(yx)p_{X}(yx)F(x)\mathrm{d}\mu(x)$ (66)
$\displaystyle=\int\mathrm{d}R_{x^{-1}y}(x)\nabla\log
p_{X}(x)p_{X}(x)F(y^{-1}x)\mathrm{d}\mu(x).$ (67)

Hence, we get that

$\nabla\log p_{Y}(y)=\int\mathrm{d}R_{x^{-1}y}(x)\nabla\log
p_{X}(x)p_{X|Y}(x|y)\mathrm{d}\mu(x).$ (68)

###  A.5 Proof of Proposition 3.4

It follows directly from (22) that

$p_{Y}(y)=\int p_{W}(y-\alpha x_{0}-(1-\alpha
x_{1}))p_{X_{0}}(x_{0})p_{X_{1}}(x_{1})\mathrm{d}x_{0}\mathrm{d}x_{1}.$ (69)

So by considering the change of variables $w=y-\alpha x_{0}-(1-\alpha)x_{1}$,
i.e. $x_{0}=\alpha^{-1}(y-(1-\alpha)x_{1}-w)$, we obtain

$p_{Y}(y)=\alpha^{-1}\int
p_{X_{0}}(\alpha^{-1}(y-(1-\alpha)x_{1}-w))p_{W}(w)p_{X_{1}}(x_{1})\mathrm{d}w\mathrm{d}x_{1}.$
(70)

so, emphasizing here which variable we differentiate with for clarity, we get

$\displaystyle\nabla_{y}\log p_{Y}(y)$ (71) $\displaystyle=$
$\displaystyle\alpha^{-1}\int\nabla_{y}\log
p_{X_{0}}(\alpha^{-1}(y-(1-\alpha)x_{1}-w))\frac{p_{X_{0}}(\alpha^{-1}(y-(1-\alpha)x_{1}-w))p_{W}(w)p_{X_{1}}(x_{1})}{p_{Y}(y)}\mathrm{d}w\mathrm{d}x_{1}$
(72) $\displaystyle=$ $\displaystyle\alpha^{-2}\int\nabla_{x_{0}}\log
p_{X_{0}}(\alpha^{-1}(y-(1-\alpha)x_{1}-w))\frac{p_{X_{0}}(\alpha^{-1}(y-(1-\alpha)x_{1}-w))p_{W}(w)p_{X_{1}}(x_{1})}{p_{Y}(y)}\mathrm{d}w\mathrm{d}x_{1}$
(73) $\displaystyle=\alpha^{-1}\int\nabla_{x_{0}}\log
p_{X_{0}}(x_{0})\frac{p_{X_{0}}(x_{0})p_{X_{1}}(x_{1})p_{W}(y-\alpha
x_{0}-(1-\alpha x_{1}))}{p_{Y}(y)}\mathrm{d}x_{0}\mathrm{d}x_{1}.$ (74)

The result (23) follows immediately and (24) is obtained similarly.

##  Appendix B Fokker–Planck derivation

So far our derivation of $\mathrm{TSI}$ (8) relies on $\nabla_{y}\log
p_{Y|X}(y|x)=-\nabla_{x}\log p_{Y|X}(y|x)$. This identity is due to the
additive nature of the noising process, i.e. $Y=X+W$ with $W$ independent from
$X$. The change of variables used in Proposition 3.1 uses implicitly the same
property in the case of additive noise. In what follows, we provide another
derivation of $\mathrm{TSI}$ leveraging instead the Fokker-Planck equation,
the time-reversal of diffusions and backward Kolmogorov equation. We consider
the case where $Y=X+W$ with $W\sim\mathcal{N}(0,\sigma^{2}I)$.

In our setting, $Y=X_{t_{0}}$ with $t_{0}=\sigma$, $X_{0}=X$ and
$\mathrm{d}X_{t}=\mathrm{d}B_{t}$, where $(B_{t})_{t\geq 0}$ is a
$d$-dimensional Brownian motion. For any $t\geq 0$, we denote by $p_{t}$ the
density of $X_{t}$. The Fokker-Planck equation shows that
$(p_{t}(x))_{t\in[0,t_{0}]}$ satisfies the heat equation

$\partial p_{t}(x)=\tfrac{1}{2}\Delta p_{t}(x),\qquad p_{0}(x)=p_{X}(x).$ (75)

Hence, we have that

$\displaystyle\partial\log p_{t}(x)$ $\displaystyle=\tfrac{1}{2}\Delta
p_{t}(x)/p_{t}(x)$ (76) $\displaystyle=\tfrac{1}{2}\|\nabla\log
p_{t}(x)\|^{2}+\tfrac{1}{2}\Delta\log p_{t}(x).$ (77)

If we denote $F_{t}(x)=\nabla\log p_{t_{0}-t}(x)$ then we obtain by
differentiating (76) w.r.t. $x$ that for any $x\in\mathbb{R}^{d}$ and
$t\in[0,t_{0}]$

$\partial F_{t}(x)+\nabla F_{t}(x)F_{t}(x)+\tfrac{1}{2}\Delta F_{t}(x)=0.$
(78)

Hence, using the backward Kolmogorov equation we get that for any
$y\in\mathbb{R}^{d}$

$F_{t}(y)=\mathbb{E}[F_{t_{0}}(Z_{t_{0}})\
|Z_{t}=y],\qquad\mathrm{d}Z_{t}=F_{t}(Z_{t})\mathrm{d}t+\mathrm{d}B_{t}=\nabla\log
p_{t_{0}-t}(Z_{t})\mathrm{d}t+\mathrm{d}B_{t}.$ (79)

As $F_{0}=\nabla\log p_{Y}$ and $F_{t_{0}}=\nabla\log p_{X}$, we get that

$\nabla\log p_{Y}(y)=F_{0}(y)=\mathbb{E}[F_{t_{0}}(Z_{t_{0}})\
|Z_{0}=y]=\mathbb{E}[\nabla\log p_{X}(Z_{t_{0}})\ |Z_{0}=y].$ (80)

Finally, we notice that
$(Z_{t})_{t\in[0,t_{0}]}=(X_{t_{0}-t})_{t\in[0,t_{0}]}$ is the time-reversal
of $\mathrm{d}X_{t}=\mathrm{d}B_{t}$ (Cattiaux et al.,, 2023). Hence, we get
that $(Z_{0},Z_{t_{0}})$ admits the same distribution as $(Y,X)$. Combining
this result and (80) we obtain

$\nabla\log p_{Y}(y)=\mathbb{E}_{X|Y}[\nabla\log p_{X}(X)],$ (81)

which corresponds to $\mathrm{TSI}$.

##  Appendix C Combinining $\mathrm{DSI}$ and $\mathrm{TSI}$ Monte Carlo
estimates

We consider a Gaussian target
$p_{X}(x)=\mathcal{N}(x;0,\sigma_{\text{tar}}^{2}I)$ as well as additive
Gaussian noise $p_{W}(x)=\mathcal{N}(w;0,\sigma^{2}I)$. For $\alpha>0$, we set
$Y=\alpha
X+W\sim\mathcal{N}(0,(\alpha^{2}\sigma_{\text{tar}}^{2}+\sigma^{2})I)$. The
posterior density appearing in $\mathrm{DSI}$ and $\mathrm{TSI}$ is given in
this case by

$p_{X|Y}(x|y)=\mathcal{N}(x;\alpha\sigma_{\text{tar}}^{2}/(\alpha^{2}\sigma_{\text{tar}}^{2}+\sigma^{2})y,\sigma^{2}\sigma_{\text{tar}}^{2}/(\alpha^{2}\sigma_{\text{tar}}^{2}+\sigma^{2})I).$
(82)

We compute the variance of the Monte Carlo estimates of the score obtaining by
averaging $\nabla\log p_{Y|X}(y|X^{i})$ ($\mathrm{DSI}$ estimate) and
$\nabla\log p_{X}(X^{i})/\alpha$ ($\mathrm{TSI}$ estimate) over samples
$X^{i}\sim p_{X|Y}(\cdot|y)$. We have that

$\displaystyle\sum_{i=1}^{d}\mathrm{Var}_{X|Y}(\nabla_{i}\log p_{Y|X}(y|X))$
$\displaystyle=\sum_{i=1}^{d}\mathrm{Var}_{X|Y}((\alpha X-y)/\sigma^{2})$ (83)
$\displaystyle=d\alpha^{2}(\sigma_{\text{tar}}/\sigma)^{2}/(\alpha^{2}\sigma_{\text{tar}}^{2}+\sigma^{2}).$
(84)

On the other hand, we have that

$\displaystyle\sum_{i=1}^{d}\mathrm{Var}_{X|Y}(\nabla_{i}\log
p_{X}(X)/\alpha)$
$\displaystyle=\sum_{i=1}^{d}(1/\alpha^{2})\mathrm{Var}_{X|Y}(-X/\sigma_{\text{tar}}^{2})$
(85)
$\displaystyle=d(\sigma/\sigma_{\text{tar}})^{2}/(\alpha^{2}\sigma_{\text{tar}}^{2}+\sigma^{2}).$
(86)

Hence we have

$\sum_{i=1}^{d}\mathrm{Var}_{X|Y}(\nabla_{i}\log
p_{Y|X}(y|X))\leq\sum_{i=1}^{d}\mathrm{Var}_{X|Y}(\nabla_{i}\log
p_{X}(X)/\alpha)$ (87)

if and only if $\sigma_{\text{tar}}^{2}\leq\sigma^{2}/\alpha$. Again, this is
aligned with our previous observations. For $\sigma\gg 1$ then the variance of
the $\mathrm{DSI}$ estimator is lower than the variance of the $\mathrm{TSI}$
estimator. For $\sigma\ll 1$ then the variance of the $\mathrm{TSI}$ estimator
is lower than the variance of the $\mathrm{DSI}$ estimator.

We now consider any convex combination of the integrands appearing in
$\mathrm{DSI}$ and $\mathrm{TSI}$

$\displaystyle Z$ $\displaystyle=\kappa\nabla\log
p_{Y|X}(Y|X)+(1-\kappa)\nabla\log p_{X}(X)/\alpha$ (88)
$\displaystyle=\kappa\alpha
X/\sigma^{2}-(1-\kappa)X/(\alpha\sigma_{\text{tar}}^{2})-\kappa Y/\sigma^{2}$
(89)
$\displaystyle=(\kappa\alpha/\sigma^{2}-(1-\kappa)/(\alpha\sigma_{\text{tar}}^{2}))X-\kappa
Y/\sigma^{2}$ (90)
$\displaystyle=(1/\alpha\sigma_{\text{tar}}^{2})(\kappa(1+\alpha^{2}\sigma_{\text{tar}}^{2}/\sigma^{2})-1)X-\kappa
Y/\sigma^{2}.$ (91)

By construction, the expectation of $Z$ under $p_{X|Y}(x|y)$ is equal to the
score $\nabla_{y}\log p_{Y}(y)$. In order to minimize the variance of $Z$
under $p_{X|Y}(x|y)$, we set
$\kappa=1/(1+\alpha^{2}\sigma_{\text{tar}}^{2}/\sigma^{2})$. Hence when
$\sigma\gg 1$ we get that $\kappa=1$ and when $\sigma\ll 1$ we get that
$\kappa=0$. In this specific Gaussian setting we have that the estimator has
actually zero variance since
$Z=-Y/(\alpha^{2}\sigma_{\text{tar}}^{2}+\sigma^{2})=\nabla_{y}\log p_{Y}(Y)$.

##  Appendix D Preconditioning the training loss

In this section, we follow the analysis (Karras et al.,, 2022, Appendix B.6)
and derive a rescaled training loss for TSM from first principles for the
additive model $Y=\alpha X+W$.

In Karras et al., (2022), the input of the network is scaled by $c_{i}$ and
the output is scaled by $c_{o}$. An additional skip-connection is considered
with weight $c_{s}$. Hence we have

$s^{\theta}_{Y}(y)=c_{o}F_{\theta}(\sigma,c_{i}y)+c_{s}y.$ (92)

The total loss is weighted by $\lambda>0$ and we have

$\displaystyle\mathcal{L}(\theta)$
$\displaystyle=\lambda\ell_{\mathrm{TSM}}(\theta)=\lambda\mathbb{E}_{X,Y}[\|c_{o}F_{\theta}(\sigma,c_{i}Y)+c_{s}Y-\alpha^{-1}\nabla\log
p_{X}(X)\|^{2}].$ (93)

In (Karras et al.,, 2022, Appendix B.6) the hyperparameters
$\lambda,c_{i},c_{o},c_{s}$ are computed in the case of the DSM loss with
$x_{0}$-prediction using the following principles

  1. (i)

the input of the network $F_{\theta}$ should have unit variance,

  2. (ii)

the target of the regression loss should have unit variance,

  3. (iii)

the effective weighting of the loss defined as $\lambda c_{0}^{2}$ should be
equal to one,

  4. (iv)

we choose $c_{s}$ to minimize $c_{o}$ so that the errors of the network are
not amplified.

For simplicity, we assume that $p_{X}=\mathcal{N}(0,\sigma_{\text{tar}}^{2}I)$
and $W\sim\mathcal{N}(0,\sigma^{2}I)$. In this case, we have that $\nabla\log
p_{X}(x)=-x/\sigma_{\text{tar}}^{2}$. We then obtain

$\displaystyle\mathcal{L}(\theta)$
$\displaystyle=\lambda\mathbb{E}_{X,Y}[\|c_{o}F_{\theta}(\sigma,c_{i}Y)+c_{s}Y+X/(\alpha\sigma_{\text{tar}}^{2})\|^{2}]$
(94)
$\displaystyle=\lambda/(\alpha\sigma_{\text{tar}})^{4}\mathbb{E}_{X,Y}[\|-(\alpha\sigma_{\text{tar}})^{2}c_{o}F_{\theta}(\sigma,c_{i}Y)-(\alpha\sigma_{\text{tar}})^{2}c_{s}Y-\alpha
X\|^{2}]$ (95)
$\displaystyle=\lambda^{\prime}\mathbb{E}_{X,Y}[\|c_{o}^{\prime}F_{\theta}(\sigma,c_{i}^{\prime}Y)+c_{s}^{\prime}Y-\alpha
X\|^{2}],$ (96)

where

$\lambda^{\prime}=\lambda/(\alpha\sigma_{\text{tar}})^{4},\quad
c_{o}^{\prime}=-(\alpha\sigma_{\text{tar}})^{2}c_{o},\quad
c_{s}^{\prime}=-(\alpha\sigma_{\text{tar}})^{2}c_{s},\quad
c_{i}^{\prime}=c_{i}.$ (97)

We emphasize that $\lambda^{\prime}(c_{o}^{\prime})^{2}=\lambda c_{o}^{2}$.
Therefore, the rescaled effective weight is the same as the effective weight
described in Karras et al., (2022).

Using (Karras et al.,, 2022, (117), (131), (138),(144)), we get

$\displaystyle\lambda^{\prime}=(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})/(\alpha\sigma\sigma_{\text{tar}})^{2},\quad
c_{o}^{\prime}=\alpha\sigma\sigma_{\text{tar}}/(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})^{1/2},$
(98) $\displaystyle
c_{s}^{\prime}=\alpha^{2}\sigma_{\text{tar}}^{2}/(\alpha^{2}\sigma_{\text{tar}}^{2}+\sigma^{2}),\quad
c_{i}=1/\sqrt{\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2}}.$ (99)

Hence, combining this result and (97) we get that

$\displaystyle\lambda=\alpha^{2}\sigma_{\text{tar}}^{2}(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})/\sigma^{2},\quad
c_{o}=-\sigma/[\alpha\sigma_{\text{tar}}(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})^{1/2}],$
(100) $\displaystyle
c_{s}=-1/(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2}),\quad
c_{i}=1/\sqrt{\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2}}.$ (101)

Using (Karras et al.,, 2022, (151)), we get that the variance of
$\lambda\|c_{o}F_{\theta}(\sigma,c_{i}y)+c_{s}y+x/(\alpha\sigma_{\text{tar}}^{2})\|^{2}$
is equal to one for every time $\sigma,\alpha>0$, at initialization for
$F_{\theta}=0$. In the general case, using the hyperparameters
$\lambda,c_{i},c_{s},c_{o}$ we get

$\displaystyle\mathcal{L}(\theta)=\sigma_{\text{tar}}^{2}(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})/\sigma^{2}$
$\displaystyle\mathbb{E}_{X,Y}[\|\sigma
F_{\theta}(\sigma,Y/(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})^{1/2})/[\sigma_{\text{tar}}(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})^{1/2}]$
(102) $\displaystyle\qquad+\alpha
Y/(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})+\nabla\log p_{X}(X)\|^{2}].$
(103)

The score network is then given by

$s_{Y}^{\theta}(y)=-(1/\alpha)[\sigma
F_{\theta}(\sigma,y/(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})^{1/2})/\sigma_{\text{tar}}(\sigma^{2}+\alpha^{2}\sigma_{\text{tar}}^{2})^{1/2}+\alpha^{2}y/(\alpha^{2}\sigma_{\text{tar}}^{2}+\sigma^{2})].$
(104)


11institutetext: Laboratoire d’Informatique et Systèmes, Marseille, France
22institutetext: Université Aix-Merseille, Marseille, France

# Distributed search on graphs using discrete time quantum walk

Mathieu Roget 1122 0000-0003-3890-1552 Giuseppe Di Molfetta 1122
0000-0002-1261-7476

###### Abstract

Searching with coined quantum walk is a problem that has interested the
community since a long time. While most results consider spatial searches on
regular lattices, some work have introduced several models of coined quantum
walks on graphs. This work introduces a distributed searching quantum walk on
graphs. Our contribution is in two parts: (i) we introduce a new mathematical
model of a coined quantum walk on graphs designed to search both nodes or
edges; (ii) we provide an anonymous distributed scheme to implement such a
model.

###### Keywords:

Quantum Walk Distributed Algorithm Searching Algorithm

## Introduction

Spatial search using discrete time quantum walks (DTQW) has been extensively
studied in the past decades. Most of them show up an asymptotic complexity
similar to the Grover algorithm (up to a logarithmic factor) [2, 10, 1].
Compared to other search algorithms, they hold all the interesting properties
of a DTQW, such as locality or translation invariance. The Grover algorithm
has a wide range of applications, including optimisation and machine learning.
Recently, they have been used, as subroutines, to solve distributed
computational tasks [4, 8, 6], including determining the diameter of a network
[8] or the triangle finding problem [6]. However, in all these examples, some
global information is supposed to be known by each node of the network, such
as the size, and usually the network is not anonymous, namely every node has a
unique label. The absence of anonymity and the quantum properties of the
algorithm effectively address the incalculability of certain problems such as
the leader election problem. However, associating global information with a
state in a quantum algorithm is generally problematic for its implementation.
Inspired by the classical anonymous distributed computing where random walks
or Markovian processes are regularly used [9], we propose in this manuscript a
model of searching DTQW for rephrasing distributed algorithm. In the
following, we prove how such an algorithm does not require global information
and that it is possible to implement it in a completely anonymous network.

#### Contribution

Section 1 of this work introduces a Quantum Walk model well suited to search
indifferently a node or an edge in arbitrary graphs. In this model the
walker’s amplitudes are defined onto the graph’s edges and ensures
2-dimensional coin everywhere. We detail how this Quantum Walk can be used to
search a node or an edge and we show examples of this Quantum Walk on several
graphs. In this first part, the Quantum Walk is introduced formally as a
discrete dynamical system.

In Section 2 we introduce a distributed implementation of the Quantum Walk 2.
For the coin and the scattering, local quantum circuits and distributed
algorithms are provided. The implementation proposed conserves the graph
locality and does not require a node or an edge to have global information
like the graph size. The nodes (and edges) do not have unique label, and no
leader is needed.

#### Related work

While several discrete time quantum walk model on graphs have been designed
[5], the closest one is the arc-reversal quantum walk [5, 7]. Indeed, in the
arc-reversal quantum walk, the scattering step corresponds to a particular
coin in our model. The Grover coin step instead is what we define as the
scattering operator. Apart from that, in our searching algorithm, the oracle
is applied onto the edges while in the arc-reversal quantum walk it is applied
on the nodes. More importantly, our model has only coins of dimension $2$,
which is far more practical, especially for noisy-robust experimental
implementations. Finally, our model can search both edges and nodes. All in
all, the model presented in this paper is designed to solve a searching
computational task in a distributed fashion. This opens up the possibility of
solving a wide range of distributed problems.

##  1 Mathematical model of quantum walk on graphs

This section introduces the quantum walk searching algorithm on undirected
connected graphs. First we introduce the quantum walk model on graph and give
an example. Then we explain how such a quantum walk can be used to search a
marked edge or a marked node. Finally, we show two examples of searching : in
the first one we look for an edge and in the second one we look for a node.
Finally we derive the asymptotic complexities.

###  1.1 Quantum walk model on graphs

We consider an undirected connected graph $G=(V,E)$, where $V$ is the set of
vertices and $E$ the set of edges. We define the walker’s position on the
graph’s edges and a coin register of dimension two (either $+$ or $-$). The
full state of the walker at step $t$ is noted :

$\ket{\Phi_{t}}=\sum_{(u,v)\in
E}\psi_{u,v}^{+}\ket{(u,v)}\ket{+}+\psi_{u,v}^{-}\ket{(u,v)}\ket{-}.$

The graph is undirected so we indifferently use $\ket{(u,v)}$ and
$\ket{(v,u)}$ to name the same state. Similarly for the complex amplitudes
$\psi_{u,v}^{+}(t)=\psi_{u,v}^{-}(t)$. We also introduce a polarity for every
edge of $G$. For each of them, we have a polarity function $\sigma$, such that
: $\forall(u,v)\in E,\;\sigma(u,v)\in\\{+,-\\}\text{ and
}\sigma(u,v)\neq\sigma(v,u).$ Figure 1 illustrates how amplitudes and polarity
are placed with edge $(u,v)$ of polarity $\sigma(u,v)=1$. The polarity is a
necessary and arbitrary choice made at the algorithm’s initialization
independently for every edge. We discuss why it is necessary and a way to make
that choice at the end of this section. In order to update the walker from
step $t$ to step $t+1$, the overall operator $S$ is composed by to unitary
operations ; a coin operation which act only on the coin register and a
scattering operation that moves the amplitudes $\alpha_{u,v}^{\pm}$ according
to the coin state.

(a) $\sigma(u,v)=+$

(b) $\sigma(u,v)=-$

Figure 1: The edge $(u,v)$ and how are placed the amplitudes for the two
possible polarities.

#### Coin

The coin is defined as a $2\times 2$ unitary operator $C$ on the coin register
:

$\forall(u,v)\in
E,\;\ket{(u,v)}\ket{\pm}\overset{\text{coin}}{\longmapsto}(I\otimes
C)\times\ket{(u,v)}\ket{\pm}=\ket{(u,v)}(C\ket{\pm}).$

#### Scattering

The scattering operator coincides instead with the standard diffusion operator
$D_{n}$ of size $n\times n$ applied locally on each node :

$\forall u\in V,\;\left(\psi_{u,v}^{\sigma(u,v)}\right)_{v\in
V}\overset{\text{scattering}}{\longmapsto}D_{\text{deg}(u)}\times\left(\psi_{u,v}^{\sigma(u,v)}\right)_{v\in
V},$

where $D_{n}=\left(\frac{2}{n}\right)_{i,j}-I_{n}$.

#### Dynamic

The full unitary evolition of the walk finally reads :
$\ket{\Phi_{t+1}}=S\times(I\otimes C)\times\ket{\Phi_{t}}$, where
$\ket{\Phi_{0}}$ is the initial state.

As an example, one can consider the path of size 3 with the nodes
$\\{u,v,w\\}$, with polarity $\sigma(u,v)=\sigma(v,w)=+$. Figure 2 shows the
unitary evolution of the walker from step $t$ to step $t+1$, when the coin
coincide the first Pauli matrix $X$.

Figure 2: Example of a walk on a path of size 3.

#### Polarity : why and how

There are two reasons we need polarity. First, if the coin operator acts
differently on $\ket{+}$ and $\ket{-}$, then different polarities lead to
different dynamics. Polarities, in fact, determines how each edge’s state
scatter within the network, analogously to other spatial searching algorithms
on regular lattice [2, 1, 10]. Moreover, polarity is used to divide the edges
in two. So, when scattering, each node can only access one amplitude of each
edge (depending on the polarity of the edges). This makes things much simpler
for a distributed implementation, since we do not have to consider the case of
two simultaneous operations on the same amplitude during scattering. One way
to initialise the polarity on the graph is to create a coloring of $G$ at the
initial time step. Afterward, edges have their $+$ pole in the direction of
the node with the higher color. This design is especially convenient for
bipartite graphs where every node sees only $+$ polarities or only $-$
polarities.

###  1.2 Searching

The above dynamical system formally describes a quantum walks on the graph’s
edges. If the quantum state is measured, we obtain one of the edges of the
graph. This scheme is thus well suited to search edges. Thus we may now
introduce an oracle marking the desired edge $(a,b)\in E$. Analogously with
the standard spatial search, the oracle is defined as follows :

$\mathcal{O}_{f}=\left(\sum_{f((u,v))=1}\ket{(u,v)}\bra{(u,v)}\right)\otimes
R+\left(\sum_{f((u,v))=0}\ket{(u,v)}\bra{(u,v)}\right)\otimes I_{2},$

where $f$ is the classical oracle equals to $1$ if and only if the edge is
marked, and $0$ otherwise. The operator $R$ is a special coin operator which
is applied only to the marked edge. Without lack of generality, in the
following we set $C=X$ and $R=-X$.

The Algorithm 1 describes the procedure to search an edge. The initial state
of the quantum walk is initialized to be diagonal and the additional oracle
operation is added to the dynamic. There are two parameters that characterize
the algorithm performances here; the probability of success $P$ of returning
the marked edge and the hitting time $T$. The number of oracle calls of this
algorithm is $O(T)$. In practice, we want to choose $T$ such that $P$ is
maximal.

Algorithm 1 Search a marked edge

1:$G=(V,E)$ undirected, connected graph

2:$f$ the classical oracle

3:$T\in\mathbb{N}$ the hitting time

4:function Search($G,f,T$)

5: $\displaystyle\ket{\Phi}\leftarrow\frac{1}{\sqrt{2|E|}}\sum_{(u,v)\in
E}\ket{(u,v)}\ket{+}+\ket{(u,v)}\ket{-}$ $\triangleright$ Diagonal initial
state.

6: for $0\leq i<T$ do

7: $\displaystyle\ket{\Phi}\leftarrow S\times(I\otimes
C)\times\mathcal{O}_{f}\times\ket{\Phi}$ $\triangleright$ One step of the walk

8: end for

9: $(u,v,\pm)\leftarrow$ Measure($\ket{\Phi}$) $\triangleright$ Measure the
quantum state.

10: return $(u,v)$

11:end function

The previous algorithm can easily be made into a Las Vegas algorithm whose
answer is always correct but the running time random. Algorithm 2 shows such
transformation. Interestingly, the expected number of times Algorithm 2 calls
Algorithm 1 is $O(1/P)$ where $P$ is the probability of success of Algorithm
1. The expected complexity of Algorithm 2 is $O\left(T/P\right).$

Algorithm 2 Search a marked edge with guaranteed success

1:$G=(V,E)$ undirected, connected graph

2:$f$, the classical oracle

3:$T\in\mathbb{N}$ the hitting time

4:function GuaranteedSearch($G,f,T$)

5: $(u,v)\leftarrow(\texttt{nil},\texttt{nil})$ $\triangleright$ Initial value

6: while $f((u,v))\neq 1$ do $\triangleright$ Until we find the marked edge …

7: $(u,v)\leftarrow$ Search($G,f,T$) $\triangleright$ … search again.

8: end while

9: return $(u,v)$

10:end function

#### Remark on the complexity

The complexity of our algorithm is limited by the optimal complexity of
$O(\sqrt{K})$ for searching problem where $K$ is the total number of elements.
This is the complexity of the Grover algorithm which has been shown to be
optimal [12]. However Grover’s algorithm assume full connectivity between
elements which is not our case. The quantum walk presented here is on the edge
of $G$, which means that the optimal complexity is $O\left(\sqrt{|E|}\right)$.

###  1.3 Example : Searching an edge in the star graph

In this example we consider the star graph with $M$ edges and $M+1$ nodes; a
graph with $M+1$ nodes where every node (other than node $u_{0}$) is connected
to node $u_{0}$. Our searching algorithm performs well on this graph, as shown
by Theorem 1.1. In this section, we show the proof of Theorem 1.1 which is
mainly spectral analysis on a simplified dynamic.

###### Theorem 1.1

Let $G$ be the star graph with $M$ edges. Algorithm 1 has an optimal hitting
time $T=O\left(\sqrt{M}\right)$ and a probability of success $O(1)$ for $G$.
Algorithm 2 has expected complexity $O\left(\sqrt{M}\right)$.

###### Proof

We consider the star graph $G=(V,E)$ of size $M+1$ such that
$V=\\{u_{0},\ldots,u_{M}\\}\text{ and }E=\\{(u_{0},u_{i})\mid 1\leq i\leq
M\\}.$ We assume without a loss of generality that the marked edge is
$(u_{0},u_{1})$ and that the polarity is $\forall
i>1,\;\sigma(u_{0},u_{i})=+$.

At any time step $t$, the state of the walk is noted

$\ket{\Phi_{t}}=\sum_{i=1}^{M}\psi_{u_{0},u_{i}}^{+}(t)\ket{(u_{0},u_{i})}\ket{+}+\psi_{u_{0},u_{i}}^{-}(t)\ket{(u_{0},u_{i})}\ket{-}.$

We first show that $\forall t\in\mathbb{N},\;\forall
i>1,\;\psi_{u_{0},u_{i}}^{+}(t)=\alpha^{+}\text{ and
}\psi_{u_{0},u_{i}}^{-}(t)=\alpha^{-}.$ This greatly simplifies the way we
describe the walk dynamic. Afterward we shall provide simple spectral analysis
to extract the optimal hitting time $T$ and probability of success $P$.

Next, we prove $(Q_{t}):\quad\forall
i>1,\;\psi_{u_{0},u_{i}}^{+}(t)=\alpha_{t}^{+}\text{ and
}\psi_{u_{0},u_{i}}^{-}(t)=\alpha_{t}^{-}$ for all $t\in\mathbb{N}$.

Initialization: The initial state $\ket{\Psi_{0}}$ is the diagonal state:

$\ket{\Phi}\leftarrow\frac{1}{\sqrt{2|E|}}\sum_{(u,v)\in
E}\ket{(u,v)}\ket{+}+\ket{(u,v)}\ket{-}.$

All $\left(\psi_{u_{0},u_{i}}^{\pm}(0)\right)_{i\geq 1}$ are equals so the
property $Q_{0}$ is satisfied with
$\alpha_{t}^{+}=\alpha_{t}^{-}=\frac{1}{\sqrt{2M}}$.

Recurrence: We assume that $Q_{t}$ is true and use the walk dynamic to show
that $Q_{t+1}$ is also true. The state $\ket{\Psi_{t+1}}$ is described in
Table 1 and show that $Q_{t+1}$ is true.

Table 1: Detailed dynamic of the quantum walk for star graphs. | Step $t$ | After oracle | After coin | After Scattering
---|---|---|---|---
$\alpha^{+}$ | $\alpha^{+}_{t}$ | $\alpha^{+}_{t}$ | $\alpha^{-}_{t}$ | $\displaystyle\alpha^{+}_{t+1}=\frac{(M-2)\alpha_{t}^{-}-2\psi_{u_{0},u_{1}}^{+}}{M}(t)$
$\alpha^{-}$ | $\alpha^{-}_{t}$ | $\alpha^{-}_{t}$ | $\alpha^{+}_{t}$ | $\alpha^{-}_{t+1}=\alpha^{+}_{t}$
$\psi_{u_{0},u_{1}}^{+}$ | $\psi_{u_{0},u_{1}}^{+}(t)$ | $-\psi_{u_{0},u_{1}}^{-}(t)$ | $-\psi_{u_{0},u_{1}}^{+}(t)$ | $\displaystyle\psi_{u_{0},u_{1}}^{+}(t+1)=\frac{(2M-2)\alpha_{t}^{-}+(M-2)\psi_{u_{0},u_{1}}^{+}}{M}(t)$
$\psi_{u_{0},u_{1}}^{-}$ | $\psi_{u_{0},u_{1}}^{-}(t)$ | $-\psi_{u_{0},u_{1}}^{+}(t)$ | $-\psi_{u_{0},u_{1}}^{-}(t)$ | $\psi_{u_{0},u_{1}}^{-}(t+1)=-\psi_{u_{0},u_{1}}^{-}(t)$

Using the recurrence in Table 1, we can put the dynamic of the walk into a
matrix form.

$\psi_{u_{0},u_{1}}^{-}(t)=\frac{(-1)}{\sqrt{2M}}\qquad\text{and}\qquad
X_{t+1}=AX_{t},$

where

$X_{t}=\left(\begin{matrix}\alpha^{+}_{t}\\\ \alpha_{t}^{-}\\\
\psi_{u_{0},u_{1}}^{+}(t)\\\ \end{matrix}\right)\qquad\text{and}\qquad
A=\left(\begin{matrix}0&\frac{M-2}{M}&\frac{-2}{M}\\\ 1&0&0\\\
0&2\frac{M-1}{M}&\frac{M-2}{M}\\\ \end{matrix}\right).$

We can diagonalize matrix $A$ in $\mathbb{C}$: $A=PDP^{-1}$. The eigenvalues
are $\\{-1,e^{i\lambda},e^{-i\lambda}\\}$ where
$e^{i\lambda}=\frac{M-1+i\sqrt{2M-1}}{M}$. Numerical resolution of this
equation gives $\psi_{u_{0},u_{1}}^{+}(t)\sim\frac{i}{2}\left(e^{-i\lambda
t}-e^{i\lambda t}\right)\sim\sin(\lambda t).$ With this we can deduce the
probability $p_{t}$ of hitting the marked edge $p_{t}\sim\sin^{2}(\lambda t),$
and the optimal hitting time by solving $p_{T}=1$:

$T=\frac{\pi}{2}\frac{1}{\lambda}\sim\frac{\pi}{2}\sqrt{\frac{M}{2}}\sim\frac{\pi}{2\sqrt{2}}\sqrt{M}=O(\sqrt{M}).$

Finally we have a probability of success $P\sim 1$ for the optimal hitting
time $\displaystyle T\sim\frac{\pi}{2\sqrt{2}}\sqrt{M}=O(\sqrt{M})$.

###  1.4 Searching nodes

The searching quantum walk presented in the previous section can search one
marked edge in a graph. In order to search a node instead, we need to
transform the graph we walk on. We call this transformation starify.

###### Definition 1

Starify
Let us consider an undirected graph $G=(V,E)$. The transformation starify on
$G$ returns a graph $\tilde{G}$ for which

  * •

every node and edge in $G$ is in $\tilde{G}$ (we call them real nodes and real
edges).

  * •

for every node $u\in V$, we add a node $\tilde{u}$ (we call these new nodes
virtual nodes).

  * •

for every node $u\in V$, we add the edge $(u,\tilde{u})$ (we call these new
edges virtual edges).

We call the resulting graph $\tilde{G}$ the starified graph of $G$.

#### Searching nodes

Now that we have this transformation, we are ready to search a marked node
instead of an edge. Let us assume that we want to search a node $u$ in $G$. We
starify it into the graph $\tilde{G}$, and we mark the virtual edge
corresponding to the node $u$ (i.e. edge $(u,\tilde{u})$). We can then use the
previous searching walk to search the virtual edge and then deduce without
ambiguity the marked node $u$ of the initial graph $G$. This procedure implies
that we have to increase the size of the graph (number of edges and nodes). In
particular, increasing the number of edges is significant since we increase
the dimension of the walker (this can be significant for sparse graphs) and
one must be careful of this when computing the complexity.

A reassuring result is that searching a node on the complete graph (the
strongest possible connectivity) has optimal complexity. As stated in Theorem
1.2, the complexity is $O(\sqrt{M})$, which is optimal when searching over the
edges. Compared to a classical algorithm in $O(M)$ (a depth first search for
instance), we do have a quadratic speedup.

###### Theorem 1.2

When using the quantum algorithm to search one marked node in the starified
graph $\tilde{G}$ of the complete graph $G$ of size $N$, the probability of
success is $P\sim 1$ and the hitting time $T\sim\frac{\pi}{4}N$.

###### Proof

The proof is very similar to the one of Theorem 1.1. We show that several
edges have the same states to reduce the dynamic to a simple recursive
equation. We then solve it numerically and derive asymptotic values for the
probability of success. A detailed proof shall be provided in the
supplementary materials.

##  2 Distributed Implementation

In this section we detail a distributed implementation of the Quantum walk
presented in the previous section. We assume that both nodes and edges have
quantum and classical registers. They can transmit classical bit and apply
controlled operations over the local neighborhood. The nodes need to know
their degree but do not require a unique label. We first explain how the
quantum registers are dispatched along the graph, then we introduce a
distributed protocol to reproduce the dynamic of the walk. In this section we
assume that we want to implement the walk on a graph $G=(V,E)$ with polarity
$\sigma$.

###  2.1 Positioning the qubits

The implementation proposed uses $W$ states[3]; i.e. states of the form
$(\ket{001}+\ket{010}+\ket{100})/\sqrt{3}$. Basically the superposition of all
the unary digits. We use two qubits per edge, one for the $+$ state and one
for the $-$ state. For every node $u\in V$, we use $\ln(\text{deg}(u))+1$
qubits to apply the diffusion across all the neighborhood.

#### Edges register

The edge register represents the position of the walker and is the one
measured at the end of the algorithm. It consists of two qubits per edge,
corresponding respectively to the + and - polarity amplitudes. To represent a
unary state in this register, we need to introduce the notation
$\delta_{k}^{n}=\underbrace{0\ldots 0}_{k-1\text{ times}}1\underbrace{0\ldots
0}_{n-k\text{ times}}$ for the string of $n$ digit with $0$ everywhere except
at position $k$ where it has $1$. For instance, $\delta_{2}^{4}=0100$. For the
sake of simplifying the notations, let us arbitrary enumerate the edges of $E$
such that we have $E=\\{e_{1},\ldots,e_{|E|}\\}$. The full state of the walk
is a linear combination of all the $\left(\ket{\delta_{k}^{2|E|}}\right)_{k\in
2|E|}$, where $\ket{\delta_{2k}^{2|E|}}$ corresponds to amplitude
$\ket{e_{k}}\ket{+}$ of the mathematical model, and
$\ket{\delta_{2k+1}^{2|E|}}$ corresponds to amplitude $\ket{e_{k}}\ket{-}$. It
is important that the global state of the register remains a superposition of
$\left(\ket{\delta_{k}^{2|E|}}\right)_{k\in 2|E|}$, as this allows us to
obtain a valid solution to the search problem during measurement. In fact, no
matter which state $\delta_{k}^{2|E|}$ is measured, all the edges are
measuring state 0 except one, which is measuring state 1.

#### Nodes register

The node register consists of $\log\deg+1$ qubits per node. This register is
made up exclusively of auxiliary qubits used for the scattering operation.
During the scattering operation, the amplitudes of the edges around a given
node $u$ are moved to $u$’s auxiliary qubits. The $u$’s qubits then holds the
local amplitudes in binary format, applies the scattering operation and
finally transfers these amplitudes back to the local edges. After every step
of the quantum walk, the nodes register is in state $\ket{0\ldots 0}$ without
any entanglement with the edges register.

###  2.2 Distributed protocol for one step of the walk

The various quantum registers, their sizes, and how they are positionned on
the network has been made clear in the previous section. We know give
distributed schemes for every operation of the mathematical model: oracle,
coin and scattering.

#### Oracle

In the mathematical model, the oracle applies $R$ on the marked edge only. A
definition on the basis of this operator when $e_{1}$ is marked would be

$\left(\begin{matrix}\ket{e_{1}}\ket{+}\\\ \ket{e_{1}}\ket{-}\\\
\end{matrix}\right)\overset{oracle}{\longmapsto}R\left(\begin{matrix}\ket{e_{1}}\ket{+}\\\
\ket{e_{1}}\ket{-}\\\ \end{matrix}\right),\qquad\text{and}\qquad\forall 2\leq
k\leq|E|\forall
s=\pm,\;\ket{e_{k}}\ket{s}\overset{oracle}{\longmapsto}\ket{e_{k}}\ket{s}.$

Translated to the distributed registers it means that

$\left(\begin{matrix}\ket{\delta_{2k}^{2|E|}}\\\ \ket{\delta_{2k+1}^{2|E|}}\\\
\end{matrix}\right)\overset{oracle}{\longmapsto}R\left(\begin{matrix}\ket{\delta_{2k}^{2|E|}}\\\
\ket{\delta_{2k+1}^{2|E|}}\\\ \end{matrix}\right),\qquad\text{and identity
everywhere else}.$

We recall that in our searching algorithm, $R=-X$. This operation can be
achieved by two $Z$ gates and a swap applied only to the marked edge.

#### Coin

The coin operation is very similar to the oracle. We want to apply $C=X$ to
all edge states. In the mathematical model, it can be written as

$\forall k,\;\left(\begin{matrix}\ket{e_{k}}\ket{+}\\\ \ket{e_{k}}\ket{-}\\\
\end{matrix}\right)\overset{oracle}{\longmapsto}C\left(\begin{matrix}\ket{e_{k}}\ket{+}\\\
\ket{e_{k}}\ket{-}\\\ \end{matrix}\right).$

This can be realized by a swap applied on all edges.

#### Scattering

The scattering operation consists in applying a diffusion operator $D$ locally
around the nodes. Since scattering around each node is affecting distinct
qubits, we only need to design the circuit for one node $u$ of degree $d$. We
note $a_{1},\ldots,a_{d}$ the local edges connected to $u$. For this section
we call $\eta_{k}$ the qubits of edge $a_{k}$ accessible to $u$ following the
polarity. Note that, while the local edges (and their qubits) are indexed,
this is completely arbitrary and is not changing the final result of the
scattering. This nice property is due to the shift invariance of diffusion
operator $D$. The $\log d$ first quibts of $u$ are storing the amplitudes of
the local edges in a binary format, while the last one is marking the states
we need to act on. We note $k^{(2)}$ the binary representation of $k$. We
define the operator $\text{Tr}_{k}$ between qubit $\eta_{k}$ and $u$ register
such that

$\ket{0}\ket{0^{(2)}}\ket{0}\overset{\text{Tr}_{k}}{\longmapsto}\ket{0}\ket{0^{(2)}}\ket{0}\qquad\text{and}\qquad\ket{1}\ket{0^{(2)}}\ket{0}\overset{\text{Tr}_{k}}{\longmapsto}\ket{0}\ket{(k-1)^{(2)}}\ket{1}.\\\
$

This operator $\text{Tr}_{k}$ admits a circuit depending of $k$ composed of X,
CNOT and multi controlled Tofolli gates. Algorithm 3 shows a distributed
scheme for realizing this circuit. It is running on the node and uses two
communication methods: RequestCnot(edge, target) which applies a Not on target
controlled by the qubit in edge accessible according to polarity; and
ApplyMCT(edge) which applies a Not on the edge’s qubit controlled by the full
node register. In the worst case, this algorithm require $\log d+1$ calls to
RequestCnot and one call to ApplyMCT. We define Tr by the successive
applications of $\text{Tr}_{k}$ between edge $a_{k}$ and node $u$ for all $k$.
It is worth to note that, while one has to apply every $\text{Tr}_{k}$
successively, the application order does not matter. On the direct basis, the
complete scattering operation we want to achieve can be written as

$\underbrace{\ket{\delta_{k}^{d}}}_{\text{edges}}\ket{0^{(2)}}\ket{0}\overset{\text{Tr}}{\longmapsto}\ket{0}\ket{k^{(2)}}\ket{1}\overset{D}{\longmapsto}\ket{0}D\ket{k^{(2)}}\ket{1}\overset{\text{Tr}^{-1}}{\longmapsto}\underbrace{\left(\sum_{j=1}^{d}D_{k,j}\ket{\delta_{j}^{d}}\right)}_{\text{walker
scattered}}\otimes\ket{0^{(2)}}\ket{0}.$

The last qubit allows us to apply $D$ only when the walker was actually on
edge $k$. Following this global idea, Algorithm 4 provides a distributed
implementation for the scattering operator. Every node needs a total of
$O(d\log d)$ controlled operations (RequestCnot or ApplyMCT) with the local
edges.

Finally, Figure 3 shows an example of this distributed design of the path
graph of five nodes. The path graph has the particularity of having the same
topology as a circuit (a qubit being connected to the preceding and following
qubit). We can therefore clearly visualize the locality of the implementation
on the circuit.

Algorithm 3 Distributed scheme for $\text{Tr}_{k}$ on node $u$ of degree $d$

1:$e$ an edge connected to $u$

2:$1\leq k\leq d$

3:function $\text{Tr}_{k}$($e,k$)

4: $r\leftarrow\lceil\log d\rceil$ $\triangleright$ The size of $u$’s register
is $r+1$

5: for $0\leq i<r\mid\left((k-1)^{(2)}\right)_{i}=1$ do

6: RequestCnot($e,q_{i}$) $\triangleright$ CNOT on the $i^{\text{th}}$qubit of
$u$

7: end for

8: RequestCnot($e,q_{r}$)

9: Apply $X$ gate on all qubits $q_{i}$ of $u$ such that
$\left((k-1)^{(2)}\right)_{i}=0$.

10: ApplyMCT(e)

11: Apply $X$ gate on all qubits $q_{i}$ of $u$ such that
$\left((k-1)^{(2)}\right)_{i}=0$.

12:end function

Algorithm 4 Distributed scheme for Tr on node $u$ of degree $d$

1:$\mathcal{N}$ the set of all edges connected to $u$

2:function Tr($u$)

3: $\\{a_{1},\ldots,a_{d}\\}=\mathcal{N}$ $\triangleright$ Random or arbitrary
enumeration of the edges

4: for $1\leq k\leq d$ do

5: $\text{Tr}_{k}$($a_{k},k$)

6: end for

7:end function

Figure 3: Circuit of one step of the quantum walk for the path graph $u-v-w$.
Dashed lines signal $\text{Tr}_{1}$ circuit and its inverse while dotted lines
$\text{Tr}_{2}$. The circuit applies successively the oracle on $(u,v)$, the
coin, Tr, $D$, $\text{Tr}^{-1}$.

## Discussion

The distributed algorithm described in Section 2 is based on the quantum
walker dynamics, without requiring a leader or unique identifiers. Once the
measure is done, all edges measure $0$ except one that measure $1$, giving us
a valid solution; i.e. a solution where every edge agrees on who is marked.
The edge who measured $1$ can then send the answer to some hypervisor or
broadcast it across the network, depending of the context. Notice that this
distributed scheme requires that the network be in a $W$ state, which is very
common and widely studied in the distributed algorithm community [3, 11].

To study the performances of this algorithm, one should use the mathematical
model from Section 1. In general, it is hard to study analytically this model,
except in some rare cases, where a lot of edges have the same state during the
whole execution. Indeed, the topology of the dynamics, which is different from
that of the graph, makes its analysis mathematically difficult. In fact, we
have to consider a topology based on edges, where two edges are neighbors if
they share a node. As such, regular lattices like a grid do not actually have
the topology of a grid. However, one could also study this algorithm
numerically; preliminary non published work on this aspect have already shown
promising results for popular lattices and random graphs.

#### Aknowledgement

This work is supported by the PEPR integrated project EPiQ ANR-22-PETQ-0007,
by the ANR JCJC DisQC ANR-22-CE47-0002-01 founded from the French National
Research Agency and with the support of the french government under the France
2030 investment plan, as part of the Initiative d’Excellence d’Aix-Marseille
Université - A*MIDEX AMX-21-RID-011.

## References

  * [1] Arrighi, P., Di Molfetta, G., Márquez-Martín, I., Pérez, A.: Dirac equation as a quantum walk over the honeycomb and triangular lattices. Physical Review A 97(6), 062111 (2018)
  * [2] Bezerra, G., Lugão, P.H., Portugal, R.: Quantum-walk-based search algorithms with multiple marked vertices. Physical Review A 103(6), 062202 (2021)
  * [3] Cruz, D., Fournier, R., Gremion, F., Jeannerot, A., Komagata, K., Tosic, T., Thiesbrummel, J., Chan, C.L., Macris, N., Dupertuis, M.A., et al.: Efficient quantum algorithms for ghz and w states, and implementation on the ibm quantum computer. Advanced Quantum Technologies 2(5-6), 1900015 (2019)
  * [4] Gall, F.L., Nishimura, H., Rosmanis, A.: Quantum advantage for the local model in distributed computing. arXiv preprint arXiv:1810.10838 (2018)
  * [5] Godsil, C., Zhan, H.: Discrete Quantum Walks on Graphs and Digraphs, vol. 484. Cambridge University Press (2023)
  * [6] Izumi, T., Gall, F.L., Magniez, F.: Quantum distributed algorithm for triangle finding in the congest model. arXiv preprint arXiv:1908.11488 (2019)
  * [7] Kendon, V.: Quantum walks on general graphs. International Journal of Quantum Information 4(05), 791–805 (2006)
  * [8] Le Gall, F., Magniez, F.: Sublinear-time quantum computation of the diameter in congest networks. In: Proceedings of the 2018 ACM Symposium on Principles of Distributed Computing. pp. 337–346 (2018)
  * [9] Masuda, N., Porter, M.A., Lambiotte, R.: Random walks and diffusion on networks. Physics reports 716, 1–58 (2017)
  * [10] Roget, M., Guillet, S., Arrighi, P., Di Molfetta, G.: Grover search as a naturally occurring phenomenon. Physical Review Letters 124(18), 180501 (2020)
  * [11] Tani, S., Kobayashi, H., Matsumoto, K.: Exact quantum algorithms for the leader election problem. ACM Transactions on Computation Theory (TOCT) 4(1), 1–24 (2012)
  * [12] Zalka, C.: Grover’s quantum searching algorithm is optimal. Physical Review A 60(4), 2746 (1999)