Pythagorean Centrality for Data Selection
Djemel Ziou

Département d’informatique

Université de Sherbrooke

Sherbrooke, QC, Canada J1K 2R1

<EMAIL_ADDRESS>
Abstract

This paper provides an overview of the Pythagorean centrality measures, which
are the arithmetic, geometric, and harmonic means. Both the evolution of their
meaning through history and their geometrical interpretation are outlined.
Relevant examples of use cases for each of them are introduced, spanning a
variety of areas of knowledge. Their differences and similarities are
explored. Finally, the issue of which mean to use in different situations in
order to make advantageous predictions is addressed.
Keywords: central tendency, Pythagorean means, arithmetic mean, geometric
mean, harmonic mean, data selection.

##  1 Introduction

The concept of the mean of a set of measurements was known in civilizations
before our era. Its meaning has evolved over time, according to needs. In
ancient India, Rituparna estimated the number of leaves on a tree based on a
single twig, which he multiplied by the estimated number of twigs on the
branches. This can be seen as an intuitive predecessor to the arithmetic mean,
as the leaf count of a single twig was taken to be representative of any other
twig [1]. During the civilization of the ancient Greeks, the Pythagoreans
found several mathematical formulations of the mean which can give different
values for the same two numbers [4]. The 9th-century mathematician Al-Kindi
introduced cryptography which used several statistical concepts including the
relative frequency and the arithmetic mean [5]. The 16th and 17th centuries
brought a shift in paradigm that led to the arithmetic mean being viewed as
the ”true” mean in the eyes of most people [4].

In our current knowledge, given the non-negative measurements
$x_{1}<\cdots<x_{n}$, the mean can be any value between $x_{1}$ and $x_{n}$
[6]. However for practical reasons, a requirement is mandatory in order to
target one value among all the values. The requirement is often expressed by
informal arguments, such as the mean value must greatly resemble other
measurements, or that it must be the closest to the others, or the one that
expresses a trend of the phenomenon, or the value that compensates for errors
such as the additive noise in a signal. No matter how the informal requirement
is expressed, the mean value is the closest to the measurements, up to some
transform, according to some criterion, e.g. the value $\mu$ minimizing the
quadratic error $\sum_{i=1}^{n}(T(x_{i})-O(\mu))^{2}$, where $T()$ and $O()$
are some transforms. One may think that in order to sidestep an arbitrary
choice of a value as a mean among all, it is enough to determine the value
that satisfies a predefined criterion. But it is not so, because there can be
an infinity of criteria. In other words, we are not escaping arbitrariness
altogether, but moving from a type of arbitrariness where the mean is in
$[x_{1},x_{n}]$ to an agreed arbitrariness where the mean value fulfils some
criterion chosen arbitrarily among many others. Despite the arbitrariness, the
mean is a concept that many among us claim to know and is the basis of
reasoning in various fields. It is widely used for purposes such as
determining the price evolution of consumer products, estimating the Human
Development Index, evaluating students learning, measuring the daily air
temperature in a city, and calculating the capacitance of several capacitors
connected in series. The goal of this paper is to examine Pythagorean
centrality measures: the arithmetic mean (AM), the geometric mean (GM), and
the harmonic mean (HM). For this purpose, we will describe these three means
as a finite set of non-negative measurements and explore the relationship
between them, their history, and their use. Moreover, we will propose a
geometric interpretation of the three means in terms of dimensions of
hyperrectangles as well as a statistical interpretation in terms of data
selection. The next section is devoted to the definition and the usefulness of
Pythagorean centrality measures. Their history and the comparison between them
are described in Section 3. In Section 4, we present the means as a mechanism
of data selection. Examples of use cases are in section 5.

##  2 Pythagorean centrality measures

The term Pythagorean centrality measures encompass the arithmetic, geometric,
and harmonic means. In all this study, we limit ourselves to the case of non-
negative measurements. Given $n$ scalar measurements $x_{1}<\cdots<x_{n}$, the
interpretation given our current knowledge would suggest that a Pythagorean
centrality measure is located between $x_{1}$ and $x_{n}$ [6] and is related
to compensation of errors, tendency, balance, and representativeness. These
qualifications can be considered synonymous because all can be couched as the
value $\mu$ minimizing a criterion based on the squared error; i.e.
$\sum_{i=1}^{n}w_{i}(T(x_{i})-O(\mu))^{2}$, where $T()$ and $O()$ are some
transforms. The non-negative weight $w_{i}$ of $x_{i}$ can be interpreted as
its relative frequency, the empirical probability of observing it, and its
relevance according to a certain criterion. In what follows, we will use
centrality, trend, and mean interchangeably and focus only on the mean of a
finite set of measurements. In reality, the geometric (resp. harmonic) mean is
also an arithmetic mean of measurements transformed by logarithmic (resp.
inverse) transform as shown in Table 2. In both these cases, the measurements
in the original space are transformed, the arithmetic mean is calculated in
the transformed space, and then the mean is transported back to the original
space. The squared-error-based criterion to be minimized is expressed in the
transformed space. As we will see later, this formalization of the Pythagorean
centrality measures in Table 2 makes it possible to explain one of the data
selection mechanisms.

| AM | GM | HM
---|---|---|---
Criterion |  $min_{a}\sum_{i=1}^{n}w_{i}(y_{i}-a)^{2}$, $w_{i}\geq 0$ $\forall i$, and $\sum_{i=1}^{n}w_{i}=1$
Data transform | y=x | $y=ln(x)$ | $y=1/x$
Parametrization | $a=\mu$ | $a=ln(\mu)$ | $a=1/\mu$
Mean ($\mu$) | $\sum_{i=1}^{n}w_{i}x_{i}$ | $\prod_{i=1}^{n}x_{i}^{w_{i}}$ | $1/(\sum_{i=1}^{n}w_{i}/x_{i})$
Table 2: Pythagorean centrality measures: the minimization problem by which
each is derived, the data transform, the parametrization used, and the mean.

Most humans default to the arithmetic mean when asked to find the average of a
sample. During an undergraduate class in computer science, 30 students were
asked to calculate the mean of three numbers. As expected, all calculated the
arithmetic mean and no one thought of using another measure of centrality.
There are countless examples of everyday use for the arithmetic mean, such as
the estimation of marks obtained by students, the daily temperature, the
position of a centre of mass [10], and the Consumer Product Index (CPI) [36].
It is easy to compute and understand, and, unlike the other two means, it does
not require a transformation of measurements. The arithmetic mean is suitable
for indicating the tendency when the dispersion of measurements is not high.
For the compensation of errors, a noisy measurement can be replaced by the
arithmetic mean estimated from the closest measurements in the case of
uncorrelated and additive noise [31]. However, it is sensitive to extreme
values. For example, let’s assume that the five equally-weighted measurements
of the daily temperature reported are 8, 13, 14, 10, and 1000 degrees.
Although 1000 degrees is clearly a mistake, it influences the average daily
temperature, resulting in an unrealistic value: 209 degrees. This sensitivity
of the arithmetic mean may explain why the geometric mean has regained
popularity in many areas of knowledge. In physics, it expresses the
relationship between the classical and the relativistic Doppler effect [13].
In image processing, a geometric filter based on the geometric mean is used
for the reduction of speckle in radar images [16]. As a filter, the geometric
mean weakens the contribution of high values heavily corrupted by
multiplicative noise. In finance, the geometric mean is used to characterize
the growth rates of an investment [17]. It is also used to assess the
effectiveness of a vaccine [11]. The CPI is estimated using the geometric mean
for certain products, as will be discussed below. The harmonic mean is used
for change detection in radar images [20], estimation of the average velocity
over a trip [21], averaging the financial multiples with price in the
numerator [23], estimation of the pollutant load that can be permitted in a
water quality limited stream [18], and estimation of the properties of a
heterogeneous system of porous media [18]. In the case of multiples, the
harmonic mean can provide more accurate estimates than the arithmetic mean.
This surprising result has been addressed several times in litigation in the
United States [19]. The geometric and harmonic means favour small values and
are therefore less sensitive than the arithmetic mean to larger outliers. In
the example presented above the arithmetic mean is 209 degrees while the
geometric and harmonic means are around 27.02 and 13.36 degrees, respectively.
While the former is far from the median, which is 13 degrees, the latter is
close to it.

Figure 1: The three means of 1 and x.

Finally, the three means can be ranked. To illustrate this, the three means of
1 and $x$ are plotted in Fig. 1 as a function of $x\in[2,15]$. It can be
observed that the means are contained in the interval $[0,x]$ and AM $>$ GM
$>$ HM. The inequality is valid for any finite number of measurements. The
reader can find the formal proof in [24].

##  3 Evolution of the Pythagorean centrality measures

Note that historical clarification is required to summarize what was defined
by Pythagoras, his disciples and what has been introduced afterwards. For
example, in [9], it is stated that the geometric mean was introduced at the
beginning of the 19th century. For this purpose, the history of the three
means and the current geometrical interpretation are the focus of this
section.

###  3.1 History

The arithmetic, geometric, and harmonic means of a sample appeared during the
last centuries BCE. The Babylonians used the arithmetic mean in astronomy
[25], classical Greek music theorists employs the geometric mean in tuning
instruments [26], and the harmonic mean is the basis of a solution of a
division of loaves of bread by 10 [28]. However, it was in ancient Greece that
the study of the mean figured prominently because Greek mathematics was
primarily concerned with proportions, as is the case with ratios and the
Pythagorean theorem. The remarkable contributions of Pythagoras in the study
of proportions probably led to the three measures of centrality attributed to
Pythagoras and his disciples, i.e. Pythagorean means [22]. Much later, Al-
Kindi introduced cryptography in the 9th century, using several statistical
concepts, including relative frequency and arithmetic mean [5]. In the 11th
century, Al-Biruni used the midrange as a mean in astronomy [4]. Until the end
of the Middle Ages, the notion of mean had several interpretations. One of
these was the midrange [29], which is the value located equidistant from the
two extreme values. This is not easily generalizable to the case of more than
two measurements, which is why the definition of the arithmetic mean of a
sample as we know it today was not proposed until the 16th century [4]. There
is no consensus on the origin of the arithmetic mean of a population (i.e. the
mathematical expectation). Some date it to the 17th century, attributing it to
Blaise Pascal. Regardless, it was only formalized in the 19th century, by
Laplace [30]. Also in the Middle Ages, the harmonic mean was used to solve
what is commonly called the cistern mathematical problem that appeared in
Europe and India [27]. Details regarding the development of the harmonic mean
after this time period seem to be lacking in the literature. The geometric
mean as we know it was not proposed until the beginning of the 19th century
[9]. Surprisingly, the population geometric mean was only introduced in 2013
[9]. Note that the geometric and arithmetic means have been combined to form
the arithmetic-geometric mean used in the calculation of integrals [32, 33].
As its name suggests, this average consists of iteratively calculating the
arithmetic and geometric mean of the measurements.

