###### Abstract

The properties of decays that take place during jet formation cannot be easily
deduced from the final distribution of particles in a detector. In this work,
we first simulate a system of particles with well defined masses, decay
channels, and decay probabilities. This presents the ,,true system” for which
we want to reproduce the decay probability distributions. Assuming we only
have the data that this system produces in the detector, we decided to employ
an iterative method which uses a neural network as a classifier between events
produced in the detector by the ,,true system” and some arbitrary ,,test
system”. In the end, we compare the distributions obtained with the iterative
method to the ,,true” distributions.

###### keywords:

quantum chromodynamics, network model, data analysis, interpretability

xx 1 5 Received: date; Accepted: date; Published: date Introduction and
analysis of a method for the investigation of QCD-like tree data Marko Jercic
1*, Ivan Jercic 1 and Nikola Poljak1 Marko Jercic, Ivan Jercic, Nikola Poljak
Correspondence<EMAIL_ADDRESS>
##  1 Introduction

The properties of particle interactions determine the evolution of a quantum
chromodynamical (QCD) system. Thorough understanding of these properties can
help answer many fundamental questions in physics, such as the origin of the
Universe or the unification of forces. This is one of the important reasons to
collect data with particle accelerators, such as the Large Hadron Collider
(LHC) at CERN. However, when collecting this data, we only register complex
signals of high dimensionality which we can later interpret as signatures of
final particles in the detectors. This interpretation stems from the fact that
we, more or less, understand the underlying processes that produce the final
particles.

In essence, from the all the particles produced in a collision at the
accelerator, only the electron, the proton, the photon and the neutrinos are
stable and can be reconstructed with certainty, given that you have the proper
detector. Other particles are sometimes also directly detected, given that
they reach the active volume of the detector without first decaying. These
include muons, neutrons, charged pions and charged kaons. On the other hand,
short-lived particles will almost surely decay before reaching the detector
and we can only register the particles they decay into.

A similar situation arises with quarks, antiquarks and gluons, the building
blocks of colliding nuclei. When a high energy collision happens, a quark
within a nucleus behaves almost as if it doesn’t interact with neighbouring
particles, because of a property called asymptotic freedom. If it is struck
with a particle from the other nucleus, it can be given sufficient momentum
pointing outwards from the parent nucleus. However, we know that there are no
free quarks in nature and that this quark needs to undergo a process called
hadronisation. This is a process in which quark-antiquark pairs are generated
such that they form hadrons. Most of the hadrons are short-lived and they
decay into other, more stable, hadrons. The end result of this process is a
jet of particles whose average momentum points in the direction of the
original outgoing quark. Unfortunately, we don’t know the exact quark, nor
gluon, decay properties, which serves as a motivation for this work.

The determination of these properties is a long standing problem in particle
physics. To determine them, we turn to already produced data and try to fit
decay models onto them. With every new set of data our understanding changes.
This is evident from the fact that we want to simulate a collision event, we
can obtain, on average, slightly different results with different versions of
the same tool pythia . Therefore, even though simulation tools are regularly
reinforced with new observations from data, we can not expect the complete
physical truth from them.

Instead of trying to perform direct fits to data, we propose the use of
machine learning methods to determine the decay properties. In fact, the onset
of these methods is already hinted in the traditional approach, since a
multivariate fit of decay models to data is already a form of a machine
learning technique. It is only natural to extend the existing methods since we
can’t rely entirely on simulated data. In this work, we develop an
interpretable model by first simulating a system of particles with well
defined masses, decay channels, and decay probabilities. We take this to be
the ,,true system”, whose decay properties we pretend not to know and want to
reproduce. Mimicking the real world, we assume to only have the data that this
system produces in the detector. Next, we employ an iterative method which
uses a neural network as a classifier between events produced in the detector
by the ,,true system” and some arbitrary ,,test system”. In the end, we
compare the distributions obtained with the iterative method to the ,,true”
distributions.

This paper is organized as follows: in the Materials and methods section we
describe the developed artificial physical system and the algorithm used to
recover underlying probability distributions of the system. Also, we present
in detail the methodology used to obtain the presented results. In the Results
section we present our findings to see whether our hypothesis holds true. We
conclude the paper with the Discussion section…

##  2 Materials and methods

The code used for the development the particle generator, the neural network
models and the calculations is written in the the Python programming language
using the Keras module with the TensorFlow2 backend keras . The calculations
were performed using a standardized PC setup equipped with an NVIDIA Quadro
p6000 graphics processing unit.

###  2.1 The physical system

In particle physics, jets are detected as collimated streams of particles. The
jet production mechanism is in essence clear: partons from the initial hard
process undergo the fragmentation and hadronization processes. In this work,
we develop a simplified physical model in which the fragmentation process is
modeled as cascaded $1\rightarrow 2$ independent decays of partons with a
constant number of decays. This way, any single jet can be represented as a
perfect binary tree of depth $N$, corresponding to $2^{N}$ particles in the
final state. Since the initial parton properties are set, jets can be
described by $2^{N}-1$ decay parameters. We represent each decay of a mother
parton of mass $M$ by four real numbers
$(\frac{m_{1}}{M},\frac{m_{2}}{M},\theta,\phi)$, where $m_{1}$ and $m_{2}$ are
the masses of the daughter particles and $\theta$ and $\phi$ are the polar and
azimuthal angle of the lighter particle, as measured from the rest frame of
the mother particle. For simplicity we make all the decays isotropic, which
isn’t necessarily true in real processes. To fully define our physical system
we set a decay probability distribution function $p(m_{1},m_{2}|M)$, the
details of which are given in the following subsection. The aim of our
proposed algorithm is to recover these underlying probability distributions,
assuming we have no information on them, using only a dataset consisting of
jets described with final particles’ four-momenta, as one would get from a
detector.

###  2.2 Particle generator

To generate the jets, we developed an algorithm where we take a particle of
known mass that undergoes three successive decays. We consider only the
possibility of discrete decays, in the sense that the decay product masses and
decay probabilities are well defined. We consider a total of 10 types of
particles, labelled A – J, which can only decay into each other. The masses
and the decay probabilities of these particles are given in Table 1. In this
scenario, the ,,decay probabilities” $p$ are given by the ratios of decay
amplitudes. Thus, the total sum of the probabilities for a given particle to
decay into others has to be one, and the probabilities describe the number of
produced daughters per $N$ decays, scaled by $1/N$.

particle | A | B | C | D | E
---|---|---|---|---|---
mass | 0.1 | 0.6 | 1.3 | 1.9 | 4.4
$p$ / channel | 1 | A | 0.7 | B | 1 | C | 0.3 | A+C | 0.6 | C+C
|  |  | 0.3 | A+A |  |  | 0.3 | A+A | 0.4 | E
|  |  |  |  |  |  | 0.4 | D |  |
particle | F | G | H | I | J
mass | 6.1 | 8.4 | 14.2 | 18.1 | 25
$p$ / channel | 0.5 | A+A | 0.9 | B+B | 0.6 | D+D | 1 | F+G | 0.5 | F+I
| 0.5 | B+C | 0.1 | A+F | 0.25 | D+E |  |  | 0.4 | G+H
|  |  |  |  | 0.15 | E+F |  |  | 0.1 | E+E
Table 1: Allowed particle decays in the discrete model. The designation
$p$/channel shows the probability that a mother particle will decay into
specific daughters.

Particles A–E are set to be long lived and can thus be detected in the
detector, which only sees the decay products after several decays. This can be
seen in table 1 as a probability for a particle to decay into itself. In this
way, we assure two things: first, that we have stable particles and second,
that each decay in the binary tree is recorded, even if it is represented by a
particle decaying into itself. Particles A and C are completely stable, since
they only have one ,,decay” channel, in which they decay back into themselves.
On the other hand, particles F–I are hidden resonances: if one of them appears
in the i-th step of the decay chain, it will surely decay into other particles
in the next, (i+1)-th step of the chain.

To create a jet, we start with particle J, which we call the mother particle,
and allow it to decay in one of the decay channels. Each of the daughter
particles then decays according to their decay channels, and this procedure
repeats a total of 3 times. In the end, we obtain a maximum of 8 particles
from the set A–E, with known momenta as measured from the rest frame of the
mother particle. An example of a generated jet is given in Fig.1.

for tree= grow=east, edge=->, parent anchor=east, child anchor=west, s
sep=1pt, l sep=1cm , [J [F [A[A]] [A[A]] ] [I [F [B] [C] ] [G [B] [B] ] ] ]

Figure 1: An example of the operation of the discrete jet generator. The
mother particle J decays into particles I and F. According to decay
probabilities, this happens in half the generated jets. The daughter particles
subsequently decay two more times, leaving only stable, detectable particles
in the final state.

