# New method for renormalon subtraction using Fourier transform

Y. Hayashia, Y. Suminoa and H. Takaurab aDepartment of Physics, Tohoku
University, Sendai, 980–8578 Japan
bTheory Center, KEK, Tsukuba, Ibaraki 305-0801, Japan

###### Abstract

To improve accuracy in calculating QCD effects, we propose a method for
renormalon subtraction in the context of the operator-product expansion. The
method enables subtracting renormalons of various powers in $\Lambda_{\rm
QCD}$ efficiently and simultaneously from single-scale observables. We apply
it to different observables and examine consistency with theoretical
expectations.

††journal: Physics Letters B

TU–1114, KEK–TH–2288

In view of an outstanding success of the Standard Model for particle physics
today, it is becoming highly important to calculate observables of particle
physics with high precision. This would be needed, among others, to probe the
next energy scale that dictates the laws of elementary particle physics. In
particular, it is still challenging to calculate various QCD effects
accurately, even though many technologies for such calculation have been
developed over the past decades.

Use of the operator-product expansion (OPE) in combination with renormalon
subtraction is a way to achieve high precision calculation of QCD effects. In
the OPE of an observable, ultraviolet (UV) and infrared (IR) contributions are
factorized into the Wilson coefficients and nonperturbative matrix elements,
respectively. The former are calculated in perturbative QCD, while the latter
are determined by nonperturbative methods. It has been recognized, however,
that (IR) renormalons are contained in the respective parts. IR renormalons
are the source of rapid growth of perturbative coefficients, which lead to bad
convergence of perturbative series. They limit achievable accuracies and
induce inevitable uncertainties even when using the OPE framework. Since a
physical observable as a whole should not contain renormalon uncertainties,
the necessary task is to separate the renormalons from the respective parts
and cancel them.

Historically cancellation of renormalons made a strong impact in heavy
quarkonium physics. The perturbative series for the quarkonium energy levels
turned out to be poorly convergent when the quark pole mass was used to
express the levels, following the conventional wisdom. When they were
reexpressed by a short-distance quark mass, the convergence of the
perturbative series improved dramatically. This was understood as due to
cancellation of the ${\cal O}(\Lambda_{\rm QCD})$ renormalons between the pole
mass and binding energy [1, 2, 3]. A similar cancellation was observed in the
$B$ meson partial widths in the semileptonic decay modes [4, 5, 6]. These
features were applied in accurate determinations of fundamental physical
constants such as the heavy quark masses [7, 8], some of the Cabbibo-
Kobayashi-Maskawa matrix elements [9, 10], and the strong coupling constant
$\alpha_{s}$ [11].

Analyses including cancellation of renormalons beyond the ${\cal
O}(\Lambda_{\rm QCD})$ renormalon of the pole mass have started only recently.
The ${\cal O}(\Lambda_{\rm QCD})$ and ${\cal O}(\Lambda_{\rm QCD}^{3})$
renormalons of the static QCD potential $V_{\rm QCD}(r)$ were subtracted and
combined with the corresponding nonperturbative matrix elements [12], which
were determined by comparing the renormalon-subtracted leading Wilson
coefficient with a lattice result. By the renormalon subtraction, the
perturbative uncertainty of the Wilson coefficient reduced considerably and
the matching range with the lattice data became significantly wider. By the
matching $\alpha_{s}$ was determined with good accuracy, which agreed with
other measurements. Also the ${\cal O}(\Lambda_{\rm QCD}^{4})$ renormalon
contained in the lattice plaquette action was subtracted and absorbed into the
local gluon condensate [13].

In this paper we generalize the method which was previously established only
for $V_{\rm QCD}(r)$ [14]. In a general case it was difficult to find a simple
integral representation of renormalons which enables separation of them from
the rest. A transparent way was found to understand the corresponding
mechanism for $V_{\rm QCD}(r)$ using properties of the Fourier transform [15].
We are now able to apply it to a general single-scale observable, and we
present the necessary formulation. The formulation is simple and easy to
compute, once the standard perturbative series is given for the observable.

We apply our method to subtract the following renormalon(s) in our test
analyses: the ${\cal O}(\Lambda_{\rm QCD}^{4})$ renormalon of the Adler
function; the ${\cal O}(\Lambda_{\rm QCD}^{2})$ renormalon of the $B\to
X_{u}\ell\bar{\nu}$ decay width; the ${\cal O}(\Lambda_{\rm QCD})$ and ${\cal
O}(\Lambda_{\rm QCD}^{2})$ renormalons simultaneously of the $B$ or $D$ meson
mass. The renormalons of ${\cal O}(\Lambda_{\rm QCD}^{2})$ in the $B$ decay
and $B$, $D$ meson masses are subtracted for the first time in this paper. We
show that our results meet theoretical expectations, e.g., good convergence
and consistent behavior with the OPE.

Consider a general dimensionless observable $X(Q)$ with a characteristic hard
scale $Q$, whose perturbative expansion is given by
$X=\sum_{n}c_{n}\alpha_{s}^{n+1}$, where $\alpha_{s}=\alpha_{s}(\mu)$ and
$\mu$ is the renormalization scale.

An ambiguity induced by renormalons is defined from the discontinuities of the
corresponding singularities in the Borel plane:

$\delta
X=\frac{1}{b_{0}}\int_{C_{+}-C_{-}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!du\,e^{-u/(b_{0}\alpha_{s})}B_{X}(u)\,,~{}~{}~{}B_{X}=\sum_{n}\frac{c_{n}}{n!}\Bigl{(}\frac{u}{b_{0}}\Bigr{)}^{n}\,.$
(1)

In the complex Borel ($u$) plane the integral contours $C_{\pm}(u)$ connect
the origin and $\infty\pm i\epsilon$ infinitesimally above/below the positive
real axis on which the discontinuities are located.
$b_{0}=(11-2n_{f}/3)/(4\pi)$ denotes the one-loop coefficient of the QCD beta
function, where $n_{f}$ is the number of quark flavors.

The location of a renormalon singularity $u_{*}$ and the form of $\delta X$
due to this singularity (apart from its overall normalization $N_{u_{*}}$) can
be determined from the OPE and renormalization-group (RG) equation [16].
Typically $\delta X(Q)\approx N_{u_{*}}(\Lambda_{\overline{\rm
MS}}/Q)^{2u_{*}}$, up to systematically calculable corrections [expressed by
the anomalous dimension and series expansion in $\alpha_{s}(Q)$], where
$\Lambda_{\overline{\rm MS}}$ denotes the dynamical scale defined in the
$\overline{{\rm MS}}$ scheme.

To regulate the divergent behavior of $X(Q)$, we adopt the usual “principal
value (PV) prescription” and take the PV of the Borel resummation integral,
that is, take the average over the contours $C_{\pm}(u)$,

$[X(Q)]_{\rm PV}=\frac{1}{b_{0}}\int_{0,\rm
PV}^{\infty}du\,e^{-u/(b_{0}\alpha_{s})}B_{X}(u)\,.$ (2)

In this regularization, the renormalons are minimally subtracted from $X(Q)$
by definition. How to evaluate the PV integral in eq. (2) is non-trivial in
practice with the limited number of known perturbative coefficients. In this
paper we propose a new method to obtain the PV integral in a systematic
approximation using a different integral representation.

To evaluate $[X(Q)]_{\rm PV}$, we extend the formulation of ref. [14], which
works for $V_{\rm QCD}(r)$. In the static QCD potential, renormalons are
located at $u_{*}=1/2$, $3/2$, $\dots$. It is known [2, 3, 17, 15] that these
renormalons are eliminated (or highly suppressed) and the perturbative series
exhibits good convergence in the momentum-space potential, i.e., the Fourier
transform of the coordinate-space potential $V_{\rm QCD}(r)$. Equivalently,
$V_{\rm QCD}(r)$ is given by the inverse Fourier integral of the momentum-
space potential, which is (largely) free from the renormalons. This indicates
that the renormalons of $V_{\rm QCD}(r)$ arise from (IR region of) the inverse
Fourier transform. Using the formulation of ref. [14], one can avoid the
renormalon uncertainties reviving from the inverse Fourier transform and give
a renormalon-subtracted prediction for $V_{\rm QCD}(r)$. This is realized by
proper deformation of the integration contour of the inverse Fourier
transform. In this method, one minimally subtracts the renormalons using
Fourier transform, and the renormalon-subtracted prediction is actually
equivalent to the PV integral (2) as long as the momentum-space potential is
free of renormalons. We propose a generalized method using an analogous
mechanism. The key to achieving this goal is to find a proper Fourier
transform such that the renormalons of an original quantity are suppressed.

For a general observable $X(Q)$, let $r=|\vec{x}|=Q^{-1/a}$ and define a
Fourier transform of $r^{2au^{\prime}}X$ (where $a$ and $u^{\prime}$ are
parameters) into “momentum ($\tau$) space” as

$\tilde{X}(\tau)=\int
d^{3}\vec{x}\,{e^{i\vec{\tau}\cdot\vec{x}}}\,{r^{2au^{\prime}}}X(r^{-a})\,,~{}~{}~{}(\tau=|\vec{\tau}|)\,,$
(3)

whose typical energy scale is now $\tau^{a}$. Since the Borel resummation and
Fourier transform mutually commute, it follows that

$\delta\tilde{X}(\tau)=\int
d^{3}\vec{x}\,{e^{i\vec{\tau}\cdot\vec{x}}}\,{r^{2au^{\prime}}}\delta
X(r^{-a})\,,$ (4)

where $\delta\tilde{X}$ is defined from $\tilde{X}$ similarly to eq. (1).
Substituting the leading form $N_{u_{*}}(\Lambda_{\overline{\rm
MS}}/Q)^{2u_{*}}$ of $\delta X$ to eq. (4) we obtain

$\delta\tilde{X}\approx-\frac{4\pi N_{u_{*}}\Lambda_{\overline{\rm
MS}}^{2u_{*}}}{\tau^{3+2a(u_{*}+u^{\prime})}}\sin(\pi
a(u_{*}+u^{\prime}))\,\Gamma(2a(u_{*}+u^{\prime})+2)\,.$ (5)

In the case of the static potential $X(Q)=rV_{\rm QCD}(1/r)$, with the choice
$a=1$, $u^{\prime}=-1/2$, $\tilde{X}(\tau)$ reduces to the standard momentum
space potential and the sine factor cancels the leading renormalon at
$u_{*}=1/2$ and simultaneously suppresses the dominant renormalons at
$u_{*}=3/2,\dots$ [15]. (In particular, in the large–$\beta_{0}$ approximation
IR renormalons are totally absent in the momentum-space potential.) In the
case of a general observable $X$, we can adjust the parameters $a$ and
$u^{\prime}$ to cancel or suppress the dominant renormalons of $\tilde{X}$.111
In addition we can vary the dimension of the Fourier transform to
$d^{n}\vec{x}$. For simplicity we set $n=3$.  The level of suppression depends
on the observable, but at least the first two renormalons closest to the
origin can always be suppressed.222 In principle we can include corrections to
eq. (5) generated by anomalous dimension and higher order terms of
$\alpha_{s}(Q)$ in $\delta X$, and we can adjust additional parameters to
cancel (or suppress more severely) renormalons including the corrections. On
the other hand, if these corrections are absent, the current procedure cancels
the corresponding renormalon exactly.

We can reconstruct $X(Q)$ by the inverse Fourier transform. After integrating
out the angular variables, naively we obtain

$X(Q)=\frac{r^{-2au^{\prime}-1}}{2\pi^{2}}\int_{0}^{\infty}d\tau\,\tau\sin(\tau
r)\tilde{X}(\tau)\,.$ (6)

The left-hand side has renormalons, while the dominant renormalons are
suppressed in $\tilde{X}$ on the right-hand side. The dominant renormalons are
generated by the $\tau$-integral of logarithms $\log(\mu^{2}/\tau^{2a})^{n}$
at small $\tau$ in the perturbative series for $\tilde{X}(\tau)$. When we
consider resummation of the logarithms by RG alternatively (as we will do in
practice), they stem from the singularity of the running coupling constant
$\alpha_{s}(\tau)$ in $\tilde{X}$ located on the positive $\tau$ axis. $\delta
X$ is generated by the integral surrounding the discontinuity of this
singularity. The expected power dependence on $\Lambda_{\overline{\rm MS}}$ is
obtained once we expand $\sin(\tau r)$ in $\tau$.

We propose to compute the renormalon-subtracted $X(Q)$ in the PV prescription,
$[X(Q)]_{\rm PV}$, in the following way. We take the principal value of the
above integral, that is, take the average over the contours $C_{\pm}(\tau)$:

$[X(Q)]_{\rm FTRS}=\frac{r^{-2au^{\prime}-1}}{2\pi^{2}}\int_{0,\rm
PV}^{\infty}d\tau\,\tau\sin(\tau r)\tilde{X}(\tau)\,.$ (7)

Here, $\tilde{X}(\tau)$ is evaluated by RG-improvement up to a certain order
(see below) and has a singularity (Landau singularity) on the positive
$\tau$-axis. In the case of $V_{\rm QCD}(r)$, this quantity coincides with the
renormalon-subtracted leading Wilson coefficient of $V_{\rm QCD}(r)$ used in
the analyses [14, 12]. Since renormalons of $\tilde{X}(\tau)$ are suppressed,
the only source of renormalons in eq. (6) is from the integral of the
singularity of $\tilde{X}(\tau)$. Then the PV prescription in Eq. (7), which
minimally regulates the singularity of the integrand (or more specifically,
that of the running coupling), corresponds to the minimally renormalon-
subtracted quantity (2). We will give an argument for equivalence of eqs. (2)
and (7) up to the N4LL approximation in [18]. (See also [15].) We note that
the equivalence holds only when $\tilde{X}(\tau)$ does not have renormalons.
If renormalons remain in $\tilde{X}(\tau)$, renormalons cannot be removed from
the $Q$-space quantity merely by the PV integral in Eq. (7), which only
regulates the Landau pole of the running coupling. Hence, the renormalon
suppression in $\tau$-space quantity [cf. Eq. (4)] is crucial for renormalon
subtraction.

$\tilde{X}(\tau)$ in the NkLL approximation is calculated in the following
manner. From the coefficients of the series up to $k$-th order perturbation
$X(Q)=\sum_{n=0}^{k}c_{n}\alpha_{s}(Q)^{n+1}$, $\tilde{X}$ is given by

$\tilde{X}(\tau)\to\tilde{X}^{(k)}(\tau)=\frac{4\pi}{\tau^{3+2au^{\prime}}}\sum_{n=0}^{k}\tilde{c}_{n}(0)\,\alpha_{s}(\tau^{a})^{n+1}\,,$
(8)

where $\tilde{c}_{n}(L_{\tau})$ is defined by the following relation

$F(\hat{H},L_{\tau})\sum_{n=0}^{\infty}c_{n}\alpha_{s}^{n+1}=\sum_{n=0}^{\infty}\tilde{c}_{n}(L_{\tau})\alpha_{s}^{n+1}\,,$
(9) $F(u,L_{\tau})=-\sin\left(\pi
a(u+u^{\prime})\right)\Gamma\left(2a(u+u^{\prime})+2\right)e^{L_{\tau}u}\,.$
(10)

Here,
$L_{\tau}=\log(\mu^{2}/\tau^{2a})\,,\,\hat{H}=-\beta(\alpha_{s})\frac{\partial}{\partial\alpha_{s}}\,$.
Thus, $\tilde{c}_{n}$ is given explicitly by the coefficients of the original
series $c_{0}\,,\,c_{1}\,,\,\cdots\,,\,c_{n}$ as

$\displaystyle\tilde{c}_{0}(L_{\tau})=F(0,0)c_{0}\,,\,$ (11)
$\displaystyle\tilde{c}_{1}(L_{\tau})=F(0,0)c_{1}+\partial_{u}F(0,L_{\tau})b_{0}c_{0}\,,\,$
(12) $\displaystyle~{}~{}~{}~{}~{}\vdots~{}~{}\,.$

The relation (9) is a straightforward consequence of the Fourier transform.
Since the renormalons in $\tilde{X}(\tau)$ are suppressed and its perturbative
series has a good convergence, it is natural to perform RG improvement in the
$\tau$ space, and a higher order (large $k$) $\tilde{X}^{(k)}(\tau)$ would be
a more accurate approximation of $\tilde{X}(\tau)$. Accuracy tests by going to
higher orders will be given in the test analyses below, where we estimate a
higher order effect by using the 5-loop coefficient of the QCD beta function.

In numerical evaluation of eq. (7) it is useful to decompose the PV integral
into two parts $X_{0}(Q)$ and $X_{\rm pow}(Q)$ by deforming the integral
contour in the complex plane:

$[X(Q)]_{\rm FTRS}=\bar{X}_{0}(Q)+\bar{X}_{\rm pow}(Q)\,,\,$ (13)
$\bar{X}_{0}(Q)=\frac{-r^{-2au^{\prime}-1}}{2\pi^{2}}\int_{0}^{\infty}dt\,t\,e^{-tr}\,{\rm
Im}\,[\tilde{X}(\tau=it)]\,,$ (14) $\bar{X}_{\rm
pow}(Q)=\frac{r^{-2au^{\prime}-1}}{4\pi^{2}i}\int_{C_{*}}d\tau\,\tau\cos(\tau
r)\tilde{X}(\tau)\,.$ (15)

The integration contour $C_{*}$ is shown in Fig. 1. $\bar{X}_{0}$ is defined
as the integral on the imaginary axis in the upper half $\tau$ plane
($\tau=it$), where $e^{-i\tau r}$ turns to a damping factor $e^{-tr}$.
$\bar{X}_{\rm pow}$ can be expanded by $r=Q^{-1/a}$ once the expansion of
$\cos(\tau r)$ is performed in $\tau$, and the coefficients of this power
series are real.333 When we take the integration contour as $C_{\pm}$ instead
of the PV integral in eq. (7), we also have power dependence with imaginary
coefficients whose sign depends on which contour is chosen. The power series
with imaginary coefficients is identified as renormalon uncertainties. They do
not appear in eq. (7), where the average over $C_{\pm}$ is taken.  It should
be expanded at least to the order of the eliminated renormalon. Then eq. (13)
gives the renormalon-subtracted prediction of a general observable $X(Q)$ with
an appropriate power accuracy of $1/Q$.444 The results scarcely change by
varying the truncation order beyond the minimum necessary order, in the tested
range of $Q$ in the examples below.

Figure 1: Integration contour $C_{*}$ in Eq. (15). The cut singularity, shown
by the wavy line, is caused by the Landau pole of the running coupling.

[In the case of the static potential, eq. (14) is equal to $U_{1}(r)$ of eq.
(61) in [14], while eq. (15) is equal to the renormalon-free part ${\cal
A}/r+{\cal C}r+\cdots$ of eq. (63) in the same paper with $Q=1/r$.]

The advantages of our method can be stated as follows. First, our formulation
to subtract renormalons works without knowing normalization constants
$N_{u_{*}}$ of the renormalons to be subtracted, following the above
calculation procedure. In other methods [19, 20, 21, 22, 23], in order to
subtract renormalons one needs normalization constants $N_{u_{*}}$ of the
corresponding renormalons. Normalization constants of renormalons far from the
origin are generally difficult to estimate. Although we certainly need to know
large order perturbative series to improve the accuracy of renormalon-
subtracted results, the above feature of our method practically facilitates
subtracting multiple renormalons even with small number of known perturbative
coefficients. Secondly, we can give predictions free from the unphysical
singularity around $Q\sim\Lambda_{\overline{\rm MS}}$ caused by the running of
the coupling, in the same way as the previous study of $V_{\rm QCD}(r)$ [14].
Since renormalons and the unphysical singularity are the main sources
destabilizing perturbative results at IR regions, the removal of these factors
is a marked feature of our method.

We will test validity of the above method (renormalon subtraction using
Fourier transform: “FTRS method”) by applying it to different observables in
the following.555 As seen in eq. (5), the Fourier transform generates
artificial UV renormalons in $\tilde{X}$. They are Borel summable, and we
perform the Borel summation whenever the induced UV renormalons are located
closer to the origin than the IR renormalons of our interest [18].

Adler function

The Adler function is defined from the photon vacuum polarization function
$\Pi(-Q^{2})$ in the Euclidean momentum region ($Q^{2}>0$):

$\displaystyle D(Q^{2})=12\pi^{2}Q^{2}\frac{d}{dQ^{2}}\Pi(-Q^{2})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}=Q^{2}\int_{0}^{\infty}\\!\\!\\!ds\,\frac{R(s)}{(s+Q^{2})^{2}}\,.$
(16)

The second line shows that it is also expressed by the $R$-ratio through
dispersion relation. The OPE is given by

$D(Q^{2})=C_{1}+2\pi^{2}\sum_{f}Q_{f}^{2}C_{GG}\,\frac{\left<{0}\,\right|\frac{\alpha_{s}}{\pi}G^{a\mu\nu}G^{a}_{\mu\nu}\left|\,{0}\,\right>}{Q^{4}}+\cdots,$
(17)

where $C_{1}$ and $C_{GG}$ denote the Wilson coefficients of the operators
$\boldmath\bf{1}$ and $G^{a\mu\nu}G^{a}_{\mu\nu}$, respectively. For
simplicity, we omit ${\cal O}(\alpha_{s})$ correction to $C_{GG}$
($C_{GG}=1$). We set $n_{f}=2$ and the quark masses to zero, hence we ignore
the quark condensate $m_{i}\langle\bar{\psi}_{i}{\psi_{i}}\rangle$.

We apply FTRS to the leading Wilson coefficient $C_{1}$ and subtract the
renormalon at $u=2$, i.e., at order $(\Lambda_{\rm QCD}/Q)^{4}$, so that the
local gluon condensate
$\langle(\alpha_{s}/\pi)G^{a\mu\nu}G^{a}_{\mu\nu}\rangle$ is also well defined
without renormalon uncertainty of the same order of magnitude. We choose
$a=1/2$, $u^{\prime}=-2$ such that the renormalons at $u=2,4,6,\cdots$ are
suppressed.
$\tilde{C}_{1}(\tau)=\sum_{n=0}^{4}\tilde{d}_{n}\alpha_{s}(\tau)^{n}$ in eq.
(7) is readily obtained from the NNNLO perturbative calculation of the Adler
function [24]. Here, $\alpha_{s}(\tau)$ denotes the 5-loop running coupling
constant in the $\overline{\rm MS}$ scheme [25]. (For simplicity we basically
use the 5-loop running constant.)

We change the scale by a factor 2 or $1/2$ from $\mu=\tau$ in
$\tilde{C}_{1}(\tau)$. Consistently with the theoretical expectation, we
verify that the scale dependence of $C_{1}^{\rm FTRS}$ is considerably smaller
than the fixed-order or RG-improved calculation without renormalon
subtraction, especially at low energy region.