###  3.2 Geometric interpretation

The geometric interpretation is only provided for equal weights measurements.
Let us consider a sample of only two values, $x_{1}$ and $x_{2}$ with equal
weights. The three means can be represented geometrically using a circle of
diameter $x_{1}+x_{2}$, as depicted in Fig. 2.

Figure 2: Representation of the three means of two numbers.

The lengths of BG and GC are $x_{1}$ and $x_{2}$, respectively. The arithmetic
mean coincides with the radius of the circle (OH), HG corresponds to the
geometric mean, and HD to the harmonic mean [2]. This relationship among the
three means is only applicable to two numbers. We suggest an alternative
geometric interpretation that can be applied to more than two numbers $(n\geq
2)$. Let us consider the hyperrectangle $R_{n}$, where $x_{1}$, …, $x_{n}$ are
its dimensions. In $R_{n}$, there is a total of $2^{(n-1)}n$ edges. As a
result, each edge is repeated $2^{(n-1)}$ times. For instance, in $R_{3}$,
each one of the 3 edges would be repeated 4 times. The hyperperimeter of
$R_{n}$ is given by:

$P_{n}=2^{(n-1)}(x_{1}+x_{2}+...+x_{n})$ (1)

According to this formula, the hyperperimeter is the arithmetic mean
multiplied by a factor whose value is contingent on the number of dimensions:

$AM=\frac{P_{n}}{2^{(n-1)}n}$ (2)

The hypervolume $V_{n}=x_{1}x_{2}...x_{n}$ of $R_{n}$ is the geometric mean
raised to the nth power:

$GM=\sqrt[n]{V_{n}}$ (3)

The hyperrectangle $R_{n}$ is made up of $n$ pairs of equal
$(n-1)$-dimensional hyperrectangles. For example, in $R_{3}$, there are 6
rectangles that are equal to one another two by two. Let us consider the
$j^{th}$ hyperrectangle $R_{n-1,j}$ with hypervolume
$V_{n-1,j}=\prod\limits_{i\neq j}^{n}x_{i}$. The arithmetic mean of the
hypervolumes of the $2n$ hyperrectangles $R_{n-1,j}$ is:

$A_{n}=\frac{1}{2n}\sum^{2n}_{j=1}V_{n-1,j}=\frac{1}{n}\sum^{n}_{j=1}V_{n-1,j}$
(4)

The harmonic mean is expressed in terms of a ratio between the hypervolume of
$R_{n}$ and the arithmetic mean of the hypervolumes of these n $R_{n-1}$:

$H_{n}=\frac{V_{n}}{A_{n}}$ (5)

In other terms, this is the ratio between the geometric mean of edges to the
power of $n$ and the arithmetic mean $A_{n}$ of the hypervolumes of $n$
hyperrectangles of dimension $n-1$. To conclude, the geometric interpretation
of the three Pythagorean centrality measures leads to translate the inequality
between them into an inequality between geometric primitives of a
hyperrectangle.

##  4 Data selection

Figure 3: Weighted Pythagorean centrality measures of 1 and $x$ as a function
of $x$ and the weight $w$ of $x$. The upper row shows the arithmetic,
geometric, and harmonic means respectively. The bottom row shows the
velocities of the means when $x=0.1$ (left) and $x=0.7$ (right).

Remember that the non-negative weight $w_{i}$ of $x_{i}$ encodes its relative
frequency, the empirical probability of observing it, or its relevance
according to a certain criterion. Let us call it, the w-weight. When
estimating the mean, the w-weights implement a selection mechanism for
measurements, since they indicate how each of them contributes to the
estimated mean. Note that the selection of a measurement is not taken in the
sense of all-or-nothing. Assigning a weight to a measurement is a selection
since the measurement’s contribution is strengthened or weakened depending on
this weight [12, 3, 34]. Fig. 3 displays the three weighted means of two
numbers 1 and $x$, by varying $x\in[0,1[$ and the weight $w\in]0,1]$ of $x$.
For the three weighted means, the more $w$ increases, the more strongly they
are attracted by $x$. The difference between the three means is explained by
the velocity at which each move towards $x$ as $w$ increases. For a given $w$,
the velocity of a mean can be set to the square of the difference between the
mean and $x$. Fig. 3 presents the velocities of the three means for $x=0.1$
and $x=0.7$. For a small $x$, the harmonic mean approaches $x$ more quickly
than the geometric mean which in turn approaches $x$ more quickly than the
arithmetic mean. This explains why the harmonic mean is the most sensitive to
small values. When $x$ increases, the difference between the velocities of the
three means is reduced.

We propose to make explicit another mechanism for selecting measurements which
does not seem to have been made explicit before. The mechanism is expressed in
the original space (the measurement space). For the sake of simplicity, first
we consider that all measurements have the same w-weight. Let us return to the
measurements $x_{1}<\cdots<x_{n}$. A centrality measure can be seen as a value
belonging to the interval $[x_{1},x_{n}]$, where each of the measurements
contributes to its identification. The closer a measure is to it, the more
influence it has on it to keep it close. This conception implicitly
establishes an analogy with the Coulomb’s law; i.e. both $x$ and the mean are
seen as electric charges whose force of attraction decreases with the square
of the distance between them. This analogy leads to basing attraction on the
Euclidean distance $d(x,\mu)$ between a measurement $x$ and a Pythagorean
centrality measure $\mu$ of the measurements. Since the influence of $x$ on
$\mu$ decreases as $d(x,\mu)$ increases, an attraction function can be any
decreasing function on $d(x,\mu)$, such as a Cauchy or Gaussian pdfs. Fig. 4
presents an example of the attraction function, which we call the Cauchy
attraction function $f_{p}(x)=\frac{d^{2}(x_{n},x_{1})}{d^{2}(x,\mu_{p})+1}$,
where $d^{2}(a,b)=(a-b)^{2}$ and $p~{}\in$ {arithmetic, geometric, harmonic}.
The squared distance $d^{2}(x_{n},x_{1})$ is the range of measurements used as
the normalization constant. In the denominator, one is added to $d^{2}()$ in
order to avoid division by zero.

Figure 4: The attraction functions over the three central tendencies obtained
by regular sampling of $[0.1,10]$.

Fig. 4 depicts the attraction functions obtained by regular sampling of
$[0.1,10]$. The arithmetic, geometric, and harmonic means are around 5.05,
3.85, and 2.13 respectively. It can be noted that the measurements close to
the means have the greatest attraction, regardless of the means used. The
attraction of high values is higher in the case of the arithmetic mean. For
example, the attraction function of the arithmetic mean at $x=10$ is greater
than that of the geometric mean, which in turn is greater than that of the
harmonic mean. Conversely, at $x=0.1$ the attraction function of the harmonic
mean is greater than that of the geometric mean which in turn is greater than
that of the arithmetic mean. Let us examine how the two selection mechanisms
(w-weights and attraction function) can be combined. The combination is
straightforward because the weighted attraction function can be either Cauchy
$f_{p}(x_{i})\propto 1/(w_{i}d^{2}(x_{i},\mu_{p})+1)$ or Gaussian with a non-
constant variance $g_{p}(x_{i})\propto exp(-w_{i}(x_{i}-\mu_{p})^{2})$. Fig. 5
shows the student grades and the Cauchy weighted attraction functions for the
three means. The observations deduced from Fig. 4 remain valid. Moreover, the
data are noisy which explains the difference between the weighted attraction
functions at the maximum points.

Figure 5: Frequencies of student grades and the three weighted attraction
functions over the three central tendencies.

The third data selection mechanism is implemented by performing
transformations to calculate geometric and harmonic means (see Table 2).
Indeed, these two transformations modify the measurements in such a way as to
enhance some and weaken the others. For example, the logarithm accentuates the
difference between small values and attenuates the difference between large
values. Someone may think that the inverse transform neutralizes the effects
of the measurements transform. However, the back-transformed arithmetic mean
calculated in the transformed space does not cancel the transformation of the
measurements because the two transformations are not linear.

To conclude, a weighted mean implements three data selection mechanisms. One
of the mechanisms is based on a frequency or any priori knowledge and the two
other on the measurement itself. Thus, choosing which Pythagorean centrality
measure to use defines how the measurements are considered and which ones are
relevant. It follows that the right choice of a mean can be the cause of
appropriate subsequent use of measurements such as accurate predictive
analysis and fair decision-making.

##  5 Applications

###  5.1 Predictions

Let’s take the example of an airline that needs to estimate the optimal number
of additional tickets to sell [35]. These seats are those purchased by
passengers, called no-shows, who do not show up for a flight. For estimation
purposes, the number of no-shows must be predicted as accurately as possible,
because prediction errors lead to financial losses for the airline. A free
seat is a loss of earnings and an oversold seat means that a person will be
compensated for not boarding. To make the prediction, the airline uses the
number of no-shows per flight over the last 200 flights. These data are
indicated in Table 3. For each number of no-shows, the table gives the number
of flights and the relative frequency of these flights, considered as an
empirical probability of the flights.

No-shows | 1 | 2 | 3 | 4 | 5 | 6 | Total
---|---|---|---|---|---|---|---
Number of flights | 70 | 40 | 10 | 20 | 20 | 40 | 200
Probability | 0.35 | 0.20 | 0.05 | 0.10 | 0.10 | 0.20 | 1
Table 3: Number of no-shows on the last 200 flights, the number of flights per
number of no-shows, and the empirical probability of the flights.

The prediction is made by a consultant whose compensation is linked to the
accuracy of the predictions. If $x$ is the expected number of ”no-shows” and
$n$ the actual number, the compensation is given by the following gain
function expressed in a transformed space:

$g(x,n)=1000-30(T(x)-T(n))^{2}$ (6)

As shown in the table 2, $T(a)=a$ gives the arithmetic mean, $T(a)=ln(a)$ the
geometric mean, and $T(a)=1/a$ the harmonic mean. Because $n$ is unknown, a
way of finding the prediction $x$ is to use the expectation of the gain
function, named the return function.

