# Black-box sampling for weakly smooth Langevin Monte Carlo using
$p$-generalized Gaussian smoothing

Anh D Doan, Xin Dang  Dao Nguyen11footnotemark: 1 Departments of Mathematics,
University of Mississippi, Oxford, MS 38655, USA<EMAIL_ADDRESS>
###### Abstract

Discretization of continuous-time diffusion processes is a widely recognized
method for sampling. However, the canonical Euler-Maruyama discretization of
the Langevin diffusion process, also named as Langevin Monte Carlo (LMC),
studied mostly in the context of smooth (gradient-Lipschitz) and strongly log-
concave densities, a significant constraint for its deployment in many
sciences, including computational statistics and statistical learning. In this
paper, we establish several theoretical contributions to the literature on
such sampling methods. Particularly, we generalize the Gaussian smoothing,
approximate the gradient using p-generalized Gaussian smoothing and take
advantage of it in the context of black-box sampling. We first present a non-
strongly concave and weakly smooth black-box LMC algorithm, ideal for
practical applicability of sampling challenges in a general setting.

###### keywords:

Langevin Monte Carlo, weakly smooth, $p$-generalized Gaussian, black-box
sampling

††journal: Journal of LaTeX Templates

##  1 Introduction

Over the last decade, Bayesian inference has become one of the most prevalent
inferring instruments for a variety of disciplines, including the
computational statistics and statistical learning [1]. In general, Bayesian
inference seeks to generate samples of the posterior distribution of the form:

$\pi(\mathrm{x})=\mathrm{e}^{-U(x)}/\int_{\mathbb{R}^{d}}\mathrm{e}^{-U(y)}\mathrm{d}y,$

where the function $U(\mathrm{x})$, also known as the potential function, is
convex. The most conventional approaches, random walks Metropolis Hasting [1],
often struggle to select a proper proposal distribution for sampling. As a
result, it has been proposed to consider continuous dynamics which inherently
leave the objective distribution $\pi$ invariant. Probably one of the most
well-known example of these continuous dynamic applications is the over-damped
Langevin diffusion [2] associated with $U$, and its Euler-Maruyama
discretization defines by

$\mathrm{x}_{k+1}=\mathrm{x}_{k}-\eta_{k}\nabla
U(\mathrm{x}_{k})+\sqrt{2\eta}\xi_{k},$ (1)

where $(\eta_{k})_{k\geq 1}$ is a sequence of step sizes which can be kept
constant or decreasing to $0$, and $\xi_{k}\sim\mathcal{N}(0,\ I_{d\times d})$
are independent Gaussian random vectors. Nonetheless, there is a critical gap
in the theory of discretization of an underlying stochastic differential
equation (SDE) to the potential broad spectrum of applications in statistical
inference. In particular, the application of techniques from SDEs
traditionally requires $U(\mathrm{x})$ to have Lipschitz-continuous gradients.
This requirement prohibits many typical utilizations [3]. [4] has recently
established an original approach to deal with weakly smooth (possibly non-
smooth) potential problems directly. Their technique rests on results obtained
from the optimization community, perturbing a gradient evaluating point by a
Gaussian. In this paper, we show that Gaussian smoothing can be generalized to
$p$-generalized Gaussian smoothing, which is novel in both Bayesian and
optimization communities. In a more general context, this novel perturbation
process provides additional features, comparable to Gaussian smoothing while
preserving the same result when a $p$-generalized Gaussian is a Gaussian. In
particular, the bias and variance are bounded for a general $p$ while these
bounds are equivalent to Gaussian smoothing for the case $p=2$. In addition,
we also improve the convergence result of [4] in Wasserstein distance by a
simpler approach. The paper is organized as follows. Section 2 sets out the
notation and context necessary to give our core theorems, generalizing the
results obtained from [4] for Gaussian smoothing. Section 3 broadens these
outcomes of stochastic approximation of the potential (negative log-density)
and composite structure of the potential, while Section 4 shows our
conclusions.

##  2 Weakly smooth LMC using $p$-generalized Gaussian smoothing

The objective is to sample from a distribution
$\pi\propto\exp(-U(\mathrm{x}))$, where $\mathrm{x}\in\mathbb{R}^{d}$. We
furnish the space $\mathbb{R}^{d}$ with the regular Euclidean norm
$\|\cdot\|=\|\cdot\|_{2}$ and use $\left\langle\ ,\ \right\rangle$ to specify
inner products. While sampling from the exact distribution $\pi(\mathrm{x})$
is generally computationally demanding, it is largely adequate to sample from
an approximated distribution $\tilde{\pi}(\mathrm{x})$ which is in the
vicinity of $\pi(\mathrm{x})$ by some distances. In this paper, we use
Wasserstein distance and briefly define it in Appendix A. Suppose the
following conditions, identical to the first two assumptions of [4]:

A.1 $U$ is convex and sub-differentiable. Specifically, for all
$\mathrm{x}\in\mathbb{R}^{d}$, there exists a sub-gradient of $U,\ \nabla
U(\mathrm{x})\in\partial U(\mathrm{x})$, ensuring that
$\forall\mathrm{y}\in\mathbb{R}^{d}:U(\mathrm{y})\geq
U(\mathrm{x})+\left\langle\nabla U(\mathrm{x}),y-x\right\rangle.$

A.2 There exist $L<\infty$ and $\alpha\in[0,1]$ so that $\forall\mathrm{x},\
\mathrm{y}\in\mathbb{R}^{d}$, we obtain

$\|\nabla U(\mathrm{x})-\nabla U(\mathrm{y})\|_{2}\leq
L\|\mathrm{x}-\mathrm{y}\|_{2}^{\alpha},$ (2)

where $\nabla U(\mathrm{x})$ represents an arbitrary sub-gradient of $U$ at
$\mathrm{x}$.

######  Remark 1.

Condition 2 is known as $(L,\alpha)$-weakly smoothness or Holder continuity of
the (sub)gradients of $U$. A feature that follows straightfowardly from Eq.(2)
is that:

$U({\mathrm{y})\leq U(\mathrm{x})+\left\langle\nabla U(\mathrm{x}),\
\mathrm{y}-\mathrm{x}\right\rangle+\frac{L}{1+\alpha}\|\mathrm{y}-\mathrm{x}\|^{1+\alpha},\
\forall\mathrm{x},\ \mathrm{y}\in\mathrm{R^{d}}}.$ (3)

When $\alpha=1$, it is equivalent to the standard smoothness (Lipschitz
continuity of the gradients), whereas at the opposite extreme, when
$\alpha=0,\ U$ is (possibly) non-smooth and Lipschitz-continuous.

[4] perturbs the gradient evaluating point by a Gaussian so that they can
utilize weakly smooth Lipschitz potentials directly without any additional
structure. Here, we generalize their approach to perturbation using
$p$-generalize Gaussian smoothing. Particularly, consider $\mu\geq 0$,
$p$-generalized Gaussian smoothing $U_{\mu}$ of $U$ is defined as:

$U_{\mu}(\mathrm{y}):=\mathrm{E}_{\xi}[U(\mathrm{y}+\mu\xi)]=\frac{1}{\kappa}\int_{\mathbb{R}^{d}}U(\mathrm{y}+\mu\xi)e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi,$

where
$\kappa\stackrel{{{}_{def}}}{{=}}\int_{\mathbb{R}^{d}}e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi=\frac{2^{d}\Gamma^{d}(\frac{1}{p})}{p^{d-\frac{d}{p}}}$
and $\xi\sim N_{p}(0,I_{d\times d})$ (the $p$-generalized Gaussian
distribution). The rationale for taking into account the $p$-generalized
Gaussian smoothing $U_{\mu}$ rather than $U$ is that it typically benefits
from superior smoothness properties. In particular, $U_{\mu}$ is smooth albeit
$U$ is not. In addition, $p$-generalized Gaussian smoothing is more
generalized than Gaussian smoothing in the sense that it contains all normal
distribution when $p=2$ and all Laplace distribution when $p=1$. This family
of distributions allows for tails that are either heavier or lighter than
normal and in the limit as well as containing all the continuous uniform
distribution. It is possible to study $p\in\mathbb{R}^{+}$ but to simplify the
proof, we are only interested in $p\in[1,2]$. Here we examine some primary
features of $U_{\mu}$ based on adapting those results of [4].

######  Lemma 1.

If $U:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a convex function which
satisfies Eq.(2) for some $L<\infty$ and $\alpha\in[0,1]$, then:

(i) $\forall\mathrm{x}\in\mathbb{R}^{d}$ :
$|U_{\mu}({\displaystyle\mathrm{x})-U(\mathrm{x})|=U_{\mu}(\mathrm{x})-U(\mathrm{x})\leq\frac{L\mu^{1+\alpha}d^{\frac{1+\alpha}{p}}}{1+\alpha}}$.

(ii) $\forall\mathrm{x},\ \mathrm{y}\in\mathbb{R}^{d}$:
${\displaystyle\|\nabla U_{\mu}(\mathrm{y})-\nabla
U_{\mu}(\mathrm{x})\|_{2}\leq\frac{Ld^{\frac{1-\alpha}{p}}}{\mu^{1-\alpha}(1+\alpha)^{1-\alpha}}\|\mathrm{y}-\mathrm{x}\|_{2}}$.

###### Proof.

See Appendix B.1 ∎

##  3 Black-box sampling for regularized potential

To study the convergence of the continuous-time process (involving strong
convexity), we work with regularized potentials that have the following
composite structure:

$\overline{U}(\mathrm{x}):=U(\mathrm{x})+\frac{\lambda}{2}\left\|x\right\|^{2},$
(4)

where $\lambda>0$. Observe that by the triangle inequality:

$\displaystyle\|\nabla\overline{U}(\mathrm{x})-\nabla\overline{U}(\mathrm{y})\|_{2}$
$\displaystyle\leq\|\nabla U(\mathrm{x})-\nabla
U(\mathrm{y})\|_{2}+\lambda\|\mathrm{x}-\mathrm{y}\|_{2}$ (5)
$\displaystyle\leq
L\|\mathrm{x}-\mathrm{y}\|_{2}^{\alpha}+\lambda\|\mathrm{x}-\mathrm{y}\|_{2}.$

From which, by employing Lemma 1, we get:

${\displaystyle\|\nabla\overline{U}_{\mu}(\mathrm{x})-\nabla\overline{U}_{\mu}(\mathrm{y})\|_{2}\leq\left(\frac{Ld^{\frac{1-\alpha}{p}}}{\mu^{1-\alpha}(1+\alpha)^{1-\alpha}}+\lambda\right)\|\mathrm{x}-\mathrm{y}\|_{2}.}$
(6)

In this case, $U(\cdot)$ is $(L,\alpha)$-weakly smooth (possibly with
$\alpha=0$, $U$ is non-smooth and Lipschitz-continuous). We now analyze LMC
where we perturb the points at which gradients of $\overline{U}$ are evaluated
by a $p$-generalized Gaussian random variable. Note that it remains unclear
whether it is possible to achieve such bounds for (LMC) without this
perturbation. Recall that LMC in terms of the potential $\overline{U}$ can be
specified as:

$\mathrm{x}_{k+1}=\mathrm{x}_{k}-\eta\nabla\overline{U}_{\mu}(\mathrm{x}_{k})+\sqrt{2\eta}\varsigma_{k},$
(7)

where $\varsigma_{k}\sim N(0,\ I_{d\times d})$ are independent Gaussian random
vectors. This method is actually the Euler-Maruyama discretization of the
Langevin diffusion. Rather than working with the algorithm specified by Eq. 7,
we define an estimate $g_{\mu,n}(\theta)$, of the gradient
$\nabla\overline{U}_{\mu}(\theta)$ as follows:

$\displaystyle
g_{\mu,n}(x)=\frac{1}{n}\sum_{i=1}^{n}\frac{\overline{U}_{\mu}(x+\mu\xi_{i})-\overline{U}_{\mu}(x)}{\mu}\left(\xi_{i}\circ\left|\xi_{i}\right|^{p-2}\right)$
(8)

where $\xi_{i}\sim N_{p}(0,I_{d})$ and are i.i.d. and $\circ$ stands for the
Hadamard product and $\left|\cdot\right|$ is used for absolute value of each
component of the vector $\xi_{i}$. The gradient estimator in Equation 8 is a
consequence of a p-generalized Gaussian identity we derive in the proof of
Lemma 1. Based on the above estimate of the gradient, we obtain the following
result.

######  Lemma 2.

For any $\mathrm{x_{k}}\in\mathbb{R}^{d}$, then
$g_{\mu,n}(x_{k})=\nabla\overline{U}_{\mu}(x_{k})+\zeta_{k}$ is an estimator
of $\nabla\overline{U}_{\mu}$ such that