Figure 2:  Comparison of Adler function by OPE in FTRS method and that by
phenomenological determination. They agree reasonably well.

As a phenomenological input in the $n_{f}=2$ case, we use the phenomenological
model [26] for $R(s)$, after projecting it to only the $n_{f}=2$ sector. It is
inserted into eq. (16) to obtain $D_{\rm pheno}$. This is compared with

$D_{\rm OPE}(Q^{2})=C_{1}^{\rm FTRS}(Q/\Lambda_{\overline{\rm
MS}})+A\,\biggl{(}\frac{\Lambda_{\overline{\rm MS}}}{Q}\biggr{)}^{4}\,.$ (18)

Here, the parameters $A=-25(12)$ and $\Lambda_{\overline{\rm MS}}=271(39)$ MeV
are determined by a fit to minimize $|D_{\rm pheno}-D_{\rm OPE}|$ in the range
$0.6~{}\text{GeV}^{2}\leq Q^{2}\leq 2~{}\text{GeV}^{2}$; see Fig. 2. (The
errors inside the brackets are estimated purely perturbatively by changing the
scale by a factor 2 or $1/2$ from $\mu=\tau$ in $\tilde{X}(\tau)$, for a
reference.) In this difference there scarcely remains a room for, e.g., a term
proportional to $(\Lambda_{\overline{\rm MS}}/Q)^{2}$. To improve accuracy by
going to higher orders, it may be important to deal with the $u=-1$ UV
renormalon of the Adler function properly. This speculation is based on our
analyses with the 5-loop perturbative coefficient in the large-$\beta_{0}$
approximation.

The above first analysis is fairly crude, with unknown uncertainties included
in the model cross section, etc.666 For instance, the UV renormalons (included
in the Adler function on its own) may give non-negligible contributions, and
we have not taken them into account.  Nevertheless, it may be informative to
compare the above result with the determination $\Lambda_{\overline{\rm
MS}}^{(2)}=310(20)$ MeV in the two-flavor lattice simulation [27], with which
we observe a rough consistency.

$B\to X_{u}\ell\bar{\nu}$ decay width

We consider the decay width of the process $B\to X_{u}\ell\bar{\nu}$. The OPE
is given within Heavy Quark Effective Theory (HQET) as an expansion in
$1/m_{b}$ [28]:

$\Gamma=\Gamma_{0}\biggl{[}\gamma_{1}-\frac{\mu_{\pi}^{2}}{2m_{b}^{2}}+\frac{\mu_{G}^{2}}{2m_{b}^{2}}+{\cal
O}(m_{b}^{-3})\biggr{]}\,,$ (19)

where $\Gamma_{0}$ denotes the partonic decay width without QCD corrections.
$\gamma_{1}$ denotes the Wilson coefficient of the operator $\boldmath\bf{1}$,
and the ${\cal O}(\Lambda_{\rm QCD}^{2})$ nonperturbative matrix elements are
denoted as

$\mu_{\pi}^{2}=\left<{B}\,\right|h_{b}^{\dagger}\vec{D}^{2}h_{b}\left|\,{B}\,\right>\,,~{}~{}~{}\mu_{G}^{2}=\left<{B}\,\right|h_{b}^{\dagger}\vec{\sigma}\cdot
g\vec{B}h_{b}\left|\,{B}\,\right>\,.$ (20)

For simplicity we omit the Wilson coefficients multiplying these matrix
elements [29, 30].

$\Gamma_{0}$ is proportional to $m_{b}^{5}$, and it is known that the $u=1/2$
renormalon is canceled when a short-distance $b$-quark mass is used instead of
the pole mass to express the decay width $\Gamma$ [5]. Thus, the leading
renormalon of $\gamma_{1}$ is at $u=1$.

We apply FTRS to $\gamma_{1}$ and subtract the dominant part of this
renormalon, such that $\mu_{\pi}^{2}$ and $\mu_{G}^{2}$ are also well defined.
We choose $a=1$, $u^{\prime}=-1$, hence renormalons at $u=1,2,3,\cdots$ are
suppressed.
$\tilde{\gamma}_{1}(\tau)=\sum_{n=0}^{2}\tilde{s}_{n}\alpha_{s}(\tau)^{n}$ is
derived from the two-loop calculation of $\Gamma$ [31], after expressing it by
the $\overline{\rm MS}$ mass $\overline{m}_{b}$.

Figure 3:  Comparison of the decay width for $B\to X_{u}\ell\bar{\nu}$ given
by $\gamma_{1}^{\rm FTRS}$ and by RG improvement without renormalon
subtraction (both at NLL and expressed by $\overline{m}_{b}$). The former is
less scale dependent than the latter at $\Lambda_{\overline{\rm
MS}}^{2}/\overline{m}_{b}^{2}\hbox{ \raise 3.0pt\hbox
to0.0pt{$>$}\raise-3.0pt\hbox{$\sim$} }0.02$.

We change the scale by a factor 2 or $1/2$ from $\mu=\tau$ in
$\tilde{\gamma}_{1}(\tau)$. For a hypothetically small value
$\overline{m}_{b}\sim 1$–1.5 GeV, the scale dependence of $\gamma_{1}^{\rm
FTRS}$ is smaller than the fixed-order or RG-improved calculation without
subtraction of the $u=1$ renormalon; see Fig. 3. However, in the relevant
region $\overline{m}_{b}\approx 4$ GeV, there are no significant differences
in the scale dependence.

The scale dependence of $\gamma_{1}^{\rm FTRS}$ at the current accuracy is
still rather large, where only the corrections up to ${\cal
O}(\alpha_{s}^{2})$ are known. We expect a healthy convergence behavior of
$\gamma_{1}^{\rm FTRS}$ since we subtracted the dominant renormalons. We
estimate higher order results by including log dependent terms, dictated by RG
equation, at the 5-loop level. We consider RG improvement for perturbative
series in $\tau$ space rather than $\bar{m}$ space, in accordance with the
concept that there are no IR renormalons in the $\tau$-space quantity. (Note
that RG improvements in two spaces are not equivalent when we do not know
exact perturbative coefficients. In this treatment, the $u=1$ renormalon is
induced in the perturbative series in $\bar{m}$ space.) Then, a large
uncertainty arises in the result of the usual RG improvement. By FTRS, we
obtain a significantly small uncertainty thanks to the renormalon subtraction,
while the analyses are base on the same perturbative series. However, in order
to have an insight into the true size of the $u=1$ renormalon, we need more
terms of the perturbative series.

$M_{B}$ and $M_{D}$

In HQET the mass of a heavy-light system $H=B^{(*)},D^{(*)}$ is given as an
expansion in $1/m_{h}$ ($h=b,c$) as [32]

$M_{H}=m_{h}+\bar{\Lambda}+\frac{\mu_{\pi}^{2}}{2m_{h}}+w(s)\frac{\mu_{G}^{2}}{2m_{h}}+{\cal
O}(m_{h}^{-2})\,.$ (21)

Here, $m_{h}$ denotes the pole mass of $h$; $\bar{\Lambda}$ represents the
contribution of ${\cal O}(\Lambda_{\rm QCD})$ from light degrees of freedom;
$s=0,1$ denotes the spin of $H$, and $w(0)=-1$, $w(1)=1/3$. For simplicity we
omit the Wilson coefficients multiplying $\mu_{\pi}^{2}$ and $\mu_{G}^{2}$. By
the heavy quark symmetry $\bar{\Lambda}$, $\mu_{\pi}^{2}$ or $\mu_{G}^{2}$ are
common for $h=b,c$. The renormalons in $m_{h}$ become manifest when we express
it by a short-distance mass (we use the $\overline{\rm MS}$ mass
$\overline{m}_{h}$).

We apply FTRS to $m_{h}$ and subtract the renormalons at $u=1/2$ and $u=1$. We
choose $a=2$, $u^{\prime}=-1/2$, hence renormalons at $u=1/2,1,3/2,2,\cdots$
are suppressed. Projecting out $\mu_{G}^{2}$, we define