$R(x)=1000-\sum\limits^{6}_{n=1}30(T(x)-T(n))^{2}pr(n)$ (7)

Because $T()$ is invertible, maximizing $R(x)$ w.r.t $x$ allows to find the
best predictor:

$x^{*}=T^{-1}(\sum\limits^{6}_{n=1}T(n)\text{pr}(n))$ (8)

Table 4 provides the three return functions, their maximum values, and the
best predictors. For the consultant, the harmonic mean is the best scenario,
since it allows her or him to earn more. In the absence of ground truth and
other economic data, it is difficult to choose the best predictor for the
airline. However, the airline could object to the use of a return function and
therefore a central tendency measure.

Return function $R(x)$ | Best predictor $x^{*}$ | Best return $R(x^{*})$
---|---|---
$1000-30\sum_{n=1}^{6}(x-n)^{2}pr(n)$ | AM = 3.00 | 883.00
$1000-30\sum_{n=1}^{6}(ln(x)-ln(n))^{2}pr(n)$ | GM=2.34 | 984.26
$1000-30\sum_{n=1}^{6}(1/x-1/n)^{2}pr(n)$ | HM=1.83 | 996.27
Table 4: Return functions expressed in the original and transformed spaces of
data, their maximum values, and their best predictors

###  5.2 Consumer Price Index

We briefly describe the case of Canada. The CPI represents price change by
comparing, over time, the cost of the same basket of goods and services. For
decades, it was calculated using different formulas, including the arithmetic
mean, but since 1995, for many products, the geometric mean has been used
[36]. The arithmetic mean is still used for three out of 691 products (rents,
tuition, and passenger vehicle insurance premiums) on the grounds that they
are unlikely to have outliers. It is also used to group products into broader
categories, such as food, housing and transport and to calculate the all-items
CPI. The idea of using the geometric mean instead of the arithmetic mean to
estimate CPI was put forward by economist W. S. Jevons in 1863 [7]. Later,
however, he proposed the use of the harmonic mean instead. The reason for this
choice was never explicitly stated in his published work [8]. Let’s examine
the effect of Pythagorean centrality on the CPI using Canadian data acquired
during the two years 2002 and 2017. Note that, the Canadian government tracks
a wide variety of products and services. Based on consumer habits, weights are
assigned to products indicating their importance to consumers. Table 5
summarizes the weights of the individual product groups and their CPIs in
2017, taking 2002 as the reference period. Table 6 shows the overall CPI as a
function of Pythagorean centrality measures. As expected, the highest CPI is
given by the arithmetic mean and the lowest by the harmonic mean. The latter
is 1.3% lower than the CPI based on the arithmetic mean. It seems that there
are not yet formal arguments to choose the central tendencies to use for the
calculation of the CPI. This could open the door to arbitrariness and
arguments of a political nature. Nevertheless, a good understanding of data
selection could guide the choice of the appropriate central tendencies.

Product group | Weight | CPI
---|---|---
Food | 0.1648 | 141.5
Shelter | 0.2736 | 137.8
Household operations, | 0.1280 | 121.4
furnishings, and equipment |  |
Clothing and footwear | 0.0517 | 91.1
Transportation | 0.1995 | 133.0
Health and personal care | 0.0479 | 123.4
Recreation, education, and reading | 0.1024 | 111.3
Alcohol, beverages, tobacco, | 0.0321 | 158.7
products, and cannabis |  |
Table 5: Weights and CPI associated with each basket category as of January 2017 [37, 38]. Arithmetic mean | Geometric mean | Harmonic mean
---|---|---
130.20 | 129.40 | 128.50
Table 6: All-items CPI given by the three Pythagorean means.

###  5.3 Ellipse fitting

Fitting ellipses to given points is a problem that arises in computer
graphics, remote sensing, geology, environment, and statistics, among many
others. Different ideas have been implemented, including least-squares error
minimizing and probability density estimation [39, 14, 15]. For 2D scatter
points, we want to estimate an enclosing ellipse. An ellipse is defined by its
centre and its principal axes. The centre can be estimated by using a
Pythagorean centrality measure of the coordinates of the points. The
directions of the principal axes are the eigenvectors of the covariance matrix
estimated in the original data space or in the transformed space and
transported back to the original space. The lengthening of the ellipse on each
of the two axes is determined by the eigenvalues of the covariance matrix.
Each of the three means is used for the calculation of both the center of mass
and the principal axes of the ellipse. We can see in Fig. 6 that the direction
of the principal axes of the ellipses is almost the same, but there is no
spatial coincidence of their centres. This is expected because the harmonic
mean and the associated ellipse are shifted toward the origin of the data
space while the arithmetic mean and the associated ellipse are the farthest.

Figure 6: Ellipses obtained by using the arithmetic mean (solid), geometric
mean (-.), and harmonic mean (..).

## Conclusion

The arithmetic, geometric, and harmonic means make it possible to estimate the
centrality of measurable phenomena. The arithmetic mean is the most used. The
use of geometric and harmonic means reduces the influence of large
measurements because they give greater influence to small measurements, which
leads to the inequality: AM $>$ GM $>$ HM. Initially, this inequality was
displayed in the case of the central tendencies of two measurements. We have
provided a generalization that gives a geometric interpretation of these means
in the case of several measurements. Through examples, we have shown the
importance of selecting the appropriate mean in the sale of airline tickets,
the ellipse estimation, and the estimation of the Consumer Price Index. We
have also explained three data selection mechanisms implemented by Pythagorean
means. It should be noted that the choice of the mean has been the subject of
several lawsuits. Ultimately, this paper is just a brief overview of the
concept of Pythagorean centrality measures. Much remains to be done to
understand and better exploit the concept of centrality for data selection.
For example, it’s hard to summarize and fully understand a set of measurements
by using a single centrality measure. Selecting measurements using different
centrality measures could provide a more complete picture of data.

## Acknowledgement

I would like to thank A. Hajji for his help in the development of this work.

## References

  * [1] A. Bakker. The Early History of Average Values and Implications for Education. Journal of Statistics Education 11, 2003.
  * [2] C. Alsina and R. Nelsen. Math Made Visual: Creating Images for Understanding Mathematics. Mathematical Association of America, 2006.
  * [3] M.W. Ayech, D. Ziou. Segmentation of Terahertz imaging using k-means clustering based on ranked set sampling. Expert Systems with Applications 42, pp. 2959-2974, 2015.
  * [4] C. Eisenhart. The Development of the Concept of the Best Mean of a Set of Measurements from Antiquity to Present Day. Annual Meeting of the American Statistical Association, 1971.
  * [5] L. D. Broemeling. An Account of Early Statistical Inference in Arab Cryptology. The American Statistician 4, pp. 255-257, 2011.
  * [6] J. M. Borwein and P. B. Borwein. Pi and the AGM: A Study in the Analytic Number Theory and Computational Complexity. Wiley-Interscience, 1987.
  * [7] W. S. Jevons. A Serious Fall in the Value of Gold Ascertained: And Its Social Effects Set Forth. E. Stanford, 1863.
  * [8] F. Coggeshall. The Arithmetic, Geometric, and Harmonic Means. The Quarterly Journal of Economics 1, p. 83-86 , 1886.
  * [9] R. M. Vogel. The geometric mean? Communications in Statistics - Theory and Methods 51, pp. 82-9, 2022.
  * [10] J. R. Taylor. Mécanique classique. De Boeck supérieur, 2012.
  * [11] K. A. Earle et al. Evidence for antibody as a protective correlate for COVID-19 vaccines. Vaccine 39, pp. 4423-4428, 2021.
  * [12] S. Boutemedjet, D. Ziou. and N. Bouguila. Unsupervised Feature Selection for Accurate Recommendation of High-Dimensional Image Data. NIPS, pp. 177-184, 2007.
  * [13] E. Baird. Special relativity considered as an average of earlier theories. Research Gate, 2020.
  * [14] K. C. Chen, N. Bouguila, D. Ziou. Quantization-free parameter space reduction in ellipse detection. Expert Systems with Applications 38, pp. 7622-7632, 2011.
  * [15] A. B. Goumeidane, D. Ziou, and N. Nacereddine. Scale Space Radon Transform for Non Overlapping Thick Ellipses Detection. $11^{th}$ Int, Conf. on Image Processing Theory, Tools and Applications (IPTA), pp. 1-6, 2022.
  * [16] N. Gasnier, L. Denis, and F. Tupin. On the use and denoising of the temporal geometric mean for SAR time series. IEEE Geoscience and Remote Sensing Letters, 2021.
  * [17] J. Hull. Options, futures and other derivatives. Pearson Education, 2015.
  * [18] J. F. Limbrunner, R. M. Vogel, and L. C. Brown. Estimation of Harmonic Mean of a Lognormal Variable. Journal of Hydrologic Engineering 5, pp. 59-66, 2000.
  * [19] G. E. Matthews. When Averaging Multiples, the Arithmetic Mean Is Inferior to the Harmonic Mean. Business Valuation Review 40, pp. 61–67, 2021.
  * [20] G. Quin, B. Pinel-Puysségur, and J.-M. Nicolas. Comparison of Harmonic, Geometric and Arithmetic means for change detection in SAR time series. European Conference on Synthetic Aperture, pp. 255-258, 2012.
  * [21] R. Falk, A. Lann, and S. Zamir. Average Speed Bumps: Four Perspectives on Averaging Speeds. CHANCE 18, pp. 25-32, 2005.
  * [22] M. K. Faradj. Which mean do you mean?: an exposition on means. LSU Master’s Theses, 2004.
  * [23] G. E. Matthews. When Averaging Multiples, Apply the Harmonic Mean. Business Valuation Update 12, 2006\.
  * [24] B. L. Burrows and R. F. Talbot. Which mean do you mean? International Journal of Mathematical Education in Science and Technology 17, pp. 275-284, 1986.
  * [25] Y. Dodge. Arithmetic Mean. The Concise Encyclopedia of Statistics, Springer, 2008.
  * [26] C. Huffman. Archytas. Stanford Encyclopedia of Philosophy, 2008.
  * [27] P. Ballew. The Harmony of the Harmonic Mean, and more Related Problems. Pat’s Blog. https://pballew.blogspot.com/2019/12/the-harmony-of-harmonic-mean-and-more.html, 2019.
  * [28] A. Spalinger. The Rhind Mathematical Papyrus as a Historical Document. Studien zur Altägyptischen Kultur. Helmut Buske Verlag. 17, pp. 295–337, 1990.
  * [29] A. Bakker and K. P. E. Gravemeijer. An Historical Phenomenology of Mean and Median. Educational Studies in Mathematics 62, pp. 149-168, 2006.
  * [30] P. S. Laplace. Essai philosophique sur les probabilités. Courcier, imprimeur-Libraire pour les mathématiques, 1814.
  * [31] D. Ziou and S. Tabbone. Edge Detection Techniques: An Overview. International Journal of Pattern Recognition and Image Analysis 8, pp. 537-559, 1998.
  * [32] G. Almkvist and B. Berndt. Gauss, Landen, Ramanujan, the Arithmetic-Geometric Mean, Ellipses, $\pi$, and the Ladies Diary. The American Mathematical Monthly 95, pp. 585-608, 1988.
  * [33] M. Abramowitz and I. A. Stegun. The Process of the Arithmetic-Geometric Mean. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, pp. 571 and 598-599, 1972.
  * [34] N Bouguila and D Ziou. A countably infinite mixture model for clustering and feature selection. Knowledge and information systems 33 pp. 351-370, 2012.
  * [35] M. Holt and S. M. Scariano. Mean, Median, and Mode from a Decision Perspective. Journal of Statistics Education 17, 2009.
  * [36] Chapter 6 – Calculation of the Consumer Price Index. Statistics Canada, https://www150.statcan.gc.ca/n1/pub/62-553-x/2014001/chap/chap-6-eng.htm. Accessed 6 January 2021.
  * [37] Basket weights of the Consumer Price Index, Canada, provinces, Whitehorse, Yellowknife and Iqaluit. Statistics Canada, https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1810000701. Accessed 13 January 2021.
  * [38] Consumer Price Index by product group, monthly, percentage change, not seasonally adjusted, Canada, provinces, Whitehorse, Yellowknife and Iqaluit. Statistics Canada, https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1810000413. Accessed 13 January 2021.
  * [39] C. Y. Wong, S. C. F. Lin, T. R. Ren, and N. M. Kwok. A survey on ellipse detection methods. IEEE Int. Symposium on Industrial Electronics, pp. 1105-1110, 2012.