$\displaystyle\|E[\zeta_{k}\mid x_{k}]\|^{2}$
$\displaystyle\leq\left(M+\lambda\right)^{2}\mu^{2}d^{\frac{2}{p}}$
$\displaystyle~{}E\left[\left\|\zeta_{k}-E[\zeta_{k}\mid
x_{k}]\right\|^{2}\right]$
$\displaystyle\leq\frac{1}{n}\left(\frac{1}{2}\left(M+\lambda\right)\mu\left(d+3\right){}^{3/p}+\sqrt{2}\left(d+2\right)^{\frac{2}{p}}\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|\right)^{2}$

###### Proof.

See Appendix B.2 ∎

Let the distribution of the $k^{th}$ iterate $\mathrm{y}_{k}$ be represented
by $\overline{\pi}_{\mu,k}$, and let
$\overline{\pi}_{\mu}\propto\exp(-\overline{U}_{\mu})$ be the distribution
with $\overline{U}_{\mu}$ as the potential. It is straightfoward that
$p$-generalized Gaussian smoothing of regularized potential remains strong
convexity so we omit the proof. Our overall tactics for showing our core
results are as follows. First, we prove that the $p$-generalized Gaussian
smoothing does not alter the objective distribution substantially in term of
the Wasserstein distance, by bounding
$W_{2}(\overline{\pi},\overline{\pi}_{\mu})$ (Lemma 3). We then deploy a
result on mixing times of Langevin diffusion with stochastic gradients, which
enables us to bound $W_{2}(\overline{\pi}_{\mu,k},\overline{\pi}_{\mu})$.

######  Lemma 3.

Assume that $\overline{\pi}\propto\exp(-\overline{U})$ and
$\overline{\pi}_{\mu}\propto\exp(-\overline{U}_{\mu})$. We deduce the
following bounds, for any $\lambda\geq 0$

$W_{2}^{2}(\overline{\pi},\
\overline{\pi}_{\mu})\leq\frac{4(d+\lambda\|x^{*}\|_{2}^{2})}{\lambda}\left(a+e^{a}-1\right).$

where $x^{*}$ is a unique minimizer of $\overline{U}$,
$a=\frac{L\mu^{1+\alpha}d^{\frac{1+\alpha}{p}}}{1+\alpha}+\frac{1}{2}\lambda\mu^{2}\left(d+1\right)^{\frac{2}{p}}$.
If, in addition, $a\leq 0.1$ and for sufficient small $\lambda$, then

$W_{2}(\overline{\pi},\ \overline{\pi}_{\mu})\leq 3\sqrt{\frac{da}{\lambda}},$

###### Proof.

See Appendix B.3 ∎

Our primary outcome is reported in the subsequent theorem.

######  Theorem 1.

Let the initial iterate $\mathrm{y}_{0}$ be drawn from a probability
distribution $\overline{\pi}_{\mu,0}$. If the step size $\eta$ satisfies
$\eta<2/(M+2\lambda)$ and for sufficiently small $\lambda$, then:

$\displaystyle
W_{2}(\overline{\pi}_{\mu,K},{\pi})\leq(1-0.5\lambda\eta)^{K/2}W_{2}(\overline{\pi}_{\mu,0},\overline{\pi}_{\mu})+1.9\frac{\left(M+\lambda\right)}{\lambda}(\eta
d)^{1/2}+2\frac{\left(M+\lambda\right)}{\lambda}\mu
d^{1/p}+3\sqrt{\frac{da}{\lambda}}$
$\displaystyle+\frac{1}{\sqrt{\lambda}}\cdot\frac{1}{\sqrt{n}}\eta^{1/2}\left(M+\lambda\right)\mu(d+3)^{3/p}+\frac{\sqrt{\left(M+\lambda\right)}}{\sqrt{\lambda}}\cdot\frac{1}{\sqrt{n}}\eta^{1/2}d^{1/2}(d+2)^{1/p}+C\lambda\left(dlogd\right)^{4},$

where ${\displaystyle
M=\frac{Ld^{\frac{1-\alpha}{p}}}{\mu^{1-\alpha}(1+\alpha)^{1-\alpha}}}$,
${\displaystyle
a=\frac{L\mu^{1+\alpha}d^{\frac{1+\alpha}{p}}}{1+\alpha}+\frac{1}{2}\mu^{2}\lambda(d+1)^{\frac{2}{p}}}$

###### Proof.

Adapting the technique of analysis of [2, 5], by a triangle inequality, the
Wasserstein distance between $\overline{\pi}_{\mu,K}$ and $\pi$ is bounded by:

$W_{2}(\overline{\pi}_{\mu,K},\pi)\leq
W_{2}(\overline{\pi}_{\mu,K},\,\overline{\pi}_{\mu})+W_{2}(\overline{\pi}_{\mu},\overline{\pi})+W_{2}(\overline{\pi},\,\pi)$.

We use [2] Theorem 9 to bound the first term
$W_{2}(\overline{\pi}_{K},\,\overline{\pi}_{\mu})$. Recall that
$\overline{U}_{\mu}$ is continuously differentiable, $(M+\lambda)$-smooth
(with
$M=\frac{Ld^{(\frac{1-\alpha}{p})}}{\mu^{1-\alpha}(1+\alpha)^{1-\alpha}}$) and
$\lambda$-strongly convex. Additionally, $\\{\mathrm{y}_{k}\\}_{k=1}^{K}$ can
be viewed as iterates of over-damped Langevin MCMC with respect to the
potential specified by $\overline{U}_{\mu}$ and is updated using biased noisy
gradients of $\overline{U}_{\mu}$. For any random vector $\boldsymbol{X}$, we
define the norm
$\|\boldsymbol{X}\|_{L_{2}}=(E[\|\boldsymbol{X}\|_{2}^{2}])^{1/2}$. Let
$L_{0}$ be a random vector drawn from $\overline{\pi}_{\mu}$ such that
$W_{2}(\nu_{k},\overline{\pi}_{\mu})=\|L_{0}-x_{k}\|_{L_{2}}$ and
$E[\zeta_{k}|x_{k},L_{0}]=E[\zeta_{k}|x_{k}]$. Let $W$ be a d-dimensional
Brownian motion independent of $(x_{k},L_{0},\zeta_{k})$, such that
$W_{\eta}=\sqrt{\eta}\,\zeta_{k+1}$. We define the stochastic process $L$ so
that
$L_{t}=L_{0}-\int_{0}^{t}\nabla\overline{U}_{\mu}(L_{s})\,ds+\sqrt{2}\,W_{t},\
\forall\,t>0$. It is clear that this equation implies that ${\displaystyle
L_{\eta}}=L_{0}-\int_{0}^{\eta}\nabla\overline{U}_{\mu}(L_{s})\,ds+\sqrt{2}\,W_{\eta}=L_{0}-\int_{0}^{\eta}\nabla\overline{U}_{\mu}(L_{s})\,ds+\sqrt{2\eta}\,{\zeta}_{k+1}$.
Furthermore, $\\{L_{t}:t\geq 0\\}$ is a diffusion process having
$\overline{\pi}_{\mu}$ as the stationary distribution. Since the initial value
$L_{0}$ is drawn from $\overline{\pi}_{\mu}$, we have
$L_{t}\sim\overline{\pi}_{\mu}$ for every $t\geq 0$.

Let us define $\Delta_{k}=L_{0}-x_{k},\;\Delta_{k+1}=L_{\eta}-x_{k+1}$ and
$V=\int_{0}^{\eta}[\nabla\overline{U}_{\mu}(L_{s})-\nabla\overline{U}_{\mu}(L_{0})]ds$
then
$\|V\|_{L_{2}}\leq\frac{1}{2}\big{(}\eta^{4}(\lambda+M)^{3}d\big{)}^{1/2}+\frac{2}{3}(2\eta^{3}d)^{1/2}(\lambda+M)\leq.95(\lambda+M)\eta^{\frac{3}{2}}d^{\frac{1}{2}}$
for $\eta$ small enough by Lemma 4 in [2]. Thus,

$\|\Delta_{k+1}\|_{L_{2}}=\|\Delta_{k}-\eta[\nabla\overline{U}_{\mu}(x_{k}+\Delta_{k})-\nabla\overline{U}_{\mu}(x_{k})]-V+\eta\zeta_{k}\|_{L_{2}}$
$\leq\\{\|\Delta_{k}-\eta\lambda\Delta_{k}\|_{L_{2}}^{2}+\eta^{2}\|\zeta_{k}-E[\zeta_{k}\mid
x_{k}]\|_{L_{2}}^{2}\\}^{1/2}+\|V\|_{L_{2}}+\eta\|E[\zeta_{k}\mid
x_{k}]\|_{L_{2}}$
$\leq\left\\{(1-\lambda\eta)^{2}\|\Delta_{k}\|_{L_{2}}^{2}+\frac{\eta^{2}}{n}\left(\frac{1}{2}\left(M+\lambda\right)\mu(d+3)^{3/p}+(\|\nabla\overline{U}_{\mu}(L_{0})\|_{L_{2}}+\left(M+\lambda\right)\|\Delta_{k}\|_{L_{2}})\sqrt{2}\left(d+2\right)^{1/p}\right)^{2}\right\\}^{1/2}$
$\qquad+.95\left(M+\lambda\right)(\eta^{3}d)^{1/2}+\left(M+\lambda\right)\mu\eta
d{}^{1/p}$
$\leq\left\\{(1-\lambda\eta)^{2}\|\Delta_{k}\|_{L_{2}}^{2}+\frac{2\eta^{2}}{n}\left(\frac{1}{2}\left(M+\lambda\right)\mu(d+3)^{3/p}+\|\nabla\overline{U}_{\mu}(L_{0})\|_{L_{2}}\sqrt{2}\left(d+2\right)^{1/p}\right)^{2}\right\\}^{1/2}$
$+\frac{2\left(M+\lambda\right)^{2}\eta^{2}\left(d+2\right)^{\frac{2}{p}}}{n(1-\lambda\eta)}\|\Delta_{k}\|_{L_{2}}+.95\left(M+\lambda\right)(\eta^{3}d)^{1/2}+\left(M+\lambda\right)\mu\eta
d^{1/p}$
$\leq\left\\{(1-\lambda\eta)^{2}\|\Delta_{k}\|_{L_{2}}^{2}+\frac{2\eta^{2}}{n}\left(\frac{1}{2}\left(M+\lambda\right)\mu(d+3)^{3/p}+\sqrt{2\left(M+\lambda\right)d}\left(d+2\right)^{1/p}\right)^{2}\right\\}^{1/2}$
$+\frac{\sqrt{2}\left(M+\lambda\right)^{2}\eta^{2}\left(d+2\right)^{2/p}}{n(1-\lambda\eta)}\|\Delta_{k}\|_{L_{2}}+.95\left(M+\lambda\right)(\eta^{3}d)^{1/2}+\left(M+\lambda\right)\mu\eta
d^{1/p}.$

Here we use the fact that $\sqrt{a^{2}+b+c}\leq\sqrt{a^{2}+b}+\frac{c}{2a}$
and $E[\|\nabla\overline{U}_{\mu}(X_{s})\|^{2}]\leq\left(M+\lambda\right)d$.
By Lemma 9 in [2], for $n$ large enough,
$\frac{\sqrt{2}\left(M+\lambda\right)^{2}\eta^{2}\left(d+2\right)^{2/p}}{n(1-\lambda\eta)}\leq
0.5\lambda\eta$, the above inequality leads to

$\displaystyle\|\Delta_{k}\|_{L_{2}}$
$\displaystyle\leq\left(1-0.5\lambda\eta\right)^{k}\|\Delta_{0}\|_{L_{2}}+1.9\frac{\left(M+\lambda\right)}{\lambda}(\eta
d)^{1/2}+2\frac{\left(M+\lambda\right)}{\lambda}\mu d^{1/p}$
$\displaystyle\qquad\frac{1}{\sqrt{\lambda}}\cdot\frac{1}{\sqrt{n}}\eta^{1/2}\left(M+\lambda\right)\mu(d+3)^{3/p}+\frac{\sqrt{\left(M+\lambda\right)}}{\sqrt{\lambda}}\cdot\frac{1}{\sqrt{n}}\eta^{1/2}d^{1/2}(d+2)^{1/p}.$

Therefore, we obtain the bound in Wasserstein distance.

$\displaystyle W_{2}(\overline{\pi}_{\mu,k},\overline{\pi}_{\mu})$
$\displaystyle\leq\left(1-0.5\lambda\eta\right)^{k}W_{2}(\overline{\pi}_{\mu,0},\overline{\pi}_{\mu})+1.9\frac{\left(M+\lambda\right)}{\lambda}(\eta
d)^{1/2}+2\frac{\left(M+\lambda\right)}{\lambda}\mu d^{1/p}$
$\displaystyle\qquad\frac{1}{\sqrt{\lambda}}\cdot\frac{1}{\sqrt{n}}\eta^{1/2}\left(M+\lambda\right)\mu(d+3)^{3/p}+\frac{\sqrt{\left(M+\lambda\right)}}{\sqrt{\lambda}}\cdot\frac{1}{\sqrt{n}}\eta^{1/2}d^{1/2}(d+2)^{1/p}.$