###  2.3 Introduction to the algorithm

Let’s assume we have two distinct datasets: one that consists of samples from
a random variable X distributed with an unknown probability density $p(x)$,
which we call the ,,real” dataset, and the other, which consists of samples
from a random variable Y distributed with a known probability density $q(x)$,
which we call the ,,test” dataset. We would like to perform a hypothesis test
between $H_{0}:p=p(x)$ and $H_{1}:p=q(x)$ using a likelihood-ratio test. The
approach we use follows earlier work employs the Neyman–Pearson lemma NNNP1 ;
NNNP2 ; NNNP3 . This lemma states that the likelihood ratio, $\Lambda$, given
by:

$\Lambda(p\mid q)\equiv\frac{{\mathcal{L}}(x\mid real)}{{\mathcal{L}}(x\mid
test)}=\frac{p(x)}{q(x)}$ (1)

is the most powerful test at the given significance level NeyPear .

We can obtain an approximate likelihood ratio $\Lambda$ by transforming the
output of a classifier used to discriminate between the two datasets. Assume
that the classifier is a neural network optimized by minimizing the
crossentropy loss. In this case, the network output gives the probability of
$x$ being a part of the real dataset $C_{nn}(x)=p(real\mid x)$ NNProbability .
If the datasets consist of the same number of samples, we can employ the
Bayes’ theorem in a simple manner:

$\displaystyle p(real\mid x)$ $\displaystyle=$ $\displaystyle\frac{p(x\mid
real)p(real)}{p(x\mid real)p(real)+p(x\mid test)p(test)}$ (2) $\displaystyle=$
$\displaystyle\frac{p(x\mid p_{\textrm{real}})}{p(x\mid real)+p(x\mid
test)}=\frac{\Lambda}{\Lambda+1}\,.$

A simple inversion of Eq.2 gives:

$\Lambda=\frac{p(x)}{q(x)}=\frac{C_{\textrm{NN}}(x)}{1-C_{\textrm{NN}}(x)},$
(3)

$p(x)=\frac{C_{\textrm{NN}}(x)}{1-C_{\textrm{NN}}(x)}q(x).$ (4)

Therefore, in ideal conditions, the unknown probability density $p(x)$
describing the real dataset can be recovered with the help of the known
probability density $q(x)$ and a classifier, using 4. It must be noted that 4
is strictly correct only for optimal classifiers, which are unattainable. In
our case, the classifier is optimized by minimizing the crossentropy loss
defined by:

$L=-\frac{1}{n}\sum_{i=1}^{n}\left[y(x_{i})\ln
C_{\textrm{NN}}(x_{i})+(1-y(x_{i}))\ln(1-C_{\textrm{NN}}(x_{i}))\right]\,,$
(5)

where $y(x_{i})$ is 1 if $x_{i}$ is a part of the real dataset, and 0 if
$x_{i}$ is a part of the test dataset. We can safely assume that the final
value of loss of the suboptimal classifier is greater than the final value of
loss of the optimal classifier:

$L_{\textrm{optimal}}<L<\ln{2}\,.$ (6)

The value of $\ln 2$ is obtained under the assumption of the worst possible
classifier. To prove our findings, in the next step we regroup the sums in the
loss function in two parts, corresponding to the real and the test
distributions:

$-\frac{1}{n}\sum_{i\in real}\ln
C_{\textrm{NN}}^{\textrm{optimal}}(x_{i})<-\frac{1}{n}\sum_{i\in real}\ln
C_{\textrm{NN}}(x_{i})<-\frac{1}{n}\sum_{i\in real}\ln\frac{1}{2},$ (7)

$-\frac{1}{n}\sum_{i\in
test}\ln\left[1-C_{\textrm{NN}}^{\textrm{optimal}}(x_{i})\right]<-\frac{1}{n}\sum_{i\in
test}\ln\left[1-C_{\textrm{NN}}(x_{i})\right]<-\frac{1}{n}\sum_{i\in
test}\ln\frac{1}{2}.$ (8)

After expanding inequality 7 we obtain:

$-\frac{1}{n}\sum_{i\in
real}\ln\left[\frac{C_{\textrm{NN}}^{\textrm{optimal}}(x_{i})}{1-C_{\textrm{NN}}^{\textrm{optimal}}(x_{i})}\right]<-\frac{1}{n}\sum_{i\in
real}\ln\left[\frac{C_{\textrm{NN}}(x_{i})}{1-C_{\textrm{NN}}(x_{i})}\right]<-\frac{1}{n}\sum_{i\in
real}\ln 1.$ (9)

According to Eq.4, we can recover the real probability density $p(x)$ when
using the optimal classifier. However, if one uses a suboptimal classifier, a
slightly different probability density $p^{\prime}(x)$ will be calculated.
Since the ratios that appear as arguments of the logarithms in Eq.9 correspond
to distribution ratios, it follows that:

$-\frac{1}{n}\sum_{i\in
real}\ln\left[\frac{p(x_{i})}{q(x_{i})}\right]<-\frac{1}{n}\sum_{i\in
real}\ln\left[\frac{p^{\prime}(x_{i})}{q(x_{i})}\right]<-\frac{1}{n}\sum_{i\in
real}\ln 1.$ (10)

After some simplification this becomes:

$\sum_{i\in real}\ln p(x_{i})>\sum_{i\in real}\ln p^{\prime}(x_{i})>\sum_{i\in
real}\ln q(x_{i}).$ (11)

If an analogous analysis is carried out for inequality 8 we get:

$\sum_{i\in test}\ln p(x_{i})<\sum_{i\in test}\ln p^{\prime}(x_{i})<\sum_{i\in
test}\ln q(x_{i}).$ (12)

From this, it can be seen that probability density $p^{\prime}(x)$ is on
average closer to the real probability density $p(x)$ than to the test
probability density $q(x)$. In a realistic case, Eq.4 can’t be used to
completely recover the real probability density $p(x)$. However, it can be
used in an iterative method; starting with a known distribution $q(x)$, we can
approach the real distribution more and more with each iteration step.

###  2.4 A simple example

Let us illustrate the recovery of an unknown probability density by using a
classifier on a simple example. We start with a set of 50 000 real numbers
generated from a random variable with a probability density given by

$p_{\textrm{real}}(x)=\frac{1}{4}\mathcal{N}(-1,1)+\frac{3}{4}\mathcal{N}(3,1)\,,$
(13)

where $\mathcal{N}(\mu,\sigma^{2})$ denotes a normal distribution. A histogram
of values in this set is shown in 2. Let’s now assume we don’t know
$p_{\textrm{real}}(x)$ and want to recover it using the procedure outlined in
the previous subsection. This set will be denoted as the ,,real” dateset and
the underlying probability density will be denoted as the ,,real” probability
density.

Figure 2: (a) The normalized probability density for the example, given by Eq.
13. (b) A histogram of values sampled from the set generated by the same
equation.

To construct the ,,test” dataset, we generate values with a uniform
probability density in the interval $\left[-10,10\right]$. Finally, we
construct a simple neural network which is used as a classifier that
distinguishes the examples from the real dataset from examples from the test
dataset. The classifier we use is a simple feed-forward neural network with
100 hidden units using a ReLU activation function. The activation function of
the final neural network output is the sigmoid function, which we use to
constrain the output values to the interval [0,1]. After the classifier is
trained to discriminate between the two datasets by minimizing the binary
crossentropy loss, we evaluate its output at 200 equidistant points between
-10 and 10. Using Eq.4, the probability distribution
$p_{\textrm{calculated}}(x)$ is calculated using the classifier outputs. The
calculated $p_{\textrm{calculated}}(x)$ is compared to the real probability
density $p_{\textrm{real}}(x)$ and is shown in 3.

Although the resulting probability density differs from the real probability
density due to the non-ideal classifier, we can conclude that the calculated
$p_{\textrm{calculated}}(x)$ is considerably closer to $p_{\textrm{real}}(x)$
than to uniform probability density $q(x)$ used to generate the test dataset.
Now, if we use the calculated $p_{\textrm{calculated}}(x)$ to construct a new
test dataset and repeat the same steps, we can improve the results even more.
This procedure can therefore iteratively improve the resemblance of
$p_{\textrm{calculated}}(x)$ to $p_{\textrm{real}}(x)$ to the point where the
datasets are so similar that the classifier cannot distinguish between them.
In this simple example convergence is reached after the 5th iteration, since
no significant improvement is observed afterwards. The calculated probability
density $p_{\textrm{calculated}}(x)$ after the final iteration is shown in 3
compared to the real distribution $p_{\textrm{real}}(x)$. It is clear that in
this case the procedure converges, and we could possibly obtain a better match
between $p_{\textrm{calculated}}(x)$ and $p_{\textrm{real}}(x)$ if we used a
more optimal classifier.