# Leveraging Self-Consistency for Data-Efficient
Amortized Bayesian Inference

Marvin Schmitt
University of Stuttgart, Germany
&Daniel Habermann
TU Dortmund, Germany
&Paul-Christian Bürkner
TU Dortmund, Germany
Ullrich Köthe
Heidelberg University, Germany
&Stefan T. Radev
Rensselaer Polytechnic Institute, USA

###### Abstract

We propose a method to improve the efficiency and accuracy of amortized
Bayesian inference (ABI) by leveraging universal symmetries in the
probabilistic joint model $p(\theta,y)$ of parameters $\theta$ and data $y$.
In a nutshell, we invert Bayes’ theorem and estimate the marginal likelihood
based on approximate representations of the joint model. Upon perfect
approximation, the marginal likelihood is constant across all parameter values
by definition. However, approximation error leads to undesirable variance in
the marginal likelihood estimates across different parameter values. We
formulate violations of this symmetry as a loss function to accelerate the
learning dynamics of conditional neural density estimators. We apply our
method to a bimodal toy problem with an explicit likelihood (likelihood-based)
and a realistic model with an implicit likelihood (simulation-based).

##  1 Introduction

Computer simulations are ubiquitous in today’s world, and their widespread
application in the sciences has heralded a new era of _simulation
intelligence_ [18]. Typically, scientific simulators describe a mapping from
latent parameters $\theta$ to observable data $y$. This _forward problem_ is
probabilistically described by the likelihood $p(y\,|\,\theta)$. The _inverse
problem_ of reasoning about the unknown parameters $\theta$ given observed
data $y$ and a prior $p(\theta)$ is described by the posterior
$p(\theta\,|\,y)=p(\theta)\,p(y\,|\,\theta)/p(y)$. The posterior represents a
coherent way to combine all available information in a probabilistic system
[12] and quantify epistemic uncertainty [17]. For complex models, the marginal
likelihood $p(y)$ in the denominator is a high-dimensional integral, rendering
the posterior analytically intractable.

Different from traditional likelihood-based inference [24, 7], simulation-
based inference (SBI) circumvents explicit likelihood evaluation and relies
purely on random samples from a simulation program [8]. In the face of
analytically intractable simulators, previous research has explored other
properties of such programs for learning surrogate likelihood functions or the
likelihood ratio [4, 6, 5]. As a general perspective on simulation
intelligence, amortized Bayesian inference (ABI) is concerned with enabling
fully probabilistic inference in real-time [27, 14, 2].

A core principle of ABI lies in tackling probabilistic problems (forward,
inverse, or both) with neural networks. By re-casting an intractable
probabilistic problem as forward passes through a trained generative neural
network, the required computational time reduces from hours (MCMC) to just a
few seconds (ABI). Yet, there is rarely a free lunch, and neural ABI
algorithms require a potentially long upfront training phase. Nevertheless,
the latter is subsequently repaid with real-time inference on new data sets,
thereby _amortizing_ the training time. To this end, past work has primarily
focused on amortized neural posterior estimation (NPE) for cases where no
explicit likelihood is available [27, 14, 2, 11, 30]. Recently, simulation-
based inference and surrogate modeling have been tackled _jointly_ by learning
an approximate neural posterior $q_{\phi}(\theta\,|\,y)$ and a neural
surrogate likelihood $q_{\eta}(y\,|\,\theta)$ in tandem, an approach coined
neural posterior and likelihood estimation (NPLE;[28, 35, 13]).

Perfect  Symmetry |
---|---
Posterior Approximation |
Figure 1: The performance of the posterior approximator is evaluated via the
variance of the corresponding marginal likelihood estimates. Top row: For the
true posterior (or a perfect approximation thereof), the estimated marginal
likelihood is constant for any parameter value
$\smash{\tilde{\theta}\sim\tilde{p}(\theta)}$. Bottom row: For an imperfect
approximate posterior, the estimated marginal likelihood varies across
different parameter values. Hence, the inherent symmetry of the probabilistic
joint model $p(\theta,y)$ is violated by its approximate representation.
Minimizing the variance of the marginal likelihood estimates pushes the
estimated marginal likelihood towards uniformity. This restores the symmetry
of the unified representation, which is equivalent to improving the
approximate posterior.

##  2 Leveraging Probabilistic Symmetries

The joint model $p(\theta,y)$ implies a fundamental symmetry between marginal
likelihood $p(y)$, prior $p(\theta)$, likelihood $p(y\,|\,\theta)$, and
posterior $p(\theta\,|\,y)$. By inverting Bayes’ theorem, we obtain the
equality

$p(y)=p(\theta)\,p(y\,|\,\theta)\big{/}p(\theta\,|\,y),$ (1)

which must still hold if any component of the joint model is represented
through a _perfect_ approximator $q(\cdot)$. However, we cannot directly use
Equation 1 as a loss function for learning the posterior because the marginal
likelihood on the LHS is notoriously difficult to approximate with high
precision [21, 20]. Instead, we exploit the fact that $p(y)$ is constant
across all parameters $\theta$ (LHS), even though its computation (RHS) hinges
on an arbitrary but fixed parameter value $\theta$. In other words, if we
choose $K$ parameter values $\tilde{\theta}_{1},\ldots,\tilde{\theta}_{K}$,
all computed marginal likelihoods must be _equal_ , regardless of the
individual parameter $\tilde{\theta}_{k}$. We refer to this phenomenon as the
_self-consistency criterion_ :

$\dfrac{p(\tilde{\theta})\,p(y\,|\,\tilde{\theta})}{p(\tilde{\theta}\,|\,y)}=\text{const}\;\forall\tilde{\theta}{\in}\,\Theta\;\,\Longleftrightarrow\,\;\dfrac{p(\tilde{\theta}_{1})\,p(y\,|\,\tilde{\theta}_{1})}{p(\tilde{\theta}_{1}\,|\,y)}=\ldots=\dfrac{p(\tilde{\theta}_{K})\,p(y\,|\,\tilde{\theta}_{K})}{p(\tilde{\theta}_{K}\,|\,y)}\;\;\tilde{\theta}_{1},\ldots,\tilde{\theta}_{K}{\in}\,\Theta$
(2)

###  2.1 Integrating self-consistency into ABI

The _self-consistency loss_ measures the expected degree of violation of the
self-consistency criterion via the variance of the estimated (log) marginal
likelihood across different parameter values $\theta$,

$\mathcal{L}_{\text{SC}}:=\mathbb{E}_{p(y)}\Big{[}\operatorname{Var}_{\tilde{\theta}\sim\tilde{p}(\theta)}\Big{(}\log
p(\tilde{\theta})+\log p(y\,|\,\tilde{\theta})-\log
q_{\phi}(\tilde{\theta}\,|\,y)\Big{)}\Big{]},$ (3)

where the analytic likelihood $p(y\,|\,\tilde{\theta})$ may be replaced with
an approximate $q_{\eta}(y\,|\,\tilde{\theta})$, as demonstrated in Experiment
2. If the variance in Equation 3 is zero, the estimated marginal likelihood is
constant across the parameter space $\Theta$ and the approximation is
consistent (cf. Figure 1). Increasing deviations from a constant estimated
marginal likelihood across the parameter space lead to a larger variance. The
self-consistency loss $\mathcal{L}_{\text{SC}}$ can be seamlessly added to
popular losses for NPE or NPLE. For instance, using the maximum likelihood
loss for NPE with conditional normalizing flows, we obtain:

$\mathcal{L}_{\text{NPE-
SC}}:=\mathbb{E}_{p(y)}\Big{[}\underbrace{\mathbb{E}_{p(\theta\,|\,y)}\big{[}{-}\log
q_{\phi}(\theta\,|\,y)\big{]}}_{\text{NPE
loss}}+\underbrace{\lambda\operatorname{Var}_{\tilde{\theta}\sim\tilde{p}(\theta)}\Big{(}\log
p(\tilde{\theta})+\log p(y\,|\,\tilde{\theta})-\log
q_{\phi}(\tilde{\theta}\,|\,y)\Big{)}}_{\text{self-consistency
loss\;}\mathcal{L}_{\text{SC}}\;\text{with
weight\;}\lambda\in\mathbb{R}_{+}}\Big{]}$ (4)

###  2.2 Monte Carlo estimation of the self-consistency loss