By [6] Propositions 1 and 2,

$W_{2}^{2}(\overline{\pi},\pi)\leq 111\lambda m_{2}(\pi)^{2},$

where
$m_{2}(\pi):=\int_{\mathbb{R}^{d}}\|x-x^{*}\|_{2}^{2}\,\pi(x)\,dx=C\left(dlogd\right)^{2}$
for some constant $C$ independent of $d.$ Combining with Lemma 3 produces the
desired result. ∎

######  Remark 2.

When taking the limit as $\mu$ goes to $0$, our gradient estimate is unbiased,
our bound is equivalent to
$O(d^{\frac{2+p}{2p}}+\frac{d^{\frac{1+\alpha+p}{2p}}}{\sqrt{\lambda}})$ for
non-regularization potential. In comparison to the previous result of [4],
which is of order $O(\frac{d}{\lambda})$ when considered Lipschitz constant
$L$ as a constant and $p=2$, our approach has better dependencies on strongly
convex parameter $\lambda$. In general, our result is strictly better whenever
$p<2$. In addition, our method is more generalized than Gaussian smoothing,
especially when we desire to explore the space by heavier tail distribution
than normal.

##  4 Conclusion

We derive polynomial-time theoretical assurances for an LMC that uses
$p$-generalized Gaussian smoothing. The algorithm we proposed is a
generalization of the recent Gaussian smoothing LMC method. It is potential to
broaden our results to sampling from distributions with non-smooth and non-
convex structure or integrate into derivative-free LMC algorithm. In addition,
we speculate that the dependence on $d$ and $\lambda$ is not optimal and can
be improved to match those obtained for the 2-Wasserstein distance using
proximal assumption.

##  Appendix A Distance measures

Define a transference plan $\zeta$, a distribution on
$(\mathbb{R}^{d}\times\mathbb{R}^{d},\
\mathcal{B}(\mathbb{R}^{d}\times\mathbb{R}^{d}))$ (where
$\mathcal{B}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ is the Borel $\sigma$-field
of ($\mathbb{R}^{d}\times\mathbb{R}^{d}$)) so that
$\zeta(A\times\mathbb{R}^{d})=P(A)$ and $\zeta(\mathbb{R}^{d}\times A)=Q(A)$
for any $A\in\mathcal{B}(\mathbb{R}^{d})$. Let $\Gamma(P,\ Q)$ designate the
set of all such transference plans. Then the 2-Wasserstein distance is
formulated as:

$W_{2}(P,Q)\
:=\left(\inf_{\zeta\in\Gamma(P,Q)}\int_{\mathrm{x},\mathrm{y}\in\mathbb{R}^{d}}\|\mathrm{x}-\mathrm{y}\|_{2}^{2}\mathrm{d}\zeta(\mathrm{x},\
\mathrm{y})\right)^{1/2}.$

##  Appendix B Proofs

###  B.1 Proof of Lemma 1

###### Proof.

(i). Since $U_{\mu}(\mathrm{x})=\mathrm{E}_{\xi}[U(\mathrm{x}+\mu\xi)]$,
$U(\mathrm{x})=\mathrm{E}_{\xi}[U(x)]$ and $E_{\xi}(\xi)=0$, which implies
$E_{\xi}\mu\left\langle\nabla U(\mathrm{x}),\ \xi\right\rangle=0$, we obtain

$U_{\mu}(\mathrm{x})-U(\mathrm{x})=E\left[U(\mathrm{x}+\mu\xi)-U(\mathrm{x})-\mu\left\langle\nabla
U(\mathrm{x}),\ \xi\right\rangle\right].$

Because $U$ is convex and $\mu>0$, from A.1,
$U(\mathrm{x}+\mu\xi)-U(\mathrm{x})-\mu\left\langle\nabla U(\mathrm{x}),\
\xi\right\rangle\geq 0$ for every $\xi$ and $x$, so $U_{\mu}(\mathrm{x})\geq
U(\mathrm{x}),\ \forall\mathrm{x}$. By the definition of the density of
$p$-generalized Gaussian distribution, we also have:

$U_{\mu}(\mathrm{x})-U(\mathrm{x})=\frac{1}{\kappa}\int_{\mathbb{R}^{d}}[U(\mathrm{x}+\mu\xi)-U(\mathrm{x})-\mu\left\langle\nabla
U(\mathrm{x}),\
\xi\right\rangle]e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi.$

Applying Eq. 3 and previous inequality:

$\displaystyle|U_{\mu}(\mathrm{x})-U(\mathrm{x})|$
$\displaystyle\leq\frac{L}{\kappa(1+\alpha)}\mu^{1+\alpha}\int_{\mathbb{R}^{d}}\left\|\xi\right\|_{2}^{(1+\alpha)}e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi$
$\displaystyle\leq\frac{L}{\kappa(1+\alpha)}\mu^{1+\alpha}\int_{\mathbb{R}^{d}}\left\|\xi\right\|_{p}^{(1+\alpha)}e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi$
$\displaystyle=\frac{L}{\kappa(1+\alpha)}\mu^{1+\alpha}\int_{\mathbb{R}^{d}}\left\|\xi\right\|_{p}^{2p(1+\alpha)/(2p)}e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi.$

The second inequality follows from
$\left\|\xi\right\|_{p}\geq\left\|\xi\right\|_{2}$ for any $p\leq 2$. Since
$1\leq p$ we have $1+\alpha\leq 2p$ and $x^{\frac{1+\alpha}{2p}}$ is concave,
it follows that

$\displaystyle U_{\mu}(\mathrm{x})-U(\mathrm{x})$
$\displaystyle\leq\frac{L}{1+\alpha}\mu^{1+\alpha}E\left[\left\|\xi\right\|_{p}^{\frac{2p(1+\alpha)}{2p}}\right]$
$\displaystyle\leq\frac{L}{1+\alpha}\mu^{1+\alpha}\left[E\left\|\xi\right\|_{p}^{2p}\right]^{\frac{1+\alpha}{2p}}$
$\displaystyle\leq\frac{L}{1+\alpha}\mu^{1+\alpha}\left[\frac{2p\Gamma\left(\frac{d}{p}+2\right)}{\Gamma\left(\frac{d}{p}\right)}\right]^{\frac{1+\alpha}{2p}}$
$\displaystyle=\frac{L}{1+\alpha}\mu^{1+\alpha}\left(\frac{2pd\left(d+p\right)}{p^{2}}\right)^{\frac{1+\alpha}{2p}}$
$\displaystyle=\frac{L}{1+\alpha}\mu^{1+\alpha}\left(\frac{2d\left(d+p\right)}{p}\right)^{\frac{1+\alpha}{2p}}.$

Finally, for sufficiently large $d$, we have

$|U_{\mu}(\mathrm{x})-U(\mathrm{x})|\leq\frac{L}{1+\alpha}\mu^{1+\alpha}d^{\frac{1+\alpha}{p}}.$

(ii). We adapt the technique of [7] to $p$-generalized Gaussian smoothing. Let
$y=\mathrm{x}+\mu\xi$, then $U_{\mu}(\mathrm{x})$ is rewritten in another form
as

$U_{\mu}(\mathrm{x})=\frac{1}{\mu^{d}\kappa}\int_{\mathbb{R}^{d}}U(y)e^{-\frac{1}{p\mu^{p}}\left\|y-x\right\|_{p}^{p}}\mathrm{d}y.$

Now taking the gradient with respect to $x$ of $U_{\mu}(\mathrm{x})$ gives

$\nabla_{x}U_{\mu}(\mathrm{x})=\frac{1}{\mu^{d}\kappa}\nabla_{x}\int_{\mathbb{R}^{d}}U(y)e^{-\frac{1}{p\mu^{p}}\left\|y-x\right\|_{p}^{p}}\mathrm{d}y.$

By Fubini Theorem with some regularity (i.e. $E|U(\mathrm{y})|<\infty$), we
can exchange the gradient and integral and get

$\displaystyle\nabla_{x}U_{\mu}(\mathrm{x})$
$\displaystyle=\frac{1}{\mu^{d}\kappa}\int_{\mathbb{R}^{d}}\nabla_{x}\left(U(y)e^{-\frac{1}{p\mu^{p}}\left\|y-x\right\|_{p}^{p}}\right)\mathrm{d}y$
$\displaystyle=\frac{1}{\mu^{d}\kappa}\int_{\mathbb{R}^{d}}U(y)\nabla_{x}\left(e^{-\frac{1}{p\mu^{p}}\left\|y-x\right\|_{p}^{p}}\right)\mathrm{d}y$
$\displaystyle=\frac{1}{\mu^{d}\kappa}\int_{\mathbb{R}^{d}}U(y)e^{-\frac{1}{p\mu^{p}}\left\|y-x\right\|_{p}^{p}}\frac{1}{\mu^{p}}\mathrm{\left\|y-x\right\|_{p}^{p-1}\nabla_{x}\left(\left\|y-x\right\|_{p}\right)d}y.$
$\displaystyle=-\frac{1}{\mu^{d}\kappa}\int_{\mathbb{R}^{d}}U(y)e^{-\frac{1}{p\mu^{p}}\left\|y-x\right\|_{p}^{p}}(y-x)\circ\left|y-x\right|^{p-2}dy.$

where $\circ$ stands for the Hadamard product and $\left|\cdot\right|$ is used
for absolute value of each component of the vector $y-x$. Thereby, by changing
variable back to $\xi$, we deduce

$\displaystyle\nabla_{x}U_{\mu}(\mathrm{x})$
$\displaystyle=-\frac{1}{\mu\kappa}\int_{\mathbb{R}^{d}}U(\mathrm{x}+\mu\xi)e^{-\frac{1}{p}\left\|\xi\right\|_{p}^{p}}\xi\circ\left|\xi\right|^{p-2}\mathrm{d}\xi.$
$\displaystyle=-E_{\xi}\left[U(\mathrm{x}+\mu\xi)\xi\circ\left|\xi\right|^{p-2}\right]$

Thus, by Jensen’s inequality:

${\displaystyle\|\nabla U_{\mu}(\mathrm{y})-\nabla
U_{\mu}(\mathrm{x})\|_{2}}\leq\frac{1}{\mu\kappa}\int_{\mathbb{R}^{d}}\left\|U(\mathrm{y}+\mu\xi)-U(\mathrm{x}+\mu\xi)\right\|_{2}\left\|\xi\circ\left|\xi\right|^{p-2}\right\|_{2}e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi.$
(9)

Further, by using generalized Holder inequality for $1\leq p\leq 2$,
$\left\|\xi\circ\left|\xi\right|^{p-2}\right\|_{2}$ can be bounded as:

$\left\|\xi\circ\left|\xi\right|^{p-2}\right\|_{2}=\left\|\xi^{p-1}\right\|_{2}\\\
\leq\left\|\xi^{p-1}\right\|_{p}\\\ =\left\|\xi^{p-1}\cdot
1_{d}\right\|_{p}\\\
\leq\left\|\xi\right\|_{p}^{p-1}\left\|1_{d}\right\|_{p}^{2-p}\\\
=\left\|\xi\right\|_{p}^{p-1}d^{\frac{2-p}{p}}.$

Thus, applying Jensen’s inequality, we derive:

${\displaystyle\|\nabla U_{\mu}(\mathrm{y})-\nabla
U_{\mu}(\mathrm{x})\|_{2}\leq\frac{1}{\mu}d^{\frac{2-p}{p}}\left(\frac{p^{1-\frac{1}{p}}}{2\Gamma(\frac{1}{p})}\right)^{d}\int_{\mathbb{R}^{d}}|U(\mathrm{x}+\mu\xi)-U(\mathrm{y}+\mu\xi)|\left\|\xi\right\|_{p}^{p-1}e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi.}$
(10)

Using Eq. 3, we obtain:

$\displaystyle|U(\mathrm{x}+\mu\xi)-U(\mathrm{y}+\mu\xi)|$
$\displaystyle\leq\min\left\\{\left\langle\nabla
U(\mathrm{y}+\mu\xi),\mathrm{x}-\mathrm{y}\right\rangle+\frac{L}{1+\alpha}\|\mathrm{y}-\mathrm{x}\|_{2}^{1+\alpha},\right.$
$\displaystyle\left.\left\langle\nabla
U(\mathrm{x}+\mu\xi),\mathrm{y}-\mathrm{x}\right\rangle+\frac{L}{1+\alpha}\|\mathrm{y}-\mathrm{x}\|_{2}^{1+\alpha}\right\\}$
$\displaystyle\leq\frac{1}{2}\left\langle\nabla U(\mathrm{y}+\mu\xi)-\nabla
U(\mathrm{x}+\mu\xi),\
\mathrm{x}-\mathrm{y}\right\rangle+\frac{L}{1+\alpha}\|\mathrm{y}-\mathrm{x}\|_{2}^{1+\alpha}$
$\displaystyle\leq\frac{L}{1+\alpha}\|\mathrm{y}-\mathrm{x}\|_{2}^{1+\alpha},$

