Automation of Beam and Jet functions at NNLO

Guido Bell, Kevin Brune, Goutam Das$\star$, and Marcel Wald

Theoretische Physik 1, Center for Particle Physics, Universität Siegen,

Walter-Flex-Strasse 3, 57068 Siegen, Germany

$\star$<EMAIL_ADDRESS>
$\dagger$ SI-HEP-2021-25, P3H-21-069

15th International Symposium on Radiative Corrections:
Applications of Quantum Field Theory to Phenomenology,
FSU, Tallahasse, FL, USA, 17-21 May 2021

## Abstract

We present a novel framework to streamline the calculation of jet and beam
functions to next-to-next-to-leading order (NNLO) in perturbation theory. By
exploiting the infrared behaviour of the collinear splitting functions, we
factorise the singularities with suitable phase-space parametrisations and
perform the observable-dependent integrations numerically. We have implemented
our approach in the publicly available code pySecDec and present first results
for sample jet and beam functions.

###### Contents

  1. 1 Introduction
  2. 2 Jet functions
    1. 2.1 NLO calculation
    2. 2.2 NNLO calculation
    3. 2.3 Renormalisation
    4. 2.4 Results
  3. 3 Beam functions
    1. 3.1 Calculational details
    2. 3.2 Results
  4. 4 Conclusions

##  1 Introduction

In recent years, Soft-Collinear Effective Theory (SCET) has been successful in
describing observables at lepton as well as hadron colliders through the
resummation of large logarithms appearing in different corners of phase space.
The resummation in the SCET framework relies on an underlying factorisation
theorem consisting of hard ($H$), beam ($B$), jet ($J$), and soft ($S$)
functions,

$\displaystyle d\sigma=H\cdot\prod_{i}B_{i}\otimes\prod_{j}J_{j}\otimes S\,.$
(1)

All these functions can be calculated in perturbation theory order-by-order in
the strong-coupling expansion. The resummation can be performed by evolving
them to a common renormalisation scale by solving the corresponding
renormalisation group equations (RGE). The calculation of these functions
often becomes challenging at higher orders. In particular, the jet, beam, and
soft quantities depend on the specific observable and need to be calculated on
a case-by-case basis. In recent years, there have been efforts in automating
the calculation of these perturbative ingredients. Whereas di-jet soft
functions are now available to next-to-next-to-leading order (NNLO) for many
SCET-1 and SCET-2 observables through the publicly available package SoftSERVE
[1, 2, 3], an extension to general N-jet soft functions is currently in
progress [4]. There have also been efforts in automating the calculation of
jet and beam functions at NLO [5, 6], and in this work we plan to bring these
efforts to the NNLO level, which is essential to achieve high-precision
resummations at collider processes.

##  2 Jet functions

The jet functions appear in SCET factorisation theorems whenever one considers
processes with coloured partons in the final state at both lepton or hadron
colliders. In this work, we focus on quark jet functions, which are defined in
terms of the collinear field operator
$\chi=W^{\dagger}_{\bar{n}}\frac{{\not{n}}{\not{\bar{n}}}}{4}\psi$ via

$\displaystyle\Big{[}\frac{{\not{n}}}{2}\Big{]}_{{}_{\beta\alpha}}J_{q}(\tau,\mu)=\frac{1}{\pi}\sum_{X}(2\pi)^{d}\delta\Big{(}Q-\sum_{i}k_{i}^{-}\Big{)}\delta^{(d-2)}\Big{(}\sum_{i}k_{i}^{\perp}\Big{)}\bra{0}\chi_{\beta}\ket{X}\bra{X}\bar{\chi}_{\alpha}\ket{0}{\cal
M}(\tau;\\{k_{i}\\})\,,$ (2)

where $W_{\bar{n}}$ denotes a collinear Wilson line, and we introduced light-
cone coordinates with $k_{i}^{-}=\bar{n}\cdot k_{i}$, $k_{i}^{+}=n\cdot k_{i}$
and a transverse component $k_{i}^{\perp,\mu}$ satisfying $n\cdot
k_{i}^{\perp}=\bar{n}\cdot k_{i}^{\perp}=0$, along with $n^{2}=\bar{n}^{2}=0$
and $n\cdot\bar{n}=2$. The sum over $X$ refers to the phase space of the
final-state particles and ${\cal M}(\tau;\\{k_{i}\\})$ denotes a generic
measurement function in Laplace space, with $\tau$ being the corresponding
Laplace variable and $\\{k_{i}\\}$ the momenta of the final-state particles.
The expansion of the bare jet functions in the strong coupling
($a_{s}=\alpha_{s}/4\pi$) in Laplace space can then be written as

$\displaystyle J_{q}^{0}(\tau)$
$\displaystyle=1+\sum_{k=1}^{\infty}\left(Z_{\alpha}a_{s}\right)^{k}\left(\mu^{2}\bar{\tau}^{2}\right)^{k\epsilon}J_{q}^{(k)}(\tau,\epsilon)\,,$
(3)

where $\bar{\tau}=\tau e^{\gamma_{E}}$ with $\gamma_{E}\simeq 0.5772$ being
Euler’s constant, and $Z_{\alpha}=1-a_{s}\beta_{0}/\epsilon$ is the strong-
coupling renormalisation constant in $d=4-2\epsilon$ dimensions in the
$\overline{\text{MS}}$-scheme.

###  2.1 NLO calculation

At NLO one encounters a two-body phase space for the calculation of the jet
functions in terms of an on-shell gluon ($k^{\mu}$) and a quark ($p^{\mu}$)
momentum. However, according to (2) the sum of their large components is
constrained by the total jet energy $Q$ and their transverse momenta must
balance each other. This leads to the following parametrisation (defining
$\bar{z}=1-z$),

$\displaystyle k_{-}=zQ,\qquad
p_{-}=\bar{z}Q,\qquad|\vec{k}_{\perp}|=|\vec{p}_{\perp}|=k_{T},\qquad\cos\theta_{k}=1-2t_{k}\,.$
(4)

At NLO the collinear matrix element is proportional to the well-known
splitting function [7],

$\displaystyle P_{q^{*}\to gq}^{(0)}(z)$
$\displaystyle=\frac{C_{F}}{z}\left[1+\bar{z}^{2}-\epsilon z^{2}\right]\,.$
(5)

We then parametrise the one-emission measurement function in the form (see
also [1, 2, 3]),

$\displaystyle{\cal M}_{1}(\tau;p,k)$ $\displaystyle=\exp\left[-\tau
k_{T}\left(\frac{k_{T}}{zQ}\right)^{n}f(z,t_{k})\right]\,,$ (6)

where the exponential is a result of the Laplace transformation of the
momentum-space measurement function, and the function $f(z,t_{k})$ is assumed
to be finite and non-zero in the singular limit of the matrix element $z\to
0$. With this ansatz for the generic measurement function, we arrive at the
following master formula for the NLO quark jet functions,

$\displaystyle J_{q}^{(1)}(\tau,\epsilon)=$ $\displaystyle\left(\tau
Q\right)^{\frac{-2n\epsilon}{1+n}}\frac{8e^{-\gamma_{E}\epsilon}}{(1+n)\sqrt{\pi}}\,\frac{\Gamma\left(-\frac{2\epsilon}{1+n}\right)}{\Gamma\left(\frac{1}{2}-\epsilon\right)}$
$\displaystyle\times\;\int_{0}^{1}\\!dz\;\,z^{-1-\frac{2n\epsilon}{1+n}}\;\left[zP_{q^{*}\to
gq}^{(0)}(z)\right]\;\int_{0}^{1}\\!d\\!t_{k}\;\,(4t_{k}\bar{t}_{k})^{-\frac{1}{2}-\epsilon}\;f(z,t_{k})^{\frac{2\epsilon}{1+n}}\,.$
(7)

Notice that the $k_{T}$-integration is already performed, and we are thus left
with two remaining integrations over the splitting variable $z$ and the
angular variable $t_{k}$, which we perform numerically. As seen above, all
singularities are nicely factorised in this representation in terms of the
gamma function and the monomial $z^{-1-\frac{2n\epsilon}{1+n}}$. The explicit
dependence on the observable then enters through the parameter $n$ defined in
(6) and the function $f(z,t_{k})$, which parametrises the dependence on the
splitting variable and the azimuthal angle $\theta_{k}$.

###  2.2 NNLO calculation

The NNLO contribution to the jet functions involves two different kinds of
contributions, $viz.$ the real-virtual (RV) and the real-real (RR)
contribution. In the RV case, the calculation follows similarly to the NLO
case since the phase space still consists of two particles. In particular, the
measurement function is again given by (6) and the matrix element is now
related to the one-loop correction to the splitting function [8, 9, 10],

$\displaystyle\mathbb{P}_{q^{*}\to gq}^{(1)}(z)=$
$\displaystyle\frac{C_{F}}{z^{1-\epsilon}\bar{z}^{-\epsilon}}\bigg{\\{}(1+\bar{z}^{2}-\epsilon
z^{2})\bigg{[}C_{F}+(C_{F}-C_{A})\,\Big{(}1-\frac{\epsilon^{2}}{1-2\epsilon}\Big{)}-C_{A}\;{}_{2}F_{1}\Big{(}1,-\epsilon;1-\epsilon;-\frac{\bar{z}}{z}\Big{)}$
$\displaystyle\hskip
19.91692pt+(C_{A}-2C_{F})\,{}_{2}F_{1}\Big{(}1,-\epsilon;1-\epsilon;-\frac{z}{\bar{z}}\Big{)}\bigg{]}+(C_{F}-C_{A})\,\frac{\epsilon^{2}\bar{z}(1+\bar{z})}{1-2\epsilon}\bigg{\\}}\,.$
(8)

In terms of this, the master formula for the RV contribution takes the form,

$\displaystyle J^{(2)}_{q,RV}(\tau,\epsilon)=$ $\displaystyle\left(\tau
Q\right)^{\frac{-4n\epsilon}{1+n}}\frac{4^{2+\epsilon}\pi~{}e^{-2\gamma_{E}\epsilon}}{1+n}\,\frac{\Gamma(-\frac{4\epsilon}{1+n})\cot(\pi\epsilon)}{\epsilon~{}\Gamma(1/2-\epsilon)^{2}}$
$\displaystyle\times\;\int_{0}^{1}\\!dz\;\,z^{-1-\frac{4n\epsilon}{1+n}}\;\left[z~{}\mathbb{P}_{q^{*}\to
gq}^{(1)}(z)\right]\;\int_{0}^{1}\\!d\\!t_{k}\;\,(4t_{k}\bar{t}_{k})^{-\frac{1}{2}-\epsilon}\;f(z,t_{k})^{\frac{4\epsilon}{1+n}}\,.$
(9)

The calculation of the RR contribution is, on the other hand, more involved
due to the complicated singularity structure of the $1\to 3$ splitting
functions [11, 12, 13]. In particular, it turns out one cannot factorise all
the overlapping singularities of the matrix elements with a single
parametrisation in this case. In order to tackle the RR contribution, we then
start from the following parametrisation of the three-body phase space,

$\displaystyle a=\frac{k_{-}l_{T}}{l_{-}k_{T}},\qquad
b=\frac{k_{T}}{l_{T}},\qquad z=\frac{k_{-}+l_{-}}{Q},\qquad
q_{T}=\sqrt{(k_{-}+l_{-})(k_{+}+l_{+})}\,,$ (10)

where $k^{\mu}$ and $l^{\mu}$ are the momenta of the emitted partons, while we
again denote the final-state momentum of the mother quark by $p^{\mu}$. Here,
$a$ is a measure of the rapidity difference of the emitted daughter particles,
$b$ is the ratio of their transverse momenta, and $z$ and $q_{T}$ parametrise
the dependence on their total light-cone momenta. In addition, we also need
three angles $\theta_{k},\theta_{l},\theta_{kl}$, which we rewrite in terms of
$t_{k},t_{l},t_{kl}$, similar to the NLO parametrisation discussed above. We
also find it convenient to remap the variables $a$ and $b$ to the unit
hypercube, which automatically introduces four sectors for each contribution.