The variance in Equation 3 must be empirically estimated based on finite
samples $\\{\tilde{\theta}_{k}\\}_{k=1}^{K}$ from some proposal distribution
$\tilde{p}(\theta)$. Ideally, the Monte Carlo samples should cover regions of
high variance in the landscape of the estimated marginal likelihood. Due to
the probabilistic symmetry of the joint distribution, high-density regions in
the approximate posterior have the potential to cause large deviations in the
estimated marginal likelihood landscape (cf. Figure 1). Consequently, we
choose the approximate posterior as the proposal, $\tilde{\theta}\sim
q_{\phi}(\theta\,|\,y)$, to render the Monte Carlo estimate efficient. The
approximate posterior has one caveat: It is spectacularly bad at very early
stages of training. Therefore, we can enable the self-consistency loss after
an initial warm-up (e.g., after 5 epochs in Experiment 1) or use a schedule
for progressively increasing its weight $\lambda$ during training.

##  3 Empirical evaluation

###  3.1 Experiment 1: Gaussian mixture model

A

B

C

Figure 2: Experiment 1. Our self-consistent estimator (SC10, SC100, SC500)
outperforms the NPE baseline using the same neural architecture trained on an
identical simulation budget. Adding the self-consistency loss (SC) improves
density estimation (A) and sampling (B), judged both visually and via MMD
between approximate and true posteriors (C). Pink star $\boldsymbol{\star}$
marks the true parameter $\theta$.

We illustrate our method using a 2-dimensional Gaussian mixture model with two
symmetrical components [11]. The symmetrical components are equally weighted
and have equal variance,

$\theta\sim\mathcal{N}(\theta\,|\,\mathbf{0},\mathbf{I}),\quad\quad y\sim
0.5\,\mathcal{N}(y\,|\,\theta,\mathbf{I}/2)+0.5\,\mathcal{N}(y\,|\,{-}\theta,\mathbf{I}/2),$
(5)

where $\mathcal{N}(\cdot\,|\,\mu,\Sigma)$ is a Gaussian distribution with
location $\mu$ and covariance matrix $\Sigma$. Each simulated data set $y$
consists of ten observations. The simulation budget is strictly limited to
$N{=}1024$ data sets for training. We train four networks for 35 epochs each:
NPE (baseline), SC10, SC100, and SC500, where SC stands for self-consistency
and 10/100/500 is the number $K$ of Monte Carlo samples for estimating the
variance in Equation 4. All networks use the same neural spline flow
architecture [9] with a DeepSet [36] as a summary network and a heavy-tailed
Student-$t_{100}$ latent distribution [1] for the posterior approximator. We
enable the self-consistency loss after five epochs (cf. Section 2.2).

#### Results

While NPE struggles to fit the posterior, our additional self-consistency loss
drastically improves the approximate posterior _for the same simulation
budget_ (see Figure 2). All self-consistent variants outperform NPE through
(i) more accurate posterior densities; and (ii) better samples, as indexed by
visual inspection and by lower maximum mean discrepancy (MMD;[16]) to samples
from the true posterior on 50 new data sets. All approximators are well-
calibrated (see Appendix C). Increasing the number of consistency samples
beyond $K{=}100$ does not visibly improve performance. We parallel this
experiment with an approximate likelihood in Appendix E and also observe
superior performance of our self-consistent approximator with respect to
density estimation and sampling.

###  3.2 Experiment 2: Oscillatory Hes1 expression model

As a more realistic example, we apply our method to an experimental data set
in biology [31]. Upon serum stimulation of various cell lines, the
transcription factor Hes1 exhibits sustained oscillatory transcription
patterns [23]. The measured concentration of Hes1 mRNA can be modeled with a
set of three differential equations which are governed by four parameters
$\theta=(p_{0},h,k_{1},\nu)$ with fixed initial conditions according to [10]
(see Appendix D for details).

We use a fixed simulation budget of $N{=}512$ data sets for neural network
training. This simulation budget enables amortized inference, yet it is orders
of magnitude smaller than the required budget for approximate Bayesian
computation algorithms [33] such as ABC-SMC [32] for a single observed data
set [31]. In order to use an approximate likelihood $q_{\eta}(y\,|\,\theta)$
in the self-consistency loss, we train a neural surrogate likelihood in tandem
with the neural posterior (NPLE;[28]). Both NPLE (baseline) and our self-
consistent approximator are trained for 70 epochs and use the same neural
spline flow architecture with spectral normalization [22] and a heavy-tailed
Student-$t_{50}$ latent distribution [1]. The self-consistent approximator
uses $K{=}500$ Monte Carlo samples, and the self-consistency loss is activated
after 10 epochs.

A  NPLE

B

C  SC (Ours)

D

Figure 3: Experiment 2. The baseline NPLE approximator shows deficient
simulation-based calibration, as indexed by ECDF lines outside the gray 95%
confidence bands (A). In contrast, our self-consistent approximator is well-
calibrated (C). Samples from the posterior predictive distribution on real
experimental data (black dots; [31]) are comparable between NPLE (B) and SC
(D).

#### Results

Our self-consistent approximator with approximate likelihood shows superior
simulation-based calibration [34, 29] compared to the NPLE baseline,
particularly with respect to the parameters $p_{0}$ and $h$ (see Figure 3A and
C). The posterior predictive distributions of both methods have a comparable
fit to the real experimental time series $y_{\text{real}}$ from [31] (see
Figure 3B and D).

##  4 Conclusion

In this paper, we proposed a new method to leverage inherent symmetry in a
joint model $p(\theta,y)$ to improve amortized Bayesian inference. In two
experiments, we illustrated that the combination of simulation-based inference
and (approximate) likelihood-based learning increases the efficiency of neural
posterior estimation. Concretely, we demonstrated that an additional self-
consistency loss leads to (i) better densities and (ii) better samples from
the approximate posterior in situations with scarce training data The latter
occurs frequently in real-world applications of simulation-based inference
(e.g., [38, 37, 3]), where the simulation program is a computational
bottleneck or other reasons strictly limit the simulation budget for neural
network training.

While this paper focused on amortized Bayesian inference with conditional
normalizing flows, our self-consistency loss can readily be applied to
sequential simulation-based inference [26, 15, 13, 35]. Likewise, other
conditional density estimators like score modeling [11, 30, 25] or flow-
matching [19] may certainly benefit from our additional loss function as well.
Finally, future research could explore variations and extensions of our
proposed method, such as different proposal distributions $\tilde{p}(\theta)$
for more efficient Monte Carlo estimates, or better loss functions altogether
that build on the principle of self-consistency. A selection of frequently
asked questions (FAQ) which a reader might have is answered in Appendix A.

## Acknowledgments

MS was supported by the Cyber Valley Research Fund (grant number: CyVy-
RF-2021-16) and the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy EXC-2075 - 390740016 (the
Stuttgart Cluster of Excellence SimTech). DH was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under project
508399956. UK was supported by the Informatics for Life initiative funded by
the Klaus Tschira Foundation.

## References

  * [1] Simon Alexanderson and Gustav Eje Henter.  Robust model training and generalisation with studentising flows, 2020\.
  * [2] Grace Avecilla, Julie N Chuong, Fangfei Li, Gavin Sherlock, David Gresham, and Yoav Ram.  Neural networks enable efficient and accurate simulation-based inference of evolutionary parameters from adaptation dynamics.  PLoS Biology, 20(5):e3001633, 2022.
  * [3] Ayush Bharti, Louis Filstroff, and Samuel Kaski.  Approximate Bayesian computation with domain expert in the loop.  In International Conference on Machine Learning, pages 1893–1905. PMLR, 2022.
  * [4] Johann Brehmer, Kyle Cranmer, Gilles Louppe, and Juan Pavez.  A guide to constraining effective field theories with machine learning.  Physical Review D, 98(5):052004, 2018.
  * [5] Johann Brehmer, Felix Kling, Irina Espejo, and Kyle Cranmer.  Madminer: Machine learning-based inference for particle physics.  Computing and Software for Big Science, 4:1–25, 2020.
  * [6] Johann Brehmer, Gilles Louppe, Juan Pavez, and Kyle Cranmer.  Mining gold from implicit models to improve likelihood-free inference.  Proceedings of the National Academy of Sciences, 117(10):5242–5249, 2020.
  * [7] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell.  Stan: A probabilistic programming language.  Journal of statistical software, 76(1), 2017.
  * [8] Kyle Cranmer, Johann Brehmer, and Gilles Louppe.  The frontier of simulation-based inference.  Proceedings of the National Academy of Sciences, 2020.
  * [9] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  Neural spline flows.  Advances in neural information processing systems, 32, 2019.
  * [10] Sarah Filippi, Chris Barnes, Julien Cornebise, and Michael P. H. Stumpf.  On optimality of kernels for approximate bayesian computation using sequential monte carlo, June 2011.
  * [11] Tomas Geffner, George Papamakarios, and Andriy Mnih.  Compositional score modeling for simulation-based inference, 2022.
  * [12] Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin.  Bayesian Data Analysis (3rd Edition).  Chapman and Hall/CRC, 2013.
  * [13] Manuel Glöckler, Michael Deistler, and Jakob H. Macke.  Variational methods for simulation-based inference.  In International Conference on Learning Representations, 2022.
  * [14] Pedro J Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, et al.  Training deep neural density estimators to identify mechanistic models of neural dynamics.  Elife, 2020.
  * [15] David Greenberg, Marcel Nonnenmacher, and Jakob Macke.  Automatic posterior transformation for likelihood-free inference.  In International Conference on Machine Learning, 2019.
  * [16] A Gretton, K. Borgwardt, Malte Rasch, Bernhard Schölkopf, and AJ Smola.  A Kernel Two-Sample Test.  The Journal of Machine Learning Research, 13:723–773, 2012.
  * [17] Eyke Hüllermeier and Willem Waegeman.  Aleatoric and Epistemic Uncertainty in Machine Learning: An Introduction to Concepts and Methods.  Machine Learning, 110(3):457–506, March 2021.
  * [18] Alexander Lavin, Hector Zenil, Brooks Paige, et al.  Simulation intelligence: Towards a new generation of scientific methods.  arXiv preprint, 2021.
  * [19] Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le.  Flow matching for generative modeling.  In The 11th International Conference on Learning Representations, 2023.
  * [20] Sanae Lotfi, Pavel Izmailov, Gregory Benton, Micah Goldblum, and Andrew Gordon Wilson.  Bayesian model selection, the marginal likelihood, and generalization.  In International Conference on Machine Learning, pages 14223–14247. PMLR, 2022.
  * [21] Xiao-Li Meng and Wing Hung Wong.  Simulating ratios of normalizing constants via a simple identity: a theoretical exploration.  Statistica Sinica, 1996.
  * [22] Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida.  Spectral normalization for generative adversarial networks.  In International Conference on Learning Representations, 2018.
  * [23] Hiroshi Momiji and Nicholas A.M. Monk.  Dissecting the dynamics of the hes1 genetic oscillator.  Journal of Theoretical Biology, 254(4):784–798, oct 2008.
  * [24] Radford M. Neal.  MCMC using Hamiltonian dynamics.  May 2011.
  * [25] Lorenzo Pacchiardi and Ritabrata Dutta.  Score matched neural exponential families for likelihood-free inference.  J. Mach. Learn. Res., 23:38–1, 2022.
  * [26] George Papamakarios, David Sterratt, and Iain Murray.  Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows.  2019\.
  * [27] Stefan T Radev, Ulf K Mertens, Andreas Voss, Lynton Ardizzone, and Ullrich Köthe.  BayesFlow: Learning complex stochastic models with invertible neural networks.  IEEE transactions on neural networks and learning systems, 2020\.
  * [28] Stefan T. Radev, Marvin Schmitt, Valentin Pratz, Umberto Picchini, Ullrich Köthe, and Paul-Christian Bürkner.  JANA: Jointly Amortized Neural Approximation of Complex Bayesian Models.  In Robin J. Evans and Ilya Shpitser, editors, Proceedings of the 39th Conference on Uncertainty in Artificial Intelligence, volume 216 of Proceedings of Machine Learning Research, pages 1695–1706. PMLR, 2023.
  * [29] Teemu Säilynoja, Paul-Christian Bürkner, and Aki Vehtari.  Graphical test for discrete uniformity and its applications in goodness-of-fit evaluation and multiple sample comparison.  Statistics and Computing, 32(2):1–21, 2022.
  * [30] Louis Sharrock, Jack Simons, Song Liu, and Mark Beaumont.  Sequential neural score estimation: Likelihood-free inference with conditional score based diffusion models, 2022.
  * [31] Daniel Silk, Paul D.W. Kirk, Chris P. Barnes, Tina Toni, Anna Rose, Simon Moon, Margaret J. Dallman, and Michael P.H. Stumpf.  Designing attractive models via automated identification of chaotic and oscillatory dynamical regimes.  Nature Communications, 2(1), 2011.
  * [32] S. A. Sisson, Y. Fan, and Mark M. Tanaka.  Sequential monte carlo without likelihoods.  Proceedings of the National Academy of Sciences, 104(6):1760–1765, February 2007.
  * [33] Scott A Sisson, Yanan Fan, and Mark Beaumont.  Handbook of approximate Bayesian computation.  CRC Press, 2018.
  * [34] Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari, and Andrew Gelman.  Validating Bayesian inference algorithms with simulation-based calibration.  arXiv preprint, 2018.
  * [35] Samuel Wiqvist, Jes Frellsen, and Umberto Picchini.  Sequential neural posterior and likelihood approximation.  arXiv preprint, 2021.
  * [36] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Ruslan Salakhutdinov, and Alexander Smola.  Deep sets, 2017.
  * [37] Jice Zeng, Michael D Todd, and Zhen Hu.  Probabilistic damage detection using a new likelihood-free Bayesian inference method.  Journal of Civil Structural Health Monitoring, 13(2-3):319–341, 2023.
  * [38] Yi Zhang and Lars Mikelsons.  Sensitivity-guided iterative parameter identification and data generation with BayesFlow and PELS-VAE for model calibration.  Advanced Modeling and Simulation in Engineering Sciences, 10(1):1–28, 2023.