where the second inequality comes from the minimum being smaller than the
mean, and the last inequality is achieved by convexity of $U$ (which implies
$\left\langle\nabla U(\mathrm{x})-\nabla U(\mathrm{y}),x-y\right\rangle\geq
0$, $\forall\mathrm{x}$, $y$). Thus, combining with Eq. 10, we induce:

$\displaystyle\|\nabla U_{\mu}(\mathrm{y})-\nabla U_{\mu}(\mathrm{x})\|_{2}$
$\displaystyle\leq\frac{L}{\mu(1+\alpha)}d^{\frac{2-p}{p}}\|\mathrm{y}-\mathrm{x}\|_{2}^{1+\alpha}\left(\frac{p^{1-\frac{1}{p}}}{2\Gamma(\frac{1}{p})}\right)^{d}\int_{\mathbb{R}^{d}}\left\|\xi\right\|_{p}^{p-1}e^{-\left\|\xi\right\|_{p}^{p}/p}\mathrm{d}\xi$
$\displaystyle=\frac{L}{\mu(1+\alpha)}d^{\frac{2-p}{p}}\|\mathrm{y}-\mathrm{x}\|_{2}^{1+\alpha}\left[\frac{p^{\frac{p-1}{p}}\Gamma(\frac{n+p-1}{p})}{\Gamma(\frac{n}{p})}\right]$
(11)
$\displaystyle\leq\frac{L}{\mu(1+\alpha)}\|\mathrm{y}-\mathrm{x}\|_{2}^{1+\alpha}d^{\frac{1}{p}}.$
(12)

Finally, combining A.2 and (12):

$\|\nabla U_{\mu}(\mathrm{y})-\nabla U_{\mu}(\mathrm{x})\|_{2}=\|\nabla
U_{\mu}(\mathrm{y})-\nabla U_{\mu}(\mathrm{x})\|_{2}^{\alpha}\cdot\|\nabla
U_{\mu}(\mathrm{y})-\nabla U_{\mu}(\mathrm{x})\|_{2}^{1-\alpha}$
$\leq\frac{Ld^{(\frac{1-\alpha}{p})}}{\mu^{1-\alpha}(1+\alpha)^{1-\alpha}}\|\mathrm{y}-\mathrm{x}\|_{2},$

as claimed. ∎

###  B.2 Proof of Lemma 2

###### Proof.

Under smooth assumption,

$\displaystyle\|E[\zeta_{k}\mid x_{k}]\|^{2}$
$\displaystyle=\|E[g_{\mu,n}(x_{k})-\nabla\overline{U}_{\mu}(x_{k})\mid
x_{k}]\|^{2}$
$\displaystyle\leq\|E[g_{\mu,1}(x_{k})-\nabla\overline{U}_{\mu}(x_{k})\mid
x_{k}]\|^{2}$
$\displaystyle=\left\|E\left[\nabla\overline{U}_{\mu}(x_{k}+\mu\xi)|x_{k}\right]-\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}$
$\displaystyle\leq
E\left[\left(M+\lambda\right)\mu\left\|\xi\right\|\right]^{2}$
$\displaystyle=\left(M+\lambda\right)^{2}\mu^{2}E^{2}\left(\left\|\xi\right\|\right)$
$\displaystyle\stackrel{{{}_{1}}}{{\leq}}\left(M+\lambda\right)^{2}\mu^{2}d^{\frac{2}{p}},$

where $1$ follows from Lemma 4 below.

For the variance, by independence of p-generalized Gaussian sample
$\xi_{i}$’s, we have the following observation.

$\displaystyle E\left[\left\|\zeta_{k}-E[\zeta_{k}\mid
x_{k}]\right\|^{2}\right]$
$\displaystyle=E\left[\left\|g_{\mu,n}(x_{k})-\overline{U}_{\mu}(x_{k})\right\|^{2}\right]$
$\displaystyle=\frac{1}{n}E\left[\left\|g_{\mu,1}(x_{k})-\overline{U}_{\mu}(x_{k})\right\|^{2}\right]$

In order to obtain the variance of $g_{\mu,1}(x_{k})$, we split the centered
error term into three parts.

$\displaystyle\zeta_{k}-E[\zeta_{k}\mid x_{k}]$
$\displaystyle=\frac{\overline{U}_{\mu}(x_{k}+\mu\xi)-\overline{U}_{\mu}(x_{k})}{\mu}\left(\xi\circ\left|\xi\right|^{p-2}\right)-\nabla\overline{U}_{\mu}(x_{k})$
$\displaystyle=\frac{\overline{U}_{\mu}(x_{k}+\mu\xi)-\overline{U}_{\mu}(x_{k})-\left\langle\nabla\overline{U}_{\mu}(x_{k}),\mu\xi\right\rangle}{\mu}\left(\xi\circ\left|\xi\right|^{p-2}\right)+\left\langle\nabla\overline{U}_{\mu}(x_{k}),\xi\right\rangle\left(\xi\circ\left|\xi\right|^{p-2}\right)-\nabla\overline{U}_{\mu}(x_{k})$
$\displaystyle\stackrel{{{}_{def}}}{{=}}A+B-C.$

We have

$\displaystyle E[\|A\|^{2}\mid x_{k}]$
$\displaystyle=E\left[\left.\left\|\frac{\overline{U}_{\mu}(x_{k}+\mu\xi)-\overline{U}_{\mu}(x_{k})-\left\langle\nabla\overline{U}_{\mu}(x_{k}),\mu\xi\right\rangle}{\mu}\xi\circ\left|\xi\right|^{p-2}\right\|^{2}\right|x_{k}\right]$
$\displaystyle\leq
E\left[\left(\tfrac{1}{2}\left(M+\lambda\right)\mu\left\|\xi\right\|^{2}\right)^{2}\left\|\xi\circ\left|\xi\right|^{p-2}\right\|^{2}\mid
x_{k}\right]$ $\displaystyle\leq
E\left[\tfrac{1}{4}\left(M+\lambda\right)^{2}\mu^{2}\left\|\xi\right\|^{4}\left\|\xi\right\|_{p}^{2p-2}d^{\frac{4-2p}{p}}\mid
x_{k}\right]$ $\displaystyle\leq
E\left[\tfrac{1}{4}\left(M+\lambda\right)^{2}\mu^{2}d^{\frac{4-2p}{p}}\left\|\xi\right\|_{p}^{2p+2}\mid
x_{k}\right]$
$\displaystyle\leq\tfrac{1}{4}M\left(M+\lambda\right)^{2}\mu^{2}d^{\frac{4-2p}{p}}E[\left\|\xi\right\|_{p}^{2p+2}]$
$\displaystyle\leq\frac{1}{4}\left(M+\lambda\right)^{2}\mu^{2}d^{\frac{4-2p}{p}}(d+p+1)^{\frac{2p+2}{p}},$
$\displaystyle\leq\frac{1}{4}\left(M+\lambda\right)^{2}\mu^{2}(d+p+1)^{\frac{6}{p}},$
$\displaystyle E[\|B\|^{2}\mid x_{k}]$
$\displaystyle=E[\|\left\langle\nabla\overline{U}_{\mu}(x_{k}),\xi\right\rangle\left(\xi\circ\left|\xi\right|^{p-2}\right)\|^{2}\mid
x_{k}]$ $\displaystyle\leq
E[\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}\left\|\xi\right\|^{2}\left\|\xi\circ\left|\xi\right|^{p-2}\right\|^{2}\mid
x_{k}]$ $\displaystyle\leq
E[\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}\left\|\xi\right\|^{2}\left\|\xi\right\|_{p}^{2p-2}d^{\frac{4-2p}{p}}\mid
x_{k}]$ $\displaystyle\leq
E\left[d^{\frac{4-2p}{p}}\left\|\xi\right\|_{p}^{2p}\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}\mid
x_{k}\right]$ $\displaystyle\leq
d^{\frac{4-2p}{p}}E\left[\left\|\xi\right\|_{p}^{2p}\right]\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}$
$\displaystyle\leq
d^{\frac{4-2p}{p}}d\left(d+1\right)\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}$
$\displaystyle\leq\left[d+1\right]^{\frac{4}{p}}\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}.$

$\displaystyle E[\|B-C\|^{2}\mid x_{k}]$ $\displaystyle\leq
2E[\left\|B\right\|^{2}\mid x_{k}]+2E[\left\|C\right\|^{2}\mid x_{k}]$
$\displaystyle=2\left[d+1\right]^{\frac{4}{p}}\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}+2\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}$
$\displaystyle\leq
2\left[d+2\right]^{\frac{4}{p}}\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|^{2}$

Thus the variance is bounded by

$\displaystyle E\left[\left\|\zeta_{k}-E[\zeta_{k}\mid
x_{k}]\right\|^{2}\right]$
$\displaystyle\leq\left(E\left[\left\|A\right\|^{2}\right]{}^{\frac{1}{2}}+E\left[\left\|B-C\right\|{}^{2}\right]^{\frac{1}{2}}\right)^{2}$
$\displaystyle\leq\left(\frac{1}{2}(M+\lambda)\mu\left(d+3\right){}^{3/p}+\sqrt{2}\left(d+2\right)^{\frac{2}{p}}\left\|\nabla\overline{U}_{\mu}(x_{k})\right\|\right)^{2}.$

∎

###  B.3 Proof of Lemma 3

###### Proof.

This proof adapts the technique of the proof of [6] Proposition 1. Without
loss of generality we may assume that
${\displaystyle\int_{\mathbb{R}^{p}}\exp(-\overline{U}(x))dx=1}$. We first
give upper and lower bounds to the normalizing constant of
$\overline{\pi}_{\mu}$, that is

$c_{\mu}\stackrel{{{}_{\triangle}}}{{=}}\int_{\mathbb{R}^{d}}\overline{\pi}(x)e^{-\left|\overline{U}_{\mu}(\mathrm{x})-\overline{U}(\mathrm{x})\right|}dx.$

The constant $c_{\mu}$ is an expectation with respect to the density
$\overline{\pi}_{\mu}$, it can be trivially upper bounded by 1. From Lemma 1,
$\left|\overline{U}_{\mu}(\mathrm{x})-\overline{U}(\mathrm{x})\right|=\left|E_{\xi}\left(U_{\mu}(\mathrm{x})-U(\mathrm{x})+\frac{\lambda}{2}\left\|x+\mu\xi\right\|^{2}-\frac{\lambda}{2}\left\|x\right\|^{2}\right)\right|=\left|E_{\xi}\left(U_{\mu}(\mathrm{x})-U(\mathrm{x})\right)\right|+\frac{\lambda}{2}\mu^{2}E\left\|\xi\right\|^{2}\leq
a=\frac{L\mu^{1+\alpha}d^{\frac{1+\alpha}{p}}}{1+\alpha}+\frac{1}{2}\lambda\mu^{2}\left(d+1\right)^{\frac{2}{p}}$,
where we have used $E\left\|\xi\right\|^{2}\leq\left(d+1\right)^{\frac{2}{p}}$
from Lemma 4\. So, we have
$e^{-\left|\overline{U}_{\mu}(\mathrm{x})-\overline{U}(\mathrm{x})\right|}\geq
e^{-a}.$ This fact yields $e^{-a}\leq c_{\mu}\leq 1.$ Now we control the
distance between densities $\overline{\pi}$ and $\overline{\pi}_{\mu}$ at any
fixed $x\in\mathbb{R}^{d}$:

$\displaystyle\left|\overline{\pi}(x)-\overline{\pi}_{\mu}(x)\right|$
$\displaystyle=\overline{\pi}(x)\left|1-\frac{e^{-\left|\overline{U}_{\mu}(\mathrm{x})-\overline{U}(\mathrm{x})\right|}}{c_{\mu}}\right|$
$\displaystyle\leq\overline{\pi}(x)\left\\{\left(1-e^{-\left|\overline{U}_{\mu}(\mathrm{x})-\overline{U}(\mathrm{x})\right|}\right)+e^{-\left|\overline{U}_{\mu}(\mathrm{x})-\overline{U}(\mathrm{x})\right|}\left(\frac{1}{c_{\mu}}-1\right)\right\\}$
$\displaystyle\leq\overline{\pi}(x)\left(\left|\overline{U}_{\mu}(\mathrm{x})-\overline{U}(\mathrm{x})\right|+e^{a}-1\right).$

The second inequality is trivial while the last inequality follows from
$1-e^{-x}\leq x$ for any $x\geq 0$. To bound $W_{2}$, we use an inequality
from [8](Theorem 6.15, page 115):

$W_{2}^{2}(\mu,\ \nu)\leq
2\int_{\mathbb{R}^{d}}\|x\|_{2}^{2}|\mu(x)-\nu(x)|dx.$