To factorise the divergences of the RR contribution, we employ a mixed
strategy of sector-decomposition steps, non-linear transformations and
selector functions (details will be given in [14]). We then introduce a
similar generic measurement function for the two-emission case,

$\displaystyle{\cal M}_{2}(\tau;p,k,l)$ $\displaystyle=\exp\left[-\tau
q_{T}\left(\frac{q_{T}}{zQ}\right)^{n}{\cal
F}(a,b,z,t_{k},t_{l},t_{kl})\right]\,,$ (11)

whose exact form depends on the specific parametrisation and the divergences
of the matrix element. While the RR contribution involves three colour
structures ($C_{F}^{2},C_{F}C_{A},C_{F}T_{F}n_{f}$), the calculation
simplifies significantly for the $C_{F}T_{F}n_{f}$ piece, where the above
parametrisation can be used to factorise all divergences. Specifically, the
singularities arise in the limits $q_{T}\to 0$, $z\to 0$ and in the overlap of
$a\to 1$ and $t_{kl}\to 0$ in this case. The overlapping divergence is a
result of the configuration where the partons with momenta $k^{\mu}$ and
$l^{\mu}$ become collinear to each other. It can be resolved with a simple
substitution $(a,t_{kl})\to(u,v)$ as shown in [1, 2, 3].

The calculation of the $C_{F}^{2}$ and $C_{F}C_{A}$ colour structures, on the
other hand, involve more complicated singularity patterns. Moreover, one needs
to ensure that the function ${\cal F}$ defined in (11) stays finite and non-
zero in the singular limits of the matrix elements, which poses additional
constraints on the phase-space parametrisations. We have developed a strategy
that satisfies these constraints – while properly factorising all
singularities – which requires, however, to introduce about a dozen different
parametrisations for the $C_{F}^{2}$ and $C_{F}C_{A}$ contributions.

After factorising all divergences into monomials, one needs to perform a
Laurent expansion to expose the divergences in the dimensional regulator
$\epsilon$, followed by numerical integrations of the associated coefficients.
To perform these steps, we use the publicly available package pySecDec [15].
The numerical integrations are then performed with the Vegas routine as
implemented in the Cuba library inside the pySecDec framework.

###  2.3 Renormalisation

In Laplace space the renormalisation of the jet functions takes a
multiplicative form, $J_{q}=Z_{J_{q}}J_{q}^{0}$, and the corresponding
renormalisation group equation reads,

$\displaystyle\frac{d}{d\ln\mu}\;J_{q}(\tau,\mu)$
$\displaystyle=\bigg{[}2g(n)\,\Gamma_{\rm cusp}(\alpha_{s})\,L+\gamma^{\rm
J}(\alpha_{s})\bigg{]}\;J_{q}(\tau,\mu)\,,$ (12)

where $g(n)=(n+1)/n$,
$L=\ln\big(\mu\bar{\tau}/\left(Q\bar{\tau}\right)^{\frac{n}{1+n}}\big{missing})$,
and $\Gamma_{\rm cusp}(\alpha_{s})$ and $\gamma^{\rm J}(\alpha_{s})$ are the
cusp and non-cusp anomalous dimensions, respectively. Expanding the anomalous
dimensions in the form $G(\alpha_{s})=\sum_{n=0}^{\infty}G_{n}a_{s}^{n+1}$,
the solution of the RGE becomes up to two-loop order,

$\displaystyle J_{q}(\tau,\mu)=1\,+\,$ $\displaystyle
a_{s}(\mu)\,\bigg{\\{}g(n)\,\Gamma_{0}L^{2}+\gamma_{0}^{J}L+c_{1}^{J}\bigg{\\}}+a_{s}^{2}(\mu)\,\bigg{\\{}g(n)^{2}\frac{\Gamma_{0}^{2}}{2}L^{4}+g(n)\,\Big{(}\gamma_{0}^{J}+\frac{2\beta_{0}}{3}\Big{)}\,\Gamma_{0}L^{3}$
$\displaystyle+\left(g(n)\,(\Gamma_{1}+\Gamma_{0}c_{1}^{J})+\gamma_{0}^{J}\,\Big{(}\frac{\gamma_{0}^{J}}{2}+\beta_{0}\Big{)}\right)L^{2}+\bigg{(}\gamma_{1}^{J}+c_{1}^{J}(\gamma_{0}^{J}+2\beta_{0})\bigg{)}\,L+c_{2}^{J}\bigg{\\}}\,.$
(13)

The renormalisation constant $Z_{J_{q}}$ satisfies a similar RGE as (12) and
the solution follows as

$\displaystyle
Z_{J_{q}}(\tau,\mu)=1\,+\,a_{s}(\mu)\,\left\\{-g(n)\frac{\Gamma_{0}}{2}\frac{1}{\epsilon^{2}}-\bigg{(}g(n)\Gamma_{0}L+\frac{\gamma_{0}^{J}}{2}\bigg{)}\,\frac{1}{\epsilon}\right\\}+a_{s}^{2}(\mu)\,\bigg{\\{}g(n)^{2}\frac{\Gamma_{0}^{2}}{8}\frac{1}{\epsilon^{4}}$
$\displaystyle\quad+\left(g(n)^{2}\frac{\Gamma_{0}^{2}}{2}L+g(n)\Gamma_{0}\Big{(}\frac{\gamma_{0}^{J}}{4}+\frac{3\beta_{0}}{8}\Big{)}\right)\frac{1}{\epsilon^{3}}+\bigg{(}g(n)^{2}\frac{\Gamma_{0}^{2}}{2}L^{2}+g(n)\frac{\Gamma_{0}}{2}\left(\gamma_{0}^{J}+\beta_{0}\right)L$
$\displaystyle\qquad\;-\frac{g(n)\Gamma_{1}}{8}+\frac{(\gamma_{0}^{J})^{2}}{8}+\frac{\beta_{0}\gamma_{0}^{J}}{4}\bigg{)}\frac{1}{\epsilon^{2}}-\left(g(n)\frac{\Gamma_{1}}{2}L+\frac{\gamma_{1}^{J}}{4}\right)\frac{1}{\epsilon}\bigg{\\}}\,.$
(14)

This form of $Z_{J_{q}}$ provides a strong check of our calculation for the
higher poles, since the observable-independent anomalous dimensions
$\Gamma_{0}$, $\Gamma_{1}$ and $\beta_{0}$ are all known (we use the
conventions from [1, 2, 3]). On the other hand, we can use (13) and (14) to
extract the non-cusp anomalous dimensions $\gamma_{0}^{J}$ and
$\gamma_{1}^{J}$ from the $1/\epsilon$ poles, and the non-logarithmic
coefficients $c_{1}^{J}$ and $c_{2}^{J}$ of the renormalised jet functions
from the finite terms of the bare NNLO calculation.

###  2.4 Results

With this setup, we computed the jet functions for the event-shape observables
thrust and angularities. As the NLO case is trivial, we focus here on the NNLO
numbers, which we present in the form
$\gamma_{1}^{J}=\gamma_{1}^{C_{F}}\,C_{F}^{2}+\gamma_{1}^{C_{A}}\,C_{F}C_{A}+\gamma_{1}^{n_{f}}\,C_{F}T_{F}n_{f}$,
and similarly for $c_{2}^{J}$. Our numbers for thrust are summarised in Table
1, which show a very good agreement with the analytic results from [16]. Our
uncertainty estimates, on the other hand, seem to be too conservative at
present and require further investigations. We also find excellent agreement
for angularities for the $\gamma_{1}^{J}$ coefficients with the literature
[1], as can be seen from the left panel of Figure 1. On the right panel, we
present our numbers for the $c_{2}^{J}$ coefficients for different angularity
values, which significantly improve the Event2 fits results from [17].

|  Analytic |  This work
---|---|---
$\gamma_{1}^{n_{f}}$ | $-26.699$ | $-26.699(8)$
$\gamma_{1}^{C_{A}}$ | $-6.520$ | $-6.522(130)$
$\gamma_{1}^{C_{F}}$ | $21.220$ | $21.219(119)$

|  Analytic |  This work
---|---|---
$c_{2}^{n_{f}}$ | $-10.787$ | $-10.787(15)$
$c_{2}^{C_{A}}$ | $-2.165$ | $-2.169(189)$
$c_{2}^{C_{F}}$ | $4.655$ | $4.654(146)$

Table 1: Two-loop non-cusp anomalous dimension (left) and finite term of the
renormalised jet function (right) for thrust. The analytic results have been
extracted from [16].

Figure 1: The two-loop non-cusp anomalous dimensions (left) are compared
against the semi-analytic results (BRT) of [1]. The new $c_{2}^{J}$
coefficients (right) are obtained for different angularity values. The
numerical uncertainties are too small to be visible on the plot.

##  3 Beam functions

The beam functions are defined via proton matrix elements of collinear field
operators. In contrast to the jet functions discussed above, the beam
functions are non-perturbative objects, which must be matched onto parton
distribution functions (PDF) to extract the perturbative information. At the
partonic level, the quark-to-quark beam function is defined e.g. as

$\displaystyle\frac{1}{2}\,\left[\frac{{\not{n}}}{2}\right]_{\beta\alpha}{\cal
B}_{qq}(x,\tau,\mu)=\sum_{X}\delta\Big{(}(1-x)P_{-}-\sum_{i}k_{i}^{-}\Big{)}\,\bra{P}\bar{\chi}_{\alpha}\ket{X}\bra{X}\chi_{\beta}\ket{P}{\cal
M}(\tau;\\{k_{i}\\})\,,$ (15)

where $\ket{P}$ is now a partonic state with momentum $P^{\mu}$. In the
general case, the matching onto the partonic PDF $f_{ij}$ then takes the form,

$\displaystyle{\cal B}_{ij}(x,\tau,\mu)$
$\displaystyle=\sum_{k}\int_{x}^{1}\frac{dz}{z}\,{\cal
I}_{ik}\Big{(}\frac{x}{z},\tau,\mu\Big{)}~{}f_{kj}(z,\mu)\,.$ (16)

While we use the same notation for the Laplace-space measurement function
${\cal M}(\tau;\\{k_{i}\\})$ as in the previous section, the beam functions
${\cal B}_{ij}$ and the matching kernels ${\cal I}_{ij}$ are distribution-
valued in the momentum fraction $x$. To avoid this, we perform an additional
Mellin transformation, and we denote the corresponding quantities in Mellin
space by $\widehat{\cal B}_{ij}$ and $\widehat{\cal I}_{ij}$. In Mellin space,
the bare quantities have a similar perturbative expansion as the jet functions
in (3).

###  3.1 Calculational details

The collinear matrix elements in the definition of the beam functions can be
extracted from the well-known splitting functions using crossing symmetry.
Notice that at NNLO, the calculation of the beam functions only involves a
two-body phase space, for which we employ similar parametrisations as in the
jet-function case discussed above. In the following, we focus on the quark-to-
quark beam function for transverse momentum resummation, which is defined in
SCET-2, and as such requires an additional rapidity regulator. To this end, we
adopt the analytic regulator from [18] in a symmetrised version (see also [1,
2, 3]), and we follow the collinear anomaly approach [19] to extract the final
matching kernels. Specifically, we write the product of soft and beam-function
kernels in Mellin space $(N_{1},N_{2})$ in the form

$\displaystyle\left[{\cal S}(\tau,\mu,\nu)\;\widehat{\cal
I}_{qq}(N_{1},\tau,\mu,\nu)\;\widehat{\cal
I}_{\bar{q}\bar{q}}(N_{2},\tau,\mu,\nu)\right]_{q^{2}}\stackrel{{\scriptstyle\alpha=0}}{{\equiv}}\left(\bar{\tau}^{2}q^{2}\right)^{-F_{q\bar{q}}(\tau,\mu)}\;\widehat{I}_{qq}(N_{1},\tau,\mu)\;\widehat{I}_{\bar{q}\bar{q}}(N_{2},\tau,\mu)\,,$
(17)