Figure 3: The calculated $p_{\textrm{calculated}}(x)$ (blue line) compared to
the real probability density $p_{\textrm{real}}$(x) (orange line). (a) The
left panel shows the comparison after one iteration of the algorithm,
alongside the starting ,,test” distribution (green line).(b) The right panel
shows the comparison after the 5th iteration.

In essence, a simple histogram could be used in this simple example to
determine the underlying probability distribution instead of using the method
described above. However, in case of multivariate probability distributions,
which can be products of unknown probability distributions, a histogram
approach would not prove useful.

###  2.5 Extension to jets

We would now like to apply the described procedure on the datasets that
contain jets. Every jet, represented by a binary tree of depth $N$, consists
of $2^{N}-1$ independent decays producing a maximum of $2^{N}$ particles in
the final state. Since all the decays are isotropic in space, a jet can be
described with a 4 $\times$ $(2^{N}-1)$–dimensional vector $\vec{x}$ and a
probability distribution function given by:

$p\left(\vec{x}\right)=\prod_{i}^{2^{N}-1}p(m_{1}^{i},m_{2}^{i}|M)p(\theta^{i})p(\phi^{i}),$
(14)

where $i$ denotes the decay index and ($m_{1}^{i}$, $m_{2}^{i}$, $\theta^{i}$,
$\phi_{i}$) are the components of the vector $\vec{x}$. Since both angles are
uniformly spatially distributed, they contribute to the probability with a
simple constant factor. Therefore, when plugging $p\left(\vec{x}\right)$ from
Eq.14 into Eq.4 we can omit angles, since the constant factors will cancel
each other out:

$\prod_{i}^{2^{N}-1}p(m_{1}^{i},m_{2}^{i}|M)=\frac{C_{NN}(\vec{x})}{1-C_{NN}(\vec{x})}\prod_{i}^{2^{N}-1}q(m_{1}^{i},m_{2}^{i}|M).$
(15)

Taking the logarithm of both sides:

$\sum_{i}^{2^{N}-1}\ln
p(m_{1}^{i},m_{2}^{i}|M)=\ln{C_{NN}(\vec{x})}-\ln({1-{C_{NN}(\vec{x})}})+\sum_{i}^{2^{N}-1}\ln
q(m_{1}^{i},m_{2}^{i}|M).$ (16)

Unfortunately, we cannot explicitly obtain the probability $p(m_{1},m_{2}\mid
M)$ directly from Eq.16 without solving a linear system of equations. This
task proves to be computationally exceptionally challenging due to the high
dimensionality of the dataset. In order to avoid this obstacle, we introduce a
neural network $f$ to approximate $\ln p(m_{1},m_{2}|M)$. We can optimize this
neural network by minimizing the mean squared error applied to the two sides
of Eq.16.

###  2.6 The 2 Neural Networks (2NN) algorithm

At this point we are ready to recover the underlying probability distributions
from an existing dataset that consists of jets described by the four-momenta
of the final particles. We denote the jets from this dataset as ,,real”. The
building blocks of the full recovery algorithm are two independent neural
networks; the aforementioned classifier $C_{NN}$ and the neural network $f$.
Based on the usage of 2 neural networks, we dubbed the algorithm 2NN. The
detailed architectures of both networks are given in Appendix A.

The workflow of the 2NN algorithm is simple: first we initialize the
parameters of both neural networks. Then, we generate a test dataset using the
neural network $f$. The test dataset and the real dataset are fed into the
classifier network, which produces a set of linear equations in the form of
Eq.16. We approximate the solution to these by fitting the neural network $f$,
which in turn produces a new test dataset. The procedure is continued
iteratively until there are no noticeable changes in the difference of the
real and test distributions. More detailed descriptions of the individual
steps are given in the next subsections.

####  2.6.1 Generating the test dataset

After the parameters of the neural network $f$ are initialized, we need to
generate a test dataset of jets with known decay probabilities $q(\vec{x})$.
The input of the neural network $f$ is a vector consisting of 3 real numbers:
$a=m_{1}/M$, $b=m_{2}/M$ and $M$. We denote the output of the neural network
with $f(a,b,M)$. Due to conservation laws, the sum $a+b$ needs to be strictly
less or equal to 1. We can assume $a\leq b$ without any loss of generality. In
order to manipulate with probabilities a partition function:

$Z(M)=\int_{\Omega}e^{f(a,b,M)}\,\mathrm{d}a\mathrm{d}b$ (17)

needs to be calculated. Here, $\Omega$ denotes the entire probability space
and is shown as the gray area in the left panel of Fig. 4. To calculate the
integral in the above expression, the probability space is discretized into
650 equal areas shown in the right panel of Fig.4. These areas are obtained by
discretizing the parameters $a$ and $b$ into equidistant segments of length
0.2. After the discretization, the partition function $Z(M)$ then becomes:

$Z(M)\approx\sum_{j}\sum_{k}e^{f(a_{j},b_{k},M)}\,.$ (18)

$a$$b$$a=b$$a+b=1$$\Omega$$a$$b$$a=b$$a+b=1$

Figure 4: (a) The left panel shows the entire allowed probability space of the
parameters $a$ and $b$, designated by $\Omega$. Due to conservation laws,
$a+b\leq 1$ needs to hold true. To describe our system, we selected the case
where $a\leq b$, which we can do without loss of generality. (b) The right
panel shows the discretized space $\Omega$, as used to evaluate the partition
function.

To generate the jets which form the test dataset, we must generate each decay
in the cascading evolution using the neural network $f$. Each of the decays is
generated by picking a particular pair of parameters $(a,b)$ from the 650
possible pairs which form the probability space for a given mass $M$. The
decay probability is then given by:

$q(m_{1},m_{2}\mid M)=\frac{e^{f(a,b,M)}}{Z(M)}\,.$ (19)

After applying this procedure we have a test dataset in which each jet is
represented as a list of $2^{N}$ particles and their four-momenta. For each
decay, we also store the pairs $(a^{i},b^{i})$ as well the corresponding decay
probabilities.

####  2.6.2 Optimizing the classifier

The classifier used in this work is a convolutional neural network. The input
to these type network are sets of images. For this purposes, all the jets are
preprocessed by transforming the list of particles’ four-momenta into jet
images. Two 32$\times$32 images are produces for a single jet. In both images
the axes correspond to the decay angles $\theta$ and $\phi$, while the pixel
values are either the energy of the momentum of the particle found in that
particular pixel. If a pixel contains two or more particles, their energy and
momenta are summed and stored as pixel values. The transformation of the jet
representations is done on both the real and the test datasets. We label the
,,real” jet images with the digit 1 and ,,test” jet images with the digit 0.
The classifier is then optimized by minimizing the binary crossentropy loss
between the real and the test datasets. The optimization is performed by ADAM
algorithm adam . It is important to note that the sizes of both datasets need
to be the same.

####  2.6.3 Optimizing the neural network $f$

After the classifier is optimized, a new jet dataset is generated by using the
neural network $f$. Just as earlier, the generated jets are first transformed
to jet images and then fed to the classifier. Since we have access to each of
the decay probabilities for each jet, the right side of Eq.16 can be easily
calculated for all the jet vectors $\vec{x}$ in the dataset. This way we can
obtain the desired log value of the total probability for each jet
$p(\vec{x})$:

$\ln
p(\vec{x})=\ln{C_{NN}(\vec{x})}-\ln({1-{C_{NN}(\vec{x})}})+\sum_{i}^{2^{N}-1}\ln
q(m_{1}^{i},m_{2}^{i}|M).$ (20)

Finally, we update the parameters of the neural network $f$ by minimizing the
expression given by:

$L=\frac{1}{n}\sum_{i}^{n}\left[\sum_{j}^{2^{N}-1}f(a_{i}^{j},b_{i}^{j},M_{j})-\ln
p_{i}(\vec{x})\right]^{2},$ (21)

where $i$ denotes the jet index and $j$ denotes the decay index in a
particular jet. After this step, the weights of the neural network are updated
in such a way that the network output values $f(a,b,M)$ are on average closer
to the real log value of $p(m_{1},m_{2}\mid M)$. The updated network $f$ is
used to generate the test dataset in the next iteration.

###  2.7 Evaluation of the 2NN algorithm

Upon completion of each iteration of the algorithm, the underlying probability
densities can be obtained from the output values of the neural network $f$
according to Eq.19. In the Results section the 2NN algorithm is evaluated in
terms of the Kullback-Leibler divergence (KL) in the following way KLD :

$KL(M)=\sum_{j}\sum_{k}p_{\textrm{real}}(m_{1}^{j},m_{2}^{k}\mid M)\left[\ln
p_{\textrm{real}}(m_{1}^{j},m_{2}^{k}\mid
M)-f(a^{j},b^{k},M)+\ln{Z(M)}\right]$ (22)