Combining this with the bound on
$\left|\overline{\pi}(x)-\overline{\pi}_{\mu}(x)\right|$ shown above, we have

$W_{2}^{2}(\overline{\pi},\ \overline{\pi}_{\mu})\leq
2\int_{\mathbb{R}^{d}}\|x\|_{2}^{2}\overline{\pi}(x)\left(a+e^{a}-1\right)dx.$

By [9], the following bound on the second moment, centered on the mode holds

$\int_{\mathbb{R}^{d}}\|x-x^{*}\|_{2}^{2}\overline{\pi}(x)dx\leq\frac{d}{\lambda}.$

In addition, by Young inequality

$\int_{\mathbb{R}^{d}}\|x\|_{2}^{2}\overline{\pi}(x)dx\leq
2\int_{\mathbb{R}^{d}}\|x-x^{*}\|_{2}^{2}\overline{\pi}(x)dx+2\int_{\mathbb{R}^{d}}\|x^{*}\|_{2}^{2}\overline{\pi}(x)dx$
$\leq\frac{2d}{\lambda}+2\|x^{*}\|_{2}^{2}.$

Therefore

$W_{2}^{2}(\overline{\pi},\
\overline{\pi}_{\mu})\leq\frac{4(d+\lambda\|x^{*}\|_{2}^{2})}{\lambda}\left(a+e^{a}-1\right),$

which is the claim of the Lemma. Finally, $a<0.1$ ensures that $e^{a}-1\leq
1.06a$. Follow from previous inequality we have,

$W_{2}(\overline{\pi},\
\overline{\pi}_{\mu})\leq\sqrt{\frac{8.24(d+\lambda\|x^{*}\|_{2}^{2})a}{\lambda}}.$

Treating $L,\,\mu,\,\|x^{*}\|_{2}^{2}$ as constants, choose $\lambda$ small
enough so that $8.24\lambda\|x^{*}\|_{2}^{2}<0.76d$, then
$W_{2}(\overline{\pi},\ \overline{\pi}_{\mu})$ is less than
$3\sqrt{\frac{da}{\lambda}}$. ∎

###  B.4 Additional lemma

######  Lemma 4.

If $\xi\sim N_{p}\left(0,I_{d}\right)$ then
$d^{\left\lfloor\frac{n}{p}\right\rfloor}\leq
E(\left\|\xi\right\|_{p}^{n})\leq\left[d+\frac{n}{2}\right]^{\frac{n}{p}}$where$\left\lfloor
x\right\rfloor$ denotes the largest integer less than or equal to $x.$ If
$n=kp,$ then $E(\left\|\xi\right\|_{p}^{n})=d..(d+k-1)$.

###### Proof.

From [10], we have
$E(\left\|\xi\right\|_{p}^{n})=p^{\frac{n}{p}}\frac{\Gamma\left(\frac{d+n}{p}\right)}{\Gamma\left(\frac{d}{p}\right)}.$

Since $\Gamma$ is an inscreasing function,
$p^{\frac{n}{p}}\frac{\Gamma\left(\frac{d+n}{p}\right)}{\Gamma\left(\frac{d}{p}\right)}\geq
p^{\frac{n}{p}}\frac{\Gamma\left(\frac{d}{p}+\left\lfloor\frac{n}{p}\right\rfloor\right)}{\Gamma\left(\frac{d}{p}\right)}=p^{\frac{n}{p}}\frac{d}{p}\ldots\left(\frac{d}{p}+k-1\right)\geq
d^{\left\lfloor\frac{n}{p}\right\rfloor}.$

If $n=kp$ for $k\in N$ then
$E(\left\|\xi\right\|_{p}^{n})=p^{\frac{n}{p}}\frac{d}{p}\ldots\left(\frac{d}{p}+k-1\right).$
If $n\neq kp$, let $\left\lfloor\frac{n}{p}\right\rfloor=k$. Since $\Gamma$ is
log-convex, by Jensen’s inequality for any $p\geq 1$, we acquire

$\displaystyle\left(1-\frac{n}{p\left\lfloor\frac{n}{p}\right\rfloor+p}\right)\log\Gamma\left(\frac{d}{p}\right)+\frac{n}{p\left\lfloor\frac{n}{p}\right\rfloor+p}\log\Gamma\left(\frac{d}{p}+\left\lfloor\frac{n}{p}\right\rfloor+1\right)$
$\displaystyle\geq\log\Gamma\left(\left(1-\frac{n}{p\left\lfloor\frac{n}{p}\right\rfloor+p}\right)\frac{d}{p}+\frac{n}{p\left\lfloor\frac{n}{p}\right\rfloor+p}\left(\frac{d}{p}+\left\lfloor\frac{n}{p}\right\rfloor+1\right)\right)$
$\displaystyle\geq\log\Gamma\left(\frac{d+n}{p}\right)>0.$

Raising $e$ to the power of both sides, we get

$\Gamma\left(\frac{d}{p}\right)^{\left(1-\frac{n}{p\left\lfloor\frac{n}{p}\right\rfloor+p}\right)}\Gamma\left(\frac{d}{p}+\left\lfloor\frac{n}{p}\right\rfloor+1\right)^{\frac{n}{p\left\lfloor\frac{n}{p}\right\rfloor+p}}\geq\Gamma\left(\frac{d+n}{p}\right),$

which implies that

$\begin{array}[]{cc}\left[\frac{\Gamma\left(\frac{d}{p}+\left\lfloor\frac{n}{p}\right\rfloor+1\right)}{\Gamma\left(\frac{d}{p}\right)}\right]^{\frac{n}{p\left\lfloor\frac{n}{p}\right\rfloor+p}}&\geq\frac{\Gamma\left(\frac{d+n}{p}\right)}{\Gamma\left(\frac{d}{p}\right)}\\\
\left[\frac{d}{p}\ldots\left(\frac{d}{p}+\left\lfloor\frac{n}{p}\right\rfloor\right)\right]^{\frac{n}{p\left\lfloor\frac{n}{p}\right\rfloor+p}}&\geq\frac{\Gamma\left(\frac{d+n}{p}\right)}{\Gamma\left(\frac{d}{p}\right)}.\end{array}$

Combining with
$E(\left\|\xi\right\|_{p}^{n})=p^{\frac{n}{p}}\frac{\Gamma\left(\frac{d+n}{p}\right)}{\Gamma\left(\frac{d}{p}\right)}$
gives the conclusion. ∎

## References

  * [1] C. Robert, G. Casella, Monte Carlo statistical methods, Springer Science & Business Media, 2013.
  * [2] A. S. Dalalyan, A. Karagulyan, User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient, Stochastic Processes and their Applications 129 (12) (2019) 5278–5311.
  * [3] A. Durmus, E. Moulines, M. Pereyra, Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau, SIAM Journal on Imaging Sciences 11 (1) (2018) 473–506.
  * [4] N. S. Chatterji, J. Diakonikolas, M. I. Jordan, P. L. Bartlett, Langevin Monte Carlo without smoothness, arXiv preprint arXiv:1905.13285.
  * [5] L. Shen, K. Balasubramanian, S. Ghadimi, Non-asymptotic Results for Langevin Monte Carlo: Coordinate-wise and Black-box Sampling, arXiv preprint arXiv:1902.01373.
  * [6] A. S. Dalalyan, L. Riou-Durand, A. Karagulyan, Bounding the error of discretized Langevin algorithms for non-strongly log-concave targets, arXiv preprint arXiv:1906.08530.
  * [7] Y. Nesterov, V. Spokoiny, Random gradient-free minimization of convex functions, Foundations of Computational Mathematics 17 (2) (2017) 527–566.
  * [8] C. Villani, Optimal transport: old and new, Vol. 338, Springer Science & Business Media, 2008.
  * [9] A. Durmus, E. Moulines, et al., High-dimensional Bayesian inference via the unadjusted Langevin algorithm, Bernoulli 25 (4A) (2019) 2854–2882.
  * [10] W.-D. Richter, Generalized spherical and simplicial coordinates, Journal of Mathematical Analysis and Applications 336 (2) (2007) 1187–1202.


# Supplemental Materials for
‘Demonstration of the Two-Fluxonium Cross-Resonance Gate’

Ebru Dogan  Department of Physics, University of Massachusetts-Amherst,
Amherst, MA, USA  Dario Rosenstock  Department of Physics, University of
Massachusetts-Amherst, Amherst, MA, USA  Loïck Le Guevel  Department of
Physics, University of Massachusetts-Amherst, Amherst, MA, USA  Department of
Electrical and Computer Engineering, University of Massachusetts-Amherst, MA,
USA  Haonan Xiong  Department of Physics, Joint Quantum Institute, and Center
for Nanophysics and Advanced Materials, University of Maryland, College Park,
MD, USA  Raymond A. Mencia  Department of Physics, Joint Quantum Institute,
and Center for Nanophysics and Advanced Materials, University of Maryland,
College Park, MD, USA  Aaron Somoroff  Department of Physics, Joint Quantum
Institute, and Center for Nanophysics and Advanced Materials, University of
Maryland, College Park, MD, USA  Konstantin N. Nesterov  Department of Physics
and Wisconsin Quantum Institute, University of Wisconsin–Madison, Madison, WI,
USA  Maxim G. Vavilov  Department of Physics and Wisconsin Quantum Institute,
University of Wisconsin–Madison, Madison, WI, USA  Vladimir E. Manucharyan
Department of Physics, Joint Quantum Institute, and Center for Nanophysics and
Advanced Materials, University of Maryland, College Park, MD, USA  Chen Wang
<EMAIL_ADDRESS>Department of Physics, University of Massachusetts-Amherst,
Amherst, MA, USA

##  1 Pulse predistortion via active reflection cancellation