where $\alpha$ is the rapidity regulator, $\nu$ is the rapidity scale and
$q^{2}$ refers to the hard scale in the problem, e.g. the invariant mass of
the Drell-Yan pair. On the right-hand-side of (17), the large rapidity
logarithms are then resummed to all orders through the anomaly coefficient
$F_{q\bar{q}}(\tau,\mu)$. Due to our specific choice of the rapidity
regulator, the collinear and anti-collinear matching kernels on the left-hand-
side are furthermore symmetric under the exchange of
$n^{\mu}\leftrightarrow\bar{n}^{\mu}$, and it is therefore sufficient to
compute only one of them explicitly. On the other hand, we also need the soft
function in the same regularisation scheme, for which we rely on SoftSERVE [1,
2, 3].

The renormalised anomaly coefficient then satisfies the RGE

$\displaystyle\frac{d}{d\ln\mu}\;F_{q\bar{q}}(\tau,\mu)=2\,\Gamma_{\rm
cusp}(\alpha_{s})\,,$ (18)

whose solution, up to two loops, takes the form

$\displaystyle
F_{q\bar{q}}(\tau,\mu)=a_{s}(\mu)\,\Big{\\{}2\Gamma_{0}L+d_{1}\Big{\\}}+a_{s}^{2}(\mu)\,\Big{\\{}2\beta_{0}\Gamma_{0}L^{2}+2\,(\Gamma_{1}+\beta_{0}d_{1})\,L+d_{2}\,\Big{\\}}$
(19)

where $L=\ln(\mu\bar{\tau})$ and $d_{i}$ are the non-logarithmic terms of the
renormalised anomaly coefficient. The matching kernels, on the other hand,
obey the following RGE in Mellin space,

$\displaystyle\frac{d}{d\ln\mu}\;\widehat{I}_{qq}(N,\tau,\mu)=2\left[\Gamma_{\rm
cusp}(\alpha_{s})\,L+\gamma^{B}(\alpha_{s})\right]\widehat{I}_{qq}(N,\tau,\mu)-2\sum_{j}\;\widehat{I}_{qj}(N,\tau,\mu)\,\widehat{P}_{jq}(N,\mu)\,,$
(20)

where $\gamma^{B}$ is the non-cusp anomalous dimension, and the sum in the
last term of the right-hand-side runs over all partons. In this notation,
$\widehat{P}_{qq}$ ($\widehat{P}_{gq}$) is simply the Mellin-transformed
splitting function $P_{q\to qg^{*}}$ ($P_{q\to gq^{*}}$). From the solution of
(20), we finally extract the non-logarithmic pieces
$\textstyle\widehat{I}_{qq,i}(N)$ by expanding in $a_{s}(\mu)$ and setting
$\mu=1/\bar{\tau}$, similar to (13).

###  3.2 Results

At NNLO we have calculated all contributions except for the $C_{F}C_{A}$
colour structure appearing in the $P_{q\to ggq^{*}}$ splitting kernel.

|  Analytic |  This work
---|---|---
$\gamma_{1}^{n_{f}}$ | $-11.395$ | $-11.392(8)$
$\gamma_{1}^{C_{F}}$ | $10.610$ | $10.594(40)$

|  Analytic |  This work
---|---|---
$d_{2}^{n_{f}}$ | $4.148$ | $4.147(4)$
$d_{2}^{C_{F}}$ | $0$ | $0.024(45)$

Table 2: Two-loop non-cusp anomalous dimension (left) and anomaly coefficient
(right) for $p_{T}$-resummation. The analytic results have been extracted from
[20, 21].

As can be seen in Table 2, our results for the remaining colour structures of
the anomaly coefficient $d_{2}$ as well as the non-cusp anomalous dimension
$\gamma_{1}^{B}$ are consistent within numerical uncertainties with the
analytic results from [20, 21]. In Figure 2, we display the non-logarithmic
contribution to the renormalised matching kernel
$\textstyle\widehat{I}_{qq,2}(N)$ as a function of the Mellin parameter $N$.
For both colour structures computed so far, we again find a very good
agreement with the known analytic results.

Figure 2: Comparison of the two-loop non-logarithmic contribution to the
quark-to-quark matching kernel for $p_{T}$ resummation against the analytic
results of [20, 21].

##  4 Conclusions

We have presented a formalism to automate the calculation of two-loop jet and
beam functions. Due to the complicated divergence structure of the underlying
collinear matrix elements, we encountered many overlapping singularities in
the RR contribution, which we disentangled with the help of a mixed strategy
based on sector decomposition, non-linear transformations and selector
functions. We furthermore validated our setup against known results for the
thrust jet function and the transverse-momentum-dependent beam function, and
we obtained a new prediction for the angularity jet function. While we have
not yet finished the implementation of the RR contribution in the beam-
function case, we plan to look into more observables soon. In the longer term,
we also envisage to provide a public code for the calculation of jet and beam
functions, similar to the SoftSERVE distribution.

#### Funding information

This research was supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under grant 396021762 - TRR 257.

## References

  * [1] G. Bell, R. Rahn and J. Talbert,  _Two-loop anomalous dimensions of generic dijet soft functions_ ,  Nucl. Phys. B 936, 520 (2018),  10.1016/j.nuclphysb.2018.09.026,  1805.12414.
  * [2] G. Bell, R. Rahn and J. Talbert,  _Generic dijet soft functions at two-loop order: correlated emissions_ ,  JHEP 07, 101 (2019),  10.1007/JHEP07(2019)101,  1812.08690.
  * [3] G. Bell, R. Rahn and J. Talbert,  _Generic dijet soft functions at two-loop order: uncorrelated emissions_ ,  JHEP 09, 015 (2020),  10.1007/JHEP09(2020)015,  2004.08396.
  * [4] G. Bell, B. Dehnadi, T. Mohrmann and R. Rahn,  _Automated Calculation of N-jet Soft Functions_ ,  PoS LL2018, 044 (2018),  10.22323/1.303.0044,  1808.07427.
  * [5] K. Brune,  _NLO calculation of jet and beam functions in Soft-Collinear Effective Theory_ ,  Master’s thesis, Universität Siegen (2019).
  * [6] A. Basdew-Sharma, F. Herzog, S. Schrijnder van Velzen and W. J. Waalewijn,  _One-loop jet functions by geometric subtraction_ ,  JHEP 10, 118 (2020),  10.1007/JHEP10(2020)118,  2006.14627.
  * [7] G. Altarelli and G. Parisi,  _Asymptotic Freedom in Parton Language_ ,  Nucl. Phys. B 126, 298 (1977),  10.1016/0550-3213(77)90384-4.
  * [8] Z. Bern and G. Chalmers,  _Factorization in one loop gauge theory_ ,  Nucl. Phys. B 447, 465 (1995),  10.1016/0550-3213(95)00226-I,  hep-ph/9503236.
  * [9] Z. Bern, V. Del Duca, W. B. Kilgore and C. R. Schmidt,  _The infrared behavior of one loop QCD amplitudes at next-to-next-to leading order_ ,  Phys. Rev. D 60, 116001 (1999),  10.1103/PhysRevD.60.116001,  hep-ph/9903516.
  * [10] D. A. Kosower and P. Uwer,  _One loop splitting amplitudes in gauge theory_ ,  Nucl. Phys. B 563, 477 (1999),  10.1016/S0550-3213(99)00583-0,  hep-ph/9903515.
  * [11] J. M. Campbell and E. W. N. Glover,  _Double unresolved approximations to multiparton scattering amplitudes_ ,  Nucl. Phys. B 527, 264 (1998),  10.1016/S0550-3213(98)00295-8,  hep-ph/9710255.
  * [12] S. Catani and M. Grazzini,  _Collinear factorization and splitting functions for next-to-next-to-leading order QCD calculations_ ,  Phys. Lett. B 446, 143 (1999),  10.1016/S0370-2693(98)01513-5,  hep-ph/9810389.
  * [13] S. Catani and M. Grazzini,  _Infrared factorization of tree level QCD amplitudes at the next-to-next-to-leading order and beyond_ ,  Nucl. Phys. B 570, 287 (2000),  10.1016/S0550-3213(99)00778-6,  hep-ph/9908523.
  * [14] G. Bell, K. Brune, G. Das and M. Wald,  to appear .
  * [15] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk and T. Zirke,  _pySecDec: a toolbox for the numerical evaluation of multi-scale integrals_ ,  Comput. Phys. Commun. 222, 313 (2018),  10.1016/j.cpc.2017.09.015,  1703.09692.
  * [16] T. Becher and M. Neubert,  _Toward a NNLO calculation of the $\bar{B}\to X_{s}\gamma$ decay rate with a cut on photon energy. II. Two-loop result for the jet function_,  Phys. Lett. B 637, 251 (2006),  10.1016/j.physletb.2006.04.046,  hep-ph/0603140.
  * [17] G. Bell, A. Hornig, C. Lee and J. Talbert,  _$e^{+}e^{-}$ angularity distributions at NNLL′ accuracy_,  JHEP 01, 147 (2019),  10.1007/JHEP01(2019)147,  1808.07867.
  * [18] T. Becher and G. Bell,  _Analytic Regularization in Soft-Collinear Effective Theory_ ,  Phys. Lett. B 713, 41 (2012),  10.1016/j.physletb.2012.05.016,  1112.3907.
  * [19] T. Becher and M. Neubert,  _Drell-Yan Production at Small $q_{T}$, Transverse Parton Distributions and the Collinear Anomaly_,  Eur. Phys. J. C 71, 1665 (2011),  10.1140/epjc/s10052-011-1665-7,  1007.4005.
  * [20] T. Gehrmann, T. Luebbert and L. L. Yang,  _Transverse parton distribution functions at next-to-next-to-leading order: the quark-to-quark case_ ,  Phys. Rev. Lett. 109, 242003 (2012),  10.1103/PhysRevLett.109.242003,  1209.0682.
  * [21] T. Gehrmann, T. Luebbert and L. L. Yang,  _Calculation of the transverse parton distribution functions at next-to-next-to-leading order_ ,  JHEP 06, 155 (2014),  10.1007/JHEP06(2014)155,  1403.6451.


# On Semi-supervised Estimation
of Discrete Distributions under $f$-divergences