APPENDIX

##  Appendix A Frequently Asked Questions (FAQ)

Q: How can I reproduce the results?
Code to reproduce all results will be released upon archival publication.

Q: How does the self-consistency loss change my current amortized Bayesian
workflow?
If you have access to an explicit likelihood, you need to implement a
likelihood.log_prob method to evaluate $p(y\,|\,\tilde{\theta})$. This is
straightforward with common frameworks such as tensorflow_probability or
scipy.stats. Equivalently, you need to implement a prior.log_prob method for
your prior distribution, regardless of whether you use NPE (explicit
likelihood) or NPLE (implicit likelihood).

Q: When is it useful to add the self-consistency loss?
When the simulation program is computationally costly or the simulation budget
is fixed, our experiments suggest that adding the self-consistency loss to
your current optimization objective might help you get more out of the data.

Q: What about learned summary statistics?
Our method is fully compatible with end-to-end learning of summary statistics
alongside the neural approximator [27, 28]. In fact, Experiment 1 uses a
DeepSet [36] to learn fixed-length summary statistics from the observables,
which are then passed to the posterior approximator.

Q: When do I activate the self-consistency loss during training?
This depends on the complexity of the problem. Your approximate posterior (and
approximate likelihood, if applicable) should be sufficiently good so that (i)
the proposals for $\tilde{\theta}$ cover relevant regions; and (ii) the log
density estimates for the posterior (and likelihood, if applicable) have an
acceptable quality for the Monte Carlo estimate in Equation 3. Future work
might address this question with analyses of the learning dynamics which
govern the neural approximator training.

Q: Why aren’t the posterior samples in Experiment 1 perfectly aligned with the
true parameter?
The simulated data sets in the Gaussian mixture model only consist of ten
observations from the Gaussian mixture model with locations $\theta$ and
${-}\theta$. Due to aleatoric uncertainty in the data-generating process [17],
the empirical information in the sample does not even suffice to inform the
true posteriorto concentrate on the true data-generating parameter $\theta$.
Instead, the goal of an approximate posterior is to match the true posterior,
which all self-consistent approximators achieve.

##  Appendix B Algorithm

The Monte Carlo estimation of the self-consistency loss in Equation 3 and
Equation 4 is outlined in Algorithm 1.

Algorithm 1 Self-consistency loss for finite training

0: $N$ training data tuples $\\{(\theta^{(i)},y^{(i)})\\}_{i=1}^{N}$

0: Number of self-consistency samples $K$

1: for $i=1,\ldots,N$ do

2: < compute other losses such as NPE/NPLE loss >

3: for $k=1,\ldots,K$ do

4: $\tilde{\theta}_{k}\sim\tilde{p}(\theta)$ {sample from proposal
$\tilde{p}(\theta)$}

5: $\log\hat{p}_{k}(y^{(i)})=\begin{cases}\log p(\tilde{\theta}_{k})+\log
p(y^{(i)}\,|\,\tilde{\theta}_{k})\;\,-\log
q_{\phi}(\tilde{\theta}_{k}\,|\,y^{(i)})\;\text{\\{explicit likelihood
(NPE)\\}}\\\ \log p(\tilde{\theta}_{k})+\log
q_{\eta}(y^{(i)}\,|\,\tilde{\theta}_{k})-\log
q_{\phi}(\tilde{\theta}_{k}\,|\,y^{(i)})\text{\\{implicit likelihood
(NPLE)\\}}\end{cases}$

6: end for

7:
$\mathcal{L}_{SC}^{(i)}=\operatorname{Var}(\\{\log\hat{p}_{k}(y^{(i)})\\}_{k=1}^{K})${variance
of logML estimates}

8: end for

9: return $\frac{1}{N}\sum_{i=1}^{N}\mathcal{L}_{SC}^{(i)}${average over
training data}

##  Appendix C Details about Experiment 1

All approximators (NPE and self-consistent ones) are well-calibrated according
to simulation-based calibration [34, 29], as illustrated in Figure 4.

NPE |  | SC10
---|---|---
|  |
SC100 |  | SC500
|  |
Figure 4: Experiment 1. All approximators are well-calibrated.

##  Appendix D Details about Experiment 2

In Experiment 2, we apply our method to an experimental data set in biology
[31]. Upon serum stimulation of various cell lines, the transcription factor
Hes1 exhibits sustained oscillatory transcription patterns [23]. The
concentration of Hes1 mRNA can be modeled by a set of three differential
equations,

$\odv{m}{t}=-k_{deg}m+\frac{1}{1+(p_{2}/p_{0})^{h}},\quad\odv{p_{1}}{t}=-k_{deg}p_{1}+\nu
m-k_{1}p_{1},\quad\odv{p_{2}}{t}=-k_{deg}p_{2}+k_{1}p_{1}$ (6)

with degradation rate $k_{deg}$, Hes1 mRNA concentration $m$, cytosolic Hes1
protein concentration $p_{1}$, and nuclear Hes1 protein concentration $p_{2}$.
The parameters $\theta{=}\\{p_{0},h,k_{1},\nu\\}$ govern the dynamics of the
differential equations, and we estimate them in the unbounded $\log$ space to
facilitate inference. $p_{0}$ corresponds to the amount of Hes1 protein in the
nucleus when the rate of transcription of Hes1 mRNA is at half of its maximum
value, $h$ is the Hill coefficient, $k_{1}$ is the rate of transport of Hes1
protein into the nucleus and $\nu$ is the rate of translation of Hes1 mRNA
[31].

In accordance to [10], we use fixed initial conditions
$m_{0}{=}2,p_{1}{=}5,p_{2}{=}3$ and set $k_{reg}{=}0.03$. In our model we
regard the observed mRNA concentrations $y_{t}$ as noisy measurements of the
true underlying mRNA concentration $m_{t}$ with unit Gaussian observation
error:

$y_{t}\sim\mathcal{N}(m_{t},1)$ (7)

Silk et al. [31] used quantitative real-time PCR to collect the real
experimental data

$y_{t}=[1.20,5.90,4.58,2.64,5.38,6.42,5.60,4.48]$

where the first observation $y_{1}$ is measured after 30 minutes, and all
subsequent values are measured in $30$ minute intervals [10]. The mRNA
measures $y_{t}$ refer to fold changes relative to a control sample. The
Bayesian model uses Gamma priors on all parameters,

$p_{0}\sim\Gamma(2,1),\quad h\sim\Gamma(10,1),\quad
k_{1}\sim\Gamma(2,50),\quad\nu\sim\Gamma(2,50),$ (8)

where $\Gamma(a,b)$ denotes the Gamma distribution with shape $a$ and rate
$b$.

##  Appendix E Gaussian mixture model with implicit likelihood

We repeat Experiment 1 with a fully simulation-based approach which does not
need an explicit likelihood to estimate the self-consistency loss. To this
end, we replace the explicit likelihood $p(y\,|\,\theta)$ in Equation 3 with
the approximate likelihood $q_{\eta}(y\,|\,\theta)$, which is represented by a
neural network and learned simultaneously with the neural posterior
approximator. The full loss follows as