Within our current fridge setup, an impedance mismatch in the signal
transmission chain severely impedes the integrity of the microwave pulses at
the two qubit frequencies: Clean microwaves pulses played into our qubit input
lines generate several reflections hitting the qubit at a later time (on the
order of 10’s of ns) with different amplitude and phase, leading to severe
deviation from the intended rotation. More critically, the pulse reflections
overlap with subsequent pulses resulting in unpredictable behavior of gate
sequences such as contrast loss in power/time Rabi and strongly distorted
chevron patterns (see Supplementary Fig. 1(a) and 1(c)). We attribute the
higher-than-usual reflections to a non-cryogenic directional coupler placed at
the mixing chamber with its coupled port connected to the cavity through a few
eccosorb filters (see Supplementary Fig. 12(b). The un-matched nature of the
cavity inputs combined with the coupled port load deviating from the ideal
50-$\Omega$ value at cryogenic temperature and the absence of attenuation on
the reflection path results in this critical situation.

Supplementary Figure 1: (a, b) Rabi oscillation of Qubit $A$ as a function of
drive time and amplitude under a single-port drive without (a) and with (b)
the predistortion: Problematic contrast loss is resolved and the correct
oscillation frequency is recovered for (b). (c, d) Rabi oscillation of Qubit
$A$ as a function of drive time and drive detuning (the “Chevron pattern”) for
a single-port drive without (c) and with the predistortion (d), driven at an
amplitude corresponding to the $\pi$ pulse used throughout the experiment. The
time axis does not include the ramp time of the pulse, which is 6 ns for each
edge. The reflection problem applies to both qubits and both drive lines:
Similar traces are observed for all four combinations.

We mitigated this issue by applying a simple predistortion model to the qubit
pulses: Our model consists of duplicating the ideal pulse into a train of
delayed pulses whose timing, amplitudes, and phases are chosen to overlap and
cancel the mismatch-induced reflections at the qubit port (Supplementary Fig.
2(a)). All corrections are applied in the background and do not affect the
effective qubit gate times.

Supplementary Figure 2: (a, b) Illustration of the generated pulse envelope
and resulting envelope at the qubit port without any correction (a) and with
predistortion (b). The phases for reflections in (b) are not the same as in
(a) as each correction implemented creates its own reflections that changes
the phase of the ones following it. (c, d) The pulse sequences for
characterizing the quadrature-phase (c) and in-phase (d) components of
reflections: We sweep over the delay time $\Delta t$ for a given number of
repetitions n. $\tau$ is a time long enough (80 ns) for reflections from each
set to not overlap with the next set or the readout. We start by n=1 and
proceed by increasing up to n=7 for further amplification of the features. For
the sequence in (d), the sign of the $\pi-$pulse changes for each iteration.
(e) An example trace from the sequence from (c) (n=1) for characterizing the
quadrature-phase reflection. Each jump that makes the qubit deviate from its
previous state can be attributed to an additional reflection. The delay times
with which these reflections follow the original pulse can be pinned down by
summing the $X_{\pi}$ pulse length with times $t_{i}$ pointed out with red
arrows in the figure. For both in-phase and quadrature-phase reflection
characterizations, we use $X_{\pi}$ pulses as fast as 8 ns. (f) The set of
envelopes that would need to be played from AWG I and Q channels in order to
have a clean $X_{\pi}$ rotation for Qubit $A$.

###  1.1 Reflection characterization

Pulse distortions due to impedance mismatch issues are not unknown to the
superconducting qubit community but mainly arise in fast flux lines driving a
low-impedance inductor and the suggested pulse predistortion schemes mostly
target the flux drive lines [1, 2]. Our case of the microwave pulse reflection
in our qubit drive lines is more relatable to an earlier work [3] from which
we make use of a pulse sequence characterizing the Q-quadrature reflection and
extend it to our own model.

We start by defining a reflection with respect to the main pulse by three
parameters: its delay time, relative amplitude and relative phase. More
conveniently in the qubit context, the reflections in the IQ plane are defined
with the relative in-phase and relative quadrature-phase components.
Equivalently, we can characterize these reflection parameters independently
with relative amplitudes for in-phase and quadrature-phase components at their
relevant delay times.

The first step of the characterization is pinning down the delays of these
reflections. An intuitive way of extracting them is by varying the delay
between two following rotations: $X_{\pi}$ and $X_{-\pi}$: This is a sequence
that isolates and accumulates the quadrature-phase component for the
reflections while cancelling the in-phase component [3]. For short delay times
where the reflections for $X_{\pi}$ arrive after the pulse $X_{-\pi}$ with no
overlaps, all reflections due to $X_{\pi}$ should be cancelled out by the
reflections of $X_{-\pi}$: The measurement should result in ground state. As
the delay increases, reflections of $X_{\pi}$ start overlapping with the pulse
$X_{-\pi}$. For sufficiently short pulses, the result is the step-like
evolution of the qubit state where we can clearly see the effects of
individual reflections (see Supp. Fig. 2(e)): It is possible to identify each
reflection from such a trace as each of the “jumps” can be attributed to an
additional reflection hitting the qubit. The points with the highest
derivatives on these traces can be pinned down to determine the delay times
for these reflections.

We then use the same experiment with our targeted gate length and adjust the
quadrature reflection component associated with each delay until the qubit
state remains in ground for all delays, a sign of a good predistortion
canceling the reflections at the qubit ports. $X_{\pi},X_{-\pi}$ pairs are
repeated up to 7 times, amplifying the deviation from ground, in a single
measurement to finely tune the quadrature components (Supplemenary Fig. 2(c)).

The in-phase component is tuned in a similar measurement consisting of varying
delay between $X_{\pi}$ and $Y_{\pi}$ rotations (Supplementary Fig. 2(d)). The
pulse sequence for the in-phase component characterization is sensitive to
both in-phase and quadrature components: It is beneficial to have a reasonably
well tuned cancellation of the quadrature component before moving on with the
in-phase component. We then need to proceed with an iterative modification of
both I and Q of each reflection, until both pulse sequences result in the
designated final states with no significant deviations.

Both qubits’ drives require this protocol to be implemented, with the case of
Qubit $B$ being more challenging than Qubit $A$ for being a combination of two
individual drives. A sample envelope played from the AWG for an $X_{\pi}$ for
the Qubit $A$ can be seen from Supplementary Fig. 2(f): Rather than a single
gaussian on the I channel for a $\pi$ rotation around the $X$-axis, a set of
two envelopes modified for reflection cancellations are played on both I and Q
channels.

We depend on this predistortion scheme for having meaningful and reliable
operations in our system. An AllXY experiment [4] as a visual indicator of the
level of improvement can be seen from Supplementary Fig. 3.

Supplementary Figure 3: AllXY experiment[4] for the Qubit $A$ without (orange)
and with (blue) the reflection cancellation scheme after the same set of tune
up procedures. AllXY experiment consists of three sets of sequences: The first
set of points land in the ground state, the second in the equator and the
third in the excited state; pulse errors will show up as deviations from this
configuration.

##  2 Calibration of initial state population with an entangling gate

Without a single-shot readout, we take extra steps to calibrate the initial
ground state population of the qubits, i.e. to measure our state
initialization fidelity. This may be accomplished by making use of higher
transitions: Comparing the amplitudes of $\ket{1}$-$\ket{2}$ oscillations with
and without population inversion [5] would provide the information about the
initial state populations. However, this process involves extra complexity for
a two-fluxonium system whose $\ket{1}$-$\ket{2}$ transitions have frequency
and readout visibility strongly dependent on the state of the other qubit. In
order to conveniently measure and track our qubit initialization fidelity over
the long course of our experiment without requiring extra control resources,
we make use of the entangling gate available to us, the $\text{CX}_{\pi}$
gate, to calibrate qubit populations. Intuitively, if the conditional Rabi
drive generates the same contrast as the unconditional Rabi drive, it proves
that the control qubit is in a pure $\ket{1}$ state. A quantitative procedure
can be developed around this principle to calibrate any deviation of the
control qubit from the pure state.

We stress that the lack of knowledge about the initial state purity does not
represent an impediment to benchmarking the two-qubit gate: The fact that the
two-qubit gate fidelity is obtained via the decay rate of the RB measurements
allows extraction of the the $\text{CX}_{\pi}$ gate fidelity as long as we can
measure the relative qubit populations renormalized to the initial states. Our
joint readout scheme, to be elaborated in Supplementary Section 4, provides
such measurements to certify the $\text{CX}_{\pi}$ gate fidelity regardless of
the initialization fidelity. Nevertheless, for definiteness we report actual
qubit populations as calibrated in this Section throughout the manuscript.

The initial state populations Qubit $A$ (control) and Qubit $B$ (target) are
calculated with two complementary sets of measurements, each yielding the
population information for one of the two qubits. The way to design such a
scheme would not be unique: As long as we have a number of independent
measurements sufficient to set up a system of equations with robust solutions,
any of them would be applicable. We describe our procedure under three
approximations and in the end we will discuss the caveats of these
approximations. These assumptions can be stated as the following: 1) The whole
population lies within the computational space
($p_{00}+p_{01}+p_{10}+p_{11}=1$); 2) Initial residual populations for the two
qubits are independent from each other; 3) The CX gate is perfect.

Supplementary Figure 4: Pulse sequences for the set of measurements to
determine the residual population of qubit A. All four measurements result in
oscillations as we sweep the phase of the $\pi/2$ pulse for Qubit $B$.

###  2.1 Measurement protocol for individual qubit populations

To measure Qubit $A$’s population, we apply four Ramsey-like measurements on
Qubit $B$ (Supplementary Fig. 4) to obtain a set of 4 equations. In the
presence of an entangling gate, the conditional flips result in different
oscillation amplitudes for Qubit $B$, depending on the residual population of
Qubit $A$. In the ideal case of no residual population for Qubit $A$ and $B$,
these 4 measurements yield oscillations with the following contrasts in
readout voltages:

$\displaystyle(a)$ $\displaystyle\xrightarrow[]{}M_{01}-M_{00}$
$\displaystyle(b)$ $\displaystyle\xrightarrow[]{}M_{01}-M_{00}$
$\displaystyle(c)$ $\displaystyle\xrightarrow[]{}M_{11}-M_{10}$
$\displaystyle(d)$ $\displaystyle\xrightarrow[]{}M_{11}-M_{10}$

where $M_{ij}$ is the demodulated complex readout voltage for the pure
$\ket{ij}$ state. For this ideal case Qubit $A$ has all the population in the
ground state so the contrasts $(a)-(b)$ and $(c)-(d)$ are same: There will be
no difference between the unconditional and conditional oscillations since
$\textsc{CX}_{\pi}$ has no effect on Qubit $B$ (target) with Qubit $A$
(control) being in ground.

Now we can take into account that Qubit $A$ has some non-zero population $e$
in its excited state:

$\displaystyle(a)$
$\displaystyle\xrightarrow[]{}(1-e)(M_{01}-M_{00})+e(M_{11}-M_{10})$
$\displaystyle(b)$
$\displaystyle\xrightarrow[]{}(1-e)(M_{01}-M_{00})-e(M_{11}-M_{10})$
$\displaystyle(c)$
$\displaystyle\xrightarrow[]{}(1-e)(M_{11}-M_{10})+e(M_{01}-M_{00})$
$\displaystyle(d)$
$\displaystyle\xrightarrow[]{}(1-e)(M_{11}-M_{10})-e(M_{01}-M_{00})$

The voltages $M_{ij}$ can be eliminated by a correct manipulation of these 4
equations:

$e=\frac{(c)-(d)}{(a)+(b)+(c)-(d)}\;\;\;\text{or}\;\;\;e=\frac{(a)-(b)}{(a)-(b)+(c)+(d)}$
(S1)

to give the residual population for the control qubit.

Qubit B’s initial excited state population $\epsilon$ should also be taken
into account to complete this approach. However, it merely brings an overall
$(1-2\epsilon)$ factor to all four equations that gets cancelled out during
the elimination (S1).

The population $\epsilon$ of Qubit B can be measured in a completely analogous
manner with these four sets of sequences if we have a $\text{CX}_{\pi}$ gate
in the reversed direction, with B as control and A as target. While this CXB→A
gate is not native to our system, it can be constructed by the native CXA→B
gate and single-qubit Hadamard gates. Given that $\textsc{CX}_{A\rightarrow
B}$ and $\textsc{CNOT}_{A\rightarrow B}$ only vary by an S-gate on Qubit $A$
(which is accounted for with a virtual-Z in the software), we apply the
decomposition (CNOT${}_{B\rightarrow
A})=\textsc{H}_{A}\textsc{H}_{B}(\textsc{CNOT}_{A\rightarrow
B})\textsc{H}_{A}\textsc{H}_{B}$ for reverting the direction of our gate.

###  2.2 Initial state populations

Using the procedure above, we have measured the residual excited-state
populations of both qubits after applying the reset protocol detailed in
Supplementary Section 3. Over the course of multiple weeks we obtain the
average residual population $e$ of Qubit $A$ to be $1.0\pm 0.3\%$, and the
average residual population $\epsilon$ of Qubit $B$ to be $5.0\pm 1.3\%$. If
we assume the initial qubit states are uncorrelated, we can report the average
initial state population after the reset protocol to be
$\vec{P}_{i}=[P_{00},P_{01},P_{10},P_{11}]\approx[0.94,0.05,0.01,0.00]$. It is
worth considering our assumptions and possible caveats to this initial state
population measurements:

First, since the $\text{CX}_{\pi}$ gate is not perfect, the measurement
outcome of $(b)$ and $(d)$ (from Supplementary Fig. 4) may be subject to a
relative error bounded by the infidelity of the CX gate, $1-\mathcal{F}$. From
Eq. (S1), recall that $(a)\approx(b)$ and $(c)\approx(d)$ (and $(a)$ and $(c)$
are similar in magnitude), and so the possible absolute error of the
calculated population $e$ and $\epsilon$ is bounded by
$\frac{1}{2}(1-\mathcal{F})$, or $0.005$ for a 99% gate fidelity.

Second, it is possible that the initial qubit states are correlated after our
side-band cooling mechanism. However, the most correlated scenario of our
initial state given the data would be
$\vec{P}_{i}\approx[0.95,0.04,0.00,0.01]$, which corresponds to a change
comparable to the measurement noise level. Essentially, when one of the qubit
state is proven to be sufficiently pure, it strongly limits the impact of
possible correlations.

Third, the initial states may have residual population in non-computational
states, e.g. in $\ket{2}$. This population is expected to be very small
considering that we managed to deplete the $\ket{1}$ state efficiently in our
initialization and that the $\ket{1}$-$\ket{2}$ transition has high frequency
and short relaxation time. In addition, if Qubit $A$ is in $\ket{2}$, the Rabi
oscillation of $B$ under CX drive would register distinct signature (i.e.
oscillations of a different period) that is entirely absent in the experiment
(to sub-1% level). We cannot place as stringent of a bound to the possibility
of Qubit $B$ starting in its $\ket{2}$ state, but this small population would
not affect the experiment in any way expect for the overall scale of readout
voltages. For all purposes, we may simply consider $\vec{P}$ as the population
vector normalized within the computation subspace.

##  3 State preparation with multi-photon cooling

Our initial state preparation protocol aims for the reset of Qubit $A$:
Assuming a subsystem composed of $A$ and the cavity $R$, we drive a two-photon
sideband transition from the state $\ket{1}_{A}\ket{0}_{R}$ to
$\ket{0}_{A}\ket{1}_{R}$ which then quickly decays into
$\ket{0}_{A}\ket{0}_{R}$ (Supplementary Fig. 5(a)). Initially without any
reset protocol, Qubit $A$ is close to a 60%-40% mixture state.

To enact this multi-photon reset protocol, we apply a 15 $\mu$s-long square
pulse at a frequency of 3.524 GHz, which is equal to half the difference
between the cavity and Qubit $A$ frequencies.