H.S.Melihcan Erol and Lizhong Zheng  Dept. EECS and RLE
Massachusetts Institute of Technology
Cambridge, MA 02139
{hsmerol<EMAIL_ADDRESS>
###### Abstract

We study the problem of estimating the joint probability mass function (pmf)
over two random variables. In particular, the estimation is based on the
observation of $m$ samples containing both variables and $n$ samples missing
one fixed variable. We adopt the minimax framework with $l^{p}_{p}$ loss
functions. Recent work established that univariate minimax estimator
combinations achieve minimax risk with the optimal first-order constant for
$p\geq 2$ in the regime $m=o(n)$, questions remained for $p\leq 2$ and various
$f$-divergences. In our study, we affirm that these composite estimators are
indeed minimax optimal for $l^{p}_{p}$ loss functions, specifically for the
range $1\leq p\leq 2$, including the critical $l_{1}$ loss. Additionally, we
ascertain their optimality for a suite of $f$-divergences, such as KL,
$\chi^{2}$, Squared Hellinger, and Le Cam divergences.

##  I Introduction

Estimating the probability mass function (pmf) is a pivotal task in
statistics, and the minimax methodology has emerged as a fundamental approach
in this context, a fact well-established in [1]. Pioneering studies adeptly
navigated the realm of minimax risk in pmf estimation under $l^{2}_{2}$ loss,
culminating in the identification of the quintessential minimax estimator [2,
3, 4]. Further research refined these insights by accurately determining the
primary constant for the minimax risk concerning KL-divergence, $l_{1}$, and
$f$-divergences [5, 6, 7].

The recent growth in dataset sizes often exceeds the capacity for thorough
labeling, leading to datasets with a significant number of samples missing
specific information. The estimators discussed in studies such as [2, 3, 4, 5,
6, 7] fall short as they fail to effectively leverage these unlabeled samples,
resulting in suboptimal performance. This challenge underscores the importance
of refined semi-supervised learning techniques that adeptly combine both
labeled and unlabeled data. This paper strives to contribute to this area by
enhancing the approach introduced in [8], extending its reach to include a
spectrum of $f$-divergences including KL, $\chi^{2}$, Squared Hellinger, and
Le Cam divergences.

Echoing the framework in [8], our analysis is concentrated on two dependent
random variables, $X$ and $Y$, with a distribution denoted by $p_{XY}$. We
scrutinize two distinct datasets: the first composed of $m$ i.i.d. samples of
$(x_{i},y_{i})$ sourced from $p_{XY}$, and the second consisting of $n$
samples of $x_{j}$ derived from the marginal distribution $p_{X}$. This work
extends and enriches the foundational principles set forth in [8] by
demonstrating that the composition of univariate minimax estimators retains
the first-order risk term’s accuracy across $l^{p}_{p}$ loss functions for the
range $1\leq p\leq 2$ and encompasses a broad array of $f$-divergences. This
spectrum includes, notably, total-variation, $\chi^{2}$, Squared Hellinger,
and Lecam divergences, thereby broadening the scope of the initial framework
which focused on $l^{p}_{p}$ loss functions for $p\geq 2$.

##  II Notations & Preliminaries

We say $a_{n}=o(b_{n})$ if $\lim\sup_{n}\frac{a_{n}}{b_{n}}=0$,
$a_{n}=O(b_{n})$ if $\lim\sup_{n}\frac{a_{n}}{b_{n}}=K<\infty$,
$\Theta(a_{n})=b_{n}$ if $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. $\delta_{x}$
indicates a unit point mass at $x\in{\mathcal{X}}$. We write for
$x,y\in\mathbb{R}$, $x\wedge y\triangleq\min(x,y)$, $x\vee
y\triangleq\max(x,y)$. We denote the space of probability distributions over
the finite set ${\mathcal{X}}$ by $\Delta_{{\mathcal{X}}}$. We reserve the
symbols $k_{x}=\left|{\mathcal{X}}\right|$ and
$k_{y}=\left|{\mathcal{Y}}\right|$. We use the upper case of a letter to
denote a random variable and the lower case to indicate the realization of
that random variable. Let
$\mathcal{L}:\Delta_{{\mathcal{X}}}\times\Delta_{{\mathcal{X}}}\rightarrow\mathbb{R}$
be a loss function.

In the minimax setting [1], we assume that nature adversarially chooses a
distribution $p_{X}$; $n$ samples $u^{n}=\\{x_{i}\\}^{n}_{i=1}$ are drawn
i.i.d. from this distribution; and the goal is to design an estimator
$\hat{q}:{\mathcal{X}}^{n}\to\Delta_{\mathcal{X}}$ based on the samples
$u^{n}$ to minimize the expected loss. We denote the associated risk by
$r^{\mathcal{L}}_{n}$, and the problem is formulated as:

$r^{\mathcal{L}}_{n}\triangleq\min_{\hat{q}_{X}}\max_{p_{X}\in\Delta_{{\mathcal{X}}}}\operatorname{\mathbb{E}}_{U^{n}}\left[{\mathcal{L}(p_{X},\hat{q}_{X}(U^{n}))}\right]$
(1)

In this work we will consider two different forms of loss function:

###  II-A $l^{p}_{p}$ losses for $1\leq p<2$

In this case ${\mathcal{L}}$ is chosen to be:

$\mathcal{L}(p_{X},q_{X})=\|p_{X}-q_{X}\|^{p}_{p}\triangleq\sum_{x\in{\mathcal{X}}}\Big{(}\left|p_{X}(x)-q_{X}(x)\right|\Big{)}^{p}$
(2)

In particular the risk becomes for $l^{p}_{p}$ losses:

$r^{p}_{n}\triangleq\min_{\hat{q}_{X}}\max_{p_{X}}\operatorname{\mathbb{E}}_{U^{n}}\left[{\|p_{X}-\hat{q}_{X}(U^{n})\|^{p}_{p}}\right]$
(3)

The choices of $l^{2}_{2}$ and $l_{1}$ are has been considered in the prior
work [2, 4, 6]. For general $p$ as we consider in this work the minimax risk
scales as $r^{p}_{n}=\Theta(n^{-\frac{p}{2}})$ [8]. Similarly we will adopt
the convention in [8] and denote the constant of this rate by
$r^{n}_{p}=C_{p}(1+o(1))n^{-\frac{p}{2}})$.

One of the main goals of this paper is to study the extension of (3) into a
multivariate semi-supervised setting. In this new setting, the nature chooses
a model $p_{XY}$ of a pair of jointly distributed random variables $X,Y$. The
estimator has access to two datasets: a collection of complete
$l^{m}\triangleq\\{(x^{\prime}_{i},y^{\prime}_{i})\\}^{m}_{i=1}$ and a
collection of incomplete samples $u^{n}\triangleq\\{x_{i}\\}^{n}_{i=1}$
generated i.i.d from $p_{XY}$ and $p_{X}$. Our goal is to design an estimator
with minimal expected risk. We formulate this new problem as:

$\displaystyle
R^{p}_{m,n}\triangleq\min_{\hat{q}_{XY}}\max_{p_{XY}}\operatorname{\mathbb{E}}_{U^{n},L^{m}}\left[{\|p_{XY}-\hat{q}_{XY}(U^{n},L^{m})\|^{p}_{p}}\right]$
(4)

Throughout our study of (4), we will consider the following auxiliary
problems:

$\displaystyle R^{p}_{m}\triangleq\min_{\hat{q}_{Y\mid
X}}\max_{p_{XY}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\|p_{XY}-p_{X}\hat{q}_{Y\mid
X}(L^{m})\|^{p}_{p}}\right]$
$\displaystyle\bar{R}^{p}_{m}\triangleq\min_{\hat{q}_{Y\mid
X}}\max_{p_{X}}\operatorname{\mathbb{E}}_{L^{m}_{X}}\left[{\max_{p_{Y\mid
X}}\operatorname{\mathbb{E}}_{L^{m}_{Y}}\left[{\|p_{XY}-p_{X}\hat{q}_{Y\mid
X}(L^{m})\|^{p}_{p}}\right]}\right]$

###  II-B $f$-divergences

$f$-divergences represent a crucial family of statistical loss functions. In
the context of discrete distributions, f-divergences assume the form:

$\mathcal{L}(p_{X},q_{X})=D_{f}(p_{X}\|q_{X})\triangleq\sum_{x\in{\mathcal{X}}}q_{X}(x)f\left(\frac{p_{X}(x)}{q_{X}(x)}\right)$
(5)

for a function strictly convex $f:{\mathcal{X}}\rightarrow\mathbb{R}$ with
$f(1)=0$. The minmax estimation problem under consideration for this loss is
speacialized to:

$r^{f}_{n}\triangleq\min_{\hat{q}_{X}}\max_{p_{X}\in\Delta_{\delta}}\operatorname{\mathbb{E}}_{U^{n}}\left[{D_{f}(p_{X}\|\hat{q}_{X}(U^{n}))}\right]$
(6)

where
$\Delta_{\delta}=\\{p_{X}\in{\mathcal{P}}({\mathcal{X}}):p_{X}(x)\geq\delta,\quad\forall
x\in{\mathcal{X}}\\}$ and we further assume that $f$ is thrice differentiable
and is subexponential, i.e.
$\limsup_{x\rightarrow\infty}\frac{\left|f(x)\right|}{e^{cx}}=0$ for all
$c\in\mathbb{R}$. Under these constraint minmax risk is
$r^{f}=C_{f}n^{-1}+o(n^{-1})$ $C_{f}\triangleq
f^{\prime\prime}(1)(\left|{\mathcal{X}}\right|-1)/2$ [7]. In addition we
impose on $f$ that $\lim_{x\rightarrow 0^{+}}\left|f(x)\right|<\infty$ and
$\lim_{x\rightarrow 0^{+}}x\left|f^{\prime}(x)\right|<\infty$. In the sequel,
we will hide the underlying set $\Delta_{\delta}$. We note that these mild
assumptions are satisfied by common f-divergences including
KL,$\chi^{2}$,Squared Hellinger and Le Cam divergences.

Similarly we define the main problem of interest for $f$-divergences:

$\displaystyle
R^{f}_{m,n}\triangleq\min_{\hat{q}_{XY}}\max_{p_{XY}}\operatorname{\mathbb{E}}_{U^{n},L^{m}}\left[{D_{f}(p_{XY}\|\hat{q}_{XY}(U^{n},L^{m}))}\right]$

and auxilary problems:

$\displaystyle R^{f}_{m}\triangleq\min_{\hat{q}_{Y\mid
X}}\max_{p_{XY}}\operatorname{\mathbb{E}}_{L^{m}}\left[{D_{f}(p_{Y\mid
X}\|\hat{q}_{Y\mid X}(L^{m})\mid p_{X})}\right]$
$\displaystyle\bar{R}^{f}_{m}\triangleq\min_{\hat{q}_{Y\mid
X}}\max_{p_{X}}\mathbb{E}_{L^{m}_{X}}{\max_{p_{Y\mid
X}}\operatorname{\mathbb{E}}_{L^{m}_{Y}}\left[{D_{f}(p_{Y\mid
X}\|\hat{q}_{Y\mid X}(L^{m})\mid p_{X})}\right]}$

where we introduce the conditional $f$-divergence:

$D_{f}(p_{Y\mid X}\|q_{Y\mid
X}|p_{X})=\sum_{x\in{\mathcal{X}}}p_{X}(x)D_{f}(p_{Y\mid X=x}\|q_{Y\mid X=x})$

##  III Results

###  III-A $l^{p}_{p}$ loss functions

######  Theorem 1 (Theorem 1 of [8])

Let $\hat{q}^{*}_{n}$ be a minimax optimal estimator for $r^{p}_{n}$. Then the
conditional composition $\hat{q}^{*,m}_{Y\mid X}$ based on $\hat{q}^{*}_{n}$
is minimax optimal for $\bar{R}^{p}_{m}$:

$\displaystyle\max_{p_{X}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\max_{p_{Y|X}}\|p_{XY}-p_{X}\hat{q}^{*,m}_{Y\mid
X}\|^{p}_{p}}\right]=\bar{R}^{p}_{m}$ (7)

###### Theorem 2

Let $1\leq p\leq 2$ and $\hat{q}^{*}_{n}$ be a first order minimax optimal
estimator for $r^{p}_{n}$. Then the conditional composition
$\hat{q}^{*,m}_{Y\mid X}$ based on $\hat{q}^{*}_{n}$ is first order minimax
optimal for $R^{p}_{m}$:

$\displaystyle\max_{p_{XY}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\|p_{XY}-p_{X}\hat{q}^{*,m}_{Y\mid
X}\|^{p}_{p}}\right]=R^{p}_{m}+o\left(R^{p}_{m}\right)$ (8)

######  Theorem 3 (Theorem 3 of [8])

Let $m=o(n)$:

$\displaystyle\left|R^{p}_{m,n}-R^{p}_{m}\right|\leq
O\left({m^{-\frac{p-1}{2}}(n)^{-1/2}}\right)$ (9)

###### Theorem 4

Let $1\leq p\leq 2$ and $\hat{q}^{*}_{n}$ be a first-order optimal estimator
for $r^{p}_{n}$ with minimax risk:

$\displaystyle\max_{p_{X}}\operatorname{\mathbb{E}}\left[{\|p_{X}-\hat{q}^{*}_{n}\|^{p}_{p}}\right]=\min_{\hat{q}}\max_{p_{X}}\operatorname{\mathbb{E}}\left[{\|p_{X}-\hat{q}_{n}\|^{p}_{p}}\right]\simeq\frac{C_{p}}{n^{\frac{p}{2}}}$

Then the joint composition $\hat{q}^{*,m,n}_{XY}$ based on $\hat{q}_{n}$ is
first order minimax optimal for $R^{p}_{m,n}$ in the regime $m=o(n).$
Furthermore,