where the sum is performed over the whole probability space. The KL-divergence
is a non-negative measure of the difference between two probability densities
defined on same probability space. If the probability densities are identical,
the KL divergence is zero.

###  2.8 Hardware and software

Code for calculations in this reasearch are is written in Python programming
langauge using Tensorflow 2 and Numpy modules. GPU unit NVIDIA Quadro p6000
obtained from the NVIDIA Grant for academic resarch is used to increase the
speed of the performed calculations.

##  3 Results

In this section we present our findings after applying the $2NN$ algorithm on
500 000 jets created using the particle generator described in 2.2. In each
iteration, the classifier is optimized using 50 000 randomly picked jets from
the ,,real” dataset and 50 000 jets generated using the neural network $f$. To
optimize the neural network $f$, we use 50 000 jets as well. The algorithm
performed 800 iterations. After the final iteration of the $2NN$ algorithm we
obtain the calculated probability densities, which can be then used to
generate samples of jets. First, we show the energy spectrum of the particles
in the final state in jets generated by the calculated probabilities in 5.
This spectrum is directly compared to the energy spectrum of particles taken
from jets belonging to the ,,real” dataset and shown on Figure 5.

Figure 5: The energy spectrum of the particles in the final state in jets
generated by the calculated probabilities, compared to the energy spectrum of
particles taken from jets belonging to the ,,real” dataset.

The plotted spectra are obtained using 10 000 jets from each dataset. The
error bars in the histogram are smaller than the marker size and are hence not
visible. A resemblance between the two spectra is notable, especially at
higher energies. This points to the fact that the calculated probabilities are
approximately correct, so we can use them to generate samples of jets that
resemble ,,real” jets. To further examine the calculated probability densities
we need to reconstruct the hidden resonances which are not found in the final
state. For this purpose, the calculated probability densities for mother
particle masses of $M=25.0$, $M=18.1$, $M=14.2$ and $M=1.9$ are analyzed and
compared to the real probability densities in the following subsections. These
masses are chosen since they match the masses of the hidden resonances, as was
introduced in table 1.

###  3.1 Mother particle with mass $M$ = 25.0

The calculated 2$d$-probability density $p(m_{1},m_{2}\mid M)$ is shown in
Figure 6, compared to the real probability density. A simple eye reveals that
3 possible decays of particle of mass $M=25.0$ are recognized by the
algorithm. After dividing the probability space as in panel (c) in Figure 6
with lines $m_{2}>16.0$ and $m_{2}<10.0$, we calculate the mean and the
variance of the data on each subspace. As a result, we obtain
$(m_{1},m_{2})=(18.1\pm 0.5,6.1\pm 0.5)$ for $m_{2}>16.0$,
$(m_{1},m_{2})=(14.0\pm 0.7,8.4\pm 0.7)$ for $16.0\leq m_{2}>10.0$ and
$(m1,m2)=(4.8\pm 0.2,4.6\pm 0.2)$ for $m_{2}\leq 10.0$. These mean values
closely agree with the masses of the resonances expected as the products of
decays of the particle with mass $M=25.0$. The calculated small variances
indicate that the algorithm is very precise. The total decay probabilities for
each of the subspaces are equal to $p_{1}=0.48$, $p_{2}=0.47$, $p_{3}=0.05$,
which approximately agree with the probabilities of decay channels of the
particle with mass $M=25.0$, as defined in table 1.

(a)

(b)

(c)

Figure 6: The calculated probability density for a decaying particle of mass
$M=25.0$. (a) The left panel shows the density evaluated on the entire
discretized probability space. (b) The probability density of ,,real” data.
(c) A division of the probability space into three subspaces, in order to
isolate particular decays.

These results show that we can safely assume that the $2NN$ algorithm
successfully recognizes all the decay modes of the particle that initiates a
jet. To quantify the difference between the calculated probabilty density and
the real probability density, we use the KL-divergence.

Figure 7: The KL-divergence between the calculated and the real probability
densities, evaluated in the case of particle of mass $M=25.0$. The presented
results are averaged over 50-iteration intervals. The error bars represent the
standard deviation calculated on the same intervals.

Figure 7 shows the dependence of the KL-divergence on the iteration of the
$2NN$ algorithm. First, we observe an initial steep decrease in the value of
the divergence. Large variations in divergence value are observed later. This
is an indicator that the approximate probability density is found relatively
quickly - after a few hundred iterations. As the algorithm decreases the width
of the peaks found in the probability distribution, the KL-divergence becomes
very sensitive to small variations in the location of these peaks and can
therefore vary by a large relative amount.

###  3.2 Mother particle with mass $M$ = 18.1

A similar analysis is performed for the particle with mass $M=18.1$. The
calculated probability density is shown in Figure 8 compared to the expected
probability density. In this case, only one decay is allowed, so a division
into probability subspaces is not necessary, as was for the case when
$M$=25.0. The calculated mean and the variance of the shown probability
density are $(m_{1},m_{2})=(5.9\pm 0.4,8.2\pm 0.6)$. In this case, just as in
the former, the calculated values closely agree with the only possible decay,
in which the mother particle decays into two particles of masses 6.1 and 8.4.
Also, just as in the previous subsection, the obtained result is very precise.
Therefore, the algorithm can successfully find hidden resonances, as well as
recognize the decay channels, without ever seeing them in the final state in
the ,,real” dataset.

Figure 8: The calculated probability density for a decaying particle of mass
$M=18.1$. (a) The calculated density evaluated on the entire discretized
probability space. (b) The probability density of ,,real” data.

The calculated KL-divergence in the case of particle with mass $M=18.1$
decreases over time in a very smooth manner, as can be seen in Figure 9. We
believe this could be due to the simpler expected probability density, which
algorithm manages to find very quickly.

Figure 9: The KL-divergence between the calculated and the real probability
densities, evaluated in the case of particle of mass $M=18.1$. The presented
results are averaged over 50-iteration intervals. The error bars represent the
standard deviation calculated on the same intervals.

###  3.3 Mother particle with mass $M$ = 14.2

Figure 10 shows the 2$d$-probability density for the decaying particle of mass
$M=14.2$. In this case, we can identify 3 possible decay channels, which are
not as clearly separated as the channels in the previous subsections. Similar
to the case of decaying particle of mass $M=25.0$, we divided the probability
space into 3 subpaces, each of which covered one of the possible decays. In
this case, the three subspaces cover areas where $m_{2}\leq 4.0$,
$4.0<m_{2}\leq 5.5$ and $m_{2}>5.5$. The mean values of the probability
density on each of the subspaces are $(m_{1},m2)=(2.4\pm 0.5,2.9\pm 0.7)$,
$(m_{1},m_{2})=(2.7\pm 0.7,4.3\pm 0.3)$ and $(m_{1},m_{2})=(4.4\pm 0.4,6.2\pm
0.3)$, respectively. The allowed decays of a mother particle with mass $M$ =
14.2 in the ,,real” data are into channels with masses $(1.9,1.9)$,
$(1.9,4.4)$ and $(4.4,6.2)$, which agree with the calculated results. However,
in this case the calculations show higher variance, especially for decays
where one of the products is a particle with mass 1.9. The total probabilities
of decay in each of the subspaces are 0.89, 0.05 and 0.06, respectively. The
relative probabilities of decay channels into particles with masses (4.4, 6.1)
and (1.9, 4.4) are approximately the same as expected. However, the algorithm
predicts more decays in the channel (1.9,1.9) than expected. The KL-divergence
shows a steady decrease with occasional spikes, as shown on Figure 11.

(a)

(b)

(c)

Figure 10: The calculated probability density for a decaying particle of mass
$M=14.2$. (a) The left panel shows the density evaluated on the entire
discretized probability space. (b) The probability density of ,,real” data.
(c) A division of the probability space into three subspaces, in order to
isolate particular decays. Figure 11: KL-divergence between calculated and
real probability density evaluated for the $M=14.2$. Results are averaged over
the intervals of 50 iteration. Error bars represent standard deviation on the
same interval

###  3.4 Mother particle with mass $M$ = 1.9