$\displaystyle\overline{M}_{H,{\rm OPE}}\equiv\frac{1}{4}\bigl{(}M_{H,{\rm
OPE}}^{s=0}+3M_{H,{\rm OPE}}^{s=1}\bigr{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\equiv m_{h}^{\rm
FTRS}+\bar{\Lambda}^{\rm FTRS}+\frac{(\mu_{\pi}^{2})^{\rm FTRS}}{2m_{h}^{\rm
FTRS}}\,.$ (22)

Here, $m_{h}^{\rm FTRS}$ denotes the principal value of the pole mass
(expressed by the $\overline{\rm MS}$ mass, calculated in the FTRS method).

First we examine the scale dependence of $m_{h}^{\rm FTRS}$ and compare with
the fixed-order and RG-improved calculation without renormalon subtraction. We
use the perturbative series up to ${\cal O}(\alpha_{s}^{4})$ for
pole-$\overline{\rm MS}$ mass relation [33, 34]. Since the contribution of the
$u=1/2$ renormalon is known to be significant for $m_{h}$, (only) in this
comparison, $dm_{h}/d\overline{m}_{h}$ is used to cancel the $u=1/2$
renormalon. There have been estimates that the contribution of the $u=1$
renormalon is small [21]. (For example, it is absent in the large-$\beta_{0}$
approximation.) Consistently with such estimates we observe no significant
difference in the scale dependence if $\overline{m}_{h}\hbox{ \raise
3.0pt\hbox to0.0pt{$>$}\raise-3.0pt\hbox{$\sim$} }1$ GeV, even after
subtracting the $u=1$ renormalon.

We determine $\bar{\Lambda}$ and $(\mu_{\pi}^{2})^{\rm FTRS}$ in eq. (22) by
comparing $\overline{M}_{H,{\rm OPE}}$ to the experimental values of
$\overline{M}_{B}$ and $\overline{M}_{D}$. In this analysis we include non-
zero charm mass effects to the bottom mass (and non-decoupling effects from
bottom to the charm mass) up to ${\cal O}(\alpha_{s}^{3})$ [35]. We use the
values of $\overline{m}_{b}$, $\overline{m}_{c}$, $\alpha_{s}$ from the
Particle Data Group as input parameters [36]. We obtain

$\displaystyle\bar{\Lambda}_{\rm
FTRS}=0.495(15)_{\mu}(49)_{\overline{m}_{b}}(12)_{\overline{m}_{c}}(13)_{\alpha_{s}}(0)_{\rm
f.m.}\text{GeV}\,,$ (23) $\displaystyle(\mu_{\pi}^{2})_{\rm
FTRS}=-0.12(13)_{\mu}(15)_{\overline{m}_{b}}(11)_{\overline{m}_{c}}(4)_{\alpha_{s}}(0)_{\rm
f.m.}\text{GeV}^{2},$ (24)

where the errors denote, respectively, that from the scale dependence for
$\mu=2\tau$ or $\mu=\tau/2$, from the errors of the input $\overline{m}_{b}$,
$\overline{m}_{c}$, $\alpha_{s}$, and from the finite $m_{c}$ corrections
(non-decoupling bottom effects) in loops for $m_{b}$ ($m_{c}$).

The error from the scale dependence is a measure of perturbative uncertainty.
From eq. (5) and our parameter choice, the renormalons at $u=1/2$ and $1$
should be removed in our result. For $\bar{\Lambda}$, about 3 per cent error
shows successful subtraction of the $u=1/2$ renormalon. On the other hand, the
scale dependence of $(\mu_{\pi}^{2})^{\rm FTRS}$ is not smaller than ${\cal
O}(\Lambda_{\rm QCD}^{2})$. We estimate that this is not because of the
contribution from the $u=1$ renormalon but due to the insufficient number of
known terms of the perturbative series. Examination using the 5-loop beta
function combined with an estimated 5-loop coefficient in the
large-$\beta_{0}$ approximation indeed indicates that this uncertainty can be
reduced at higher order [$\bar{\Lambda}^{\rm(est)}_{\rm
FTRS}=0.488(4)_{\mu}\,\text{GeV}$ and $(\mu_{\pi}^{2})^{\rm(est)}_{\rm
FTRS}=-0.09(7)_{\mu}\text{GeV}^{2}$].

In summary, in all the above analyses we observed good consistency with
theoretical expectations. In particular, from the $B,\,D$ meson masses we
obtained

$\displaystyle\bar{\Lambda}_{\rm FTRS}=0.495\pm
0.053~{}\text{GeV}\,,~{}~{}~{}~{}~{}(\mu_{\pi}^{2})_{\rm FTRS}=-0.12\pm
0.23~{}\text{GeV}^{2}\,,$ (25)

by subtracting ${\cal O}(\Lambda_{\rm QCD})$ and ${\cal O}(\Lambda_{\rm
QCD}^{2})$ renormalons simultaneously for the first time. In general, it is
desirable to know more terms of the relevant perturbative series in order to
make conclusive statements about the effects of subtracting renormalons beyond
the ${\cal O}(\Lambda_{\rm QCD})$ renormalon. The above analyses also show
that the 5-loop QCD beta function is a crucial ingredient in improving
accuracies. We anticipate that the FTRS method can be a useful theoretical
tool for precision QCD calculations in the near future.

## Acknowledgements

Y.H. acknowledges support from GP-PU at Tohoku University. The work of Y.H.
was also supported in part by Grant-in-Aid for JSPS Fellows (No. 21J10226)
from MEXT, Japan. The works of Y.S. and H.T., respectively, were supported in
part by Grant-in-Aid for scientific research (Nos. 20K03923 and 19K14711) from
MEXT, Japan.
Note added: After completion of this work, we learned that the ${\cal
O}(\alpha_{s}^{3})$ corrections to the $B$ meson semileptonic decay width has
recently been computed in expansion in $(1-m_{c}/m_{b})$ [37]. We will present
the study including this effect in [18].

## References

  * [1] A. Pineda, Ph.D. Thesis (1998).
  * [2] A. H. Hoang, M. C. Smith, T. Stelzer and S. Willenbrock, Phys. Rev. D 59, 114014 (1999) [arXiv:hep-ph/9804227 [hep-ph]].
  * [3] M. Beneke, Phys. Lett. B 434, 115-125 (1998) [arXiv:hep-ph/9804241 [hep-ph]].
  * [4] P. Ball, M. Beneke and V. M. Braun, Phys. Rev. D 52 (1995), 3929-3948 [arXiv:hep-ph/9503492 [hep-ph]].
  * [5] I. I. Y. Bigi, M. A. Shifman, N. G. Uraltsev and A. I. Vainshtein, Phys. Rev. D 50 (1994) 2234 [hep-ph/9402360].
  * [6] M. Neubert and C. T. Sachrajda, Nucl. Phys. B 438 (1995), 235-260 [arXiv:hep-ph/9407394 [hep-ph]].
  * [7] Y. Kiyo, G. Mishima and Y. Sumino, Phys. Lett. B 752 (2016) 122 Erratum: [Phys. Lett. B 772 (2017) 878] [arXiv:1510.07072 [hep-ph]].
  * [8] A. Bazavov et al. [Fermilab Lattice, MILC and TUMQCD], Phys. Rev. D 98 (2018) no.5, 054517 [arXiv:1802.04248 [hep-lat]].
  * [9] A. H. Hoang, Z. Ligeti and A. V. Manohar, Phys. Rev. D 59 (1999) 074017 [hep-ph/9811239].
  * [10] A. Alberti, P. Gambino, K. J. Healey and S. Nandi, Phys. Rev. Lett. 114 (2015) no. 6, 061802 [arXiv:1411.6560 [hep-ph]].
  * [11] A. Bazavov, N. Brambilla, X. Garcia i Tormo, P. Petreczky, J. Soto and A. Vairo, Phys. Rev. D 86 (2012), 114031 [arXiv:1205.6155 [hep-ph]].
  * [12] H. Takaura, T. Kaneko, Y. Kiyo and Y. Sumino, Phys. Lett. B 789, 598-602 (2019) [arXiv:1808.01632 [hep-ph]].
  * [13] C. Ayala, X. Lobregat and A. Pineda, [arXiv:2009.01285 [hep-ph]].
  * [14] Y. Sumino, Phys. Rev. D 76, 114009 (2007) [arXiv:hep-ph/0505034 [hep-ph]].
  * [15] Y. Sumino and H. Takaura, JHEP 05 (2020), 116 [arXiv:2001.00770 [hep-ph]].
  * [16] M. Beneke, Phys. Rept. 317, 1 (1999). [hep-ph/9807443].
  * [17] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Phys. Rev. D 60, 091502 (1999) [arXiv:hep-ph/9903355 [hep-ph]].
  * [18] Y. Hayashi, Y. Sumino and H. Takaura, in preparation.
  * [19] T. Lee, Phys. Rev. D 67, 014020 (2003) [arXiv:hep-ph/0210032 [hep-ph]].
  * [20] C. Ayala, X. Lobregat and A. Pineda, Phys. Rev. D 99 (2019) no.7, 074019 [arXiv:1902.07736 [hep-th]].
  * [21] C. Ayala, X. Lobregat and A. Pineda, Phys. Rev. D 101 (2020) no.3, 034002 [arXiv:1909.01370 [hep-ph]].
  * [22] C. Ayala, X. Lobregat and A. Pineda, [arXiv:2005.12301 [hep-ph]].
  * [23] H. Takaura, JHEP 10, 039 (2020) [arXiv:2002.00428 [hep-ph]].
  * [24] P. A. Baikov, K. G. Chetyrkin, J. H. Kuhn and J. Rittinger, Phys. Lett. B 714, 62-65 (2012) [arXiv:1206.1288 [hep-ph]].
  * [25] P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, Acta Phys. Polon. B 48 (2017), 2135 [arXiv:1711.05592 [hep-ph]].
  * [26] D. Bernecker and H. B. Meyer, Eur. Phys. J. A 47 (2011) 148 [arXiv:1107.4388 [hep-lat]].
  * [27] S. Aoki et al. [Flavour Lattice Averaging Group], Eur. Phys. J. C 80 (2020) no.2, 113 [arXiv:1902.08191 [hep-lat]].
  * [28] A. V. Manohar and M. B. Wise, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 10 (2000) 1.
  * [29] T. Becher, H. Boos and E. Lunghi, JHEP 12 (2007), 062 [arXiv:0708.0855 [hep-ph]].
  * [30] A. Alberti, P. Gambino and S. Nandi, JHEP 01 (2014), 147 [arXiv:1311.7381 [hep-ph]].
  * [31] A. Pak and A. Czarnecki, Phys. Rev. D 78 (2008) 114015 [arXiv:0808.3509 [hep-ph]].
  * [32] A. F. Falk and M. Neubert, Phys. Rev. D 47 (1993), 2965-2981 [arXiv:hep-ph/9209268 [hep-ph]].
  * [33] P. Marquard, A. V. Smirnov, V. A. Smirnov and M. Steinhauser, Phys. Rev. Lett. 114 (2015) no.14, 142002 [arXiv:1502.01030 [hep-ph]].
  * [34] P. Marquard, A. V. Smirnov, V. A. Smirnov, M. Steinhauser and D. Wellmann, Phys. Rev. D 94 (2016) no.7, 074025 [arXiv:1606.06754 [hep-ph]].
  * [35] M. Fael, K. Schönwald and M. Steinhauser, JHEP 10 (2020), 087 [arXiv:2008.01102 [hep-ph]].
  * [36] P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020)
  * [37] M. Fael, K. Schönwald and M. Steinhauser, “Third Order Corrections to the Semi-Leptonic $b\to c$ and the Muon Decays,” [arXiv:2011.13654 [hep-ph]].


# Math-KG: Construction and Applications of Mathematical
Knowledge Graph

Jianing Wang
East China Normal University
<EMAIL_ADDRESS>Corresponding Author.

###### Abstract

Recently, the explosion of online education platforms makes a success in
encouraging us to easily access online education resources. However, most of
them ignore the integration of massive unstructured information, which
inevitably brings the problem of information overload and knowledge trek. In
this paper, we proposed a mathematical knowledge graph named Math-KG, which
automatically constructed by the pipeline method with the natural language
processing technology to integrate the resources of the mathematics. It is
built from the corpora of Baidu Baike, Wikipedia. We implement a simple
application system to validate the proposed Math-KG can make contributions on
a series of scenes, including faults analysis and semantic search. The system
is publicly available at GitHub 111https://github.com/wjn1996/Mathematical-
Knowledge-Entity-Recognition..

##  1 Introduction

Online education (OE) is a novel teaching and learning concept, which makes
full use of Internet technology to realize education intelligence. By applying
artificial intelligence (AI) technology in the teaching process, it can
provide students with personalized learning guidelines and improve teaching
efficiency. In the context of the rapid development of AI technology, the
concept of intellectual adaptation has begun to be integrated into the field
of education.

However, the number of online education platforms is becoming bigger and
bigger, most of them ignore the integration of unstructured information, which
limit to a simple classification and display of resources. In this case, it
inevitably brings the problem of information overload and knowledge trek. For
example, in mathematics, massive knowledge points provided by different
platforms are independent of each other, which results in the inconsistency of
knowledge and causes the confusion of students’ learning process. In addition,
little work concentrates on the integration of educational resources, because
the field of education contains a large number of subjects, and each of them
has great differences.

Figure 1: Construction and applications of Math-KG

The knowledge Graph was originally proposed by Google Inc. in 2012 to enable
precise question and answer and semantic search engines by allowing computers
to understand the semantics of text. The knowledge Graph aims to constructs
the knowledge base by extracting a series of triples from unstructured or
semi-structured texts. At present, many advanced knowledge extraction tools
Manning et al. (2014); Gardner et al. (2018); Li et al. (2020) are proposed to
construct knowledge graph, which has been widely used in the fields of
medicine Goodwin and Harabagiu (2013), financial Fu et al. (2018),
agricultural Sun et al. (2017), etc., but it is still a blank in mathematical
education.

In order to further promote the development of online mathematical education,
we propose a knowledge graph-based Math-KG to integrate knowledge resources.
Specifically, we define two entity types and six relation types at first, and
then propose a pipeline method for Math-KG construction. We obtain the
unstructured text corpus through the web crawler and created the dataset via
distant supervision strategy. We utilize the natural language processing
technology tool to realize the named entity recognition, entity relation
classification and knowledge fusion, and then save them in the graph database.
The specific construction process as described in section 2. We implement a
simple application system to validate the feasibility of knowledge graph,
including faults analysis and semantic search. The main contributions of this
paper are as follows:

  * •

We propose a pipeline method of mathematical knowledge graph construction, and
build a large-scale knowledge base named Math-KG.

  * •

We conduct some experiments to evaluate each part of our pipeline method, and
implement a simple application system to make validation of some educational
scenes, such as faults analysis and semantic search.

##  2 Method of Math-KG Construction

As illustrated in Fig 1, the pipeline of Math-KG construction contains two
components: Corpus Preprocessing and Information Extraction.

###  2.1 Corpus Preprocessing

In order to construct the Math-KG, it is necessary to complete the acquisition
and preprocessing of the mathematical corpus. The preprocessing stages mainly
include the following parts:

Resources Acquisition. We crawl semi-structured and unstructured data 222For
the pictures and videos, we only extract the hyperlinks as an attribution from
web resources includes Baidu Baike333https://baike.baidu.com/.,
Wikipedia444https://www.wikidata.org., etc. Specifically, for one entity
token, we capture five resources consists of categories, abstracts, main
bodies, infobox blocks and corresponding indexes.

Resources Process. For the categories, infobox blocks and corresponding
indexes, we direct view them as attributes of the entity token. The other
resources contain raw text and formulas. The texts are split into sentences
and processed by Spacy tool555https://spacy.io/. which consists of word
tokenization, part-of-speech (POS) tagging, dependency tree parse (DTP), etc.
The mathematical formulas are transformed into latex characters by
MathOCR666https://sourceforge.net/projects/mathocr/.. We release our
preprocessing code on GitHub
777https://github.com/wjn1996/scrapy_for_zh_wiki..

relation |  instruction |  example |  entity pairs
---|---|---|---
dependencies |  two knowledge points have sequential dependence in the learning process |  Learning addition${}_{e_{1}}$ helps us understand the definition of subtraction${}_{e_{2}}$. |  Dep($e_{1}$,$e_{2}$)
affiliation |  two knowledge points have containing and be contained relations |  An isosceles triangle${}_{e_{1}}$ is a special kind of triangle${}_{e_{2}}$. |  Aff($e_{1}$,$e_{2}$)
equivalence |  two entities in term of different name have the same meaning |  orthogon${}_{e_{1}}$ belongs to plane geometry which called rectangle${}_{e_{2}}$. |  Equ($e_{1}$,$e_{2}$)
opposite |  two knowledge points are conceptually opposite |  if a number is even${}_{e_{1}}$, it must not be odd${}_{e_{2}}$. |  Ant($e_{1}$,$e_{2}$)
synonyms |  two knowledge points have similar meanings |  Circles are both axial symmetric${}_{e_{1}}$ and centrally symmetric figure${}_{e_{2}}$. |  Syn($e_{1}$,$e_{2}$)
has properties |  one knowledge point is a property of another |  the area${}_{e_{1}}$ of a parallelogram${}_{e_{2}}$ is base times height. |  Pro($e_{1}$,$e_{2}$)
Table 1: Instruction and example of six relations

###  2.2 Information Extraction

Information extraction aims to extract structure triple and construct a large-
scale graph. It consists of knowledge entity recognition, entity relation
classification and knowledge fusion.

Knowledge Entity Recognition. Knowledge entity recognition can be seen as a
sequence labeling task. We define two types of knowledge entities, Conceptual
Entity and Theorem Entity. The Concept Entity indicates that the knowledge
point is an existence concept, which is represented by CON, such as triangle,
trigonometric function, etc. The Theorem Entity represents the theorem, method
or rule based on the Concept Entity , which labeled as LEG, for example, the
entity pythagorean theorem, which can be used to solve the hypotenuse of right
triangle, belongs to Theorem Entity. In order to improve the accuracy of
identification, we prefix the first character of each entity with ’B’ and mark
the rest with ’I’, other non-entities are tagged as ’O’. We first use the
entries in the encyclopedia as seeds and recall some high-quality entities by
heuristic rules, such as tf-idf, pattern matching or word2vec Mikolov et al.
(2013) similarity. Then, we construct a small dataset via distant supervision
Mintz et al. (2009) method by aligning the corpus with the recall entities set
and fine-tune by the BERT modelDevlin et al. (2018).

Entity Relation Classification. We consider relation extraction as a
classification task, which is a prerequisite for building triples. We define
six classes of relation for knowledge points, namely: dependencies,
affiliation, equivalence, antisense, synonyms, has properties. (Table 1 gives
the instruction and example of each relation.) These relations have two
directions so that Math-KG will be composed of 12 kinds of edge with an NA
888NA denotes that it has no relation between two entities.. To recognize the
relation between entities, we first use the explicit triples as seeds, which
from categories, infobox blocks and corresponding indexes that directly
provided by Baidu Baike and Wikipedia. Then, we defined some pattern rules,
and use the self-supervision method to realize pattern discovery to recall
some triples. In order to improve the quality of the seeds, we manually
performed a rough inspection of them. Although the rule-based method can
obtain high accuracy, it still lacks generalization ability. Therefore, we
follow distant supervision strategy Mintz et al. (2009) to build a large-scale
dataset by heuristically aligning with raw text, and then train by BERT
modelDevlin et al. (2018) to make the classification. However, distant
supervision assumes that all the text with the same entity pair can express
the same relation, which brings many noisy labeling data Mintz et al. (2009).
To alleviate this problem, we add a sentence-level attention layer on the top
of the BERT encoder and use the embeddings of special [CLS] token to make the
classification.

Knowledge Fusion. In order to improve the accuracy of the knowledge graph, we
used entity linker TagMe 999https://sobigdata.d4science.org/group/tagme. and
reference resolution technique Neuralcoref
101010https://github.com/huggingface/neuralcoref. to align the constructed
triples with some existing knowledge bases to improve confidence. In addition,
we unify these resources manually to ensure the consistency of knowledge.

Data set | train | dev | test | all
---|---|---|---|---
KER | 6265 | 1300 | 1385 | 8950
ERC | 9356 | 2000 | 2009 | 13365
Table 2: The statistics of knowledge entity recognition (KER) and entity relation classification (ERC) data set. Data set | precision % | recall % | F1 %
---|---|---|---
KER | 92.40 | 92.85 | 92.62
ERC | 84.30 | 85.90 | 85.09
Table 3: The evaluation results of knowledge entity recognition (KER) and
entity relation classification (ERC) data set.

(a) Faults Analysis

(b) Semantic Search

Figure 2: Applications demonstration

##  3 Evaluations of Math-KG

###  3.1 Experimental Evaluations

We release two core data set and corresponding Tensorflow codes 111111the KER
and ERC data set and corresponding codes released on the github:
https://github.com/wjn1996/Mathematical-Knowledge-Entity-Recognition. to
support the pipeline of construction, including knowledge entity recognition
(KER) and entity relation classification (ERC). The statistics of each data
set as shown in Table 2. We use BERT+CRF model to train on the KER data set.
We use BERT+Attention+Softmax to train on the ERC data set. We use large-based
BERT for both two experiments, and the batch size is set as 64, the epoch is
set as 30. We choose cross-entropy loss function to train. The F1 values of
KER and ERC are 92.62% and 85.09%, respectively. The results indicates that
the current algorithm still has room for improvement. At last, we obtain 1905
mathematical knowledge point entities and 8,019 triples. Among them, there are
1337 concept entities and 568 theorem entities. The number of triples with
relation dependencies, affiliation, equivalence, antisense, synonyms and has
properties are 2016, 2609, 76, 301, 493 and 2524, respectively.

###  3.2 Application Evaluations

In order to validate the usefulness of Math-KG, we design two application
scenes, such as faults analysis and semantic search.

Faults Analysis. Faults analysis aims to speculate the knowledge points that
student fails to master by analyzing the historical question answer results.
We use the baidu, google search engines to extract some mathematics practices.
We design a simple method to make faults analysis based on historical answer
results. For one student, we count the number of correct and incorrect
questions for each knowledge point. If the number of incorrect questions is
bigger than the correct, this student fails to master this knowledge.
Therefore, we can obtain a sub-graph for one practice, which contains other
related knowledge points. If one knowledge point often appears in the sub-
graph, it can be regarded as the source of the faults, and then be recommended
to the student. As shown in Fig 2(a), the question labeled as two knowledge
consists of Quadratic function and Quadratic equation of one variable. In
order to show more intuition, we transform the sub-graph of these knowledge
points into a tree, where the root denotes the current knowledge point of the
practice, and the leaf node denotes the related knowledge points. We calculate
the score based on the statistics at each leaf to represent the relevance to
the root knowledge point. The few red paths represent the inference evidences.

Semantic Search. We also provide semantic search, which is also a research
scope called knowledge base question answering (KBQA) Tan et al. (2016). As we
all known, knowledge graph can be viewed as a semantic network Chen et al.
(2019), which can be answered the question by multi-hop inferences.
Especially, given a question $s$, we use the Spacy tool to recognize the
corresponding topic entity and recall the multi-hop sub-graph from Math-KG. We
pre-trained our Math-KG by TransE Fan et al. (2014) to obtain the embeddings
of each knowledge entity. Then, we calculate the similarity based on a simple
neural network between the question and each candidate entities, and obtain
the top results. We display the path information from the topic entity to the
resulting entity. As shown in Fig 2(b), when querying the circumscribed circle
radius of a triangle, the system recognizes the topic entity triangle and find
a multi-hop path (the triangle, has, the circumscribed circle radius).

##  4 Conclusion

In this paper, in order to alleivate the problem of information overload and
knowledge trek in mathematics, we propose a pipeline method of Math-KG
construction to integrate the mathematical knowledge resources. We conduct
experiments to evaluate the performance of construction, and implement two
applications to validate the usefulness of Math-KG, which provides students
with faults analysis and semantic search. In future researches, we will put
our efforts into improving the accuracy and efficiency of information
extraction, and apply it to other disciplines, including physics, chemistry,
etc.

## References

  * Chen et al. (2019) Yongrui Chen, Huiying Li, and Zejian Xu. 2019.  _Convolutional Neural Network-Based Question Answering Over Knowledge Base with Type Constraint: Third China Conference, CCKS 2018, Tianjin, China, August 14–17, 2018, Revised Selected Papers_, pages 28–39.
  * Devlin et al. (2018) Jacob Devlin, Ming Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.
  * Fan et al. (2014) M. Fan, Q. Zhou, E. Chang, and Fang Zheng. 2014.  Transition-based knowledge graph embedding with relational mapping properties.  _PACLIC_ , pages 328–337.
  * Fu et al. (2018) Xiaoyi Fu, Xinqi Ren, Ole J. Mengshoel, and Xindong Wu. 2018.  Stochastic optimization for market return prediction using financial knowledge graph.  In _2018 IEEE International Conference on Big Knowledge (ICBK)_.
  * Gardner et al. (2018) Matt Gardner, Joel Grus, Mark Neumann, Oyvind Tafjord, Pradeep Dasigi, Nelson F. Liu, Matthew Peters, Michael Schmitz, and Luke Zettlemoyer. 2018.  AllenNLP: A deep semantic natural language processing platform.  In _Proceedings of Workshop for NLP Open Source Software (NLP-OSS)_ , pages 1–6, Melbourne, Australia. Association for Computational Linguistics.
  * Goodwin and Harabagiu (2013) Travis Goodwin and Sanda M. Harabagiu. 2013.  Automatic generation of a qualified medical knowledge graph and its usage for retrieving patient cohorts from electronic medical records.  In _IEEE Seventh International Conference on Semantic Computing_.
  * Li et al. (2020) Manling Li, Alireza Zareian, Ying Lin, Xiaoman Pan, Spencer Whitehead, Brian Chen, Bo Wu, Heng Ji, Shih-Fu Chang, Clare Voss, Daniel Napierski, and Marjorie Freedman. 2020.  GAIA: A fine-grained multimedia knowledge extraction system.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics: System Demonstrations_ , pages 77–86, Online. Association for Computational Linguistics.
  * Manning et al. (2014) Christopher Manning, Mihai Surdeanu, John Bauer, Jenny Finkel, Steven Bethard, and David McClosky. 2014.  The Stanford CoreNLP natural language processing toolkit.  In _Proceedings of 52nd Annual Meeting of the Association for Computational Linguistics: System Demonstrations_ , pages 55–60, Baltimore, Maryland. Association for Computational Linguistics.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeffrey Dean. 2013.  Distributed representations of words and phrases and their compositionality.  _Advances in Neural Information Processing Systems_ , 26:3111–3119.
  * Mintz et al. (2009) M. Mintz, S. Bills, R. Snow, and D. Jurafsky. 2009.  Distant supervision for relation extraction without labeled data.  In _ACL_ , pages 1003–1011. ACL.
  * Sun et al. (2017) Xiang Sun, Huarui Wu, Peng Hao, and Qingxue Li. 2017.  Research on information integration method of agricultural products producing and managing based on knowledge graph.  In _International Conference on Computer and Computing Technologies in Agriculture_.
  * Tan et al. (2016) Ming Tan, Cicero dos Santos, Bing Xiang, and Bowen Zhou. 2016.  Improved representation learning for question answer matching.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 464–473, Berlin, Germany. Association for Computational Linguistics.


# Do Counterfactual Examples Complicate Adversarial Training?

Eric Yeats1 Cameron Darwin Eduardo Ortega2 Frank Liu3 Hai Li1
1Duke University 2Arizona State University 3Old Dominion University

###### Abstract

We leverage diffusion models to study the robustness-performance tradeoff of
robust classifiers. Our approach introduces a simple, pretrained diffusion
method to generate low-norm counterfactual examples (CEs): semantically
altered data which results in different true class membership. We report that
the confidence and accuracy of robust models on their clean training data are
associated with the proximity of the data to their CEs. Moreover, robust
models perform very poorly when evaluated on the CEs directly, as they become
increasingly invariant to the low-norm, semantic changes brought by CEs. The
results indicate a significant overlap between non-robust and semantic
features, countering the common assumption that non-robust features are not
interpretable.

Figure 1: Conceptual depiction of the relationship between accuracy of
robustly trained models with proximity to counterfactual examples (CEs).
Stronger adversarial training inevitably leads to misclassification of some
clean training data, incurring downstream test performance loss. We
hypothesize that adversarially trained models are forced to become invariant
to some semantic features due to the nearby presence of true CEs.

##  1 Introduction

Leading theory by Ilyas et al. [10] asserts that adversarial vulnerability
arises from reliance of DNN classifiers on non-robust features: well-
generalizing, yet brittle features which are not comprehensible to humans. As
robust models must become invariant to this subset of predictive features,
they inevitably lose performance [20]. However, humans appear to break this
trend by simultaneously maintaining good performance and robustness on
adversarial data [16].

The existence of adversarial perturbations reflects the stark differences
between the perception of DNNs and our own. DNNs are notoriously
uninterpretable [17, 23], and there is a need for interpretable decision-
making in high-stakes scenarios. One model-agnostic way to interpret
classifier decision-making is to generate counterfactual examples (CEs):
subtle, semantic changes to an input datum which would result in a classifier
predicting a target class [21]. Following their definition from psychology,
CEs for an image-label pair $(x,y)$ can be described by the statement “if
$y^{\prime}$ were the true class of image $x$ (instead of $y$), then $x$ would
look like $x^{\prime}$.” While CEs can help users understand the decision-
making of DNN classifiers, practical methods to generate CEs are
problematically similar to those for adversarial attacks [5, 15]. Hence,
methods to generate CEs for DNNs often require robust models in some form [2,
1].

We develop new understanding of the robustness-performance tradeoff through
our study of the semantic feature distributions learned by robust classifiers.
Our study leverages independently generated CE datasets which we create with a
simple, pretrained diffusion model technique. We report that robust models,
unlike standard (non-robust) models, are more likely to lose confidence on and
misclassify clean training data which have nearby CEs, and robust models
become invariant to the semantic changes brought by CEs. Contrary to common
assumptions, our findings suggest significant overlap between non-robust and
semantically meaningful features. This conflict speaks to the complexity of
the adversarial problem and motivates alternative approaches to robustness
which can resolve or avert it. ††Corresponding Author Email:
<EMAIL_ADDRESS>
##  2 Background

#### Adversarial Examples and Adversarial Training

Adversarial examples are minute, malicious data perturbations which alter
classifier predictions [19]. For a classifier
$f_{\theta}:\mathcal{X}\rightarrow\mathcal{Y}$ where
$\mathcal{X}\in\mathbb{R}^{n}$ and $\mathcal{Y}$ is the space of categorical
distributions of support cardinality $k$, $\tilde{y}$-targeted adversarial
examples $\tilde{x}$ are often computed by altering an input $x$ to descend
the negative log loss [6]
$\mathcal{L}:\Theta\times\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}^{+}$:

$\tilde{x}=\operatorname*{arg\,min}_{x^{\prime}\in\mathcal{B}_{\varepsilon}(x)}\mathcal{L}(\theta,x^{\prime},\tilde{y}),$
(1)

where the neighborhood $\mathcal{B}_{\varepsilon}(x)$ is commonly defined with
the $L^{\infty}$ or $L^{2}$ norm. In this work, we focus on $L^{2}$
adversarial examples. Adversarial training consists of finding the optimal
model parameters $\theta^{*}$ for the robust objective [12]:

$\theta^{*}=\operatorname*{arg\,min}_{\theta}\mathbb{E}_{x,y}\max_{x^{\prime}\in\mathcal{B}_{\varepsilon}(x)}\mathcal{L}(\theta,x^{\prime},y).$
(2)

There are other robustness methods such as gradient regularization [25, 4] and
randomized smoothing [3], but adversarial training is the most well-known and
successful.

#### Diffusion Models and Classifier-free Guidance

Diffusion models learn a sequence of functions $\epsilon_{\theta}(x_{t},t)$
which predict the noise added to data by a diffusion process at time
$t\in[0,1]$ [9]. To generate new data, $\epsilon_{\theta}(x_{t},t)$ is used to
iteratively denoise samples $x_{t}$ initially drawn from a known noise prior.
The noise prediction function $\epsilon_{\theta}(x_{t},t)$ is known to
approximate the negative score at diffusion step $t$ [18, 24] (i.e.,
$\epsilon_{\theta}(x_{t},t)\approx-\sigma_{t}\nabla_{x}\log p_{t}(x_{t})$).
For class-conditional diffusion models, this property can be used to reproduce
the effect of adding targeted classifier gradients [8]:

$\epsilon^{w}_{\theta}(x_{t},y,t)\coloneq(w+1)\epsilon_{\theta}(x_{t},y,t)-w\epsilon_{\theta}(x_{t},t)\\\
\approx\epsilon_{\theta}(x_{t},y,t)+w\sigma_{t}\nabla_{x}\mathcal{L}(\theta_{t},x_{t},y),$
(3)

where $\theta_{t}$ are the equivalent parameters for a time-conditional
classifier. Coined “classifier-free guidance” with weight $w$, sampling with
$\epsilon^{w}_{\theta}(x_{t},y,t)$ as the noise prediction function generates
high-quality, conditional data.

##  3 Methodology

#### CE Dataset Generation

Our approach to CE generation for datum $x$ is to recast it as sampling from a
sequence of un-normalized distributions defined by the product of the data
distribution (represented by the class-conditional diffusion model) with an
independently diffused CE “neighborhood” distribution of scale $\sigma_{CE}$
centered on $\mu_{CE}=x$. The score of this distribution at each time $t$ is
the sum of scores of the two components, so we may simply add the diffused
neighborhood score to the diffusion model score.

We present two variants of the noise prediction function corresponding to
different choices for the neighborhood distribution: Gaussian and Boltzmann-
inspired (see appendix).

$\epsilon^{G}_{\theta}(x_{t},y,\mu_{CE},t)\coloneq\epsilon^{w}_{\theta}(x_{t},y,t)-\frac{\alpha_{t}\mu_{CE}-x_{t}}{\alpha_{t}^{2}\sigma_{CE}^{2}+\sigma_{t}^{2}}$
(4) $\epsilon^{B}_{\theta}(x_{t},y,\mu_{CE},t)\coloneq\\\
\epsilon^{w}_{\theta}(x_{t},y,t)-\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}\text{hardtanh}(\gamma_{t}\,(x_{t}-\alpha_{t}\mu_{CE}))$
(5)

where $\alpha_{t}$ is a (diffusion) time-dependent scalar decreasing in
$t\in[0,1]$, $\sigma_{t}=\sqrt{1-\alpha_{t}^{2}}$, and $\gamma_{t}$ is a time-
dependent scalar derived from the first-order Maclaurin series of the
Boltzmann-inspired scores, and is defined as:

$\gamma_{t}\coloneq\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}-\frac{\sqrt{2}}{\sigma_{t}\sqrt{\pi}}\left(\exp{\left(\frac{\sigma_{t}^{2}}{\alpha_{t}^{2}\sigma^{2}_{CE}}\right)}\,\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}\right)\right)^{-1}.$
(6)

A CE $x_{CE}$ is generated for a datum $x$ by assigning $\mu_{CE}\leftarrow x$
and sampling with $\epsilon^{G}_{\theta}(\cdot)$ or
$\epsilon^{B}_{\theta}(\cdot)$ with guidance towards target class $y_{CE}$.
Guidance $w$ and neighborhood scale $\sigma_{CE}$ are hyperparameters. The
Boltzmann-inspired distribution has a sharper mode than the Gaussian
distribution, encouraging $x-x_{CE}$ to be lower norm ($L^{2}$) and more
sparse.

(a) CIFAR10 CE Dist. ($L^{2}$)

(b) CIFAR10 CE Dist. ($L^{0}$)

(c) SVHN CE Dist. ($L^{2}$)

(d) SVHN CE Dist. ($L^{0}$)

Figure 2: CE distribution comparison. Boltzmann variant CEs produce lower-
norm, sparser changes. Best viewed in color.

##  4 Experiments

(a) CIFAR10 $\varepsilon=0$
$r^{2}=0.02$

(b) CIFAR10 $\varepsilon=1$
$r^{2}=0.33$

(c) CIFAR10 $\varepsilon=1.5$
$r^{2}=0.39$

(d) SVHN $\varepsilon=0$
$r^{2}=0.01$

(e) SVHN $\varepsilon=0.2$
$r^{2}=0.01$

(f) SVHN $\varepsilon=0.4$
$r^{2}=0.24$

Figure 3: Scatter plots of classifier confidence and average CE distance of
10000 clean training samples as adversarial training norm is increased. Robust
models are more likely to misclassify and lose confidence on data which have
closer CEs. Best viewed in color.

(a) CIFAR10 Accuracy

(b) CIFAR10 Prob. Source Predicted

(c) SVHN Accuracy

(d) SVHN Prob. Source Predicted

Figure 4: Classifier performance on 10000 training data and 200000 CE data
generated from the training samples. Best viewed in color.

We run experiments in the PyTorch framework [14] with the CIFAR10 [11] and
SVHN [13] image classification benchmarks . We employ conditional score-based
models (SBMs) [18] for all diffusion models and $L^{2}$ PGD($8$) [12]
$2\times$WideResNet-40 models [26] for all robust classifiers.