$\displaystyle
R^{p}_{n,m}=\left(k_{x}\right)^{1-\frac{p}{2}}\frac{C_{p}}{m^{\frac{p}{2}}}+o\left(\frac{1}{m^{\frac{p}{2}}}\right)$

###  III-B f-divergences

###### Theorem 5

Let $\hat{q}^{*}_{n}$ be a minimax optimal estimator for $r^{f}_{n}$. Then the
conditional composition $\hat{q}^{*,m}_{Y\mid X}$ based on $\hat{q}^{*}_{n}$
is minimax optimal for $\bar{R}^{f}_{m}$:

$\displaystyle\max_{p_{X}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\max_{p_{Y|X}}D_{f}(\check{q}^{*,m}_{Y\mid
X}\|\check{q}^{*}_{n}\mid p_{X})}\right]=\bar{R}^{f}_{m}$ (10)

###### Corollary 1

Let $\hat{q}^{*}_{n}$ be a first order minimax optimal estimator for
$r^{f}_{n}$. Then the conditional composition $\hat{q}^{*,m}_{Y\mid X}$ based
on $\hat{q}^{*}_{n}$ is first order minimax optimal for $R^{p}_{m}$:

$\displaystyle\max_{p_{XY}}\operatorname{\mathbb{E}}_{L^{m}}\left[{D_{f}(p_{Y\mid
X}\|\hat{q}^{*,m}_{Y\mid X}\mid
p_{X})}\right]=R^{f}_{m}+o\left(R^{f}_{m}\right)$ (11)

###### Theorem 6

Let $m=o(n)$:

$\left|R^{f}_{m,n}-R^{f}_{m}\right|=o(R^{f}_{m})$

###### Theorem 7

Let $\hat{q}^{*}_{n}$ be a first-order optimal estimator for $r^{p}_{n}$ with
minimax risk:

$\displaystyle\max_{p_{X}}\operatorname{\mathbb{E}}\left[{D_{f}(p_{X}\|\hat{q}^{*}_{n})}\right]=\min_{\hat{q}}\max_{p_{X}}\operatorname{\mathbb{E}}\left[{D_{f}(p_{X}\|\hat{q}_{n})}\right]\simeq
C_{f}n^{-1}$

Then the joint composition $\hat{q}^{*,m,n}_{XY}$ based on $\hat{q}_{n}$ is
first order minimax optimal for $R^{p}_{m,n}$ in the regime $m=o(n).$
Furthermore,

$\displaystyle
R^{f}_{n,m}=\left|{\mathcal{X}}\right|\frac{C_{f}}{m}+o\left(\frac{1}{{m}}\right)$

##  IV Proofs for Theorems

###  IV-A $l^{p}_{p}$ loss functions

###### Proof:

By (12) in the proof for Theorem 1,

$\displaystyle\bar{R}^{p}_{m}$ $\displaystyle=$
$\displaystyle\max_{p_{X}}\sum_{x\in{\mathcal{X}}}\sum^{m}_{i=0}{m\choose
i}\left(p_{X}(x)\right)^{i+p}\left(1-p_{X}(x)\right)^{m-i}r^{p}_{i}$ (12)
$\displaystyle=$
$\displaystyle\max_{p_{X}}\sum_{x\in{\mathcal{X}}}C_{p}\left(\frac{p_{X}(x)}{m}\right)^{\frac{p}{2}}+o\left(\frac{1}{m^{\frac{p}{2}}}\right)$
(13) $\displaystyle=$
$\displaystyle\sum_{x\in{\mathcal{X}}}C_{p}\left(\frac{\frac{1}{\left|{\mathcal{X}}\right|}}{m}\right)^{\frac{p}{2}}+o\left(\frac{1}{m^{\frac{p}{2}}}\right)$
$\displaystyle=$
$\displaystyle{\left|{\mathcal{X}}\right|}^{1-\frac{p}{2}}\frac{C_{p}}{m^{\frac{p}{2}}}+o\left(\frac{1}{m^{\frac{p}{2}}}\right)$

(13) follows from Lemma 2. In (12), the objective function is symmetric in
variables $\\{p_{X}(x):x\in{\mathcal{X}}\\}$ and for $p\leq 2$ each summand is
a concave function. Therefore the optimal solution should be the uniform
distribution, i.e. $p^{*}(x)=\frac{1}{\left|{\mathcal{X}}\right|}$ which leads
to (IV-A).

On the other hand,

$\displaystyle R^{p}_{m}$ $\displaystyle\triangleq$
$\displaystyle\min_{\hat{q}_{Y\mid
X}}\max_{p_{XY}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\|p_{XY}-p_{X}\hat{q}_{Y\mid
X}(L^{m})\|^{p}_{p}}\right]$ (15) $\displaystyle\geq$
$\displaystyle\min_{\hat{q}_{Y\mid X}}\max_{p_{Y\mid
X}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\|U_{{\mathcal{X}}}p_{Y\mid
X}-U_{{\mathcal{X}}}\hat{q}_{Y\mid X}(L^{m})\|^{p}_{p}}\right]$ (16)
$\displaystyle=$ $\displaystyle\min_{\hat{q}_{Y\mid X}}\max_{p_{Y\mid
X}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\sum_{x\in{\mathcal{X}}}\frac{1}{\left|{\mathcal{X}}\right|^{p}}\|p_{Y\mid
X=x}-\hat{q}_{Y\mid X=x}\|^{p}_{p}}\right]$ (17) $\displaystyle=$
$\displaystyle\min_{\hat{q}_{Y\mid
X}}\sum_{x\in{\mathcal{X}}}\frac{1}{\left|{\mathcal{X}}\right|^{p}}\max_{p_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\|p_{Y\mid X=x}-\hat{q}_{Y\mid
X=x}\|^{p}_{p}}\right]$ $\displaystyle=$
$\displaystyle\sum_{x\in{\mathcal{X}}}\frac{1}{\left|{\mathcal{X}}\right|^{p}}\min_{\hat{q}_{Y\mid
X=x}}\max_{p_{Y\mid X=x}}\operatorname{\mathbb{E}}_{L^{m}}\left[{\|p_{Y\mid
X=x}-\hat{q}_{Y\mid X=x}\|^{p}_{p}}\right]$ $\displaystyle=$
$\displaystyle\sum_{x\in{\mathcal{X}}}\frac{1}{\left|{\mathcal{X}}\right|^{p}}(1+o(1))\frac{C_{p}}{(\frac{m}{\left|{\mathcal{X}}\right|})^{\frac{p}{2}}}$
(20) $\displaystyle=$
$\displaystyle(1+o(1))\left|{\mathcal{X}}\right|^{1-\frac{p}{2}}\frac{C_{p}}{m^{\frac{p}{2}}}$
(21)

In (15) we substitute uniform distribution as $p_{X}$ in the inner
maximization problem to arrive at the lower bound in (16). In (IV-A) we expand
the $l^{p}_{p}$ norm. In $\eqref{eq:lem:bar_red_lp:2:4}$ we note that
optmization variables $\\{p_{Y\mid X=x}:x\in{\mathcal{X}}\\}$ and
$\\{\check{q}_{Y\mid X=x}:x\in{\mathcal{X}}\\}$ are independent in the order.
(20) follows from Lemma 1. Finally, lemma follows by the relation
$R^{p}_{m}\leq\bar{R}^{p}_{m}$ along with (IV-A) and (21) ∎

###### Proof:

By Theorem-3 of [8] $\left|R^{p}_{m,n}-R^{p}_{m}\right|\leq
O(m^{-\frac{p-1}{2}}n^{-\frac{1}{2}})$. Furthermore, by Theorem-1 the
composition estimator is minimax optimal estimator for $\bar{R}^{p}_{m}$.
However by Theorem 2 $\bar{R}^{p}_{m}-R^{p}_{m}\leq o\left(R^{p}_{m}\right)$
for $1\leq p<2$ as well. Therefore, the composition estimator is first order
minimax for $R^{p}_{m,n}$ for $1\leq p\leq 2$. ∎

###  IV-B $f$-divergences

###### Proof:

Let us define:

$\displaystyle h(p_{X},\hat{q}_{Y\mid
X})\triangleq\operatorname{\mathbb{E}}_{L^{m}_{X}}\left[{\max_{p_{Y\mid
X}}\operatorname{\mathbb{E}}_{L^{m}_{Y}}\left[{D_{f}(p_{Y\mid
X}\|\hat{q}_{Y\mid X}\mid p_{X})}\right]}\right]$ (22)

and let $\hat{q}^{**}_{Y\mid
X}:({\mathcal{X}}\times{\mathcal{Y}})^{m}\rightarrow(\Delta_{{\mathcal{Y}}})^{\left|{\mathcal{X}}\right|}$
be an estimator for the conditional distribution $p_{Y\mid X}$. Because of the
inequalities $\max_{p_{X}}h(p_{X},\hat{q}^{**}_{Y\mid
X})\geq\bar{R}^{p}_{m}\triangleq\min_{\hat{q}}\max_{p_{X}}\
h(p_{X},\hat{q})\geq\max_{p_{X}}\min_{\hat{q}}h(p_{X},\hat{q})$ it is
sufficient to show that for all $p_{X}$, $\min_{\hat{q}_{Y\mid
X}}h(p_{X},\hat{q}_{Y\mid X})=h(p_{X},\hat{q}^{**}_{Y\mid X})$. Now let us
show the latter is satisfied by the composition estimator
$\hat{q}^{*,m}_{Y\mid X}$. Set a $p_{X}$ and we define
$p_{i,x}\triangleq\operatorname{\mathbb{P}}\left({T_{x}(L_{X})=i}\right)$.
Then the left-hand side becomes:

$\displaystyle\min_{\hat{q}_{Y\mid X}}h(p_{X},\hat{q}_{Y\mid X})$
$\displaystyle=$ $\displaystyle\min_{\hat{q}_{Y\mid
X}}\operatorname{\mathbb{E}}_{L^{m}_{Y}}\left[{D_{f}(p_{Y\mid
X}\|\hat{q}_{Y\mid X}\mid p_{X})}\right]$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\min_{\hat{q}_{Y\mid
X}}\mathbb{E}_{L^{m}_{X}}\Big{[}{\max_{p_{Y\mid
X}}\mathbb{E}_{L^{m}_{Y}}\Big{[}{\sum_{x\in{\mathcal{X}}}p_{x}D_{f}(p_{Y|X=x}\|\hat{q}_{Y\mid
X=x})}\Big{]}}\Big{]}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\sum_{x\in{\mathcal{X}}}p_{x}\min_{\hat{q}_{Y|X=x}}\mathbb{E}_{L^{m}_{X}}\Big{[}{\max_{p_{Y|X=x}}\operatorname{\mathbb{E}}_{L^{m}_{Y}}\left[{D_{f}(p_{Y|X=x}\|\hat{q}_{Y|X=x})}\right]}\Big{]}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{}$ (23)
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\sum_{x\in{\mathcal{X}}}p_{X}(x)\Big{(}\sum^{m}_{i=0}p_{i,x}r^{f}_{i}\Big{)}}$
(24)

In (23) we note that the optimization variables are independent. Steps leading
(23) to (24) are given below and these steps demonstrate that
$\min_{\hat{q}_{Y\mid X}}h(p_{X},\hat{q}_{Y\mid
X})=h(p_{X},\hat{q}^{**}_{Y\mid X})$ holds for the estimator
$\hat{q}^{*,m}_{Y\mid X}$:

$\displaystyle\IEEEeqnarraymulticol{3}{l}{\min_{\hat{q}_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{m}_{X}}\left[{\max_{p_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{m}_{Y}}\left[{D_{f}(p_{Y\mid
X=x}\|\hat{q}_{Y\mid X=x})}\right]}\right]}$ $\displaystyle=$
$\displaystyle\min_{\hat{q}_{Y\mid X=x}}\sum^{m}_{i=0}p_{i,x}\max_{p_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{i}_{Y}}\left[{D_{f}(p_{Y\mid
X=x}\|\hat{q}^{i}_{Y\mid X=x})}\right]$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\sum^{m}_{i=0}p_{i,x}\min_{\hat{q}^{i}_{Y\mid
X=x}}\max_{p_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{i}_{Y}}\left[{D_{f}(p_{Y\mid
X=x}\|\hat{q}^{i}_{Y\mid X=x})}\right]}$ (25)
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\sum^{m}_{i=0}p_{i,x}\max_{p_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{i}_{Y}}\left[{D_{f}(p_{Y\mid
X=x}\|\hat{q}^{*,i}_{Y\mid X=x})}\right]}$ (26)