The last probability density we analyze is the probability density for the
mother particle with mass $M$ = 1.9. Figure 12 shows the calculated
probability density. It can be seen that one of the decay modes present in the
,,real” data, namely when the particle decays in the $(0.1,0.1)$ channel, is
not recognized by the algorithm, but the decay mode when the particle decays
in the $(0.1,1.3)$ channel is visible. If we isolate given decay as shown in
the right panel of Figure 12, we get a mean value of $(m_{1},m_{2})=(0.14\pm
0.09,1.27\pm 0.09)$, which agrees with the expected decay. We also observe
significant decay probabilities along the line $m_{1}+m_{2}=1.9$. The decays
that correspond to the points on this line in effect create particles with
zero momentum in the rest frame of the mother particle. In the lab frame this
corresponds to the daughter particles flying off in the same direction as the
mother particle. Since they reach the detector in the same time, they are
registered as one particle of total mass $M=1.9$. Thus, we can conclude that
the probabilities on this line have to add up to the total probability of the
mother particle not decaying. The calculated probabilities in the case of no
decay and in the case when decaying into particles with masses $(0.1,1.3)$ are
0.71 and 0.29, respectively. We note that relative probabilities are not
correct, but 2 of the 3 decay modes are still recognized by the algorithm. The
KL-divergence in this case can’t produce reasonable results, simply because of
multiple points in the $(m_{1},m_{2})$ phase space which produce the same
decay and is therefore omitted from the analysis.

(a)

(b)

(c)

Figure 12: The calculated probability density for a decaying particle of mass
$M=1.9$. (a) The left panel shows the density evaluated on the entire
discretized probability space. (b) The probability density of ,,real” data.
(c) A division of the probability space into three subspaces, in order to
isolate particular decays.

###  3.5 The accuracy of the classifier

The accuracy of the classifier is defined as the fraction of correctly
,,guessed” samples on a given dataset. The criterion used for guessing is
checking whether the output of the classifier, $C_{NN}$, is greater than 0.5.
The accuracy can indirectly indicate how distinguishable are some two
datasets. In our algorithm, after starting from a test probability density, we
approach the real probabilility density with increasing iteration number, so
we can expect that the two jet datasets, the ,,real” and the ,,test” dataset,
are less and less distinguishable over time. In Figure 13 we show the accuracy
of the classifier in dependence on the iteration number.

Figure 13: The calculated accuracy of the classifier in dependence on the
iteration number.

After an initially high value, the accuracy decreases with growing iteration
number, which demonstrates that the test dataset becomes more and more similar
to the real dataset. Ideally, the datasets are no longer distinguishable by a
given classifier if the evaluated accuracy reaches 0.5. Therefore, we can use
the evaluated accuracy of the classifier as a criterion for stopping the
algorithm. Other measures can also be used as the stopping criterion such as
the loss value of the classifer or the area under reciever operating
characteristic (ROC) curve of the classifier. In this work, the algorithm is
stopped after the accuracy reaches a value of 0.65, because we didn’t see any
significant decreasy in the accuracy once it reached this value. An accuracy
value of 0.65 clearly shows that the classifier is capable of further
discriminating between the two datasets. This is explained by the fact that
the neural network $f$ and its hyperparameters are not fully optimized. For
the algorithm to perform better, we need to optimize the neural network $f$
and possibly improve the architecture for the selected task.

##  4 Discussion

In this work we propose a method for calculating underlying probability
distributions in particle decays, using only the data that can be collected in
a real-world physical system. First, we developed an artificial physical
system based on the QCD fragmentation process. Next, we present the core part
of the method: the $2NN$ algorithm, which we described in detail. The
algorithm performs very well when tested on the developed physical system. It
accurately predicts most of the hidden resonant particles, as well as their
decay channels, which can occur in the evolution of jets. The energy spectra
of the particles in the final state can also be accurately reproduced.

Although tested only on the developed artificial physical system, we believe
that the method is general enough to be applicable to real-world physical
systems, such as collisions of high-energy particles, with a few possible
modifications. For example, we hope that this method can in the future prove
helpful in measuring the fragmentation functions of quarks and gluons. Also,
one could employ such a method in the search for supersymmetric particles of
unknown masses, or in measuring the branching ratios of known decays.

The $2NN$ algorithm does not specify the exact architecture of the neural
networks, nor the representation of the data used. Furthermore, the classifier
does not need to be a neural network - it can be any machine learning
technique which maximizes likelihood. Although the algorithm has a Generative
Adversarial Netowrk (GAN)-like structure, it converges readily and does not
show usual issues associated with GANs, such as mode collapse or vanishing
gradients. The downside of the presented algorithm are high computational
requirements. Continuous probability distributions, which we expect to occur
in nature, are approximated by discrete probability distributions. In quest
for higher precision and a better description of reality, one always aims to
increase the resolution of discrete steps, but this carries a high
computational cost. Also, the used neural networks are not fully optimized,
which slows down the convergence of the algorithm. In conclusion, in order to
cut down computational costs, a more thorough analysis of convergence is
needed to achieve better performance.

In future work we hope to make the method even more general and thus even more
applicable to real-world physical systems. In particular, we want to introduce
angle dependent probability distributions, which can be retrieved from some
detector data. We would also like to investigate the possibility of including
other decay modes, such as $1\rightarrow 3$ type decays.

yes

##  Appendix A Description of the neural networks

The classifier used to recover the underlying probability distributions is a
feed forward convolutional neural network (CNN), already used in NNNP3 . In
that work, the classifier was used for a different purpose, but its
architecture proved to be suitable for the task at hand.

The neural network $f$ is a feed-forward neural network consisting of 5
independent completely connected layers. Each of the layers consist of 100
neurons apart from the last one, which is a single neuron. The activation
function used in all the layers, apart from the last one, is a ReLu activation
function. The last layer has no activation function. The network was optimized
using the ADAM algorithm.

The authors contributed equally to all segments of this paper.

This research was funded by the Croatian science foundation grant
IP-2018-01-4108 “Demystifying Two Particle Correlations in pp collisions with
the upgraded Time Projection Chamber”.

###### Acknowledgements.

We gratefully acknowledge the support of NVIDIA Corporation with the donation
of the Quadro P6000 graphics processing unit used for this research. The
authors declare no conflict of interest. The following abbreviations are used
in this manuscript:
QCD | Quantum Chromodynamics
---|---
LHC | Large Hadron Collider
CERN | Conseil Européen pour la Recherche Nucléaire
2NN | 2 Neural Networks
ROC | Reciever Operating Characteristic
GAN | Generative Adversarial Network
CNN | Convolutional Neural Network
References

## References

  * (1) Sjostrand T, Mrenna S, Skands P. PYTHIA 6.4 Physics and Manual. arXiv:hep-ph/0603175, 2006
  * (2) Chollet F et al. Keras. https://keras.io, 2015
  * (3) Streit R L.A neural network for optimum Neyman-Pearson classification. IJCNN International Joint Conference on Neural Networks, 1990, 1, 685–690
  * (4) Tong X, Feng Y, Li J J. Neyman-Pearson classification algorithms and NP receiver operating characteristics, Sci. Adv.,4-2, 2018
  * (5) Jercic M, Poljak N. Exploring the Possibility of a Recovery of Physics Process Properties from a Neural Network Model, Entropy,22(9), 994, 2020
  * (6) Neyman J, Pearson E S. On the problem of the most efficient tests of statistical hypotheses. Phil. Trans. R. Soc. Lond. A. 1933, 231, 694–706
  * (7) Bishop C M. Neural networks for pattern recognition. Oxford university press, 1995
  * (8) Kingma D, Ba J. A Method for Stochastic Optimization. International Conference on Learning Representations, 2014
  * (9) Kullback S, Leibler R. A. On Information and Sufficiency. Ann. Math. Statist. 22(1), 79-86, 1951


11institutetext: Insightface 22institutetext: Imperial College London
22email: {guojia, jackyu961127}<EMAIL_ADDRESS>
# Perspective Reconstruction of Human Faces by Joint Mesh and Landmark
Regression

Jia Guo 11 Jinke Yu 11 Alexandros Lattas 22 Jiankang Deng 1122

###### Abstract

Even though 3D face reconstruction has achieved impressive progress, most
orthogonal projection-based face reconstruction methods can not achieve
accurate and consistent reconstruction results when the face is very close to
the camera due to the distortion under the perspective projection. In this
paper, we propose to simultaneously reconstruct 3D face mesh in the world
space and predict 2D face landmarks on the image plane to address the problem
of perspective 3D face reconstruction. Based on the predicted 3D vertices and
2D landmarks, the 6DoF (6 Degrees of Freedom) face pose can be easily
estimated by the PnP solver to represent perspective projection. Our approach
achieves 1st place on the leader-board of the ECCV 2022 WCPA challenge and our
model is visually robust under different identities, expressions and poses.
The training code and models are released to facilitate future research.
https://github.com/deepinsight/insightface/tree/master/reconstruction/jmlr

###### Keywords:

Monocular 3D Face Reconstruction, Perspective Projection, Face Pose
Estimation.

##  1 Introduction