$\displaystyle\mathcal{L}_{\text{NPLE-SC}}=\mathbb{E}_{p(\theta,y)}\Big{[}$
$\displaystyle\underbrace{-\log q_{\phi}(\theta\,|\,y)-\log
q_{\eta}(y\,|\,\theta)}_{\text{NPLE loss}}$ (9)
$\displaystyle+\underbrace{\lambda\operatorname{Var}_{\tilde{\theta}\sim\tilde{p}(\theta)}\Big{(}\log
p(\tilde{\theta})+\log q_{\eta}(y\,|\,\tilde{\theta})-\log
q_{\phi}(\tilde{\theta}\,|\,y)\Big{)}}_{\text{self-consistency loss
$\mathcal{L}_{\text{SC}}$ with approximate likelihood\;}}\Big{]},$

which is a combination of the NPLE loss and the self-consistency loss with
approximate likelihood.

The simulation budget is fixed to $N{=}1024$ data sets. Both NPLE and our
self-consistent approximator with $K{=}100$ Monte Carlo samples are trained
for 35 epochs, have an identical neural spline flow architecture, have a
heavy-tailed Student-$t_{100}$ latent space [1], and use an identical DeepSet
[36] to learn summary statistics of the data $y$ for the posterior
approximator. Since the Monte Carlo approximation in the self-consistency loss
now depends on both an approximate posterior and an approximate likelihood, we
only activate the self-consistency loss after 20 epochs (as opposed to 5
epochs in Experiment 1 with an explicit likelihood).

We can benchmark the methods’ performance against the true posterior since the
explicit likelihood of this simulator is known (albeit inaccessible for the
approximators). We confirm the results of Experiment 1 in the fully
simulation-based setting with only an implicit likelihood: Our self-consistent
approximator consistently outperforms the baseline NPLE approximator with
respect to posterior density and sampling (see Figure 5).

A

B

C

D

E

Figure 5: Our self-consistent posterior estimator (SC100) outperforms the NPLE
baseline using the same neural architecture trained on an identical simulation
budget. Adding the self-consistency loss leads to improved density estimation
(A) and sampling (B), judged both visually and via MMD between approximate and
true posteriors (C). Both approximators are well-calibrated (D, E). Pink star
$\boldsymbol{\star}$ marks the true parameter $\theta$.


$\frac{1}{r_{1}}+\frac{1}{r_{2}}+\frac{1}{p_{3}}+\frac{1}{r_{4}}\leq 1$ and
$A(p_{i}\to r_{i})\lesssim 1$ for all $i=1,2,4.$ Then we have

$Y(p_{1},p_{2},p_{3},p_{4})\lesssim 1.$

###### Proof.

By combining the inequality (10.3) with our assumptions on the averaging
estimates, it follows that for all functions $f_{i},i=1,2,3,4,$ on
$\mathbb{F}_{q}^{d},$

$Y(f_{1},f_{2},f_{3},f_{4})\lesssim||f_{1}||_{p_{1}}||f_{2}||_{p_{2}}||f_{3}||_{p_{3}}||f_{4}||_{p_{4}}.$

This completes the proof. ∎

The following theorem provides lots of sufficient conditions for the
boundedness of the $Y$-shaped operator.

######  Theorem 10.4.

Let $1\leq p_{1},p_{2},p_{3},p_{4}\leq\infty,$ and $Y$ be the $Y$-shaped
operator on $\mathbb{F}_{q}^{d}.$ Then $Y(p_{1},p_{2},p_{3},p_{4})\lesssim 1$
provided that one of the following conditions is satisfied:

  1. (i)

$0\leq\frac{1}{p_{1}},\frac{1}{p_{2}},\frac{1}{p_{4}}\leq\frac{d}{d+1}$ and
$\frac{1}{p_{1}}+\frac{1}{p_{2}}+\frac{d}{p_{3}}+\frac{1}{p_{4}}\leq d.$

  2. (ii)

$0\leq\frac{1}{p_{1}},\frac{1}{p_{2}}\leq\frac{d}{d+1}\leq\frac{1}{p_{4}}$ and
$\frac{1}{dp_{1}}+\frac{1}{dp_{2}}+\frac{1}{p_{3}}+\frac{d}{p_{4}}\leq d.$

  3. (iii)

$0\leq\frac{1}{p_{1}},\frac{1}{p_{4}}\leq\frac{d}{d+1}\leq\frac{1}{p_{2}}\leq
1$ and $\frac{1}{dp_{1}}+\frac{d}{p_{2}}+\frac{1}{p_{3}}+\frac{1}{dp_{4}}\leq
d.$

  4. (iv)

$0\leq\frac{1}{p_{2}},\frac{1}{p_{4}}\leq\frac{d}{d+1}\leq\frac{1}{p_{1}}\leq
1$ and $\frac{d}{p_{1}}+\frac{1}{dp_{2}}+\frac{1}{p_{3}}+\frac{1}{dp_{4}}\leq
d.$

  5. (v)

$0\leq\frac{1}{p_{1}}\leq\frac{d}{d+1}\leq\frac{1}{p_{2}},\frac{1}{p_{4}}\leq
1$ and $\frac{1}{dp_{1}}+\frac{d}{p_{2}}+\frac{1}{p_{3}}+\frac{d}{p_{4}}\leq
2d-1.$

  6. (vi)

$0\leq\frac{1}{p_{2}}\leq\frac{d}{d+1}\leq\frac{1}{p_{1}},\frac{1}{p_{4}}\leq
1$ and $\frac{d}{p_{1}}+\frac{1}{dp_{2}}+\frac{1}{p_{3}}+\frac{d}{p_{4}}\leq
2d-1.$

  7. (vii)

$0\leq\frac{1}{p_{4}}\leq\frac{d}{d+1}\leq\frac{1}{p_{1}},\frac{1}{p_{2}}\leq
1$ and $\frac{d}{p_{1}}+\frac{d}{p_{2}}+\frac{1}{p_{3}}+\frac{1}{dp_{4}}\leq
2d-1.$

  8. (viii)

$\frac{d}{d+1}\leq\frac{1}{p_{1}},\frac{1}{p_{2}},\frac{1}{p_{4}}\leq 1$ and
$\frac{d}{p_{1}}+\frac{d}{p_{2}}+\frac{1}{p_{3}}+\frac{d}{p_{4}}\leq 3d-2.$

###### Proof.

The proof uses Proposition 10.3 and the sharp averaging estimates in Lemma
2.5. The proof of this theorem is similar to that of Theorem 5.4. Therefore,
we leave the detail of the proof to readers. ∎

Conjecture 1.5 is also supported by the following theorem.

######  Theorem 10.5.

Let ${\diamondsuit}_{t}$ and $Y$ be the operators acting on the functions on
$\mathbb{F}_{q}^{2}.$ If ${\diamondsuit}_{t}(p_{1},p_{2},p_{3},p_{4})\lesssim
1$ with $1\leq p_{1},p_{2},p_{3},p_{4}\leq\infty,$ then
$Y(p_{1},p_{2},p_{3},p_{4})\lesssim 1.$

###### Proof.

Assume that ${\diamondsuit}_{t}(p_{1},p_{2},p_{3},p_{4})\lesssim 1.$ Then, by
Lemma 6.2, $(1/p_{1},1/p_{2},1/p_{3},1/p_{4})$ is contained in the convex hull
of the points $(0,0,1,0),$ $(0,1,0,0),$ $(0,0,0,1),$ $(1/2,0,1/2,1/2),$
$(2/3,2/3,0,0),$ $(1,0,0,0),$ $(2/3,0,2/3,0),$ $(1/2,1/2,1/2,0),$
$(2/3,0,0,2/3),$ $(0,2/3,2/3,0),$ $(0,0,0,0),$ $(0,0,2/3,2/3).$ By
interpolating those critical points, it suffices to check that each critical
point above satisfies one of the 8 hypotheses of Theorem 10.4 with $d=2$.
However, this can be easily shown by a direct computation. For example, for
the critical point $(1/p_{1},1/p_{2},1/p_{3},1/p_{4})=(1/2,1/2,1/2,0)$, we can
invoke the hypothesis (i) of Theorem 10.4 with $d=2$ and obtain that
$Y(2,2,2,\infty)\lesssim 1.$ In a same way, it can be easily proven for other
critical points. ∎

######  Remark 10.6.

The reverse statement of Theorem 10.5 is not true. To find a counterexample,
we can take $p_{1}=p_{3}=\infty$, $p_{2}=p_{4}=3/2.$ Indeed, by the hypothesis
(5) of Theorem 10.4 with $d=2$, we see that $Y(\infty,3/2,\infty,3/2)\lesssim
1.$ However, ${\diamondsuit}_{t}(\infty,3/2,\infty,3/2)$ is not bounded, which
follows from Lemma 6.2 with $d=2.$

The following corollary proposes some possibility that the assumption of the
subgraph in Conjecture 1.5 can be dropped.

######  Corollary 10.7.

Let $\Diamond$ and $Y$ be the operators acting on the functions on
$\mathbb{F}_{q}^{2}.$ If $\Diamond(p_{1},p_{2},p_{3},p_{4})\lesssim 1$ with
$1\leq p_{1},p_{2},p_{3},p_{4}\leq\infty,$ then
$Y(p_{1},p_{2},p_{3},p_{4})\lesssim 1.$

###### Proof.

The statement of the corollary follows immediately by combining Theorem 7.7
and Theorem 10.5. ∎

Combining all the results obtained so far, we get the following theorem:

######  Theorem 10.8.

When $d=2$ and $n=3,4$, Conjecture 1.5 is true, where we accept boundedness
results up to endpoints in the case when $G$ is the $C_{4}$ \+ diagonal and
its subgraph $G^{\prime}$ is the $K_{3}+tail.$

###### Proof.

By Corollary 5.8 for $n=3,$ and by Corollaries 8.11, 8.14, 9.13, 10.7 for
$n=4,$ we have proven that for $d=2$ and $n=3,4,$ there is the required
inclusive boundedness relationship between any two operators corresponding to
arbitrary connected ordered graph $G$ and its subgraph $G^{\prime}$ except for
the following three cases:

  * (I)

$G=C_{4}$ \+ diagonal and $G^{\prime}=C_{4}.$

  * (II)