To arrive in (26), we notice in (25) that $\min_{\hat{q}^{i}_{Y\mid
X=x}}\max_{p_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{i}_{Y}}\left[{D_{f}(p_{Y\mid
X=x}\|\hat{q}_{Y\mid X=x})}\right]$ is the problem $r^{f}_{i}$ and hence is
achieved by $\hat{q}^{*,i}_{Y\mid X}$. ∎

###### Proof:

We proceed similar to the proof for Theorem 2. We have by (24):

$\displaystyle\bar{R}^{f}_{m}=$
$\displaystyle\sum_{x}p_{x}\sum^{m}_{i=0}{m\choose
i}p^{i}_{x}(1-p_{x})^{m-i}r^{f}_{i}$ $\displaystyle=$
$\displaystyle\sum_{x}\frac{1}{p_{x}}\sum^{m}_{i=0}{m\choose
i}p^{2+i}_{x}(1-p_{x})^{m-i}r^{f}_{i}$ (27) $\displaystyle\simeq$
$\displaystyle\sum_{x}\frac{1}{p_{x}}C_{f}\frac{p_{x}}{m}=\left|{\mathcal{X}}\right|\frac{C_{f}}{m}$
(28)

In (28) we evaluate the sum (27) via Lemma 2.Similar to the proof for Theorem
2, we have the lower bound for $R^{f}_{m}$:

$\displaystyle R^{f}_{m}=$ $\displaystyle\min_{\check{q}_{Y\mid
X}}\max_{p_{XY}}\mathbb{E}_{L^{m}}[\sum_{x}p_{x}D_{f}(p_{Y\mid
X=x}\|\check{q}_{Y\mid X=x})]$ $\displaystyle\geq$
$\displaystyle\min_{\check{q}_{Y\mid X}}\max_{p_{Y\mid
X}}\mathbb{E}_{L^{m}}[\sum_{x}\frac{1}{\left|{\mathcal{X}}\right|}D_{f}(p_{Y\mid
X=x}\|\check{q}_{Y\mid X=x})]$ (29) $\displaystyle=$
$\displaystyle\sum_{x\in{\mathcal{X}}}\frac{1}{\left|{\mathcal{X}}\right|}\min_{\check{q}_{Y\mid
X=x}}\max_{p_{Y\mid X=x}}\mathbb{E}_{L^{m}}[D_{f}(p_{Y\mid
X=x}\|\check{q}_{Y\mid X=x})]$ $\displaystyle=$
$\displaystyle\sum_{x\in{\mathcal{X}}}\frac{1}{\left|{\mathcal{X}}\right|}\frac{C_{f}}{\frac{m}{\left|{\mathcal{X}}\right|}}(1+o(1))$
(30) $\displaystyle=$
$\displaystyle\left|{\mathcal{X}}\right|\frac{C_{f}}{m}(1+o(1))$ (31)

in (29) we substitute $p_{X}={\mathcal{U}}_{{\mathcal{X}}}$. In (30) we use
Lemma 1 with $p=2$ as $r^{f}_{i}=\Theta(i^{-1})$. We omit a version of this
lemma for f-divergences as the same proof applies. Finally we note that
$R^{f}_{m}\leq\bar{R}^{f}_{m}$. ∎

###### Proof:

We have the lower bound:

$R^{f}_{n+1,m}\leq R^{f}_{n,m}\implies
R^{f}_{m}=\lim_{n\rightarrow\infty}R^{f}_{m,n}\leq R^{f}_{m,n}$

Now let us prove the upperbound:

$\displaystyle R^{f}_{n,m}$ $\displaystyle=$
$\displaystyle\min_{\hat{q}_{XY}}\max_{p_{XY}}\operatorname{\mathbb{E}}_{U_{n},L_{m}}\left[{D_{f}(p_{XY}\|\check{q}_{XY}(U_{n},L_{m}))}\right]$
$\displaystyle=$
$\displaystyle\min_{\hat{q}_{XY}}\max_{p_{XY}}\operatorname{\mathbb{E}}_{U_{n},L_{m}}\left[{\sum_{x,y}\hat{q}_{XY}(x,y)f\left(\frac{p_{XY}(x,y)}{\check{q}_{XY}(x,y)}\right)}\right]$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}_{U_{n},L_{m}}\Big{[}\sum_{x,y}p_{X}(x)\hat{q}_{Y\mid
X}(y\mid x)f\left(\frac{p_{Y\mid X}(y\mid x)}{\hat{q}_{Y\mid X}(y\mid
x)}\right)}$ (32)
$\displaystyle+(p_{X}(x)-\hat{q}_{X}(x))\Gamma^{U_{n},L_{m}}_{x,y}\Big{]}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}_{U_{n},L_{m}}\Big{[}\sum_{x,y}p_{X}(x)\hat{q}_{Y\mid
X}(y\mid x)f\left(\frac{p_{Y\mid X}(y\mid x)}{\hat{q}_{Y\mid X}(y\mid
x)}\right)\Big{]}}$ $\displaystyle+$
$\displaystyle\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}_{U_{n},L_{m}}\Bigg{[}\sum_{x,y}(p_{X}(x)-\hat{q}_{X}(x))\hat{\Gamma}_{x,y}\Bigg{]}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=R^{f}_{m}+\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}_{U_{n},L_{m}}\bigg{[}\sum_{x,y}(p_{X}(x)-\hat{q}_{X}(x))\Gamma^{U_{n},L_{m}}_{x,y}\bigg{]}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq
R^{f}_{m}+\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}\Big{[}{\Big{(}{{\sum_{x}(p_{x}-\check{q}_{x})^{2}}}\Big{)}^{\frac{1}{2}}\cdot\Big{(}{{\sum_{x}(\sum_{y}\hat{\Gamma}_{x,y})^{2}}}\Big{)}^{\frac{1}{2}}}\Big{]}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq
R^{f}_{m}+\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}\Big{[}{\sum_{x}(p_{x}-\check{q}_{x})^{2}}\Big{]}\mathbb{E}\Big{[}{{{\sum_{x}(\sum_{y}\hat{\Gamma}_{x,y})^{2}}}}\Big{]}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq R^{f}_{m}+}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}\Big{[}{\sum_{x}(p_{x}-\check{q}_{x})^{2}}\Big{]}\cdot\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}\Big{[}{{{\sum_{x}(\sum_{y}\hat{\Gamma}_{x,y})^{2}}}}\Big{]}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=R^{f}_{m}+r^{2}_{n}\cdot\min_{\hat{q}_{XY}}\max_{p_{XY}}\mathbb{E}\Big{[}{{{\sum_{x}(\sum_{y}\hat{\Gamma}_{x,y})^{2}}}}\Big{]}}$
(33) $\displaystyle\leq$ $\displaystyle
R^{f}_{m}+r^{2}_{n}\cdot\max_{p_{XY}}\mathbb{E}\Big{[}{{{\sum_{x}(\sum_{y}\Gamma_{x,y})^{2}}}}\Big{]}$

In (32), we apply mean value theorem pointwise to $p_{X}\rightarrow
p_{XY}f\left(\frac{p_{XY}(x,y)}{\hat{q}_{XY}(x,y)}\right)$ where the remainder
is evaluated at a point
$\xi^{u_{n},l_{m}}_{x,y}\in(p_{X}(x),\check{q}_{X}(x))$. For the compactness
of the notation we abbreviate gradient terms with:

$\Gamma^{U_{n},L_{m}}_{x,y}\triangleq\nabla_{q}\left(q(x)\check{q}_{Y\mid
X}(y\mid x)f\left(\frac{p_{XY}(y\mid
x)}{\xi^{U_{n},L_{m}}_{y,x}\check{q}_{Y\mid X}(y\mid x)}\right)\right)$

and in (33) we chose the estimator
$\hat{q}_{XY}(x,y)=\frac{1}{\left|{\mathcal{X}}\right|\cdot\left|{\mathcal{Y}}\right|}$
for the outer minimization problem which we denote by $\Gamma_{x,y}$. Observe
that the term $\sqrt{r^{2}_{n}}=\Theta(n^{-\frac{1}{2}})$. Therefore we are
done if its multiplier
$\min_{\hat{q}_{XY}}\max_{p_{XY}}\operatorname{\mathbb{E}}\left[{{{\sum_{x}(\sum_{y}\hat{\Gamma}_{x,y})^{2}}}}\right]$
is bounded. In the sequel we show this via the pointwise boundedness of
$\hat{\Gamma}_{x,y}$ which is the case under our assumptions on $f$. First
notice that by the convexity of the perspective function the map
$p_{X}\rightarrow p_{XY}f\left(\frac{p_{XY}(x,y)}{\hat{q}_{XY}(x,y)}\right)$
is convex and thereby its gradients are monotone. Hence it can be bounded by
its values at its end points where either $\xi^{U_{n},L_{m}}_{y,x}=p_{X}(x)$
or
$\xi^{U_{n},L_{m}}_{y,x}=\check{q}_{X}(x)=\frac{1}{\left|{\mathcal{X}}\right|}$,
therefore we analyze these two cases seperately:

#### $\xi^{U_{n},L_{m}}_{y,x}=p_{X}(x)$

$\displaystyle\IEEEeqnarraymulticol{3}{l}{\left|\hat{\Gamma}_{x,y}\right|=\Big{|}\frac{1}{\left|{\mathcal{Y}}\right|}f\Big{(}\frac{p_{Y|X}(y|x)}{\frac{1}{\left|{\mathcal{Y}}\right|}}\Big{)}-\frac{p_{Y|X}(y|x)}{\frac{1}{\left|{\mathcal{Y}}\right|}}f^{\prime}\Big{(}\frac{p_{Y|X}(y|x)}{\frac{1}{\left|{\mathcal{Y}}\right|}}\Big{)}\Big{|}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq\Big{|}\frac{1}{\left|{\mathcal{Y}}\right|}f\Big{(}\frac{p_{Y|X}(y|x)}{\frac{1}{\left|{\mathcal{Y}}\right|}}\Big{)}\Big{|}+\Big{|}\frac{p_{Y|X}(y|x)}{\frac{1}{\left|{\mathcal{Y}}\right|}}f^{\prime}\Big{(}\frac{p_{Y|X}(yx)}{\frac{1}{\left|{\mathcal{Y}}\right|}}\Big{)}\Big{|}}$

In (LABEL:eq:unbounded_case:1) can be unbounded only when $p_{Y\mid
X}\rightarrow 0$, which creates the cases $f(0^{+})$ and
$0^{+}f^{\prime}(0^{+})$ which by assumption are bounded.

####
$\xi^{U_{n},L_{m}}_{y,x}=\check{q}_{X}(x)=\frac{1}{\left|{\mathcal{X}}\right|}$

Similar to (a), after substituting the estimator we arrive at:

$\displaystyle\IEEEeqnarraymulticol{3}{l}{\left|\hat{\Gamma}_{x,y}\right|=\bigg{|}\frac{1}{\left|{\mathcal{Y}}\right|}f\bigg{(}\frac{p_{XY}(x,y)}{\frac{1}{\left|{\mathcal{X}}\right|\left|{\mathcal{Y}}\right|}}\bigg{)}-\frac{p_{XY}(x,y)}{\frac{1}{\left|{\mathcal{X}}\right|}}f^{\prime}\bigg{(}\frac{p_{XY}(x,y)}{\frac{1}{\left|{\mathcal{X}}\right|\left|{\mathcal{Y}}\right|}}\bigg{)}\bigg{|}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq\bigg{|}\frac{1}{\left|{\mathcal{Y}}\right|}f\bigg{(}\frac{p_{XY}(x,y)}{\frac{1}{\left|{\mathcal{X}}\right|\left|{\mathcal{Y}}\right|}}\bigg{)}\bigg{|}+\bigg{|}\frac{p_{XY}(x,y)}{\frac{1}{\left|{\mathcal{X}}\right|}}f^{\prime}\bigg{(}\frac{p_{XY}(x,y)}{\frac{1}{\left|{\mathcal{X}}\right|\left|{\mathcal{Y}}\right|}}\bigg{)}\bigg{|}}$