Monocular 3D face reconstruction has been widely applied in many fields such
as VR/AR applications (e.g., movies, sports, games), video editing, and
virtual avatars. Reconstructing human faces from monocular RGB data is a well-
explored field and most of the approaches can be categorized into
optimization-based methods [2, 3, 14] or regression-based methods [28, 13, 7,
11].

For the optimization-based methods, a prior of face shape and appearance [2,
3] is used. The pioneer work 3D Morphable Model (3DMM) [2] represents the face
shape and appearance in a PCA-based compact space and the fitting is then
based on the principle of analysis-by-synthesis. In [3], an in-the-wild
texture model is employed to greatly simplify the fitting procedure without
optimization on the illumination parameters. To model textures in high
fidelity, GANFit [14] harnesses Generative Adversarial Networks (GANs) to
train a very powerful generator of facial texture in UV space and constrains
the latent parameter by state-of-the-art face recognition model [8, 6, 9, 1,
30, 29]. However, the iterative optimization in these methods is not efficient
for real-time inference.

For the regression-based methods, a series of methods are based on synthetic
renderings of human faces [15, 26] or 3DMM fitted data [28] to perform a
supervised training of a regressor that predicts the latent representation of
a prior face model (e.g., 3DMM [24], GCN [27], CNN [25]) or 3D vertices in
different representation formats [19, 17, 13, 7]. Genova et al. [15] propose a
3DMM parameter regression technique that is based on synthetic renderings and
Tran et al. [24] directly regress 3DMM parameters using a CNN trained on
fitted 3DMM data. Zhu et al. [28] propose 3D Dense Face Alignment (3DDFA) by
taking advantage of Projected Normalized Coordinate Code (PNCC). In [27],
joint shape and texture auto-encoder using direct mesh convolutions is
proposed based on Graph Convolutional Network (GCN). In [25], CNN-based shape
and texture decoders are trained on unwrapped UV space for non-linear 3D
morphable face modelling. Jackson et al. [19] propose a model-free approach
that reconstructs a voxel-based representation of the human face. DenseReg
[17] regresses horizontal and vertical tessellation which is obtained by
unwrapping the template shape and transferring it to the image domain. PRN
[13] predicts a position map in the UV space of a template mesh. RetinaFace
[7] directly predicts projected vertices on the image plane within the face
detector. Besides these supervised regression methods, there are also self-
supervised regression approaches. Deng et al. [11] train a 3DMM parameter
regressor based on photometric reconstruction loss with skin attention masks,
a perception loss based on FaceNet [23], and multi-image consistency losses.
DECA [12] robustly produces a UV displacement map from a low-dimensional
latent representation.

Although the above studies have achieved good face reconstruction results,
existing methods mainly use orthogonal projection to simplify the projection
process of the face. When the face is close to the camera, the distortion
caused by perspective projection can not be ignored [20, 31]. Due to the use
of orthogonal projection, existing methods can not well explain the facial
distortion caused by perspective projection, leading to poor performance.

To this end, we propose a joint 3D mesh and 2D landmark regression method in
this paper. Monocular 3D face reconstruction usually includes two tasks: 3D
face geometric reconstruction and face pose estimation. However, we avoid
explicit face pose estimation as in RetinaFace [7] and DenseLandmark [26]. The
insight behind this strategy is that the evaluation metric for face
reconstruction is usually in the camera space and the regression error on 6DoF
parameters will repeat to every vertex. Even though dense vertex or landmark
regression seems redundant, the regression error of each point is individual.
Most important, explicit pose parameter regression can result in drift-
alignment problem [5], which will make the 2D visualization of face
reconstruction unsatisfying. By contrast, the 6 Degrees of Freedom (6DoF) face
pose can be easily estimated by the PnP solver based on our predicted 3D
vertices and 2D landmarks.

To summarize, the main contributions of this work are:

  * •

We propose a Joint Mesh and Landmark Regression (JMLR) method to reconstruct
3D face shape under perspective projection, and the 6DoF face pose can be
further estimated by PnP.

  * •

The proposed JMLR achieves first place on the leader-board of the ECCV 2022
WCPA challenge.

  * •

The visualization results show that the proposed JMLR is robust under
different identities, exaggerated expressions and extreme poses.

##  2 Our Method

###  2.1 3D Face Geometric Reconstruction

Figure 1: Straightforward solutions for perspective face reconstruction. (a)
direct 3D vertices regression, (b) two-stream 3D vertices regression and 6DoF
prediction, and (c) our proposed Joint Mesh and Landmark Regression (JMLR).

In Fig. 1, we show three straightforward solutions for perspective face
reconstruction. The simplest solution as shown in Fig. 1 (a) is to directly
regress the 3D vertices and 6DoF (i.e., Euler angles and translation vector)
from the original image (i.e., $800\times 800$), however, the performance of
this method is very limited as mentioned by the challenge organizer. Another
improved solution as illustrated in Fig. 1 (b) is that the face shape should
be predicted from the local facial region while the 6DoF can be obtained from
the global image. Therefore, the local facial region is cropped and resized
into $256\times 256$, and then this face patch is fed into the ResNet [18] to
predict the 1,220 vertices. To predict 6DOF information, the region outside
the face is blackened and the original $800\times 800$ image is then resized
into $256\times 256$ as the input of another ResNet.

In RetinaFace [7], explicit pose estimation is avoided by direct mesh
regression on the image plane as direct pose parameter regression can result
in misalignment under challenging scenarios. However, RetinaFace only
considered orthographic face reconstruction. In this paper, we slightly change
the regression target of RetinaFace, but still employ the insight behind
RetinaFace, that is, avoiding direct pose parameter regression. More
specifically, we directly regress 3D facial mesh in the world space as well as
projected 2D facial landmarks on the image plane as illustrated in Fig. 1 (c).

For 3D facial mesh regression in the world space, we predict a fixed number of
$N=1,220$ vertices
($\mathbf{V}=[v_{x}^{0},v_{y}^{0},v_{z}^{0};\cdots;v_{x}^{N-1},v_{y}^{N-1},v_{x}^{N-1}]$)
on a pre-defined topological triangle context (i.e., 2,304 triangles). These
corresponding vertices share the same semantic meaning across different faces.
With the fixed triangle topology, every pixel on the face can be indexed by
the barycentric coordinates and the triangle index, thus there exists a pixel-
wise correspondence with the 3D face. In [26], 703 dense landmarks covering
the entire head, including ears, eyes, and teeth, are proved to be effective
and efficient to encode facial identity and subtle expressions for 3D face
reconstruction. Therefore, 1K-level dense vertices/landmarks are a good
balance between accuracy and efficiency for face reconstruction.

As each 3D face is represented by concatenating its $N$ vertex coordinates, we
employ the following vertex loss to constrain the location of vertices:

$\displaystyle\mathcal{L}_{vert}{=}\frac{1}{N}\sum_{i=1}^{N}||v_{i}(x,y,z)-v^{*}_{i}(x,y,z)||_{1},\vspace{-4mm}$
(1)

where $N=1,220$ is the number of vertices, $v$ is the prediction of our model
and $v^{*}$ is the ground-truth. By taking advantage of the 3D triangulation
topology, we consider the edge length loss [7]:

$\displaystyle\mathcal{L}_{edge}{=}\frac{1}{3M}\sum_{i=1}^{M}||e_{i}-e^{*}_{i}||_{1},\vspace{-4mm}$
(2)

where $M=2,304$ is the number of triangles, $e$ is the edge length calculated
from the prediction and $e^{*}$ is the edge length calculated from the ground
truth. The edge graph is a fixed topology as shown in Fig. 1 (c).

For projected 2D landmark regression, we also employ distance loss to
constrain predicted landmarks close to the projected landmarks from ground-
truth:

$\displaystyle\mathcal{L}_{land}{=}\frac{1}{N}\sum_{i=1}^{N}||p_{i}(x,y)-p^{*}_{i}(x,y)||_{1},\vspace{-4mm}$
(3)

where $N=1,220$ is the number of vertices, $p$ is the prediction of our model
and $p^{*}$ is the ground-truth generated by perspective projection.

By combining the vertex loss and the edge loss for 3D mesh regression and the
landmark loss for 2D projected landmark regression, we define the following
prospective face reconstruction loss:

$\displaystyle\mathcal{L}{=}\mathcal{L}_{vert}+\lambda_{0}\mathcal{L}_{edge}+\lambda_{1}\mathcal{L}_{land},\vspace{-4mm}$
(4)

where $\lambda_{0}$ is set to $0.25$ and $\lambda_{1}$ is set to $2$ according
to our experimental experience.

###  2.2 6DoF Estimation