#### CE Dataset Evaluation

In each experiment we generate 2 CEs for each class for at least 1000 samples
from the training set, resulting in at least 20000 CEs. We compare the two CE
generation variants with CEs generated by a robust classifier in Fig. 2. The
Boltzmann variant (type $B$) produces CEs with lower-norm ($L^{2}$) and
sparser ($L^{0}$) changes than the Gaussian variant does. The Boltzmann
variant produces CEs of norm less than or equal to the norm of CEs produced by
the robust classifier. Standard classifiers (no adversarial training) achieve
high accuracy on the CEs (i.e., they predict $y_{CE}$ when provided $x_{CE}$),
indicating good semantic quality. The remaining experiments use 200000
Boltzmann variant CEs generated from 10000 training samples (CIFAR10 or SVHN)
with $w=15$ and $\sigma_{CE}=0.2$. Please see the appendix for more
information.

#### Robust Model Evaluation

Our experience is that performance loss of robust models begins on the clean
training data, incurring performance loss on clean test data downstream. We
plot the average $L^{2}$ distance of clean training samples with the
confidence and accuracy of robust models on the clean samples in Fig. 3.
Standard models ($\varepsilon=0$) achieve very high clean training accuracy,
as expected. As robust training budget $\varepsilon$ increases, robust models
are more likely to misclassify clean training data which are closer to their
CEs. Moreover, the confidence of robust models on their training data becomes
correlated with the proximity of the data to their CEs ($r^{2}=0.39$ for
CIFAR10 and $r^{2}=0.24$ for SVHN).

Fig. 4(a) and Fig. 4(c) depict the accuracy of robust models on clean train
data and on CEs generated from the train data. Same-class CEs resemble the
original data; hence their robust accuracy trends are similar. However, robust
models perform very poorly on different-class CEs. Fig. 4(b) and Fig. 4(d)
depict the probability that the original label $y$ is predicted by a robust
model for $(x_{CE},y_{CE})$, given that the robust model was correct on
$(x,y)$. Robust models are much more likely to misclassify CEs as having the
source label, indicating that they become invariant to the low-norm, semantic
changes brought by CEs.

Fig. 5 depicts the $L^{2}$ distance distribution of Boltzmann-type CEs
generated from the original CIFAR10 training data and Boltzmann-type CEs
generated from robust CEs ($\varepsilon=1$, $\varepsilon=2$). On average,
Boltzmann-type CEs are farther away from their source data samples when the
data samples are robust model CEs (compared to using the original training
data as the source data). Since robust model CEs are generated by following
input gradient to maximize the confidence of a target class, this indicates
that robust model gradients orient towards data regions which are farther away
from CEs. Future work may investigate a link between the perceptually aligned
gradient of robust classifiers [20] and the proximity of data to CEs.

Figure 5: Distance of different-class CEs generated by the Boltzmann method
($w=15\ \sigma_{CE}=0.2$) when the input data is the original CIFAR10 train
data, CEs generated by a robust $\varepsilon=1$ model from the CIFAR10 train
data, and CEs generated by a robust $\varepsilon=2$ model from the CIFAR10
train data. Robust model CEs tend to be in data regions farther away from our
diffusion-generated CEs. Best viewed in color.

The method presented in this work can be leveraged to convert a conditional
diffusion model into an interpretable classifier. Using the diffusion CE
method with a “class with lowest average CE distance” decision rule achieves
$\sim 85\%$ accuracy on $1000$ input samples of the CIFAR10 test set.
Critically, input data for classification are never provided directly to the
diffusion model - the input data merely guides CE sampling with analytic
scores.

#### Limitations

Our study considers one combination of classifier architecture, diffusion
model, and adversarial attack. This study is limited to the $L^{2}$ norm
constraint and two image classification benchmarks. Investigation with
additional datasets, attacks, and models is left to future work.

##  5 Conclusion

We present a simple, classifier-free diffusion method to generate
counterfactual examples (CEs) which enables novel investigation of the
performance loss of robust classifiers. Our results indicate a significant
overlap between non-robust and semantically meaningful features, countering
the common assumption that non-robust features are not interpretable. Hence,
robust models must become invariant to this subset of semantic features along
with non-semantic adversarial perturbations. Our study motivates new
approaches to robust training which can resolve this issue.

## Acknowledgements

Part of this research was conducted while Eric Yeats was supported by a DOE-
SCGSR Award at Oak Ridge National Laboratory (ORNL). Hence, this material is
based upon work supported by the U.S. Department of Energy, Office of Science,
Office of Workforce Development for Teachers and Scientists, Office of Science
Graduate Student Research (SCGSR) program. The SCGSR program is administered
by the Oak Ridge Institute for Science and Education (ORISE) for the DOE.
ORISE is managed by ORAU under contract number DE-SC0014664. All opinions
expressed in this paper are the author’s and do not necessarily reflect the
policies and views of DOE, ORAU, or ORISE. This research was also supported in
part by the U.S. Department of Energy, through the Office of Advanced
Scientific Computing Research’s “Data-Driven Decision Control for Complex
Systems (DnC2S)” project.

This research is supported in part by the U.S. National Science Foundation
funding CNS-1822085 and U.S. Air Force Research Lab funding FA8750-21-1-1015.

This research used resources of the Experimental Computing Laboratory (ExCL)
at the Oak Ridge National Laboratory, which is supported by the Office of
Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

## References

  * Augustin et al. [2022] Maximilian Augustin, Valentyn Boreiko, Francesco Croce, and Matthias Hein.  Diffusion visual counterfactual explanations.  _Advances in Neural Information Processing Systems_ , 35:364–377, 2022.
  * Boreiko et al. [2022] Valentyn Boreiko, Maximilian Augustin, Francesco Croce, Philipp Berens, and Matthias Hein.  Sparse visual counterfactual explanations in image space.  In _DAGM German Conf. on Pattern Recognition_ , pages 133–148. Springer, 2022.
  * Cohen et al. [2019] Jeremy Cohen, Elan Rosenfeld, and Zico Kolter.  Certified adversarial robustness via randomized smoothing.  In _ICML_ , pages 1310–1320. PMLR, 2019.
  * Finlay and Oberman [2021] Chris Finlay and Adam M Oberman.  Scaleable input gradient regularization for adversarial robustness.  _Machine Learning with Applications_ , 3:100017, 2021.
  * Freiesleben [2022] Timo Freiesleben.  The intriguing relation between counterfactual explanations and adversarial examples.  _Minds and Machines_ , 32(1):77–109, 2022.
  * Goodfellow et al. [2014] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  _arXiv preprint arXiv:1412.6572_ , 2014.
  * Hendrycks and Gimpel [2016] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus).  _arXiv preprint arXiv:1606.08415_ , 2016.
  * Ho and Salimans [2022] Jonathan Ho and Tim Salimans.  Classifier-free diffusion guidance.  _arXiv preprint arXiv:2207.12598_ , 2022.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _Advances in neural information processing systems_ , 33:6840–6851, 2020.
  * Ilyas et al. [2019] Andrew Ilyas, Shibani Santurkar, Dimitris Tsipras, Logan Engstrom, Brandon Tran, and Aleksander Madry.  Adversarial examples are not bugs, they are features.  _Advances in neural information processing systems_ , 32, 2019.
  * Krizhevsky et al. [2009] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * Madry et al. [2017] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu.  Towards deep learning models resistant to adversarial attacks.  _arXiv preprint arXiv:1706.06083_ , 2017.
  * Netzer et al. [2011] Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, and Andrew Y Ng.  Reading digits in natural images with unsupervised feature learning.  2011\.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Pawelczyk et al. [2022] Martin Pawelczyk, Chirag Agarwal, Shalmali Joshi, Sohini Upadhyay, and Himabindu Lakkaraju.  Exploring counterfactual explanations through the lens of adversarial examples: A theoretical and empirical analysis.  In _International Conference on Artificial Intelligence and Statistics_ , pages 4574–4594. PMLR, 2022.
  * Ross and Doshi-Velez [2018] Andrew Ross and Finale Doshi-Velez.  Improving the adversarial robustness and interpretability of deep neural networks by regularizing their input gradients.  In _Proceedings of the AAAI conference on artificial intelligence_ , 2018.
  * Rudin [2019] Cynthia Rudin.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nature machine intell._ , 1(5):206–215, 2019.
  * Song et al. [2020] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2020.
  * Szegedy et al. [2013] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  _arXiv preprint arXiv:1312.6199_ , 2013.
  * Tsipras et al. [2018] Dimitris Tsipras, Shibani Santurkar, Logan Engstrom, Alexander Turner, and Aleksander Madry.  Robustness may be at odds with accuracy.  _arXiv:1805.12152_ , 2018.
  * Verma et al. [2020] Sahil Verma, John Dickerson, and Keegan Hines.  Counterfactual explanations for machine learning: A review.  _arXiv preprint arXiv:2010.10596_ , 2, 2020.
  * Wu and He [2018] Yuxin Wu and Kaiming He.  Group normalization.  In _Proceedings of the European conference on computer vision (ECCV)_ , pages 3–19, 2018.
  * Yeats et al. [2022] Eric Yeats, Frank Liu, David Womble, and Hai Li.  Nashae: Disentangling representations through adversarial covariance minimization.  In _European Conference on Computer Vision_ , pages 36–51. Springer, 2022.
  * Yeats et al. [2023] Eric Yeats, Cameron Darwin, Frank Liu, and Hai Li.  Adversarial estimation of topological dimension with harmonic score maps.  _arXiv preprint arXiv:2312.06869_ , 2023.
  * Yeats et al. [2021] Eric C Yeats, Yiran Chen, and Hai Li.  Improving gradient regularization using complex-valued neural networks.  In _International Conference on Machine Learning_ , pages 11953–11963. PMLR, 2021.
  * Zagoruyko and Komodakis [2016] Sergey Zagoruyko and Nikos Komodakis.  Wide residual networks.  _arXiv preprint arXiv:1605.07146_ , 2016.

Supplementary Material

##  6 Reproducibility Information

$L^{2}$ distance for the CEs is defined as $\|x-x_{CE}\|_{2}$, where $x$ is
the original sample from which the CE $x_{CE}$ was generated. $L^{0}$
distances ($\|x-x_{CE}\|_{0}$) employ an element-wise threshold of $0.02$ for
each pixel difference value, and the reported value is normalized (the
thresholded $L^{0}$ distance is divided by the dimension of $x$). CE quality
is measured by the accuracy of a standard classifier in predicting $y_{CE}$
given $x_{CE}$ (depicted in Fig. 6). On SVHN there appears to be sizeable drop
from robust CE accuracy to the Boltzmann CEs. This may be due to several
reasons. First, the $\epsilon=0.4$ robust classifier may be making adversarial
changes (rather than semantic), biasing its quality measure higher at that
$L^{2}$ distance. Second, diffusion models tend to ignore class conditioning
information when their inputs are simple (like those of SVHN), and they
operate instead as unconditional denoising functions. This would cause the CE
generation method to fail. Visual inspection of the CEs for SVHN indicates
that the $\varepsilon=0.4$ classifier is making some adversarial (non-
semantic) changes, and the diffusion method fails in some cases. Please refer
to the end of the supplementary material for visual depictions of CEs.

(a) CIFAR10 Standard Accuracy

(b) SVHN Standard Accuracy

Figure 6: Comparison of the accuracy attained by standard models on the CEs.
Best viewed in color.

### Robust WideResNet Experiments