(LABEL:eq:unbounded_case:2) can be unbounded only when $p_{XY}(x,y)\rightarrow
0$, which again creates the cases $f(0^{+})$ and $0^{+}f^{\prime}(0^{+})$
which by assumption are bounded. ∎

###### Proof:

By Theorem 6 $\left|R^{f}_{m,n}-R^{f}_{m}\right|\leq o(R^{f}_{m})$.
Furthermore, by Theorem 5 the composition estimator is minimax optimal
estimator for $\bar{R}^{f}_{m}$. On the other hand by Theorem 1
$\bar{R}^{f}_{m}-R^{f}_{m}\leq o\left(R^{p}_{m}\right)$. Therefore, the
composition estimator is first order minimax for $R^{f}_{m,n}$. ∎

##  V Supplementary Results

###### Lemma 1

For all $x\in{\mathcal{X}}$:

$\min_{\hat{q}}\max_{p_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{m}\sim{\mathcal{U}}_{{\mathcal{X}}}p_{Y\mid
X}}\left[{\|p_{Y\mid X=x}-\hat{q}\|^{p}_{p}}\right]\simeq
C_{p}\cdot\Big{(}\frac{\left|{\mathcal{X}}\right|}{m}\Big{)}^{\frac{p}{2}}$

###### Proof:

Let us fix constants $C>0$ and $\frac{1}{2}<\alpha<1$.Let us denote the event
${\mathcal{A}}\triangleq\\{\frac{m}{\left|{\mathcal{X}}\right|}-C\cdot
m^{\alpha}\leq T_{x}\leq\frac{m}{\left|{\mathcal{X}}\right|}+C\cdot
m^{\alpha}\\}$ and we let $s_{m}=m\vee(C\cdot
m^{\alpha}+m/\left|{\mathcal{X}}\right|)$ Observe that:

$\displaystyle\IEEEeqnarraymulticol{3}{l}{\operatorname{\mathbb{E}}_{L^{m}}\left[{\|p_{Y\mid
X=x}-\hat{q}\|^{p}_{p}}\right]\geq\operatorname{\mathbb{P}}\left({{\mathcal{A}}}\right)\cdot\operatorname{\mathbb{E}}\left[{\|p_{Y\mid
X=x}-\hat{q}\|^{p}_{p}\mid{\mathcal{A}}}\right]}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\geq\operatorname{\mathbb{P}}\left({{\mathcal{A}}}\right)\operatorname{\mathbb{E}}\left[{\|p_{Y\mid
X=x}-\hat{q}_{Y\mid X=x}\|^{p}_{p}\mid T_{x}=s_{m}}\right]}$ (36)
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\geq(1-2e^{-C^{2}m^{2\alpha-1}})\cdot\operatorname{\mathbb{E}}\left[{\|p_{Y\mid
X=x}-\hat{q}_{Y\mid X=x}\|^{p}_{p}\mid T_{x}=s_{m}}\right]}$ (37)

In (36) we lower bound the conditional expectation by noting that the risk is
smaller if more samples are from the category $x$. In (37) we use bounded
differences inequality to bound the probability of the event
$\operatorname{\mathbb{P}}\left({{\mathcal{A}}}\right)$. Finally, taking the
$\min\max$ of both sides we establish that:

$\displaystyle\min_{\hat{q}}\max_{p_{Y\mid
X=x}}\operatorname{\mathbb{E}}_{L^{m}\sim{\mathcal{U}}_{{\mathcal{X}}}p_{Y\mid
X}}\left[{\|p_{Y\mid X=x}-\hat{q}\|}\right]$
$\displaystyle\geq\left(1-2e^{-C^{2}m^{2\alpha-1}}\right)\cdot
r^{p}_{m\vee\frac{1}{\left|{\mathcal{X}}\right|}+C\cdot m^{\alpha}}\simeq
C_{p}\cdot\Big{(}\frac{\left|{\mathcal{X}}\right|}{m}\Big{)}^{\frac{p}{2}}$

∎

We note that Lemma 2, Lemma 3 are identical to Lemma 1 and Lemma 2 of [8] and
except with adjustments tailored to make them work for the case $1\leq p\leq
2$, which requires a more delicate analysis. For Lemma 2 and Lemma 3 we
introduce:

$\displaystyle H^{n}_{p}(x)$ $\displaystyle\triangleq$
$\displaystyle\sum^{n}_{i=0}{n\choose i}r^{p}_{i}x^{i+p}(1-x)^{n-i}$ (38)
$\displaystyle G^{n}_{p}(x)$ $\displaystyle\triangleq$
$\displaystyle\sum^{n}_{i=0}\binom{n}{i}r^{p}_{i}x^{i}\left(\frac{i}{n}\right)^{p}(1-x)^{n-i}$
(39)

###### Lemma 2

For $1\leq p<2$:

$H^{n}_{p}(x)=C_{p}\
\left(\frac{x}{n}\right)^{\frac{p}{2}}+o(n^{-\frac{p}{2}})$

###### Proof:

Let $c$ be the constant indicated in the statement of Lemma 3. There are two
cases to consider:

  1. 1.

$x\geq c\frac{log^{2}(n)}{n}$:

$\displaystyle H^{n}_{p}(x)$ $\displaystyle\simeq$
$\displaystyle\sum^{n}_{i=0}\binom{n}{i}\frac{C_{p}}{i^{\frac{p}{2}}+1}x^{i+p}(1-x)^{n-i}$
(40)
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\frac{C_{p}}{n^{\frac{p}{2}}}\sum^{n}_{i=0}{n\choose
i}\left(\frac{i}{n}\right)^{\frac{p}{2}}x^{i}(1-x)^{n-i}+O\left(\frac{H^{n}_{p}(x)}{\sqrt{\log{n}}}\right)}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\frac{C_{p}}{n^{\frac{p}{2}}}\left(x^{\frac{p}{2}}+O\left(n^{-\frac{p}{2}}\right)\right)+O\left(\frac{H^{n}_{p}(x)}{\sqrt{\log{n}}}\right)}$
(42)
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=C_{p}\left(\frac{x}{n}\right)^{\frac{p}{2}}+o\left(n^{-\frac{p}{2}}\right)}$
(43)

(40) holds since $r^{p}_{n}\simeq{C_{p}}{n^{-\frac{p}{2}}}$ whereas in
(LABEL:eq:lem1_2) we use Lemma 3. To obtain (42), we use Theorem 8 where we
let $f(x)=x^{\frac{p}{2}}$ and bound the error of $n$th order Bernstein
polynomial approximation $B_{n}$ as:

$\displaystyle\IEEEeqnarraymulticol{3}{C}{\left|B_{n}(x;f)-f(x)\right|=\frac{1}{n}\left|x(1-x)\frac{f^{\prime\prime}(x)}{2}\right|+o(n^{-1})}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{=\frac{1}{n}x(1-x)\frac{p}{4}(1-\frac{p}{2})x^{\frac{p}{2}-2}+o(n^{-1})}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq\frac{1}{n^{\frac{p}{2}}}\frac{p}{4}\Big{(}1-\frac{p}{2}\Big{)}\left({c\log^{2}(n)}\right)^{\frac{p}{2}-1}+o(n^{-1})}$
(44)

where in (44) we alleviate the bounds $1-x\leq 1$ and $x\geq
c\frac{log^{2}(n)}{n}$. Finally by (42) we have
$H^{n}_{p}(x)=O({n^{-\frac{p}{2}}})$ and we obtain (43) by substituting it.

  2. 2.

$x\leq c\frac{log^{2}(n)}{n}$:

$\displaystyle H^{n}_{p}(x)$ $\displaystyle=$
$\displaystyle\sum^{n}_{i=0}{n\choose i}r^{p}_{i}x^{i+p}(1-x)^{n-i}$ (45)
$\displaystyle\leq$ $\displaystyle C_{p}\ x^{p}\sum^{n}_{i=0}{n\choose
i}\frac{1}{i^{\frac{p}{2}}+1}x^{i}(1-x)^{n-i}$ (46) $\displaystyle\leq$
$\displaystyle C_{p}\ x^{p}=O\left(\frac{\log^{2p}(n)}{n^{p}}\right)$ (47)

Similarly
$C_{p}\left(\frac{x}{n}\right)^{\frac{p}{2}}=O\left(\frac{\log^{2p}(n)}{n^{p}}\right)$
when $x\leq c\frac{\log^{2}{n}}{n}$, hence
$\left|C_{p}\left(\frac{x}{n}\right)^{\frac{p}{2}}-H^{n}_{p}(x)\right|=o\left(n^{-\frac{p}{2}}\right)$.

∎

###### Lemma 3

There exists a $c>0$ such that for $x\geq c\frac{log^{2}n}{n}$:

$\left|H^{n}_{p}(x)-G^{n}_{p}(x)\right|=O\left(\frac{H^{n}_{p}(x)}{\sqrt{\log{n}}}\right)$
(48)

###### Proof:

The proof is deferred to the appendix. ∎

##  VI Conclusion

In conclusion, our study contributes to the ongoing discourse on joint
probability mass function estimation by addressing certain aspects within the
minimax framework, particularly for $l^{p}_{p}$ loss functions where $1\leq
p\leq 2$, including the important $l_{1}$ loss. Additionally, our findings
extend the understanding of minimax optimality to encompass a range of
$f$-divergences, such as KL, $\chi^{2}$, Squared Hellinger, and Le Cam
divergences. By doing so, our work not only builds upon and extends existing
methodologies in crucial directions.

## Acknowledgment

This work was supported in part by the ONR under Grant N00014-19-1-2621

## References

  * [1] A. Wald, “Statistical decision functions,” _The Annals of Mathematical Statistics_ , vol. 20, no. 2, pp. 165–205, 1949. [Online]. Available: http://www.jstor.org/stable/2236853
  * [2] S. Trybula, “Some problems of simultaneous minimax estimation,” _The Annals of Mathematical Statistics_ , vol. 29, no. 1, pp. 245–253, 1958.
  * [3] I. Olkin and M. Sobel, “Admissible and minimax estimation for the multinomial distribution and for k independent binomial distributions,” _The Annals of Statistics_ , vol. 7, no. 2, Mar. 1979. [Online]. Available: https://doi.org/10.1214/aos/1176344613
  * [4] M. Wilczyński, “Minimax estimation for the multinomial and multivariate hypergeometric distributions,” _Sankhyā: The Indian Journal of Statistics, Series A_ , pp. 128–132, 1985.
  * [5] D. Braess and T. Sauer, “Bernstein polynomials and learning theory,” _Journal of Approximation Theory_ , vol. 128, no. 2, pp. 187–206, 2004.
  * [6] Y. Han, J. Jiao, and T. Weissman, “Minimax estimation of discrete distributions under $l_{1}$ loss,” _CoRR_ , vol. abs/1411.1467, 2014. [Online]. Available: http://arxiv.org/abs/1411.1467
  * [7] S. Kamath, A. Orlitsky, D. Pichapati, and A. T. Suresh, “On learning distributions from their samples,” in _Proceedings of The 28th Conference on Learning Theory_ , ser. Proceedings of Machine Learning Research, P. Grünwald, E. Hazan, and S. Kale, Eds., vol. 40. Paris, France: PMLR, 03–06 Jul 2015, pp. 1066–1100. [Online]. Available: https://proceedings.mlr.press/v40/Kamath15.html
  * [8] H. S. Melihcan Erol, E. Sula, and L. Zheng, “On semi-supervised estimation of distributions,” in _2023 IEEE International Symposium on Information Theory (ISIT)_. IEEE, Jun. 2023.
  * [9] R. A. DeVore and G. G. Lorentz, _Constructive approximation_. Springer Science & Business Media, 1993, vol. 303.
  * [10] F. Chung and L. Lu, “Connected components in random graphs with given expected degree sequences,” _Annals of Combinatorics_ , vol. 6, no. 2, pp. 125–145, Nov. 2002. [Online]. Available: https://doi.org/10.1007/pl00012580