$K_{3}$ \+ tail and $G^{\prime}=Y$-shape.

  * (III)

$G=K_{3}$ \+ tail and $G^{\prime}=P_{3}.$

However, since $\delta(G)=\delta(G^{\prime})$ for each case of (I), (II),
(III), they do not satisfy the main hypothesis (1.6) of Conjecture 1.5. Hence,
they cannot be counterexamples contradicting Conjecture 1.5 and so there is no
counterexample against Conjecture 1.5, as required. ∎

##  11 Acknowledgements

A. Iosevich and P. Bhowmik were supported in part by the National Science
Foundation grant no. HDR TRIPODS–1934962 and the National Science Foundation
grant DMS–2154232. D. Koh was supported by Basic Science Research Programs
through National Research Foundation of Korea (NRF) funded by the Ministry of
Education (NRF-2018R1D1A1B07044469). T. Pham would like to thank to the VIASM
for the hospitality and for the excellent working conditions.

##  12 Appendix

In this appendix, we introduce the number of intersection points of two
spheres in $\mathbb{F}_{q}^{d}.$ Let $\eta$ denote the quadratic character of
$\mathbb{F}_{q}^{*}$, namely, $\eta(s)=1$ for a square number $s$ in
$\mathbb{F}_{q}^{*}$, and $\eta(s)=-1$ otherwise.

######  Definition 12.1.

Given a non-zero vector $m$ in $\mathbb{F}_{q}^{d},$ and
$t,b\in\mathbb{F}_{q},$ we define $N(m,t,b)$ to be the number of common
solutions $x\in\mathbb{F}_{q}^{d}$ of the following equations:
$||x||=t,~{}~{}m\cdot x=b.$

Notice that the value of $N(m,t,b)$ is the number of all intersection points
between the sphere $S_{t}$ and the plane $\\{x\in\mathbb{F}_{q}^{d}:m\cdot
x=b\\}.$ The explicit value of it is well known as follows.

######  Lemma 12.2.

Let $b,t\in\mathbb{F}_{q},$ and let $m$ be a non-zero element in
$\mathbb{F}_{q}^{d},d\geq 2.$ Then the following statements hold:

  1. (i)

If $||m||\neq 0$ and $b^{2}-t||m||=0,$ then

$N(m,t,b)=\left\\{\begin{array}[]{ll}q^{d-2}&\mbox{if}~{}~{}d~{}\mbox{is
even,}\\\
q^{d-2}+q^{\frac{d-3}{2}}(q-1)\eta\left((-1)^{\frac{d-1}{2}}||m||\right)&\mbox{if}~{}~{}d~{}\mbox{is
odd}.\end{array}\right.$

  2. (ii)

If $||m||\neq 0$ and $b^{2}-t||m||\neq 0,$ then

$N(m,t,b)=\left\\{\begin{array}[]{ll}q^{d-2}+q^{\frac{d-2}{2}}\eta\left((-1)^{\frac{d}{2}}(b^{2}-t||m||)\right)&\mbox{if}~{}~{}d~{}\mbox{is
even,}\\\
q^{d-2}-q^{\frac{d-3}{2}}\eta\left((-1)^{\frac{d-1}{2}}||m||\right)&\mbox{if}~{}~{}d~{}\mbox{is
odd}.\end{array}\right.$

  3. (iii)

If $||m||=0=b^{2}-t||m||,$ then

$N(m,t,b)=\left\\{\begin{array}[]{ll}q^{d-2}+\nu(t)q^{\frac{d-2}{2}}\eta\left((-1)^{\frac{d}{2}})\right)&\mbox{if}~{}~{}d~{}\mbox{is
even,}\\\
q^{d-2}-q^{\frac{d-1}{2}}\eta\left((-1)^{\frac{d-1}{2}}t\right)&\mbox{if}~{}~{}d~{}\mbox{is
odd},\end{array}\right.$

where $\nu(t)=-1$ if $t\in\mathbb{F}_{q}^{*}$ and $\nu(0)=q-1.$

  4. (iv)

If $||m||=0$ and $b^{2}-t||m||\neq 0,$ then $N(m,t,b)=q^{d-2}.$

###### Proof.

See Exercises 6.31–6.34 in [16], or one can prove it by using the discrete
Fourier analysis with the explicit value of the Gauss sum. ∎

By a direct application of Lemma 12.2, one can find the explicit number of the
intersections of two spheres over finite fields. Precisely we have the
following result.

######  Theorem 12.3.

Given a non-zero vector $m\in\mathbb{F}_{q}^{d}$ and $t,j\in\mathbb{F}_{q},$
let

$\Theta(m,t,j):=\\{x\in S_{t}:||x-m||=j\\}|.$

If $m\in S_{\ell},$ then
$|\Theta(m,t,j)|=N\left(m,t,\frac{t+\ell-j}{2}\right).$

###### Proof.

Since $||x-m||=t+\ell-2m\cdot x$ for $x\in S_{t},m\in S_{\ell},$ it is clear
that $\Theta(m,t,j)$ is the number of common solutions $x$ of the following
equations:

$||x||=t,\quad m\cdot x=\frac{t+\ell-j}{2}.$

Hence, by the definition of $N$, we obtain the required conclusion. ∎

######  Corollary 12.4.

Let $t\in\mathbb{F}_{q}^{*}$ and $\ell\in\mathbb{F}_{q}.$ Then, for every non-
zero vector $m\in S_{\ell},$ we have

$\sum_{x\in S_{t}:||x-m||=t}1\sim q^{d-2}$

excepting for the following three cases:

1) $d=2,\ell\neq 0,\eta(t\ell-\ell^{2}/4)=-1.$ 2) $d=2,\ell=0,\eta(-1)=1.$
---
3) $d=3,\ell=0,\eta(-t)=1.$

For each of those three cases, the value in the above sum takes zero. On the
other hand, if $d=2,\ell\neq 0,$ and $\eta(t\ell-\ell^{2}/4)=1$, the value in
the above sum is exactly two.

###### Proof.

It follows from Theorem 12.3 that for any $||m||=\ell,$

$\sum_{x\in S_{t}:||x-m||=t}1=N\left(m,t,\frac{\ell}{2}\right),$

and so the corollary is a direct consequence of Lemma 12.2 (i), (ii), (iii). ∎

## References

  * [1] B. Assarf, E. Gawrilow, K. Herr, M. Joswig, B. Lorenz, A. Paffenholz and T. Rehn, _Computing convex hulls and counting integer points with polymake,_ Math. Program. Comput. 9 (2017), no.1, 1-38.
  * [2] M. Bennett, D. Hart, A. Iosevich, J. Pakianathan, and M. Rudnev, Group actions and geometric combinatorics in $\mathbb{F}_{q}^{d},$ Forum Math. 29 (2017), no.1, 91-110.
  * [3] A. Carbery, B. Stones, and J. Wright, _Averages in vector spaces over finite fields_ , Math. Proc. Cambridge Philos. Soc. 144 (2008), no.1, 13-27.
  * [4] J. Chapman, M. Erdog̃an, D. Hart, A. Iosevich, and D. Koh, Pinned distance sets, Wolff’s exponent in finite fields and sum-product estimates, Math.Z. 271, (2012), no.1-2, 63-93.
  * [5] X. Du, A. Iosevich, Y. Ou, H. Wang and R. Zhang, An improved result for Falconer’s distance set problem in even dimensions, Math. Ann. 380 (2021), no.3-4, 1215-1231.
  * [6] E. Gawrilow and M. Joswig, _polymake: a framework for analyzing convex polytopes,_ Polytopes-combinatorics and computation (Oberwolfach, 1997), 43-73.
  * [7] A. Greenleaf, A. Iosevich, and K. Taylor, On $k$‐point configuration sets with nonempty interior, Mathematika, 68 (2022), no.1, 163-190.
  * [8] A. Greenleaf, A. Iosevich, and K. Taylor, Nonempty interior of configuration sets via microlocal partition optimization, arXiv:2209.02084 (2022).
  * [9] L. Guth, A. Iosevich, Y. Ou and H. Wang, On Falconer’s distance set problem in the plane, Invent. Math. 219 (2020), no.3, 779-830.
  * [10] B. Hanson, B. Lund, and O. Roche-Newton, On distinct perpendicular bisectors and pinned distances in finite fields, Finite Fields Appl. 37 (2016), 240-264.
  * [11] A. Iosevich, M. Rudnev, _Erdős distance problem in vector spaces over finite fields_ , Trans. Amer. Math. Soc. 359 (2007), no.12, 6127-6142.
  * [12] D. Koh and S. Lee, _Averages and maximal averages over product $j$-varieties in finite fields,_ Finite Fields Appl. 79 (2022), Paper No. 101984, 13 pp.
  * [13] D. Koh, T. Pham and L. A. Vinh, Extension theorems and a connection to the Erdos-Falconer distance problem over finite fields, J. Funct. Anal. 281 (2021), no.8, Paper No. 109137, 54 pp.
  * [14] D. Koh and C.Y. Shen, _Extension and averaging operators for finite fields_ , Proc. Edinb. Math. Soc. (2) 56 (2013), no.2, 599-614.
  * [15] D. Koh, C.Y. Shen, and I. Shparlinski, _Averaging operators over homogeneous varieties over finite fields_ , J. Geom. Anal. 26 (2016), no.2, 1415-1441.
  * [16] R. Lidl and H. Niederreiter, _Finite fields,_ Cambridge University Press, (1997).
  * [17] N. Lyall, A. Magyar, and H. Parshals, Spherical configurations over finite fields, Amer. J. Math. 142 (2020), no.2, 373-404.
  * [18] B. Murphy, G. Petridis, T. Pham, M. Rudnev, and S. Stevenson, _On the pinned distances problem in positive characteristic,_ J. Lond. Math. Soc. (2) 105 (2022), no.1, 469-499.
  * [19] Y. Ou and K. Taylor, Finite point configurations and the regular value theorem in a fractal setting, arXiv:2003.06218, accepted in Indiana Journal of Mathematics (2021).
  * [20] M. Rudnev, Note on the number of hinges defined by a point set in $\mathbb{R}^{2}$, Combinatorica, 40 (2020), no.5, 749–757.
  * [21] E. M. Stein, _Harmonic analysis: real-variable methods, orthogonality, and oscillatory integrals,_ With the assistance of Timothy S. Murphy. Princeton Mathematical Series, 43. Monographs in Harmonic Analysis, III. Princeton University Press, Princeton, NJ, 1993. xiv+695 pp.