To calibrate the optimal cooling tone parameters, we sweep over the frequency
of the cooling tone for various powers, comparing the on/off contrast of Qubit
$A$. Using the qubit population measurement presented in Supplementary Section
2, we can observe the effect of the cooling tone with respect to different
cooling tone lengths, reducing the residual excited state population of $A$ to
$\sim$1% and saturating after about 10 $\mu$s. For reasons yet to be
understood, this reset protocol simultaneously also resets Qubit $B$ to
$\sim$5% level. We briefly attempted further cooling of Qubit $B$ with a
$\ket{1}_{B}\ket{0}_{R}\leftrightarrow\ket{0}_{B}\ket{1}_{R}$ two-photon
drive, though it did not bring additional improvement to the initialization.

Supplementary Figure 5: (a)Multi-photon cooling scheme used for initial state
preparation (b) Effect of cooling tones with various lengths on qubit residual
populations for the tone power and frequency specified in Supplementary
Section 4 text. Each point corresponds to the average of 10 measurements (with
5000 averages each) done with same cooling configuration and the error bars
are the standard error for the distribution of 10 measurements.

##  4 Joint readout with pre-pulses

We implement a joint readout method that makes use of a number of averaged
measurements to extract the population distribution of the two qubit system
[6, 7]. Despite increasing the time and number of measurements, this simple
scheme (with no additional experimental complexity to the readout setup) makes
it possible to implement two-qubit experiments with no single-shot readout.
This used to be a common practice before Josephson parametric amplifiers
became widely available.

At the end of any measurement, the readout gives us a complex demodulated
voltage of the cavity transmission that involves contributions from all states
weighted with their corresponding populations:

$V=p_{00}M_{00}+p_{01}M_{01}+p_{10}M_{10}+p_{11}M_{11}$

where $\vec{p}=\text{diag}[\rho]=[p_{00},p_{01},p_{10},p_{11}]^{T}$ is the
population distribution for any state of the two-qubit system. This population
4-vector is particularly essential for benchmarking purposes where we need to
know how much of the final state is in $\ket{00}$. This single complex voltage
value by itself is not sufficient to pin down the population vector $\vec{p}$:
In order to obtain $\vec{p}$ at any moment, we first need to know about the
4-vector $\vec{M}$ where the components are complex demodulated voltages of
the cavity transmission for the computational states
$[M_{00},M_{01},M_{10},M_{11}]$.

We also need to have a sufficient number of independent measurements that will
let us solve for the unknowns $p_{00},p_{01},p_{10},p_{11}$. For that, any
prepared state will be measured with four different readout configurations,
varied by pre-readout single qubit rotations:

1) no rotation (II)

2) $\pi$-pulse on Qubit $B$ (IX)

3) $\pi$-pulse on Qubit $A$ (XI)

4) $\pi$-pulse on both qubits (XX)

With this scheme, each measurement will eventually yield a set of 4 measured
voltages: $\vec{V}=[V_{II},V_{IX},V_{XI},V_{XX}]^{T}$. The correct mapping
relates the four populations to the four measured voltages:

$\begin{pmatrix}V_{II}\\\ V_{IX}\\\ V_{XI}\\\ V_{XX}\\\
\end{pmatrix}=\begin{pmatrix}M_{00}&M_{01}&M_{10}&M_{11}\\\
M_{01}&M_{00}&M_{11}&M_{10}\\\ M_{10}&M_{11}&M_{00}&M_{01}\\\
M_{11}&M_{10}&M_{01}&M_{00}\\\ \end{pmatrix}\begin{pmatrix}p_{00}\\\ p_{01}\\\
p_{10}\\\ p_{11}\\\ \end{pmatrix}$

The matrix M is nothing but a permutation of the components of $\vec{M}$. The
vector $\vec{M}$ will allow us to reconstruct the inverse of the mapping
relation that will give any state’s population distribution, given the
measured voltages:

$\vec{p}=M^{-1}\vec{V}$ (S2)

In the case of perfect state preparation, obtaining components of the
$\vec{M}$ would be straightforward. Although it is not a realistic assumption,
it is intuitive to look at the case with $\vec{P_{i}}=[1,0,0,0]^{T}$ where
$\vec{P_{i}}$ is the 4-vector for the initial states of the system:

$\begin{pmatrix}V_{II}\\\ V_{IX}\\\ V_{XI}\\\ V_{XX}\\\
\end{pmatrix}=\begin{pmatrix}M_{00}&M_{01}&M_{10}&M_{11}\\\
M_{01}&M_{00}&M_{11}&M_{10}\\\ M_{10}&M_{11}&M_{00}&M_{01}\\\
M_{11}&M_{10}&M_{01}&M_{00}\\\ \end{pmatrix}\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\\\
\end{pmatrix}=\begin{pmatrix}M_{00}\\\ M_{01}\\\ M_{10}\\\ M_{11}\\\
\end{pmatrix}$

For a more realistic case, we know that the initial population vector should
have non-zero populations for the remaining computational states as well. With
the initial state populations $\vec{P_{i}}$ (obtained using the method from
Supplementary Section 2), the mapping can be written in the following way:

$\begin{pmatrix}M_{00}\\\ M_{01}\\\ M_{10}\\\ M_{11}\\\
\end{pmatrix}=\begin{pmatrix}P_{00}&P_{01}&P_{10}&P_{11}\\\
P_{01}&P_{00}&P_{11}&P_{10}\\\ P_{10}&P_{11}&P_{00}&P_{01}\\\
P_{11}&P_{10}&P_{01}&P_{00}\\\ \end{pmatrix}^{-1}\begin{pmatrix}V_{II}\\\
V_{IX}\\\ V_{XI}\\\ V_{XX}\\\ \end{pmatrix}$

which gives us the required complex demodulated voltages. With $\vec{M}$’s
components at hand, we can now construct the full mapping from equation (S2)
that gives us the population distribution for the computational states at any
point of a measurement.

##  5 Calibration and the benchmarking of Single Qubit Gates

Supplementary Figure 6: Effective circuit diagram of the capacitively coupled
two-fluxonium system with asymmetric coupling to two input ports. The cavity
mode is omitted for simplicity.

The fact that we are working in a 3D geometry with no local drives and with a
system that has a strong CR component brings up the challenge of addressing
the qubits individually: Given that both input ports are coupled to both
qubits (see Supplementary Fig. 6), any qubit $B$ drive from one of the
individual ports (at frequency $\omega_{B}$) will also be seen by Qubit $A$,
leading to conditional oscillations of $B$ (Supplementary Fig. 8(a) and (b)).
In order to realize high-fidelity single qubit gates, we need to be able to
address each qubit independently from the state of the other one. As mentioned
in the main text, the frequency of the Qubit $A$ oscillations driven from
input 1 does not have a significant dependence on the Qubit $B$ state and can
be considered as a “clean” single qubit drive by itself. On the other hand, in
order to realize a “clean” single qubit gate for B, the two individual drives
need to be combined with a scheme similar to the CR gate, only this time to
fulfill the requirement:
$\matrixelement{00}{\hat{H}_{drive}}{01}=\matrixelement{10}{\hat{H}_{drive}}{11}$
with the drive Hamiltonian from Eq. (5) of the main text. With a similar
procedure, one can find the right relative amplitude and phase between the two
drives to satisfy this condition. This vector compensation scheme is
illustrated in more detail in Supplementary Fig. 7. We can consider the drive
from each port as having two components: a direct drive, which is independent
of the state of Qubit $A$, and a CR drive, which is turned on whenever Qubit
$A$ is in $\ket{1}$. We tune the complex ratio $\eta^{\prime}$ of the drives
on each port such that the CR components are of equal amplitude and opposite
phase, thereby leaving only the direct drives. With the proper ratio, Qubit
$B$’s Rabi frequency no longer depends on the state of Qubit $A$
(Supplementary Fig. 8(c)).

Supplementary Figure 7: (a) Effective drive on qubit B, shown in the IQ plane,
when qubit A is in $\ket{0}$. The qubit is driven with a different amplitude
and phase depending on which driving port is used. (b) Effective drive on
qubit B, shown in the IQ plane, when qubit A is in $\ket{1}$. The qubit is
driven with a different effective amplitude and phase due to the CR component
of the drive now being turned on. (c) Relative amplitude and phase are applied
to the two driving ports to exactly cancel the CR drives and achieve
unconditional single qubit rotations. Supplementary Figure 8: (a) Rabi
oscillations of Qubit B driven from Port 1 (b) Rabi oscillations of Qubit B
driven from Port 2 (c) Rabi oscillations of Qubit B after applying the complex
ratio $\eta^{\prime}$ between the two driving ports.

We benchmark our single qubit gates using simultaneous single qubit randomized
benchmarking: We generate random sequences of N pairs of Cliffords, one for
each qubit, with N ranging from 1 to 200, and apply them to each qubit
simultaneously in time. At the conclusion of each sequence, we measure two
different quantities: the fidelity of the $\ket{0}$ state of qubit A and the
fidelity of the $\ket{0}$ state of qubit B. This measurement was repeated 456
times over the course of several weeks, with different random sequences of N
Clifford pairs being generated each time. The results were averaged together,
separately for each qubit, and each fit to a single exponential. The resulting
traces are shown in Supplementary Figure 9. From these traces, we can extract
the average error per Clifford on each individual qubit. Taking into account
the construction of our single qubit Clifford gate set, which utilizes virtual
Z gates and contains 1.167 physical single qubit gates per Clifford, we report
the following average fidelities per physical gate:
$\mathcal{F}_{A}=0.99835(4)$ and $\mathcal{F}_{B}=0.99734(6)$.

Supplementary Figure 9: Simultaneous single qubit randomized benchmarking. The
two traces shown are measured fidelities to the $\ket{0}$ state of qubit A and
the $\ket{0}$ state of qubit B at the conclusion of N random pairs of single-
qubit Clifford gates. Each trace is the averaged result of 456 randomized
benchmarking experiments performed over the course of several weeks. The
errorbars shown are the standard errors taken across these 456 experiments.

##  6 Calibration and benchmarking of the CR gate

The preliminary $\textsc{CX}_{\pi}$ gate tune up involves the calibration of
the relative phase, relative amplitude and the overall amplitude, as
elaborated in the main text. But we are in need of a more sophisticated
calibration procedure in order to realize a high-fidelity CX gate: Assuming
that we already have a gate with the above three parameters calibrated
reasonably well, we can now proceed with the next level of calibration. We
should also note down that we do not have any particular gate parameter
optimization procedure and these well designed tune-ups are what we rely on
for realizing high-fidelity gates.

From the observed dynamics of our CR gate, we know that the process has a
block-diagonal structure. The remaining degrees of freedom that would not be
distinguishable for the drives from Fig. 3 of main text are the phases for the
two qubits. Taking them into account, we can write our experimental CR matrix
to have the following form:

$\textsc{CR}_{exp}=\begin{pmatrix}1&0&0&0\\\ 0&e^{i\theta_{B}}&0&0\\\
0&0&0&-ie^{i\theta_{A}}\\\ 0&0&-ie^{i\theta_{A}}&0\end{pmatrix}\ $ (S3)

where $\theta_{A}$ and $\theta_{B}$ are phases picked by Qubit $A$ and $B$
during our cross-resonance gate operation. Brought in a superposition state,
the phases $\theta_{A}$ ($\theta_{B}$) would be what we would observe to be
the additional phase accumulated on the Qubits $A$ ($B$) under the CR drive.
It is possible to measure these phases using the same pulse sequences from
Supplementary Fig. 4: The phase difference between the oscillations resulting
from the measurements of (a) and (b) on qubit A (qubit B) will give us
$\theta_{A}$ ($\theta_{B})$.

###  6.1 Integrating single-qubit Z rotations

Given that our experimental cross-resonance gate realization brings along
extra phase accumulations on both Qubit $A$ and $B$, we first need to show the
equivalence of our experimental gate to a $\textsc{CX}_{\pi}$ gate.

The relevant Z-rotation operations that Qubit $A$ and $B$ undergo are
expressed as:

$\displaystyle\textsc{Z}^{A}(\theta_{A})=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&0&e^{i\theta_{A}}\\\ 0&0&e^{i\theta_{A}}&0\end{pmatrix}\ $
$\displaystyle\textsc{Z}^{B}(\theta_{B})=\begin{pmatrix}1&0&0&0\\\
0&e^{i\theta_{B}}&0&0\\\ 0&0&1&0\\\ 0&0&0&e^{i\theta_{B}}\end{pmatrix}\ $

It can be shown that our experimental CR gate is equivalent to a CX gate up to
a number of Z-rotations:

$\left[\text{CX}_{\pi}\right]=\left[\text{Z}^{B}({-\frac{\theta_{B}}{2}})\right]\left[\text{CR}_{exp}\right]\left[\text{Z}^{B}(-\frac{\theta_{B}}{2})\right]\left[\text{Z}^{A}(\frac{\theta_{B}}{2}-\theta_{A})\right]$
(S4)