blabla

##  Appendix A Supplementary Results (Cont’d)

######  Theorem 8 ([9],Theorem-3.1)

If $f$ is bounded on A, differentiable in some neighborhood of $x$, and has
second derivative $f^{\prime\prime}(x)$ for some $x\in A$, then

$\lim_{n\rightarrow\infty}n\left[B_{n}(f,x)-f(x)\right]=\frac{x(1-x)}{2}f^{\prime\prime}(x).$

######  Theorem 9 ([10], lemma-2)

Let $X_{1},\ldots,X_{n}$ be independent random variables with

$\operatorname{Pr}\left(X_{i}=1\right)=p_{i},\quad\operatorname{Pr}\left(X_{i}=0\right)=1-p_{i}.$

We consider the sum $X=\sum_{i=1}^{n}X_{i}$, with expectation
$\mathrm{E}(X)=\sum_{i=1}^{n}p_{i}$. Then, we have:

$\displaystyle\operatorname{Pr}(X\leq\mathrm{E}(X)-\lambda)\leq
e^{-\lambda^{2}/2\mathrm{E}(X)}$
$\displaystyle\operatorname{Pr}(X\geq\mathrm{E}(X)+\lambda)\leq
e^{-\frac{\lambda^{2}}{2(\mathbb{E}(X)+\lambda/3)}}$

##  Appendix B Omitted Proofs

###  B-A Proof for Lemma 3

###### Proof:

Fix $c>0$, let $\delta_{1},\delta_{2}>0$,that will be determined later, we
define
$\Delta_{p}(x,i)\triangleq\left|x^{p}-\Big{(}\frac{i}{n}\Big{)}^{p}\right|$.
As a result of triangular inequality:

$\displaystyle\left|H^{n}_{p}(x)-G^{n}_{p}(x)\right|\leq\sum^{n}_{i=0}{n\choose
i}x^{i}(1-x)^{n-i}r^{p}_{i}\ \Delta_{p}(x,i)$ (49)

Before analyzing this sum, we note that by the mean value theorem, there
exists $\xi\in(x\wedge\frac{i}{n},x\vee\frac{i}{n})$ hence
$x^{p}-\left(\frac{i}{n}\right)^{p}=\left(x-\frac{i}{n}\right)\xi^{p-1}$,
thus:

$\displaystyle\Delta_{p}(x,i)=\Big{|}x^{p}-\Big{(}\frac{i}{n}\Big{)}^{p}\Big{|}\leq
p\Big{|}x-\frac{i}{n}\Big{|}\Big{|}x\vee\frac{i}{n}\Big{|}^{p-1}$ (50)

Now we analyze the sum in (49) over the ranges $i\leq nx$, $nx>i$ separately.
For the case $i\leq nx$:

$\displaystyle\sum_{i\leq nx}{n\choose i}x^{i}(1-x)^{n-i}r^{p}_{i}\
\Delta_{p}(x,i)$ $\displaystyle=$ $\displaystyle\sum_{i\leq
n(x-\delta_{1})}{n\choose i}x^{i}(1-x)^{n-i}r^{p}_{i}\Delta_{p}(x,i)$
$\displaystyle\ +$ $\displaystyle\sum_{n(x-\delta_{1})<i<nx}{n\choose
i}x^{i}(1-x)^{n-i}r^{p}_{i}\Delta_{p}(x,i)$ $\displaystyle\leq$
$\displaystyle\sum_{i\leq n(x-\delta_{1})}{n\choose i}x^{i}(1-x)^{n-i}2$
$\displaystyle\ +$ $\displaystyle\sum_{n(x-\delta_{1})\leq i\leq nx}{n\choose
i}x^{i}(1-x)^{n-i}r^{p}_{i}\ p\left|x-\frac{i}{n}\right|x^{p-1}$
$\displaystyle\leq$ $\displaystyle
2e^{-\frac{n}{x}\delta^{2}_{1}}+\sum_{n(x-\delta_{1})\leq i\leq nx}p{n\choose
i}x^{i}(1-x)^{n-i}r^{p}_{i}\delta_{1}x^{p-1}\quad$ $\displaystyle=$
$\displaystyle 2e^{-\frac{n}{x}\delta^{2}_{1}}+\frac{p\
\delta_{1}}{x}H^{n}_{p}(x)$ (53)

In (B-A) we use (50). In (B-A), we bound the lower tail of the binomial via
Theorem 9 and observe that $\left|x-\frac{i}{n}\right|\leq\delta_{1}$ in the
range $n\cdot(x-\delta_{1})\leq i\leq n\cdot x$. We also note that in (B-A),
$r^{p}_{i}\leq 2$ and $\Delta_{p}(x,i)\leq 1$. Now we choose
$\delta_{1}=c_{1}\sqrt{-\frac{x}{n}\log\left(\frac{1}{x}H^{n}_{p}(x)\right)}$
and obtain:

$\displaystyle\eqref{eq:bound_game}$ $\displaystyle\leq$ $\displaystyle 2\
\left(\frac{H^{n}_{p}(x)}{x}\right)^{c^{2}_{1}}+p\
c_{1}\sqrt{-\frac{1}{nx}\log\left(\frac{1}{x}H^{n}_{p}(x)\right)}H^{n}_{p}(x)$

Hence to establish the lemma, we first show that
$\sqrt{-\frac{1}{nx}\log{\frac{1}{x}H^{n}_{p}(x)}}=O(\frac{1}{\sqrt{\log{n}}})$.

$\displaystyle\frac{1}{x}H^{n}_{p}(x)$ $\displaystyle\simeq$ $\displaystyle
C_{p}\sum^{n}_{i=0}\frac{x^{p-1}}{i^{\frac{p}{2}}+1}{n\choose
i}x^{i}(1-x)^{n-i}$ (55) $\displaystyle\geq$ $\displaystyle
C_{p}\frac{1}{n^{\frac{p}{2}}+1}\left(\frac{c\log^{2}{n}}{n}\right)^{{}^{p-1}}\sum^{n}_{i=0}{n\choose
i}x^{i}(1-x)^{n-i}$ $\displaystyle=$ $\displaystyle
C_{p}\frac{1}{n^{\frac{p}{2}}+1}\left(\frac{c\log^{2}{n}}{n}\right)^{{}^{p-1}}\geq\frac{k^{\prime}}{n^{\frac{3}{2}p}}$

In (LABEL:eq:lower_bound), we notice $x\geq\frac{c\log^{2}{n}}{n}$. On the
right-hand side of (55) we collect the constants in $k^{\prime}$. Therefore we
establish that:

$\displaystyle\sqrt{-\frac{1}{nx}\log{\frac{1}{x}H^{n}_{p}(x)}}\leq\sqrt{k^{\prime\prime}\frac{1}{nx}\log{n}}\leq
O\left(\frac{1}{\log(n)}\right)\quad$ (56)

where in the first inequality we use that $x\geq c\frac{\log{n}}{n}$ and
collect the constants in $k^{\prime\prime}$. Secondly, we need to show that
$\frac{H^{n}_{p}(x)}{x}$ decays sufficiently fast. To this end, we let
$\delta^{\prime}$ be a parameter to be determined later:

$\displaystyle\frac{H^{n}_{p}(x)}{x}$
$\displaystyle=\sum_{i}\frac{C_{p}}{i^{\frac{p}{2}}+1}{n\choose
i}x^{i+p-1}(1-x)^{n-i}$ $\displaystyle\leq
C_{p}x^{p-1}\operatorname{\mathbb{E}}_{I\sim\operatorname{Bin}(n,x)}\left[{\frac{1}{\sqrt{I}+1}}\right]$
$\displaystyle\leq
C_{p}x^{p-1}\left(\operatorname{\mathbb{P}}\left({\left|I-nx\right|\geq\delta^{\prime}}\right)+\frac{1}{\sqrt{nx-\delta^{\prime}}+1}\right)$
$\displaystyle\leq C_{p}x^{p-1}\left(e^{-\frac{\delta^{\prime
2}}{n}}+\frac{1}{\sqrt{nx-\delta^{\prime}}+1}\right)$ hence setting
$\delta^{\prime}=x\cdot n^{0.75}$: $\displaystyle\leq
C_{p}x^{p-1}\left(e^{-\sqrt{n}x^{2}}+\frac{1}{\sqrt{n\cdot
x}+1}(1+o(1))\right)$ $\displaystyle\leq
C_{p}x^{p-1}\left(\frac{1}{\sqrt{n\cdot x}+1}(1+o(1))\right)$
$\displaystyle\leq C_{p}\frac{1}{\sqrt{n}}x^{p-\frac{1}{2}}(1+o(1))$
$\displaystyle=O\left(n^{-\frac{1}{2}}\right)$ (57)

(57) follows from that $p\geq 1$. Therefore by choosing $c_{1}$ large enough
we ensure that
$B\left(\frac{H^{n}_{p}(x)}{x}\right)^{c^{2}_{1}}=o(n^{-\frac{p}{2}})$.
Whereas in the second case $i>nx$:

$\displaystyle\IEEEeqnarraymulticol{3}{l}{\sum_{nx<i\leq
nx+\delta_{2}}{n\choose i}x^{i}(1-x)^{n-i}r^{p}_{i}\Delta_{p}(x,i)}$
$\displaystyle\quad+$ $\displaystyle\sum_{i>nx+\delta_{2}}{n\choose
i}x^{i}(1-x)^{n-i}r^{p}_{i}\Delta_{p}(x,i)$ $\displaystyle\leq$
$\displaystyle\sum_{nx<i\leq nx+\delta_{2}}{n\choose
i}x^{i}(1-x)^{n-i}\left|x-\frac{i}{n}\right|\left(\frac{i}{n}\right)^{p-1}$
$\displaystyle\quad+$ $\displaystyle 2\sum_{i>nx+\delta_{2}}{n\choose
i}x^{i}(1-x)^{n-i}$ $\displaystyle\leq$ $\displaystyle\delta_{2}\sum_{nx<i\leq
nx+\delta_{2}}{n\choose
i}x^{i}(1-x)^{n-i}\left(x+\delta_{2}\right)^{p-1}+2e^{-\frac{n\delta^{2}_{2}}{2(x+\frac{\delta_{2}}{3})}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq\delta_{2}2^{p-1}\sum_{nx<i\leq
nx+\delta_{2}}{n\choose
i}x^{i}(1-x)^{n-i}\left(x\vee\delta_{2}\right)^{p-1}+e^{-\frac{n\delta^{2}_{2}}{\frac{2}{3}(x\vee\delta_{2})}}}$
$\displaystyle\IEEEeqnarraymulticol{3}{l}{\leq\delta_{2}\
2^{p-1}\sum_{nx<i\leq nx+\delta_{2}}{n\choose
i}x^{i}(1-x)^{n-i}\left(x\right)^{p-1}+2e^{-\frac{n\delta^{2}_{2}}{\frac{2}{3}x}}}$
(60)

Each step is justified in the corresponding step in the analysis for the range
$i\leq nx$, except now we are using the upper tail in Theorem 9. In (60), we
see that the problem is identical to (53) except for the constants. Hence we
choose
$\delta_{2}=c_{2}\sqrt{-\frac{x}{n}\log\left(\frac{1}{x}H^{n}_{p}(x)\right)}$
and by (55) we have $\delta_{2}\leq c_{3}\sqrt{x\frac{\log{n}}{n}}$. This
ensures that $\delta_{2}\leq x$ when $x\geq c\frac{\log^{2}{n}}{n}$ and the
step (60) is justified. ∎