Based on the predicted 3D vertices in the world space, projected 2D landmarks
on the image plane, and the camera intrinsic parameters, we can easily employ
the Perspective-n-Point (PnP) algorithm [21] to compute the 6D facial pose
parameters (i.e., the rotation of roll, pitch, and yaw as well as the 3D
translation of the camera with respect to the world). Even though directly
regressing 6DoF pose parameters from a single image by CNN (Fig. 1 (a)) is
also feasible, it achieves much worse performance than our method due to the
non-linearity of the rotation space.

In perspective projection, 3D face shape $V_{world}$ is transformed from the
world coordinate system to the camera coordinate system by using 6DoF face
pose (i.e., the rotation matrix $R\in\mathbb{R}^{3\times 3}$ and the
translation vector $T\in\mathbb{R}^{1\times 3}$) with known intrinsic camera
parameters $K$,

$\displaystyle V_{camera}=K(V_{world}R+T),\vspace{-4mm}$ (5)

where $K$ is related to (1) the coordinates of the principal point (the
intersection of the optical axes with the image plane) and (2) the ratio
between the focal length and the size of the pixel. The intrinsic parameters
$K$ can be easily obtained through an off-line calibration step. Knowing 2D-3D
point correspondences as well as the intrinsic parameters, pose estimation is
straightforward by calling _cv.solvePnP()._

##  3 Experimental Results

###  3.1 Dataset

In the perspective face reconstruction challenge, 250 volunteers were invited
to record the training and test dataset. These volunteers sit in a random
environment, and the 3D acquisition equipment (i.e., iPhone 11) is fixed in
front of them, with a distance ranging from about 0.3 to 0.9 meters. Each
subject is asked to perform 33 specific expressions with two head movements
(from looking left to looking right / from looking up to looking down). The
triangle mesh and head pose information of the RGB image is obtained by the
built-in ARKit toolbox. Then, the original data are pre-processed to unify the
image size and camera intrinsic, as well as eliminate the shifting of the
principal point of the camera.

In this challenge, 200 subjects with 356,640 instances are used as the
training set, and 50 subjects with 90,545 instances are used as the test set.
Note that there is only one face in each image and the location of each face
is provided by the face-alignment toolkit [4]. The ground truth of 3D vertices
and pose transform matrix of the training set is provided. The unit of 3D
ground-truth is in meters. The facial triangle mesh is made up of 1,220
vertices and 2,304 triangles. The indices of 68 landmarks [10] from 1,220
vertices are also provided.

###  3.2 Evaluation Metrics

For each test image, the challenger should predict the 1,220 3D vertices in
world space (i.e., $V_{world}\in\mathbb{R}^{1220\times 3}$) and the pose
transform matrix (i.e., the rotation matrix $R\in\mathbb{R}^{3\times 3}$ and
the translation vector $T\in\mathbb{R}^{1\times 3}$) from world space to
camera space.

$V_{world}=\begin{bmatrix}v_{x}^{0}&v_{y}^{0}&v_{z}^{0}\\\
v_{x}^{1}&v_{y}^{1}&v_{z}^{1}\\\ \vdots&\vdots&\vdots\\\
v_{x}^{1219}&v_{y}^{1219}&v_{z}^{1219}\end{bmatrix}$,
$R=\begin{bmatrix}r_{00}&r_{01}&r_{02}\\\ r_{10}&r_{11}&r_{12}\\\
r_{20}&r_{21}&r_{22}\end{bmatrix}$,
$T=\begin{bmatrix}t_{x}&t_{y}&t_{z}\end{bmatrix}$.

Then, we can compute the transformed vertices in camera space by
$V_{camera}=V_{world}R+T$. In this challenge, the error of 3D vertices and
face pose is measured in camera space. Four transformed sets of vertices are
computed as follows:

$\displaystyle V_{1}$
$\displaystyle=V^{gt}R^{gt}+T^{gt},V_{2}=V^{pred}R^{pred}+T^{pred},$ (6)
$\displaystyle V_{3}$
$\displaystyle=V^{gt}R^{pred}+T^{pred},V_{4}=V^{pred}R^{gt}+T^{gt}.$

Finally, $L_{2}$ distance between pair
$(V_{1},V_{2}),(V_{1},V_{3}),(V_{1},V_{4})$ are calculated and combined into
the final distance error:

$L_{error}=\left\|V_{1}-V_{2}\right\|_{2}+\left\|V_{1}-V_{3}\right\|_{2}+10\left\|V_{1}-V_{4}\right\|_{2},$
(7)

where geometry accuracy across different identities and expressions is
emphasized by $\times 10$. On the challenge leader-board, the distance error
is multiplied by 1,000, thus the distance error is in millimeters instead of
meters.

###  3.3 Implementation Details

Input. We employ five facial landmarks predicted by RetinaFace [7] for
normalized face cropping at the resolution of $256\times 256$. Colour jitter
and flip data augmentation are also used during our training.

(a) ResNet 34

(b) Searched ResNet

Figure 2: Computation redistribution for the backbone. (a) the default
ResNet-34 design. (b) our searched ResNet structure for perspective face
reconstruction.

Backbone. We employ the ResNet [18] as our backbone. In addition, we refer to
SCRFD [16] to optimize the computation distribution over the backbone. For the
ResNet design, there are four stages (i.e., C2, C3, C4 and C5) operating at
progressively reduced resolution, with each stage consisting of a sequence of
identical blocks. For each stage $i$, the degrees of freedom include the
number of blocks $d_{i}$ (i.e., network depth) and the block width $w_{i}$
(i.e, number of channels). Therefore, the backbone search space has 8 degrees
of freedom as there are 4 stages and each stage $i$ has 2 parameters: the
number of blocks $d_{i}$ and block width $w_{i}$. Following RegNet [22], we
perform uniform sampling of $d_{i}\leq 24$ and $w_{i}\leq 512$ ($w_{i}$ is
divisible by 8). As state-of-the-art backbones have increasing widths, we also
constrain the search space, according to the principle of $w_{i+1}\geq w_{i}$.
The searched ResNet structure for perspective face reconstruction is
illustrated in Fig 2(b).

Loss optimization. We train the proposed joint mesh and landmark regression
method for 40 epochs with the SGD optimizer. We set the learning rate as $0.1$
and the batch size is set to $64\times 8$. The training is warmed up by 1,000
steps and then the Poly learning scheduler is used.

###  3.4 Ablation Study

In Table 1, we run several experiments to validate the effectiveness of the
proposed augmentation, network structure, and loss. We split the training data
into 80% subjects for training and 20% subjects for validation. From Table 1,
we have the following observations: (1) flipping augmentation can decrease the
reconstruction error by 0.45 as we check experiments 1 and 2, (2) the capacity
of ResNet-34 is well matched to the dataset as we compare the performance
between different network structures (e.g., ResNet-18, ResNet-34, and
ResNet-50), (3) searched ResNet with similar computation cost as ResNet-34 can
slightly decrease the reconstruction error by 0.18, (4) without the topology
regularization by the proposed edge loss, the reconstruction error increase by
0.29, confirming the effectiveness of mesh regression proposed in RetinaFace
[7], (5) by the model ensemble, the reconstruction error significantly reduces
by 0.92.

Table 1: Ablation study regarding augmentation, network structure, and loss. The results are reported on 40 subjects from the training set and the rest 160 subjects are used for training. Aug | Network | Loss | Flops | Score
---|---|---|---|---
Flip | ResNet34 | $\mathcal{L}$ | 5.12G | 34.31
No Flip | ResNet34 | $\mathcal{L}$ | 5.12G | 34.76
Flip | ResNet18 | $\mathcal{L}$ | 3.42G | 34.95
Flip | ResNet50 | $\mathcal{L}$ | 5.71G | 34.70
Flip | Searched ResNet | $\mathcal{L}$ | 5.13G | 34.13
Flip | ResNet34 |  $\mathcal{L}_{vert}$+$\mathcal{L}_{land}$ | 5.12G | 34.60
Flip | ResNet34 & Searched ResNet | $\mathcal{L}$ | 10.25 G | 33.39

###  3.5 Benchmark Results

For the submission of the challenge, we employ all of the 200 training
subjects and combine all the above-mentioned training tricks (e.g. flip
augmentation, searched ResNet, joint mesh and landmark regression, and model
ensemble). As shown in Table 2, we rank first on the challenge leader-board
and our method outperforms the runner-up by 0.07.

Table 2: Leaderboard results. The proposed JMLR achieves best overall score compared to other methods. Rank | Team | Score | $\left\|V_{1}-V_{2}\right\|_{2}$ | $\left\|V_{1}-V_{3}\right\|_{2}$ | $10\left\|V_{1}-V_{4}\right\|_{2}$
---|---|---|---|---|---
1st (JMLR) | EldenRing | 32.811318 | 8.596190 | 8.662822 | 15.552306
2nd | faceavata | 32.884849 | 8.833569 | 8.876343 | 15.174936
3rd | raccoon&bird | 33.841992 | 8.880053 | 8.937624 | 16.024315