2$\times$WideResNet-40 models [26] are adversarially trained with PGD(8) [12]
in $L^{2}$ on CIFAR10 [11] and SVHN [13] for 100 epochs and learning rate
$1e-3$ using the Adam optimizer $\beta=(0.9,0.999)$. For standard models
($\varepsilon=0$), dropout with rate 0.3 is used. The data augmentations for
CIFAR10 are random horizontal flips and random crops with padding 4. No
augmentations are used for SVHN. The clean accuracies of the trained models
are listed in tables (1) and (2). All adversarial attacks are PGD(8).

We observed that robust training on SVHN would collapse for $\varepsilon\geq
0.5$, so only $\varepsilon\leq 0.4$ experiments are reported. Coincidentally,
$\varepsilon\sim 0.5$ is half the distance of a large amount of different-
class CEs to their original data on SVHN. We hypothesize that the presence of
many CEs at the $L^{2}$ distance $1$ may be related to this training collapse
phenomenon, as an adversarial budget of $\varepsilon\sim 0.5$ could make the
source class of perturbed training data highly ambiguous.

Table 1: Accuracy of CIFAR10 Models on the CIFAR10 test set
$\varepsilon$ | Clean Accuracy (%) |  | PGD(8, 0.5) Accuracy (%)
---|---|---|---
0 | 93.08 |  | 5.75
0.1 | 91.33 |  | 45.06
0.5 | 84.87 |  | 54.73
1.0 | 76.23 |  | 54.53
1.5 | 66.21 |  | 52.34
2.0 | 58.37 |  | 47.53
Table 2: Accuracy of SVHN Models on the SVHN test set
$\varepsilon$ | Clean Accuracy (%) |  | PGD(8, 0.2) Accuracy (%)
---|---|---|---
0 | 96.75 |  | 77.46
0.1 | 95.63 |  | 86.42
0.2 | 95.41 |  | 86.84
0.4 | 93.09 |  | 85.43

### Diffusion Models

In all experiments using diffusion models, we employ the variance-preserving
(VP) score-based models (SBM) of Song et al. [18]. The SBM architecture
follows a U-net structure with four symmetric stages of two ResNet blocks in
each encoding or decoding stage. Downsampling (and upsampling, respectively)
occurs in the innermost three stages (i.e., stages 2, 3, 4). 128 channels
(features) are used, and the number of features used is doubled to 256 for
stages 2, 3, and 4. Attention is applied at the center of the U-net and after
the first downsampling stage and before the last upsampling stage. The SBMs
use the SiLU activation function [7] and GroupNorm [22] with a group size of
32. Training on CIFAR10 or SVHN occurs for 1 million iterations of batch size
128 with a learning rate of $2e-4$ and Adam optimizer $\beta=(0.9,0.999)$. We
use a learning rate warmup for 5000 iterations and gradient clipping with norm
1. Class conditioning is provided as a learnable embedding which is added to
the time condition embedding. An additional learnable null embedding signifies
that no class information is provided. During conditional SBM training, class
conditions are dropped and replaced with this null embedding at a rate of
$30\%$.

##  7 CE Generation Method

CEs were generated by sampling from a sequence of un-normalized distributions
represented by the product of the data distribution (rep. by diffusion model)
and an independently diffused neighborhood distribution. The diffusion model
used is a variance-preserving (VP) score-based model of Song et al. [18] and
the sampling strategy for generation used an Euler-Maruyama predictor with
1000 discrete steps.

Since the density at each step in the sequence is represented as the product
of the diffused data distribution and an independently diffused neighborhood
distribution, sampling amounts to adding the conditional diffusion score (and
classifier free guidance) with the analytic neighborhood diffusion score.
Sampling then proceeds as normal (see [18]) with the augmented score.

#### Robust Model CEs

For CE generation with robust models, we push data in the direction of
targeted gradients with a step size of $0.05$ until a targeted class
confidence of $0.9$ or a maximum of 200 steps is reached. At each step, we
clip the pixels of the image to remain in [0, 1].

### Boltzmann-Inspired Distribution

The density of the 1D Boltzmann-inspired distribution proposed in this work is
given by:

$b(x)=\frac{1}{\sqrt{2}\,\sigma_{CE}}\exp{\left(-\frac{\sqrt{2}\,|x-\mu_{CE}|}{\sigma_{CE}}\right)},$
(7)

where $\mu_{CE}$ is its mean and $\sigma_{CE}$ is its standard deviation.
Scaling samples of this distribution with $\alpha_{t}>0$ amounts to sampling
from a new distribution with mean $\alpha_{t}\mu_{CE}$ and standard deviation
$\alpha_{t}\sigma_{CE}$. Defining $y_{t}=x-\alpha_{t}\mu_{CE}$, we convolve
this distribution with a Gaussian distribution of mean $\mu=0$ and scale
$\sigma_{t}$ to yield the (diffusion) time-dependent distribution:

$b_{t}(x)=\\\
\frac{\exp\left(\frac{\sigma_{t}^{2}}{\alpha_{t}^{2}\sigma_{CE}^{2}}-\frac{y_{t}}{\sqrt{2}\alpha_{t}\sigma_{CE}}\right)\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}-\frac{y}{\sqrt{2}\sigma_{t}}\right)}{2\sqrt{2}\alpha_{t}\sigma_{CE}}\\\
+\frac{\exp\left(\frac{\sigma_{t}^{2}}{\alpha_{t}^{2}\sigma_{CE}^{2}}+\frac{y_{t}}{\sqrt{2}\alpha_{t}\sigma_{CE}}\right)\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}+\frac{y}{\sqrt{2}\sigma_{t}}\right)}{2\sqrt{2}\alpha_{t}\sigma_{CE}}.$
(8)

Taking the logarithm of $b_{t}(x)$ and differentiating with respect to $x$, we
have:

$\nabla_{x}\log b_{t}(x)=\\\
\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}\Bigg{(}e^{\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}y_{t}}\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}+\frac{y_{t}}{\sqrt{2}\sigma_{t}}\right)\\\
-e^{-\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}y_{t}}\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}-\frac{y_{t}}{\sqrt{2}\sigma_{t}}\right)\Bigg{)}\Bigg{/}\\\
\Bigg{(}e^{\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}y_{t}}\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}+\frac{y_{t}}{\sqrt{2}\sigma_{t}}\right)\\\
+e^{-\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}y_{t}}\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}-\frac{y_{t}}{\sqrt{2}\sigma_{t}}\right)\Bigg{)}.$
(9)

Although this expression yields the exact scores for the diffused Boltzmann-
inspired distribution in 1D, it is numerically unstable. We note that
$\nabla_{x}\log b_{t}(x)$ is sigmoidal, and we elect to approximate it using
the hardtanh function. One may consider the
$\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}$ expression to define the range of the
score values, and the remainder of the expression defines a Gaussian-like
score near the mean. Hence, we define our approximate scores as:

$\nabla_{x}\log
b_{t}(x)\approx\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}\text{hardtanh}(\gamma_{t}\,y_{t}),$
(10)

where $\gamma_{t}$ is the first term of a Maclaurin series estimate of
$\nabla_{x}\log b_{t}(x)$ and is given by:

$\gamma_{t}=\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}-\frac{\sqrt{2}}{\sigma_{t}\sqrt{\pi}}\left(\exp{\left(\frac{\sigma_{t}^{2}}{\alpha_{t}^{2}\sigma^{2}_{CE}}\right)}\,\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}\right)\right)^{-1}.$
(11)

With $u=\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}$, the expression
$\frac{\exp(-u^{2})}{\text{erfc}(u)}$ is numerically unstable for large values
of $u$. However, beyond a certain point (e.g., $u\geq 20$), the function
behaves as a linear function with slope $\sqrt{\pi}$. We avoid the numerical
instability by switching to a linear approximation at $u=20$.

Figure 7: Comparison of true 1D Boltzmann-inspired scores with the proposed
hardtanh approximation (black). The mean $\mu_{CE}$ is selected to be 0. Best
viewed in color.

The exact score function $\nabla_{x}\log b_{t}(x)$ and its approximation is
displayed with $\alpha_{t}=1$ for various values of $\sigma_{t}$ and in figure
7. In our experiments, this 1D score function is applied element-wise to
vector inputs, providing a sparsifying effect on CE generation.

For completeness, we include $\nabla^{2}_{x}\log b_{t}(x)$, which was
evaluated at $y_{t}=0$ to yield the first term of the Maclaurin series of
$\nabla_{x}\log b_{t}(x)$:

$\nabla^{2}_{x}\log b_{t}(x)=\\\
\Bigg{[}\frac{2}{\sigma^{2}_{CE}}\Bigg{(}e^{\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}y_{t}}\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}+\frac{y_{t}}{\sqrt{2}\sigma_{t}}\right)\\\
+e^{-\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}y_{t}}\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}-\frac{y_{t}}{\sqrt{2}\sigma_{t}}\right)\Bigg{)}\\\
-\frac{4\,e^{-\left(\frac{\sigma_{t}^{2}}{\alpha_{t}^{2}\sigma_{CE}^{2}}+\frac{y_{t}^{2}}{2\sigma_{t}^{2}}\right)}}{\sqrt{\pi}\alpha_{t}\sigma_{CE}\sigma_{t}}\Bigg{]}\quad\Bigg{/}\\\
\Bigg{[}e^{\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}y_{t}}\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}+\frac{y_{t}}{\sqrt{2}\sigma_{t}}\right)\\\
+e^{-\frac{\sqrt{2}}{\alpha_{t}\sigma_{CE}}y_{t}}\text{erfc}\left(\frac{\sigma_{t}}{\alpha_{t}\sigma_{CE}}-\frac{y_{t}}{\sqrt{2}\sigma_{t}}\right)\Bigg{]}.$
(12)

Figure 8: Boltzmann ($w=15$ $\sigma_{CE}=0.2$) CEs from the CIFAR10 training
set. Leftmost column: original samples and class label. Other columns: CE and
$L^{2}$ distance to the original sample. Figure 9: Boltzmann ($w=15$
$\sigma_{CE}=0.2$) CEs from the CIFAR10 training set. Leftmost column:
original samples and class label. Other columns: CE and $L^{2}$ distance to
the original sample. Figure 10: Robust model ($\varepsilon=1$, confidence
threshold 0.9) CEs from the CIFAR10 training set. Leftmost column: original
samples and class label. Other columns: CE and $L^{2}$ distance to the
original sample. Figure 11: Robust model ($\varepsilon=1$, confidence
threshold 0.9) CEs from the CIFAR10 training set. Leftmost column: original
samples and class label. Other columns: CE and $L^{2}$ distance to the
original sample. Figure 12: Boltzmann ($w=15$ $\sigma_{CE}=0.2$) CEs from the
SVHN training set. Leftmost column: original samples and class label. Other
columns: CE and $L^{2}$ distance to the original sample. Figure 13: Boltzmann
($w=15$ $\sigma_{CE}=0.2$) CEs from the SVHN training set. Leftmost column:
original samples and class label. Other columns: CE and $L^{2}$ distance to
the original sample. Figure 14: Robust model ($\varepsilon=0.4$, confidence
threshold 0.9) CEs from the SVHN training set. Leftmost column: original
samples and class label. Other columns: CE and $L^{2}$ distance to the
original sample. Figure 15: Robust model ($\varepsilon=0.4$, confidence
threshold 0.9) CEs from the SVHN training set. Leftmost column: original
samples and class label. Other columns: CE and $L^{2}$ distance to the
original sample.