Since these Z-rotations can simply be applied via virtual Z’s [8] in the
software, there is a complete equivalence between our experimental cross-
resonance gate and a $\text{CX}_{\pi}$ gate as long the correct virtual Z
operations are applied. We can also trivially change $\theta_{A}$ by $\pi/2$
to convert the gate to a CNOT. For the rest of the section and for all our
benchmarking purposes, our gate is composed of the form as in Eq. (S4).

###  6.2 Fine tuning the gate for high-fidelity

Supplementary Figure 10: CR tune up sequences: These 7 experiments check the
following criteria: 1) No drive on qubit B when qubit A is off: The
calibration experiment (a) makes sure that qubit B does not rotate with
respect to X-axis when qubit A is in state 0. Experiment (b) does a similar
check for the Y-axis. 2) Resonant-$\pi$ condition when qubit A is on:
Experiment (c) checks whether qubit B over-rotates or not, given qubit A is in
state 1. Experiment (d) does the complementary check by making sure that qubit
B is driven along the big circle with qubit A being on. 3) Alignment check:
Experiment (e) checks if the CR rotation aligns with the X-axis of qubit B. 4)
Correction of single qubit frames: Experiments (f) and (g) make sure that the
extra phases on qubit A and qubit B are compensated in the correct way.

The fine tune-up procedure consists of the calibration of the following
parameters:

  1. 1.

Relative amplitude of the two drives

  2. 2.

Relative phase of the two drives

  3. 3.

Overall amplitude of the drive

  4. 4.

Detuning of the CR drive with respect to the qubit B drive (CR detuning:
$\Delta\omega_{\text{CR}}$)

  5. 5.

Misalignment of the CR drive with respect to the qubit B drive (CR angle:
$\Delta\phi_{\text{CR}}$)

  6. 6.

Phase picked up by Qubit $B$ during the CR drive

  7. 7.

Phase picked up by Qubit $A$ during the CR drive

Out of these parameters, the requirement for CR detuning (4) and CR angle (5)
has been determined heuristically: The optimal values for these two parameters
are small but non-zero ($\Delta\omega_{CR}\sim$ order of 100 kHz,
$\Delta\phi_{\text{CR}}\sim$ order of $10^{-2}$ radians) and we suspect that
their presence is related to the remaining underlying reflection issue. Our
procedure is similar to the calibration procedure from the work [9] and has
been tailored to fine tune the seven degrees of freedom of our
$\text{CX}_{\pi}$ gate. The detailed procedure and the pulse sequences can be
seen from Supplementary Fig. 10: Sweeping over the calibration parameter, the
two measurements that differ from each other by an inversion of one of their
pulses result in a crossing point that gives the correct value for the each of
these seven parameters. In order to increase the sensitivity of these
calibrations, we can increase the number of repetitions of the CR gate in the
sequence (an even number of repetitions is required for the $\theta_{A}$
calibration.)

This calibration procedure is run at regular intervals during gate
benchmarking and is completed with a visual syndrome check that includes all
the parameter checks simultaneously.

###  6.3 Quantum process tomography with joint readout

The quantum process matrix of the CXπ gate is estimated over a set of 1044
measurements with joint state readout. We first apply the two-qubit gate on 36
initial states prepared by applying all combinations of simultaneous target
and control single qubit gate: $I$, $X_{\pi/2}$, $X_{-\pi/2}$, $X_{\pi}$,
$Y_{\pi/2}$, and $Y_{-\pi/2}$. State tomography is then applied on the 36 new
states following [10], with 29 pre-read-out pulses and measurements from which
the final density matrix is found via least-square optimization. The CXπ
process matrix is estimated via a simultaneous least-square optimization over
the 36 final density matrices.

###  6.4 CXπ Randomized Benchmarking

To measure the fidelity of our $\text{CX}_{\pi}$ gate, we perform interleaved
two-qubit randomized benchmarking. To construct our two-qubit Clifford gate
set, we follow the procedures in [11] and [12], modifying it to create the
two-qubit Clifford group using our $\text{CX}_{\pi}$ gate as the generator.
The gate decomposition can be seen in Supplementary Fig. 11. Each two-qubit
Clifford contains on average 6.483 physical single qubit gates and 1.5 CX
gates.

We extract fidelity by measuring the sequence fidelity at at least 12
different N’s with N typically going from 1 to 100. To measure sequence
fidelity at each N, we generate at least 40 different random Clifford
sequences, and for each random sequence we average 5000 times to reduce noise.
The average sequence fidelity is fit according to the model $F=Ap^{m}+B$ (with
no constraints) to obtain the decay constants and the error per clifford is
$\text{EPC}=\frac{3}{4}(1-p)$. The error bars on Fig.3(a) are standard errors
$\left(\text{stdev}/\sqrt{N_{t}}\right)$ where “stdev” is the standard
deviation from the mean and $N_{t}$ is the number of runs. The uncertainties
for the EPC’s and fidelity numbers are calculated from the fit errors. The
three different kind of EPC’s ($\text{EPC}_{\text{ref}}$ for reference two
qubit-RB; $\text{EPC}_{\text{CX}}$ for interleaved RB with $\text{CX}_{\pi}$;
$\text{EPC}_{\text{idle}}$ for interleaved RB with I) are calculated
separately and average fidelities for these three are obtained via the
following relations:

$\text{Average 2Q Clifford fidelity}=(1-\text{EPC}_{\text{ref}})$

$\text{Average $\text{CX}_{\pi}$
fidelity}=(1-\text{EPC}_{\text{CX}})/(1-\text{EPC}_{\text{ref}})$

$\text{Average Idle
fidelity}=(1-\text{EPC}_{\text{idle}})/(1-\text{EPC}_{\text{ref}})$

Supplementary Figure 11: Two-qubit Clifford group composition with
$\text{CX}_{\pi}$: $\text{C}_{1}$ refers to the single-qubit Clifford group,
while the $\text{S}_{1}$ and $\text{S}_{1}^{X/2,-Y/2}$ are each three-element
subsets of $\text{C}_{1}$.

##  7 Experimental Setup

Supplementary Figure 12: (a) Room temperature setup (b) Cryogenic setup

The most fundamental microwave setup requirement for our experimental scheme
is the phase locking between the two drive tones played into the two input
ports. We realize this by using a single microwave generator that resources
all individual qubit drives: The same local oscillator feeds the two (home-
made) IQ mixers (one per control line) to maintain phase coherence during
combined drives. Its frequency is set at 840 MHz based on hardware limitations
and to minimize wrong sideband and qubit frequency overlaps. Each IQ mixer is
made with discrete 3-port mixers and power combiners/splitters carefully
chosen to cover both qubit frequencies after up-mixing (corresponding to
intermediate frequencies of 164 MHz for Qubit B and -285MHz for Qubit A). The
room temperature microwave setup is shown in Supplementary Fig. 12(a). In
addition to what is shown in that figure, we also heavily filter our lines
using low and band pass filters which are omitted in the diagram for
simplicity. The fridge setup is shown in Supplementary Fig. 12(b): The non-
cryogenic directional coupler to which we attribute the impedance mismatch can
be seen below the mixing chamber. All the Eccosorb filters are made in-house
and have been chained up (4x for each input line and 3x for the output line)
between the cavity ports and the closest microwave component. The experiment
is carried out in an Oxford Triton 500 Dilution Refrigerator with a based
temperature below 20 mK.

##  8 Fabrication of the two-fluxonium device

The fabrication procedure for the two-fluxonium device is identical to the one
in [13]: The device is fabricated on a 430 $\mu$m thickness sapphire chip with
dimensions $9\times 4$ mm. Following the cleaning of the chip, the double-
layer resist’s first layer is prepared by spinning MMA EL-13 at 5000 RPM and
baking at $180^{\circ}$C for 1 minute, followed by the second layer of the
resist: 950 PMMA A3 spun at 4000 RPM for 1 minute and baked at the same
temperature for half an hour. As the anti-charging layer, 11 nm of Al is
deposited using a Plassys deposition system. The electron beam lithography is
realized using a current of 1 nA on a 100 kV Elionix EBL system. After the
electron beam lithography, the anti-charging layer is removed by keeping the
chip in a 0.1 M potassium hydroxide (KOH) aqueous solution. The development is
done in a 3:1 IPA:DI solution at a temperature of $6^{\circ}$C for 2 minutes.
A Plassys deposition system is used for the double angle Al deposition (with a
Al deposition thickness of 20 nm for the first and 40 nm for the second layer)
with the same process steps detailed in [13]. Liftoff is realized by keeping
the chip in acetone at a temperature of $60^{\circ}$C for 3 hours, followed by
brief sonications in acetone and IPA (5 and 10 seconds). As the last step, the
chip is dried with $N_{2}$.

## References

  * Jerger  _et al._ [2019] M. Jerger, A. Kulikov, Z. Vasselin, and A. Fedorov, In Situ Characterization of Qubit Control Lines: A Qubit as a Vector Network Analyzer, Physical Review Letters 123, 150501 (2019).
  * Rol  _et al._ [2020] M. A. Rol, L. Ciorciaro, F. K. Malinowski, B. M. Tarasinski, R. E. Sagastizabal, C. C. Bultink, Y. Salathe, N. Haandbaek, J. Sedivy, and L. DiCarlo, Time-domain characterization and correction of on-chip distortion of control pulses in a quantum processor, Applied Physics Letters 116, 054001 (2020).
  * Gustavsson  _et al._ [2013] S. Gustavsson, O. Zwier, J. Bylander, F. Yan, F. Yoshihara, Y. Nakamura, T. P. Orlando, and W. D. Oliver, Improving Quantum Gate Fidelities by Using a Qubit to Measure Microwave Pulse Distortions, Physical Review Letters 110, 040502 (2013).
  * Reed [2013] M. Reed, Entanglement and Quantum Error Correction with Superconducting Qubits, arXiv:1311.6759 [cond-mat, physics:quant-ph] (2013), arXiv: 1311.6759.
  * Geerlings  _et al._ [2013] K. Geerlings, Z. Leghtas, I. M. Pop, S. Shankar, L. Frunzio, R. J. Schoelkopf, M. Mirrahimi, and M. H. Devoret, Demonstrating a Driven Reset Protocol for a Superconducting Qubit, Physical Review Letters 110, 120501 (2013).
  * Filipp  _et al._ [2009] S. Filipp, P. Maurer, P. J. Leek, M. Baur, R. Bianchetti, J. M. Fink, M. Göppl, L. Steffen, J. M. Gambetta, A. Blais, and A. Wallraff, Two-Qubit State Tomography Using a Joint Dispersive Readout, Physical Review Letters 102, 200402 (2009).
  * Chow  _et al._ [2010] J. M. Chow, L. DiCarlo, J. M. Gambetta, A. Nunnenkamp, L. S. Bishop, L. Frunzio, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Detecting highly entangled states with a joint qubit readout, Physical Review A 81, 062325 (2010).
  * McKay  _et al._ [2017] D. C. McKay, C. J. Wood, S. Sheldon, J. M. Chow, and J. M. Gambetta, Efficient $Z$ gates for quantum computing, Physical Review A 96, 022330 (2017).
  * Wei  _et al._ [2021] K. X. Wei, E. Magesan, I. Lauer, S. Srinivasan, D. F. Bogorin, S. Carnevale, G. A. Keefe, Y. Kim, D. Klaus, W. Landers, N. Sundaresan, C. Wang, E. J. Zhang, M. Steffen, O. E. Dial, D. C. McKay, and A. Kandala, Quantum crosstalk cancellation for fast entangling gates and improved multi-qubit performance, arXiv:2106.00675 [quant-ph] (2021).
  * Chow  _et al._ [2012] J. M. Chow, J. M. Gambetta, A. D. Córcoles, S. T. Merkel, J. A. Smolin, C. Rigetti, S. Poletto, G. A. Keefe, M. B. Rothwell, J. R. Rozen, M. B. Ketchen, and M. Steffen, Universal Quantum Gate Set Approaching Fault-Tolerant Thresholds with Superconducting Qubits, Physical Review Letters 109, 060501 (2012).
  * Córcoles  _et al._ [2013] A. D. Córcoles, J. M. Gambetta, J. M. Chow, J. A. Smolin, M. Ware, J. Strand, B. L. T. Plourde, and M. Steffen, Process verification of two-qubit quantum gates by randomized benchmarking, Physical Review A 87, 030301 (2013).
  * Barends  _et al._ [2014] R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C. White, J. Mutus, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, C. Neill, P. O’Malley, P. Roushan, A. Vainsencher, J. Wenner, A. N. Korotkov, A. N. Cleland, and J. M. Martinis, Superconducting quantum circuits at the surface code threshold for fault tolerance, Nature 508, 500 (2014).
  * Somoroff  _et al._ [2021] A. Somoroff, Q. Ficheux, R. A. Mencia, H. Xiong, R. V. Kuzmin, and V. E. Manucharyan, Millisecond coherence in a superconducting qubit, arXiv:2103.08578 [quant-ph] (2021).