###  3.6 Result Visualization

In this section, we display qualitative results for perspective face
reconstruction on the test dataset. As we care about the accuracy of geometric
face reconstruction and 6DoF estimation, we visualize the projected 3D
vertices on the image space both in the formats of vertex mesh and rendered
mesh. As shown in Fig. 3, the proposed method shows accurate mesh regression
results across 50 different test subjects (i.e., genders and ages). In Fig. 4
and Fig. 5, we select two subjects under 33 different expressions. The
proposed method demonstrates precise mesh regression results across different
expressions (i.e., blinking and mouth opening). In Fig. 6, Fig. 7, and Fig. 8,
we show the mesh regression results under extreme poses (e.g., large yaw and
pitch variations). The proposed method can easily handle profile cases as well
as large pitch angles. From these visualization results, we can see that our
method is effective under different identities, expressions and large poses.

Figure 3: Predicted meshes on 50 test subjects. 68 landmarks are also indexed
for better visualization. The proposed method shows stable and accurate mesh
regression results across different subjects (i.e., genders and ages). Figure
4: Predicted meshes of ID-6159 under 33 different expressions. 68 landmarks
are also indexed for better visualization. The proposed method shows stable
and accurate mesh regression results across different expressions (i.e.,
blinking and mouth opening). Figure 5: Predicted meshes of ID-239749 under 33
different expressions. 68 landmarks are also indexed for better visualization.
The proposed method shows stable and accurate mesh regression results across
different expressions (i.e., blinking and mouth opening). Figure 6: Predicted
meshes of different identities under different poses. 68 landmarks are also
indexed for better visualization. The proposed method shows stable and
accurate mesh regression results across different facial poses (i.e., yaw and
pitch). Figure 7: Predicted meshes of different identities under different
poses. 68 landmarks are also indexed for better visualization. The proposed
method shows stable and accurate mesh regression results across different
facial poses (i.e., yaw and pitch). Figure 8: Predicted meshes of different
identities under different poses. 68 landmarks are also indexed for better
visualization. The proposed method shows stable and accurate mesh regression
results across different facial poses (i.e., yaw and pitch).

##  4 Conclusion

In this paper, we explore 3D face reconstruction under perspective projection
from a single RGB image. We implement a straightforward algorithm, in which
joint face 3D mesh and 2D landmark regression are proposed for perspective 3D
face reconstruction. We avoid explicit 6DoF prediction but employ a PnP solver
for 6DoF face pose estimation given the predicted 3D vertices in the world
space and predicted 2D landmarks on the image space. Both quantitative and
qualitative experimental results demonstrate the effectiveness of our approach
for perspective 3D face reconstruction and 6DoF pose estimation. Our
submission to the ECCV 2022 WCPA challenge ranks first on the leader-board and
the training code and pre-trained models are released to facilitate future
research in this direction.

## References

  * [1] An, X., Deng, J., Guo, J., Feng, Z., Zhu, X., Yang, J., Liu, T.: Killing two birds with one stone: Efficient and robust training of face recognition cnns by partial fc. In: CVPR (2022)
  * [2] Blanz, V., Vetter, T.: A morphable model for the synthesis of 3D faces. In: SIGGRAPH (1999)
  * [3] Booth, J., Roussos, A., Ponniah, A., Dunaway, D., Zafeiriou, S.: Large scale 3d morphable models. IJCV (2018)
  * [4] Bulat, A., Tzimiropoulos, G.: How far are we from solving the 2d & 3d face alignment problem? (and a dataset of 230,000 3d facial landmarks). In: ICCV (2017)
  * [5] Chaudhuri, B., Vesdapunt, N., Wang, B.: Joint face detection and facial motion retargeting for multiple faces. In: CVPR (2019)
  * [6] Deng, J., Guo, J., Liu, T., Gong, M., Zafeiriou, S.: Sub-center arcface: Boosting face recognition by large-scale noisy web faces. In: ECCV (2020)
  * [7] Deng, J., Guo, J., Ververas, E., Kotsia, I., Zafeiriou, S.: Retinaface: Single-shot multi-level face localisation in the wild. In: CVPR (2020)
  * [8] Deng, J., Guo, J., Xue, N., Zafeiriou, S.: Arcface: Additive angular margin loss for deep face recognition. In: CVPR (2019)
  * [9] Deng, J., Guo, J., Yang, J., Lattas, A., Zafeiriou, S.: Variational prototype learning for deep face recognition. In: CVPR (2021)
  * [10] Deng, J., Roussos, A., Chrysos, G., Ververas, E., Kotsia, I., Shen, J., Zafeiriou, S.: The menpo benchmark for multi-pose 2d and 3d facial landmark localisation and tracking. IJCV (2019)
  * [11] Deng, Y., Yang, J., Xu, S., Chen, D., Jia, Y., Tong, X.: Accurate 3d face reconstruction with weakly-supervised learning: From single image to image set. In: CVPR Workshops (2019)
  * [12] Feng, Y., Feng, H., Black, M.J., Bolkart, T.: Learning an animatable detailed 3D face model from in-the-wild images. In: SIGGRAPH (2021)
  * [13] Feng, Y., Wu, F., Shao, X., Wang, Y., Zhou, X.: Joint 3D face reconstruction and dense alignment with position map regression network. In: ECCV (2018)
  * [14] Gecer, B., Ploumpis, S., Kotsia, I., Zafeiriou, S.: Ganfit: Generative adversarial network fitting for high fidelity 3d face reconstruction. In: CVPR (2019)
  * [15] Genova, K., Cole, F., Maschinot, A., Sarna, A., Vlasic, D., Freeman, W.T.: Unsupervised training for 3d morphable model regression. In: CVPR (2018)
  * [16] Guo, J., Deng, J., Lattas, A., Zafeiriou, S.: Sample and computation redistribution for efficient face detection. In: ICLR (2022)
  * [17] Güler, R.A., Trigeorgis, G., Antonakos, E., Snape, P., Zafeiriou, S., Kokkinos, I.: Densereg: Fully convolutional dense shape regression in-the-wild. In: CVPR (2017)
  * [18] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: CVPR (2016)
  * [19] Jackson, A.S., Bulat, A., Argyriou, V., Tzimiropoulos, G.: Large pose 3d face reconstruction from a single image via direct volumetric cnn regression. In: ICCV (2017)
  * [20] Kao, Y., Pan, B., Xu, M., Lyu, J., Zhu, X., Chang, Y., Li, X., Lei, Z., Qin, Z.: Single-image 3d face reconstruction under perspective projection. arXiv:2205.04126 (2022)
  * [21] Marchand, E., Uchiyama, H., Spindler, F.: Pose estimation for augmented reality: a hands-on survey. TVCG (2015)
  * [22] Radosavovic, I., Kosaraju, R.P., Girshick, R., He, K., Dollár, P.: Designing network design spaces. In: CVPR (2020)
  * [23] Schroff, F., Kalenichenko, D., Philbin, J.: Facenet: A unified embedding for face recognition and clustering. In: CVPR (2015)
  * [24] Tran, A.T., Hassner, T., Masi, I., Medioni, G.: Regressing robust and discriminative 3D morphable models with a very deep neural network. In: CVPR (2017)
  * [25] Tran, L., Liu, X.: Nonlinear 3d face morphable model. In: CVPR (2018)
  * [26] Wood, E., Baltrusaitis, T., Hewitt, C., Johnson, M., Shen, J., Milosavljevic, N., Wilde, D., Garbin, S., Sharp, T., Stojiljkovic, I., Cashman, T., Valentin, J.: 3d face reconstruction with dense landmarks. In: ECCV (2022)
  * [27] Zhou, Y., Deng, J., Kotsia, I., Zafeiriou, S.: Dense 3d face decoding over 2500fps: Joint texture & shape convolutional mesh decoders. In: CVPR (2019)
  * [28] Zhu, X., Lei, Z., Liu, X., Shi, H., Li, S.Z.: Face alignment across large poses: A 3D solution. In: CVPR (2016)
  * [29] Zhu, Z., Huang, G., Deng, J., Ye, Y., Huang, J., Chen, X., Zhu, J., Yang, T., Du, D., Lu, J., Zhou, J.: Webface260m: A benchmark for million-scale deep face recognition. TPAMI (2022)
  * [30] Zhu, Z., Huang, G., Deng, J., Ye, Y., Huang, J., Chen, X., Zhu, J., Yang, T., Lu, J., Du, D., Zhou, J.: Webface260m: A benchmark unveiling the power of million-scale deep face recognition. In: CVPR (2021)
  * [31] Zielonka, W., Bolkart, T., Thies, J.: Towards metrical reconstruction of human faces. In: ECCV (2022)