# FULL-REFERENCE POINT CLOUD QUALITY ASSESSMENT
USING SPECTRAL GRAPH WAVELETS

###### Abstract

Point clouds in 3D applications frequently experience quality degradation
during processing, e.g., scanning and compression. Reliable point cloud
quality assessment (PCQA) is important for developing compression algorithms
with good bitrate-quality trade-offs and techniques for quality improvement
(e.g., denoising). This paper introduces a full-reference (FR) PCQA method
utilizing spectral graph wavelets (SGWs). First, we propose novel SGW-based
PCQA metrics that compare SGW coefficients of coordinate and color signals
between reference and distorted point clouds. Second, we achieve accurate PCQA
by integrating several conventional FR metrics and our SGW-based metrics using
support vector regression. To our knowledge, this is the first study to
introduce SGWs for PCQA. Experimental results demonstrate the proposed PCQA
metric is more accurately correlated with subjective quality scores compared
to conventional PCQA metrics.

Index Terms—  point cloud quality assessment, full reference metric, graph
signal processing, spectral graph wavelet, support vector regression

##  1 Introduction

Point clouds are a general 3D format representing realistic 3D objects in
diverse 3D applications such as telepresence, monitoring, and holographic
display. However, the perceptual quality of point clouds often deteriorates
during scanning, compression, and transmission. Point cloud compression
methods have been proposed to achieve a good trade-off between quality and
bitrate. For instance, the MPEG committee has standardized two approaches:
geometry-based point cloud compression (G-PCC) and video-based point cloud
compression (V-PCC) [1]. Moreover, many methods aiming at improving point
cloud quality have been proposed to mitigate the impact of distortions (e.g.,
point cloud denoising [2], upsampling [3], and inpainting [4]). Accurate point
cloud quality assessment (PCQA) methods are essential to evaluate compression
and quality enhancement methods.

Methods for PCQA can be classified into full-reference point cloud quality
assessment (FR-PCQA) and no-reference point cloud quality assessment (NR-PCQA)
approaches. Unlike NR-PCQA metrics [5, 6, 7], which cannot access a reference
(noise-free) point cloud, FR-PCQA methods show reliable and stable assessment
results when a good reference is available. For example, some of the FR-PCQA
metrics [8, 9] have been utilized as criteria for evaluating point cloud
compression methods in the MPEG standardization [1]. While these methods [8,
9] primarily focus on point-wise errors, alternative approaches that consider
(i) more complex features (e.g., structural similarity [10, 11, 12] or (ii)
graph similarity [13, 14]) have been proposed to improve the correlation with
subjective assessment scores such as mean opinion score (MOS).

In recent years, learning-based FR-PCQA methods [15, 16, 17] have been shown
to obtain better correlation with subjective assessment scores. Our previous
study, full-reference quality assessment using support vector regression
(FRSVR) [17], achieves high correlation with subjective assessment scores by
integrating five different FR-PCQA metrics using support vector regression
(SVR) and achieved first place in the FR broad-range quality estimation track
in the ICIP 2023 point cloud visual quality assessment grand challenge (ICIP
2023 PCVQA grand challenge) [18]. Since it was designed for a grand challenge
focused on PCQA for compression distortion only, FRSVR achieves less
correlation with subjective assessment scores for other types of noise (e.g.,
Gaussian noise, down-sampling). In particular, the five types of scores
adopted by FRSVR do not consider multi-scale features, which, according to
previous studies [12, 14], are effective in assessing various types of
distortion. Specifically, FRSVR utilizes two point-wise error metrics to
evaluate geometric distortion but does not carry out region-to-region
comparisons. Similarly, although region-to-region comparison is performed for
color distortion, this is a single-scale feature calculated from a specified
number of neighbors.

As an alternative, in this paper, we propose a more accurate FR-PCQA method
using a spectral graph wavelet (SGW), which can provide multi-scale features
as graph signal frequencies. In addition, the proposed method is faster than
conventional multi-scale FR-PCQA methods [12, 14] because a fast calculation
technique has been proposed for spectral graph wavelet transform (SGWT) s such
as polynomial approximation [19].

Our contributions are as follows:

  1. 1.

we propose a novel spectral graph wavelet-based point cloud quality assessment
(SGW-PCQA) metric that compares SGW s calculated from the coordinate and color
signals of both reference and distorted point clouds. The previous studies on
point cloud denoising based on SGW s [20, 21, 22] demonstrate that noise in
geometric and color information can be eliminated by suppressing the high-
frequency components of SGW s. Thus, comparing SGW s calculated from reference
and distorted point clouds is useful to evaluate the amount of noise in the
distorted point cloud. As far as we know, this is the first study to introduce
SGW s for PCQA.

  2. 2.

To improve the correlation coefficients with subjective assessment scores, we
use SVR to integrate the SGW-PCQA metrics with two out of the five metrics
introduced in our previous SVR-based method, FRSVR [17]. This integrated
version, SGW-PCQA +, leads to better correlation coefficients for both
compression distortion and other types of distortions, such as Gaussian noise
and downsampling. In addition, SGW-PCQA \+ outperformed the results of our
previous work (FRSVR [17]) recorded in the ICIP 2023 PCVQA grand challenge in
terms of the correlation coefficients with MOS and processing time.

##  2 Related work

Fig. 1: Calculation of spectral graph wavelet transforms (SGWTs) for a point
cloud. Fig. 2: Overall calculation flow of the proposed PCQA metrics using
SGWs.

Since this paper focuses on FR-PCQA, we introduce conventional FR-PCQA methods
in this section.

###  2.1 3D-to-2D projection-based methods

Conventional image quality assessment techniques can be used to assess the
quality of 3D point clouds by projecting them into multiple 2D planes. In
[23], multiple projected 2D images are assessed by popular 2D image quality
assessment metrics (e.g., VIFP [24], SSIM [25]). An attention-guided PCQA
inspired by the structural similarity measure of [26] has also been proposed
[27]. However, these methods suffer from a low correlation with subjective
assessment scores since there is a loss of 3D information during projection.

###  2.2 3D model-based methods

Geometric distortions, e.g., point-to-point error [8], point-to-plane error
[9] and plane-to-plane error [28], have also been introduced, but they cannot
evaluate color distortion, so their applicability for PCQA is limited.

To evaluate color distortion, the Peak Signal-to-Noise Ratio (PSNR) [8] and
PointSSIM [10] have been proposed. However, their correlation with subjective
assessment scores is still limited since these methods compare only one
feature (e.g., color [8] or local color variance [10]). To improve the
reliability of PCQA metrics, point cloud quality metric (PCQM) [11] introduces
a weighted linear combination of geometric and color distortions. Graph-based
FR-PCQA distortions methods[13] and multi-scale methods have also been
proposed to further improve the accuracy [12, 14]. However, the improved
accuracy of [11, 12, 13, 14] comes with a large increase in computational
complexity.

###  2.3 Learning-based methods

Recently, learning-based techniques have been introduced to obtain better
correlation with subjective assessment scores. Unlike previously mentioned
approaches, these methods require training a model using subjective assessment
scores such as MOS. Examples of these approaches include [15], an end-to-end
deep-learning framework for accurate FR-PCQA and [16], a learning-based method
utilizing various PCA-based features. While [15, 16] demonstrate a high
correlation with subjective assessment scores and combine geometry and color
information, they suffer from large processing time for prediction. By using
the SVR-based method [17] to combine simple FR-PCQA metrics by SVR, we achieve
high correlation with subjective assessment scores for compression distortion
with relatively low complexity. However, this approach is less reliable for
other types of noises (e.g., Gaussian noise, down-sampling), as discussed in
Section 1.

Fig. 3: Changes of SGWs through distortions caused by (a) Gaussian noise
(geometry), (b) downsampling, (c) Gaussian noise (color), and (d) compression
error (G-PCC RAHT [1]).

##  3 Preliminaries

###  3.1 Graph and point cloud notations

We begin by introducing mathematical notations and definitions relevant to
graphs. An undirected and weighted graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ is constructed with a point cloud
$\mathrm{P}=\\{p_{i}\\},i=1,...,N$, where $\mathcal{V}$ and $\mathcal{E}$
denote the set of nodes and edges on a graph, respectively. In this paper, a
point cloud $P$ is characterized by the coordinate signals
$\bm{f}_{g}=(\bm{f}_{x},\bm{f}_{y},\bm{f}_{z})\in\mathbb{R}^{N\times 3}$ and
color (lightness) signals $\bm{f}_{L}\in\mathbb{R}^{N}$. As in [11], we
represent color using lightness in the LAB color space [29].

The adjacency matrix $\mathbf{W}=\\{w_{ij}\\}$ of the graph corresponding to a
point cloud $\mathrm{P}$ defines the edge weights between nodes (points) $i$
and $j$ using the Gaussian kernel function as:

$w_{ij}=\left\\{\begin{array}[]{ll}\mathrm{exp}\left(\frac{-||\bm{f}_{g,i}-\bm{f}_{g,j}||_{2}^{2}}{\theta^{2}}\right)&(p_{j}\in\mathcal{N}(p_{i})\textnormal{
or }p_{i}\in\mathcal{N}(p_{j}))\\\ 0&(\mathrm{otherwise}),\end{array}\right.$
(1)

where $\bm{f}_{g,i}\in\mathbb{R}^{3}$ shows the coordinate signals at point
$p_{i}$ and $\theta$ denotes the average of all pairwise distances. We have
that $p_{i}\in\mathcal{N}(p_{j})$ if $p_{i}$ is one of the neighbors of
$p_{j}$. The definition of neighbors is determined and parameterized by a
graph construction algorithm (e.g., K nearest neighbors (KNN)).

The combinatorial graph Laplacian $\mathbf{L}$ is computed by
$\mathbf{L}=\mathbf{D}-\mathbf{W}$, where $\mathbf{D}$ is the degree matrix of
$\mathbf{W}$ [30]. $\mathbf{L}$ characterizes the global smoothness of a graph
signal $\bm{f}=\\{f_{i}\\}\in\mathbb{R}^{N}$ because the graph Laplacian
quadratic form is derived from $|\nabla\bm{f}|^{2}=\sum_{i\sim
j}{w_{ij}(f_{i}-f_{j})^{2}}=\bm{f}^{\top}\mathbf{L}\bm{f},$ where the sum is
over all pairs of connected nodes $i$ and $j$, denoted by $i\sim j$. The graph
Fourier transform $\hat{\bm{f}}$ is defined as
$\hat{f}(\lambda_{l})=\sum_{i}f_{i}\phi_{l,i}$, where $\lambda_{l}$ and
$\phi_{l,i}$ denote the $l$-th eigenvalue and eigenvector of $\mathbf{L}$,
respectively.

###  3.2 Spectral graph wavelet transform (SGWT)

Fig. 1 illustrates the calculation process of SGWT [19, 31] for a point cloud
$\mathrm{P}$. SGWT s are constructed using a kernel operator
$T_{g}=g(\mathbf{L})$, which acts on a graph signal $\bm{f}$ by modulating
each graph Fourier mode: $\hat{T_{g}f(l)}=g(\lambda_{l})\hat{f}(l)$. A scaled
operator $T^{m}_{g}=g(S_{m}\mathbf{L})$ shows the scaling in the spectral
domain at scale $S_{m}$, where $m$ is the index of the band-pass filter, as
shown in Fig. 1. At this time, the wavelets are calculated by applying
$T^{m}_{g}$ to a single vertex operator: $\psi^{i}_{s}=T^{m}_{g}\delta_{i}$
where $\delta_{i}$ is the impulse on vertex $i$. Then, the SGW coefficients
$\mathbf{\Psi}_{f}\in\mathbb{R}^{M\times N}$, where $M$ indicates the number
of band-pass filters, are calculated by

$\Psi_{f}(m,i)=\sum^{N}_{l=1}g(m\lambda_{l})\hat{f}(\lambda_{l})\phi_{l,i}.$
(2)

##  4 Proposed method

###  4.1 Overview of the proposed method

Reference and distorted point clouds are represented as
$\mathrm{P_{R}}=\\{p^{R}_{i}\\},$ $i=1,...,|\mathrm{P_{R}}|$ and
$\mathrm{P_{D}}=\\{p^{D}_{j}\\},$ $j=1,...,|\mathrm{P_{D}}|$, respectively,
where $|\mathrm{P_{R}}|$ denotes the number of points of $\mathrm{P_{R}}$.
Point clouds $\mathrm{P_{R}}$ and $\mathrm{P_{D}}$ have three coordinate
signals and one color signal
$(\bm{f}^{R}_{x},\bm{f}^{R}_{y},\bm{f}^{R}_{z},\bm{f}^{R}_{L})$ and
$(\bm{f}^{D}_{x},\bm{f}^{D}_{y},\bm{f}^{D}_{z},\bm{f}^{D}_{L})$, respectively.
Fig. 2 shows the overall calculation flow of the proposed method called SGW-
PCQA, which includes the following five steps.

  1. 1.

A graph $\mathcal{G}_{R}$ with adjacency matrix $\mathbf{W_{R}}$ is
constructed from $\mathrm{P_{R}}$ by KNN.

  2. 2.

Point-to-point correspondence between $\mathrm{P_{R}}$ and $\mathrm{P_{D}}$ is
calculated. The nearest point in $\mathrm{P_{D}}$ from $p^{R}_{i}$ is selected
as $\tilde{p}^{D}_{i}$ by the nearest neighbor search using Euclidean
distance. The associated distorted point cloud
$\mathrm{\tilde{P}_{D}}=\\{\tilde{p}^{D}_{i}\\}$ has the signals
$(\tilde{\bm{f}}^{D}_{x},\tilde{\bm{f}}^{D}_{y},\tilde{\bm{f}}^{D}_{z},\tilde{\bm{f}}^{D}_{L})$.

  3. 3.

The SGW coefficients with graph $\mathcal{G}_{R}$ are calculated for each
graph signal $\bm{f}^{R}_{x}$, $\bm{f}^{R}_{y}$, $\bm{f}^{R}_{z}$,
$\bm{f}^{R}_{L}$, $\tilde{\bm{f}}^{D}_{x}$, $\tilde{\bm{f}}^{D}_{y}$,
$\tilde{\bm{f}}^{D}_{z}$ and $\tilde{\bm{f}}^{D}_{L}$. The resulting SGW
coefficients are $\bm{\Psi}^{R}_{f_{x}}$, $\bm{\Psi}^{R}_{f_{y}}$
,$\bm{\Psi}^{R}_{f_{z}}$, $\bm{\Psi}^{R}_{f_{L}}$,
$\tilde{\bm{\Psi}}^{D}_{f_{x}}$,
$\tilde{\bm{\Psi}}^{D}_{f_{y}},\tilde{\bm{\Psi}}^{D}_{f_{z}}$, and
$\tilde{\bm{\Psi}}^{D}_{f_{L}}$.

  4. 4.

The geometry PCQA score $S^{G}_{m}$ and color PCQA score $S^{G}_{m}$ are
computed by comparing the SGW coefficients at the $m$-th subband.

  5. 5.

A trained SVR model predicts a final assessment score from the scores
$S^{G}_{m}$ and $S^{C}_{m}$. Note that the parameters of the SVR model are
optimized using a training set, including the scores calculated in step 4 and
the MOS of the corresponding distorted point clouds.

We explain steps 1, 2, and 3 in Section 4.2. Steps 4 and 5 are described in
Section 4.3 and Section 4.4, respectively. Furthermore, SGW-PCQA +, which is
an extension of SGW-PCQA, is described in Section 4.5.

###  4.2 Calculation of spectral graph wavelets

This section explains how to calculate SGWs for PCQA. In graph signal
processing, the definition of graph frequency depends on the structure of the
graph [30]. Hence, if graphs are independently constructed from reference and
distorted point clouds, the definitions of their respective graph frequencies
will be different. Therefore, to perform PCQA, if we wish to compare the graph
frequencies obtained from $P_{R}$ and $P_{D}$, we first must select a common
graph for the two point clouds.

In the proposed method, a graph $\mathcal{G}_{R}$ with the adjacency matrix
$\mathbf{W_{R}}$ is calculated from the reference point cloud $\mathrm{P_{R}}$
by KNN, using (1). Then, the nearest point
$\tilde{p}^{D}_{i}\in\mathrm{P_{D}}$ from a query point $p^{R}_{i}$ is
extracted by the nearest neighbor search based on the 3D Euclidean distance to
project $\mathrm{P_{D}}$ into $\mathrm{P_{R}}$. In this way, we obtain a new
distorted point cloud $\mathrm{\tilde{P}_{D}}=\\{\tilde{p}^{D}_{i}\\}$ that
has the same node indices but different graph signals. Finally, the SGW
coefficients are calculated from coordinate and color signals of
$\mathrm{P_{R}}$ and $\mathrm{\tilde{P}_{D}}$ based on (2). Since SGWT s are
computed for each graph signal
($\bm{f}^{R}_{x},\bm{f}^{R}_{y},\bm{f}^{R}_{z},\bm{f}^{R}_{L},\tilde{\bm{f}}^{D}_{x},\tilde{\bm{f}}^{D}_{y},\tilde{\bm{f}}^{D}_{z}$
and $\tilde{\bm{f}}^{D}_{L}$), SGWT coefficients
($\bm{\Psi}^{R}_{f_{x}},\bm{\Psi}^{R}_{f_{y}},\bm{\Psi}^{R}_{f_{z}},\bm{\Psi}^{R}_{f_{L}},\tilde{\bm{\Psi}}^{D}_{f_{x}},\tilde{\bm{\Psi}}^{D}_{f_{y}},\tilde{\bm{\Psi}}^{D}_{f_{z}}$,
and $\tilde{\bm{\Psi}}^{D}_{f_{L}}$) are output.

###  4.3 SGW-based FR-PCQA metrics

Fig. 3 shows the results of our preliminary experiment, which confirmed that
SGWs change under several types of distortions. As shown in Fig. 3, since SGW
s at various subbands are changed when distortions are introduced, the quality
can be assessed by observing the changes of SGW s.

In the proposed SGW-PCQA method, to assess geometric distortion, the
differences in SGWs of $x,y$, and $z$ coordinates for each $m$-th subband are
calculated as follows:

$G^{x}_{m}=\frac{1}{|P_{R}|}\sum_{p^{R}_{i}\in\mathrm{P_{R}}}(\Psi^{R}_{f_{x}}(m,i)-\tilde{\Psi}^{D}_{f_{x}}(m,i))^{2},$
(3)
$G^{y}_{m}=\frac{1}{|P_{R}|}\sum_{p^{R}_{i}\in\mathrm{P_{R}}}(\Psi^{R}_{f_{y}}(m,i)-\tilde{\Psi}^{D}_{f_{y}}(m,i))^{2},$
(4)
$G^{z}_{m}=\frac{1}{|P_{R}|}\sum_{p^{R}_{i}\in\mathrm{P_{R}}}(\Psi^{R}_{f_{z}}(m,i)-\tilde{\Psi}^{D}_{f_{z}}(m,i))^{2},$
(5) $G_{m}=\frac{G^{x}_{m}+G^{y}_{m}+G^{z}_{m}}{3}.$ (6)

Geometric assessment scores $S^{G}_{m}$ are calculated as

$S^{G}_{m}=\frac{1}{1+G_{m}}.$ (7)

Likewise, color assessment scores $S^{C}_{m}$ are given as

$C_{m}=\frac{1}{|P_{R}|}\sum_{p^{R}_{i}\in\mathrm{P_{R}}}(\Psi^{R}_{f_{L}}(m,i)-\tilde{\Psi}^{D}_{f_{L}}(m,i))^{2},$
(8) $S^{C}_{m}=\frac{1}{1+C_{m}}.$ (9)

###  4.4 Score integration by support vector regression (SVR)

We use SVR to predict a final PCQA score. Multiple scores, $S^{G}_{m}$ and
$S^{C}_{m}$, are integrated by SVR. During the training phase, an SVR model is
trained using the MOS of the corresponding distorted point clouds and the
multiple PCQA scores calculated from training data. The Gaussian radial basis
function (RBF) kernel and the sequential minimal optimization [32] are chosen
as a kernel and solver of SVR, respectively.

###  4.5 Extension of SGW-PCQA (SGW-PCQA +)

In our preliminary experiments, the performance of SGW-PCQA can be further
improved by adding some of the FR-PCQA scores utilized in our previous work,
FRSVR [17]. Thus, we selected two types of scores utilized in FRSVR, (1)
point-to-point score ($S_{\mathrm{p2p}}$) and (2) graph total variation
(GTV)-based score ($S_{\mathrm{gtv}}$), whose addition significantly improved
the correlation with subjective scores. The point-to-point score
$S_{\mathrm{p2p}}$ is based on point-to-point error [8]. The GTV-based score
$S_{\mathrm{gtv}}$ is calculated by comparing GTV s, which quantify the
smoothness of graph signals (Please see [17] to access the detailed formulas
due to space limitation).

Although SGW s reflect graph signal variation at various scales, they cannot
consider too local and global information. Since the point-to-point and GTV-
based scores focus on very local (point-wise) and global information,
respectively, adding these scores leads to improved correlation with
subjective scores. In addition, although graph construction is required to
calculate these two scores, a graph utilized for calculating SGW s can be
reused. Thus, the additional processing time is not very large. The extension
that combines these two scores ($S_{\mathrm{p2p}}$ and $S_{\mathrm{gtv}}$)
with the SGW-based scores ($S^{G}_{m}$ and $S^{C}_{m}$) using SVR is named as
SGW-PCQA +.

##  5 EXPERIMENTS

Fig. 4: SROCCs of each SGW-PCQA score, (a) geometry-based score ($S^{G}_{m}$), and (b) color-based score ($S^{C}_{m}$), with the BASICS dataset (test set). Table 1: The correlation coefficients of the proposed SGW-PCQA method. Some of the scores were selected to train an SVR model. The blue characters mean the best result. # | Geometry score set | Color score set |  BASICS [33] |  WPC [34] | Average
---|---|---|---|---|---
PLCC | SROCC | PLCC | SROCC | PLCC | SROCC
#1 | $S^{G}_{1}$ | $S^{C}_{1}$ | 0.844 | 0.770 | 0.787 | 0.780 | 0.816 | 0.775
#2 |  $S^{G}_{1}$, $S^{G}_{2}$ |  $S^{C}_{2}$,$S^{C}_{2}$ | 0.880 | 0.815 | 0.835 | 0.815 | 0.858 | 0.815
#3 |  $S^{G}_{1}$, $S^{G}_{2}$, $S^{G}_{3}$ |  $S^{C}_{1}$, $S^{C}_{2}$, $S^{C}_{3}$ | 0.874 | 0.837 | 0.854 | 0.842 | 0.864 | 0.840
#4 |  $S^{G}_{1}$, $S^{G}_{2}$, $S^{G}_{3}$, $S^{G}_{4}$ |  $S^{C}_{1}$, $S^{C}_{2}$, $S^{C}_{3}$, $S^{C}_{4}$ | 0.870 | 0.837 | 0.863 | 0.854 | 0.867 | 0.846
#5 |  $S^{G}_{1}$, $S^{G}_{2}$, $S^{G}_{3}$, $S^{G}_{4}$, $S^{G}_{5}$ |  $S^{C}_{1}$, $S^{C}_{2}$, $S^{C}_{3}$, $S^{C}_{4}$, $S^{C}_{5}$ | 0.875 | 0.835 | 0.847 | 0.834 | 0.861 | 0.835
#6 |  $S^{G}_{1}$, $S^{G}_{2}$, $S^{G}_{3}$, $S^{G}_{4}$, $S^{G}_{5}$, $S^{G}_{6}$ |  $S^{C}_{1}$, $S^{C}_{2}$, $S^{C}_{3}$, $S^{C}_{4}$, $S^{C}_{5}$, $S^{C}_{6}$ | 0.889 | 0.843 | 0.826 | 0.814 | 0.858 | 0.829
#7 |  $S^{G}_{1}$, $S^{G}_{2}$, $S^{G}_{3}$, $S^{G}_{4}$, $S^{G}_{5}$, $S^{G}_{6}$ |  $S^{C}_{1}$, $S^{C}_{2}$, $S^{C}_{3}$ | 0.892 | 0.861 | 0.845 | 0.840 | 0.869 | 0.851
#8 |  $S^{G}_{1}$, $S^{G}_{2}$, $S^{G}_{3}$, $S^{G}_{4}$, $S^{G}_{5}$, $S^{G}_{6}$ |  $S^{C}_{1}$, $S^{C}_{2}$, $S^{C}_{3}$, $S^{C}_{4}$ | 0.888 | 0.859 | 0.842 | 0.837 | 0.865 | 0.848
#9 |  $S^{G}_{1}$, $S^{G}_{2}$, $S^{G}_{3}$, $S^{G}_{4}$, $S^{G}_{5}$, $S^{G}_{6}$ |  $S^{C}_{1}$, $S^{C}_{2}$, $S^{C}_{3}$, $S^{C}_{4}$, $S^{C}_{5}$ | 0.879 | 0.835 | 0.828 | 0.823 | 0.854 | 0.829
Table 2: The correlation coefficients and processing time [s] of the proposed and conventional methods. The “Average” indicates the average of the two datasets. The blue and red characters mean the first and second best results, respectively. Metric |  BASICS [33] |  WPC [34] | Average
---|---|---|---
PLCC | SROCC | Time [s] | PLCC | SROCC | Time [s] | PLCC | SROCC | Time [s]
Point-to-point (MSE) [8] | 0.041 | 0.735 | 9.701 | 0.401 | 0.566 | 6.832 | 0.221 | 0.651 | 8.267
Point-to-plane (MSE) [8] | 0.005 | 0.799 | 49.122 | 0.368 | 0.481 | 32.847 | 0.187 | 0.640 | 40.985
Angular similarity [28] | 0.328 | 0.306 | 39.155 | 0.296 | 0.319 | 27.336 | 0.312 | 0.313 | 33.246
Y-MSE [8] | 0.512 | 0.550 | 9.896 | 0.469 | 0.591 | 6.928 | 0.491 | 0.571 | 8.412
PointSSIM [10] | 0.605 | 0.620 | 23.682 | 0.469 | 0.471 | 15.281 | 0.537 | 0.546 | 19.482
PCQM [11] | 0.786 | 0.739 | 256.463 | 0.511 | 0.550 | 214.197 | 0.649 | 0.645 | 235.330
GraphSIM [13] | 0.801 | 0.773 | 613.652 | 0.688 | 0.691 | 382.798 | 0.745 | 0.732 | 498.225
MSGraphSIM [14] | 0.807 | 0.773 | 645.589 | 0.712 | 0.724 | 348.183 | 0.760 | 0.749 | 496.886
FRSVR [17] | 0.914 | 0.878 | 32.722 | 0.818 | 0.803 | 20.608 | 0.866 | 0.841 | 26.665
SGW-PCQA (Proposed, #7 in Table 1) | 0.892 | 0.861 | 27.387 | 0.845 | 0.840 | 15.526 | 0.869 | 0.851 | 21.457
SGW-PCQA+ (Proposed) | 0.924 | 0.888 | 30.082 | 0.873 | 0.869 | 17.266 | 0.899 | 0.879 | 23.674

###  5.1 Experimental conditions

Evaluation Criteria: To quantify the correlation to the subjective assessment
score, we employed Pearson’s linear correlation coefficient (PLCC) and
Spearman’s rank-order correlation coefficient (SROCC).

Dataset: We used the two datasets, the broad quality assessment of static
point clouds in compression scenario (BASICS) [33], consisting of 1498 point
clouds distorted by compression errors, and the Waterloo point cloud (WPC)
[34] datasets, consisting of 740 distorted point clouds contaminated by
compression, Gaussian noise, or down-sampling. The two datasets include the
MOS of the corresponding distorted point clouds to be utilized to measure PLCC
and SROCC. Since the BASICS dataset was explicitly divided into training and
test sets [33], we conducted training and evaluation with the training and
test sets, respectively. As for the WPC dataset, distorted point clouds were
segmented into five parts (20% data $\times$ 5 parts). Each segment served as
an independent test set, while the remaining data were utilized as training
data. After calculating PLCC and SROCC for each test set, the average was
shown as the evaluation result.

Machine specifications: The processing time was measured with a computer
equipped with an AMD Ryzen Threadripper 2970WX 24-Core processor, NVIDIA GTX
1080 Ti, and 128GB Random Access Memory.

Implementation details: The graph construction parameter $K$ and the number of
wavelet decomposition levels $M$ were set to 8 and 6, respectively. Since many
of the conventional methods compared in this paper do not employ GPU computing
for acceleration, GPU computing is not utilized for the implementation of the
proposed method for a fair comparison.

###  5.2 Experiment 1: Performance of a single score

We evaluated the performance of a single SGW-PCQA score, $S^{G}_{m}$ and
$S^{C}_{m}$ (see Fig.4). Since the SROCCs were calculated from only one score
in this experiment, SVR was not utilized for integrating multiple scores. As
shown in Fig. 4, the SROCC of color-based score $S^{C}_{m}$ was lower for the
high-frequency subbands. For the high frequencies of color information, there
were many cases where the differences between reference and distorted point
clouds were small. We observed that the color signals of a reference point
cloud had very few high-frequency components, even for general point clouds
other than a simple point cloud shown in Fig.3. Since small high-frequency
components are quantized to zero, the compression process does not add high-
frequency components. Comparing high-frequency components is not effective in
compression noise of color signals, which leads to low SROCCs. Thus, the SROCC
of the BASICS dataset consisting only of compression errors was particularly
low compared to that of the WPC dataset.

###  5.3 Experiment 2: Performance of the combination of multiple scores

We investigated the performance achieved by combining multiple scores using
SVR. Table 1 shows the PLCC and SROCC of the two datasets. Condition #7 in
Table 1, where all geometry scores and three color scores ($S^{C}_{1}$,
$S^{C}_{2}$, and $S^{C}_{3}$) were selected, recorded the best average PLCC
and SROCC. Since the SROCCs of color-based scores at high-frequency bands were
relatively low as shown in Fig. 4, using the color-based scores $S^{C}_{m}$ at
high-frequency bands harmed the correlation with subjective assessment scores.

###  5.4 Experiment 3: comparison with SOTA methods

We compared the proposed method with conventional FR-PCQA methods [8, 9, 10,
11, 13, 14, 17, 28]. Table 2 presents the results of the comparison experiment
with PLCC, SROCC, and average processing time across all the point clouds.
SGW-PCQA and SGW-PCQA \+ in Table 2 adopted condition #7 in Table 1 because it
achieved the best results in Section 5.3. The results demonstrate that the
proposed SGW-PCQA \+ achieved the highest PLCC and SROCC. In particular, the
PLCC and SROCC of the WPC dataset were improved against our previous study,
FRSVR [17]. The PCQA using SGW s were effective with various noise, such as
Gaussian noise and down-sampling. For compression noise, the five metrics
adopted in FRSVR [17] provided sufficiently high accuracy; however, those
metrics can be insufficient for diverse types of noise. In addition, since
SGWT s are calculated relatively fast by polynomial approximation [19], SGW-
PCQA performed faster computation compared to conventional multiscale metric
(e.g., MSGraphSIM [14]).

##  6 CONCLUSION

This paper introduced an accurate FR-PCQA method using SGW s. The proposed
metric demonstrates a great correlation by referring to the SGW s calculated
from the coordinate and color signals. In addition, we made the method more
accurate by integrating our previous method called FRSVR with the SGW-PCQA
metrics by using SVR. In the future, we will consider a more suitable graph
construction method for PCQA instead of KNN, which is adopted in the proposed
method. Since the previous denoising study [22] has shown that the values of
SGW s highly depend on the graph structure, there may be a more suitable graph
construction for PCQA.

## ACKNOWLEDGEMENTS

This work was supported by the Ministry of Internal Affairs and Communications
(MIC) of Japan (Grant no. JPJ000595).

## References

  * [1] S. Schwarz, M. Preda, V. Baroncini, M. Budagavi, P. Cesar, P. A. Chou, R. A. Cohen, M. Krivokuća, S. Lasserre, Z. Li, J. Llach, K. Mammou, R. Mekuria, O. Nakagami, E. Siahaan, A. Tabatabai, A. M. Tourapis, and V. Zakharchenko,  “Emerging MPEG standards for point cloud compression,”  IEEE Journal on Emerging and Selected Topics in Circuits and Systems, vol. 9, no. 1, pp. 133–148, 2019.
  * [2] Y. Zhou, R. Chen, Y. Zhao, X. Ai, and G. Zhou,  “Point cloud denoising using non-local collaborative projections,”  Pattern Recognition, vol. 120, pp. 108128, 2021.
  * [3] A. Akhtar, Z. Li, G. V. Auwera, L. Li, and J. Chen,  “PU-Dense: Sparse tensor-based point cloud geometry upsampling,”  IEEE Transactions on Image Processing, vol. 31, pp. 4133–4148, 2022.
  * [4] J. He, Z. Fu, W. Hu, and Z. Guo,  “Point cloud attribute inpainting in graph spectral domain,”  in 2019 IEEE International Conference on Image Processing, 2019, pp. 4385–4389.
  * [5] Q. Liu, H. Yuan, H. Su, H. Liu, Y. Wang, H. Yang, and J. Hou,  “PQA-Net: Deep no reference point cloud quality assessment via multi-view projection,”  IEEE Transactions on Circuits and Systems for Video Technology, vol. 31, no. 12, pp. 4645–4660, 2021.
  * [6] Z. Zhang, W. Sun, X. Min, T. Wang, W. Lu, and G. Zhai,  “No-reference quality assessment for 3D colored point cloud and mesh models,”  IEEE Transactions on Circuits and Systems for Video Technology, vol. 32, no. 11, pp. 7618–7631, 2022.
  * [7] A. Chetouani, M. Quach, G. Valenzise, and F. Dufaux,  “Deep learning-based quality assessment of 3D point clouds without reference,”  in 2021 IEEE International Conference on Multimedia and Expo Workshops, 2021, pp. 1–6.
  * [8] R. Mekuria, Z. Li, C. Tulvan, and P. A. Chou,  “Evaluation criteria for PCC (point cloud compression),”  in ISO/IEC JTC 1/SC29/WG11 Doc., 2016,  N16332.
  * [9] D. Tian, H. Ochimizu, C. Feng, R. Cohen, and A. Vetro,  “Geometric distortion metrics for point cloud compression,”  in 2017 IEEE International Conference on Image Processing, 2017, pp. 3460–3464.
  * [10] E. Alexiou and T. Ebrahimi,  “Towards a point cloud structural similarity metric,”  in 2020 IEEE International Conference on Multimedia Expo Workshops, 2020, pp. 1–6.
  * [11] G. Meynet, Y. Nehmé, J. Digne, and G. Lavoué,  “PCQM: A full-reference quality metric for colored 3D point clouds,”  in 2020 Twelfth International Conference on Quality of Multimedia Experience, 2020, pp. 1–6.
  * [12] D. Lazzarotto and T. Ebrahimi,  “Towards a multiscale point cloud structural similarity metric,”  in 25th International Workshop On Multimedia Signal Processing, 2023.
  * [13] Q. Yang, Z. Ma, Y. Xu, Z. Li, and J. Sun,  “Inferring point cloud quality via graph similarity,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 44, no. 6, pp. 3015–3029, 2022.
  * [14] Y. Zhang, Q. Yang, and Y. Xu,  “MS-GraphSIM: Inferring point cloud quality via multiscale graph similarity,”  in 29th ACM International Conference on Multimedia, 2021, p. 1230–1238.
  * [15] M. Tliba, A. Chetouani, G. Valenzise, and F. Dufaux,  “Point cloud quality assessment using cross-correlation of deep features,”  in Proceedings of the 2nd Workshop on Quality of Experience in Visual Multimedia Applications, 2022, p. 63–68.
  * [16] X. Zhou, E. Alexiou, I. Viola, and P. Cesar,  “PointPCA+: Extending PointPCA objective quality assessment metric,”  in 2023 IEEE International Conference on Image Processing Challenges and Workshops, 2023, pp. 1–5.
  * [17] R. Watanabe, S. N. Sridhara, H. Hong, E. Pavez, and A. Ortega,  “ICIP 2023 challenge: Full-reference and non-reference point cloud quality assessment methods with support vector regression,”  in 2023 IEEE International Conference on Image Processing Challenges and Workshops, 2023, pp. 3654–3658.
  * [18] A. Chetouani, G. Valenzise, A. Ak, E. Zerman, M. Quach, M. Tliba, M. A. Kerkouri, and P. L. Callet,  “ICIP 2023: Point cloud visual quality assessment grand challenge,” 2023,  https://sites.google.com/view/icip2023-pcvqa-grand-challenge/.
  * [19] D. K. Hammond, P. Vandergheynst, and R. Gribonval,  “Wavelets on graphs via spectral graph theory,”  Applied and Computational Harmonic Analysis, vol. 30, no. 2, pp. 129–150, 2011.
  * [20] S. Deutsch, A. Ortega, and G. Medioni,  “Manifold denoising based on spectral graph wavelets,”  in 2016 IEEE International Conference on Acoustics, Speech and Signal Processing, 2016, pp. 4673–4677.
  * [21] M. A. Irfan and E. Magli,  “Joint geometry and color point cloud denoising based on graph wavelets,”  IEEE Access, vol. 9, pp. 21149–21166, 2021.
  * [22] R. Watanabe, K. Nonaka, E. Pavez, T. Kobayashi, and A. Ortega,  “Graph wavelet-based point cloud geometric denoising with surface-consistent non-negative kernel regression,”  in 2023 IEEE International Conference on Acoustics, Speech and Signal Processing, 2023, pp. 1–5.
  * [23] E. Alexiou and T. Ebrahimi,  “Exploiting user interactivity in quality assessment of point cloud imaging,”  in 2019 Eleventh International Conference on Quality of Multimedia Experience, 2019, pp. 1–6.
  * [24] H. R. Sheikh and A. C. Bovik,  “Image information and visual quality,”  IEEE Transactions on Image Processing, vol. 15, no. 2, pp. 430–444, 2006.
  * [25] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli,  “Image quality assessment: from error visibility to structural similarity,”  IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.
  * [26] Z. Wang and Q. Li,  “Information content weighting for perceptual image quality assessment,”  IEEE Transactions on Image Processing, vol. 20, no. 5, pp. 1185–1198, 2011.
  * [27] Q. Liu, H. Su, Z. Duanmu, W. Liu, and Z. Wang,  “Perceptual quality assessment of colored 3D point clouds,”  IEEE Transactions on Visualization and Computer Graphics, vol. 29, no. 8, pp. 3642–3655, 2023.
  * [28] E. Alexiou and T. Ebrahimi,  “Point cloud quality assessment metric based on angular similarity,”  in 2018 IEEE International Conference on Multimedia and Expo, 2018, pp. 1–6.
  * [29] I. Lissner and P. Urban,  “Toward a unified color space for perception-based image processing,”  IEEE Transactions on Image Processing, vol. 21, no. 3, pp. 1153–1168, 2012.
  * [30] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst,  “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,”  IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 83–98, 2013.
  * [31] N. Leonardi and D. Van De Ville,  “Tight wavelet frames on multislice graphs,”  IEEE Transactions on Signal Processing, vol. 61, no. 13, pp. 3357–3367, 2013.
  * [32] J. Platt,  “Sequential minimal optimization: A fast algorithm for training support vector machines,”  Tech. Rep. MSR-TR-98-14, Microsoft, 1998.
  * [33] A. Ak, E. Zerman, M. Quach, A. Chetouani, A. Smolic, G. Valenzise, and P. L. Callet,  “Basics: Broad quality assessment of static point clouds in a compression scenario,”  IEEE Transactions on Multimedia, pp. 1–13, 2024.
  * [34] H. Su, Z. Duanmu, W. Liu, Q. Liu, and Z. Wang,  “Perceptual quality assessment of 3D point clouds,”  in 2019 IEEE International Conference on Image Processing, 2019, pp. 3182–3186.


# A Model Problem for First Order Mean Field Games with Discrete Initial Data

P. Jameson Graber  J. Graber: Baylor University, Department of Mathematics;
Sid Richardson Building
1410 S. 4th Street
Waco, TX 76706
Tel.: +1-254-710-
Fax: +1-254-710-3569<EMAIL_ADDRESS>and  Brady Zimmerman  B.
Zimmerman: Baylor University,
Sid Richardson Building
1410 S. 4th Street
Waco, TX 76706<EMAIL_ADDRESS>
###### Abstract.

In this article, we study a simplified version of a density-dependent first-
order mean field game, in which the players face a penalization equal to the
population density at their final position. We consider the problem of finding
an equilibrium when the initial distribution is a discrete measure. We show
that the problem becomes finite-dimensional: the final piecewise smooth
density is completely determined by the weights and positions of the initial
measure. We establish existence and uniqueness of a solution using classical
fixed point theorems. Finally, we show that Newton’s method provides an
effective way to compute the solution. Our numerical simulations provide an
illustration of how density penalization in a mean field game tends to the
smoothen the initial distribution.

###### Key words and phrases:

mean field games, optimal control

###### 2020 Mathematics Subject Classification:

49N80, 35Q89

This research was supported in part by National Science Foundation thought NSF
Grant DMS-2045027.

##  1\. Introduction

Mean field games were introduced by Lasry and Lions [LL06a, LL06b, LL07] and
Caines, Huang, and Malhamé [HMC06] in order to model the interactions between
large numbers of agents playing a differential game; see for instance [BFY13,
CD17, G+14, ACD+20] for an overview of the field. In this article, we are
interested in first-order, or deterministic, mean field games. Typically,
these are modeled using a first-order coupled system of partial differential
equations, as follows:

(1.1) $\begin{cases}-\partial_{t}u+H(x,\nabla_{x}u)=F(x,m),\quad
x\in{\mathbb{R}}^{d},t\in(0,T),\\\
\partial_{t}m-\nabla_{x}\cdot\mathinner{\left(D_{p}H(x,\nabla_{x}u)m\right)}=0,\quad
x\in{\mathbb{R}}^{d},t\in(0,T),\\\ m(x,0)=m_{0}(x),\quad
x\in{\mathbb{R}}^{n},\\\ u(x,T)=G\mathinner{\left(x,m(x,T)\right)},\quad
x\in{\mathbb{R}}^{n}.\end{cases}$

Here $H$ is the Hamiltonian for an optimal control problem, which a
representative agent solves for a given population density $m(x,t)$. The
functions $F$ and $G$ model the cost this player must pay as a function of the
density $m$. In Nash equilibrium, the density will evolve according to the
continuity equation driven by the optimal feedback vector field
$-D_{p}H(x,\nabla_{x}u)$.

There are now many results on the existence and uniqueness of suitably defined
_weak_ solutions to System (1.1) [Car15, Gra14, CG15, CGPT15] as well as
regularity of solutions [PS17, GM18, Mun22]. As far as we know, all of these
results apply only when the initial distribution $m_{0}$ is a density, i.e. it
has no singular part with respect to Lebesgue measure. It remains an open
question how to analyze weak solutions in the case where $m_{0}$ is replaced
with a more general measure, which could have singular parts.

In this article, we contribute to the theory of mean field games by analyzing
an example in which the initial measure is entirely discrete, i.e. it is a
weighted sum of Dirac measures. We simplify the problem by choosing $F=0$,
$H(x,p)=\frac{1}{2}p^{2}$, and $G(x,m(x,T))=m(x,T)$. We will also set the
space dimension $d=1$. In this case the problem can be viewed as a static,
“one-shot” game in which a player starting at a point $x$ only has to choose a
point $y$ to minimize the cost

$\frac{\mathinner{\left(x-y\right)}^{2}}{2T}+m(y,T).$

There is nothing particularly important about the value of $T$, so we will set
$T=1/2$, and we will denote by $f(y)$ the final density $m(y,T)$. Hence the
cost to a player starting at point $x$ and choosing to move to point $y$ is
given by

(1.2) $J(x,y,f)=(x-y)^{2}+f(y),$

where $f$ can in principle be any probability density function, i.e. a non-
negative measurable function such that
$\int_{-\infty}^{\infty}f(x)\operatorname{d\\!}x=1$. The initial distribution
will be given by $\displaystyle{m=\sum_{j=1}^{N}a_{j}\delta_{x_{j}}}$ with
$a_{1}+\cdots+a_{N}=1,a_{j}\geq 0$ and $x_{1},\ldots,x_{N}$ distinct points in
${\mathbb{R}}$. Thus $m$ is an empirical measure in which $a_{j}$ is the
proportion of players initially located at each $x_{j}$.

In this context, a _Nash equilibrium_ is a measure $\pi$ on
$\mathbb{R}\times\mathbb{R}$ such that:

  1. (1)

$\pi(A\times\mathbb{R})=m(A)=\sum_{j=1}^{N}a_{j}\delta_{x_{j}}(A),\text{where
}\delta_{x_{j}}(A)=$ $\begin{cases}1\text{ if }x_{j}\notin A\\\ 0\text{ if
}x_{j}\in A\end{cases}$, i.e. first marginal of $\pi$ is $m$,

  2. (2)

$\pi{(\mathbb{R}\times B)}=\int_{B}f(x)\operatorname{d\\!}{x}$ for a density
$f,$ i.e. the second marginal of $\pi$ is the density $f$, and

  3. (3)

for $\pi$ \- a.e. $(x,y)$, $y\in\operatorname{argmin}{J(x,\cdot,f)},$ i.e.
$\pi$ couples points $x$ and $y$ such that $y$ is an optimal move starting
from $x$.

In other words, given the cost function $J$, a density function $f$ is chosen
to possibly be an equilibrium measure. With the knowledge of this function
$f$, the players choose the optimal strategy for themselves, resulting in some
final distribution. If this final distribution matches $f$, then $f$ is in
fact an equilibrium. The idea of a “Self-Fulfilling Prophecy” can be helpful
to understand this: $f$ is “prophesied” to be an equilibrium measure, and, if
it turns out to result in an equilibrium, then $f$ in a sense “fulfilled its
prophecy.” The below diagram demonstrates this phenomenon.

Figure 1. Mean Field Games Diagram

In this article, our purpose is two-fold:

  1. (1)

prove there exists a unique Nash equilibrium, and

  2. (2)

compute solutions using a numerical method.

As far as we know, this constitutes the first result on existence and
uniqueness of solutions to first-order mean field game with singular measures.
Our approach relies on the fact that we take an entirely discrete measure,
which allows us to rewrite the problem in an equivalent finite-dimensional
formulation. Intuitively, all we need to compute is how each Dirac mass
$a_{j}\delta_{x_{j}}$ in the initial distribution will “fan out” into a
density of total mass $a_{j}$ centered around point $x_{j}$. As explained
below, the precise shape of this density is known through a priori
considerations. This approach makes the problem amenable to classical methods.
We believe it also provides geometric insight for what is going on more
generally in density-penalized mean field games of first order.

######  Remark 1.1.

The main results of this article were first announced in the second author’s
undergraduate thesis [Zim24].

##  2\. Reformulation of the equilibrium problem

Let $\pi$ be an equilibrium, and let $f$ be the corresponding final density.

######  Definition 2.1.

The support of $\pi$, supp $\pi$, is the set of all $(x,y)$ such that, for all
$\epsilon>0$, $\pi(B_{\epsilon}(x,y))>0,$ where
$B_{\epsilon}(x,y)=\mathinner{\left\\{(z,w)|\sqrt{(x-z)^{2}+(y-w)^{2}}<\epsilon\right\\}}$,
or, equivalently, $\pi(U)>0$ for every open set containing $(x,y)$.

The support of $f$, $\operatorname{supp}f$, is the closure of the set of all
$y$ such that $f(y)>0$.

######  Lemma 2.2.

If $(x,y)\in\operatorname{supp}\pi$, then $x=x_{j}$ for some $j$ and
$y\in\operatorname{supp}f$. Conversely, if $y\in\operatorname{supp}f$, then
there exists $x_{j}$ such that $(x_{j},y)\in\operatorname{supp}\pi$.

###### Proof.

If $(x,y)\in\operatorname{supp}\pi$, then for all $\epsilon>0$ we have
$m(B_{\epsilon}(x))=\pi(B_{\epsilon}(x)\times{\mathbb{R}})>0$ and
$\int_{B_{\epsilon}(y)}f(z)\operatorname{d\\!}z=\pi({\mathbb{R}}\times
B_{\epsilon}(y))>0$. It follows that $x=x_{j}$ for some $j$, and
$y\in\operatorname{supp}f$ since for every $\epsilon>0$ there exists $z\in
B_{\epsilon}(y)$ such that $f(z)>0$.

Conversely, suppose $y\in\operatorname{supp}f$. For every $\epsilon>0$,
$\pi({\mathbb{R}}\times B_{\epsilon}(y))=\int_{B_{\epsilon}(y)}f>0$. Since
$\pi({\mathbb{R}}\times B_{\epsilon}(y))=\sum_{j=1}^{n}\pi(\\{x_{j}\\}\times
B_{\epsilon}(y))$, for each $\epsilon$ there is a $j$ such that
$\pi(\\{x_{j}\\}\times B_{\epsilon}(y))>0$. We can then find a sequence
$\epsilon_{k}\downarrow 0$ and a fixed $j$ such that $\pi(\\{x_{j}\\}\times
B_{\epsilon_{k}}(y))>0$ for all $k$, which implies
$(x_{j},y)\in\operatorname{supp}\pi$. ∎

Now if $(x_{j},y)\in\operatorname{supp}\pi$, the definition of equilibrium
implies $y\in E_{j}$ where

$E_{j}=\operatorname{argmin}\mathinner{\left\\{y\mathrel{\mathop{\ordinarycolon}}\displaystyle{(x_{j}-y)^{2}}+f(y)\right\\}},$

i.e. $E_{j}$ is the set of minimizers $y$ for the cost $J(x_{j},y,f)$. By
Lemma 2.2, $\operatorname{supp}f\subset\bigcup_{j=1}^{N}E_{j}$; the intuitive
meaning is that every player must move into one of the sets $E_{j}$. On the
other hand, for $y\in E_{j}$ we have $f(y)=C_{j}-{(x_{j}-y)^{2}}$, where

$C_{j}=\min\mathinner{\left\\{y\mathrel{\mathop{\ordinarycolon}}\displaystyle{(x_{j}-y)^{2}}+f(y)\right\\}}.$

Thus, $f$ is completely determined by $C_{j}$ and $E_{j}$, and these in turn
are coupled by the definition of $E_{j}$. So how do we determine $C_{j}$ and
$E_{j}$? A first clue is the following proposition, whose proof is elementary:

######  Proposition 2.3.

$E_{j}=\mathinner{\left\\{y\mathrel{\mathop{\ordinarycolon}}C_{j}-{(x_{j}-y)^{2}}\geq\mathinner{\left(C_{k}-{(x_{k}-y)^{2}}\right)}_{+}\
\forall k\right\\}}$, where
$x_{+}\mathrel{\mathop{\ordinarycolon}}=\max\\{x,0\\}$.

###### Proof.

If $y\in E_{j}$, then ${(x_{j}-y)^{2}}+f(y)=C_{j}$, and at the same time
${(x_{k}-y)^{2}}+f(y)\geq C_{k}$ for all $k$ by definition of $C_{k}$.
Recalling that $f\geq 0$, we see that
$C_{j}-{(x_{j}-y)^{2}}=f(y)\geq\mathinner{\left(C_{k}-{(x_{k}-y)^{2}}\right)}_{+}$
for all $k$.

Conversely, suppose $y\notin E_{j}$. If $y\in E_{k}$ for some other $k$, then
$C_{j}-{(x_{j}-y)^{2}}<f(y)=C_{k}-{(x_{k}-y)^{2}}$. If
$y\notin\bigcup_{1}^{N}E_{k}$, then $f(y)=0$, hence $C_{j}-{(x_{j}-y)^{2}}<0$.
It follows that
$C_{j}-{(x_{j}-y)^{2}}<\mathinner{\left({C_{k}}-(x_{k}-y)^{2}\right)}_{+}$ for
at least one $k$. ∎

######  Corollary 2.4.

Let $f_{j}(y)=C_{j}-{(x_{j}-y)^{2}}$. Then
$f(y)=\max\mathinner{\left\\{0,f_{1}(y),\ldots,f_{N}(y)\right\\}}$ and
$E_{j}=\mathinner{\left\\{y\mathrel{\mathop{\ordinarycolon}}f(y)=f_{j}(y)\right\\}}$.

Intuitively, all the players initially concentrated at $x_{j}$ should spread
out according to the density function $f_{j}$ over the set $E_{j}$. We now
make this intuition rigorous.

######  Proposition 2.5.

$E_{j}=\mathinner{\left\\{y\mathrel{\mathop{\ordinarycolon}}(x_{j},y)\in\operatorname{supp}\pi\right\\}}$,
$\int_{E_{j}\setminus\hat{E}_{j}}f=0$, and $\int_{E_{j}}f=a_{j}$.

###### Proof.

If $(x_{j},y)\in\operatorname{supp}\pi$, then by definition of equilibrium, it
follows that $y\in E_{j}$. Conversely, suppose $y\in E_{j}$. As explained
below in Lemma 3.1, $E_{j}$ is an interval, and if $y$ is in the interior of
that interval, then $y\notin\cup_{k\neq j}E_{k}$. By Corollary 2.4 and Lemma
2.2 we deduce that $(x_{j},y)\in\operatorname{supp}\pi$. If $y$ is on the
boundary of $E_{j}$, there is a sequence $y_{k}\to y$ with $y_{k}$ in the
interior of $E_{j}$, so $(x_{j},y_{k})\in\operatorname{supp}\pi$ and thus
$(x_{j},y)\in\operatorname{supp}\pi$.

To prove the remaining statements, start with the following inequality:

(2.1)
$a_{j}=m(\\{x_{j}\\})=\pi\mathinner{\left(\\{x_{j}\\}\times{\mathbb{R}}\right)}=\pi\mathinner{\left(\\{x_{j}\\}\times
E_{j}\right)}\leq\pi\mathinner{\left({\mathbb{R}}\times\hat{E}_{j}\right)}=\int_{E_{j}}f.$

Summing over $j$, we get

(2.2) $1=\sum_{j=1}^{n}a_{j}\leq\sum_{j=1}^{N}\int_{E_{j}}f=\int f=1.$

It follows that $\int_{E_{j}}f=a_{j}$, as desired. ∎

With these results, we can now reformulate the problem using only the
parameters $C_{1},\ldots,C_{n}$.

######  Definition 2.6.

For a given vector $C=(C_{1},\ldots,C_{n})\in{\mathbb{R}}^{n}$,

  * •

let $f_{j}^{(C)}(x)=C_{j}-(x-x_{j})^{2}$ for $j=1,\ldots,n$;

  * •

let $f^{(C)}(x)=\max\\{f_{1}^{(C)}(x),\ldots,f_{n}^{(C)}(x),0\\}$;

  * •

let
$E_{j}^{(C)}=\mathinner{\left\\{x\mathrel{\mathop{\ordinarycolon}}f_{j}^{(C)}(x)=f^{(C)}(x)\right\\}}$
for $j=1,\ldots,n$;

  * •

let $F_{j}(C)=\int_{E_{j}}f(x)\operatorname{d\\!}{x}$ for $j=1,\ldots,n$; and

  * •

let $F(C)=\mathinner{\left(F_{1}(C),\ldots,F_{n}(C)\right)}$.

To better understand these definitions, consider the following image:

Figure 2. Visualization of the density $f(x)$

The graph of $f(x)$ is the solid-color part (i.e. the maximum). Each
individual curve $f_{j}(x)$ can be thought of the density after the population
sitting at point $x_{j}$ “fans out” and forms a “bubble” surrounding that
point.

With these definitions, we now state the problem we wish to solve as follows:

######  Problem 2.7.

Given a vector $a=(a_{1},\ldots,a_{n})$, whose components are the weights in
the discrete measure $m=\sum_{j=1}a_{j}\delta_{x_{j}}$, we want to find a
vector $C$ such that $F(C)=a$.

By the arguments we have just given, if $f(x)$ is an equilibrium density, then
$f(x)=f^{(C)}(x)$ where $C$ is the solution to Problem 2.7. Conversely,
suppose $C$ is the solution to Problem 2.7 and let $f(x)=f^{(C)}(x)$. To see
that $f(x)$ is an equilibrium density, first define a function
$T\mathrel{\mathop{\ordinarycolon}}{\mathbb{R}}\to{\mathbb{R}}$ that sends
every point in the interior of $E_{j}=E_{j}^{(C)}$ to the point $x_{j}$.
Notice that $m$ is now the _push-forward_ of the density through $T$. (We say
a measure $\nu$ is the push-forward of a measure $\mu$ through a function $T$
if
$\nu(B)=T\sharp\mu(B)\mathrel{\mathop{\ordinarycolon}}=\mu\mathinner{\left(T^{-1}(B)\right)}$;
in this case $\nu=m$ and $\operatorname{d\\!}\mu=f(x)\operatorname{d\\!}x$.)
Set $\pi$ to be the push-forward of the density $f$ by the function
$y\mapsto(T(y),y)$. It follows that $\pi$ is an equilibrium.

We now state our main result:

######  Theorem 2.8 (Existence and uniqueness).

For each $a$ such that $a_{j}>0$ for $j=1,\dots,n$, there exists a unique $C$
such that $F(C)=a$.

Theorem 2.8 immediately implies the existence and uniqueness of a Nash
equilibrium for any game where the initial distribution is discrete. We will
prove Theorem 2.8 in two sections: Section 4 shows existence, and Section 5
proves uniqueness. These proofs rely on our analysis of $F(C)$ in Section 3,
in which we establish certain formulas for $F_{j}(C)$ and its derivatives,
prove continuity and coercivity, and the convexity of the effective domain of
$F$ (namely, the set of all $C$ for which $F_{j}(C)>0$ for all $j$).

##  3\. Structure of $F(C)$

The proof of Theorem 2.8 is based on various properties of $F(C)$, as given in
Definition 2.6. Whenever $C$ is fixed, we will usually suppress the argument
$C$ and simply write $f_{j}(x)$, $f(x)$, and $E_{j}$ instead of
$f_{j}^{(C)}(x)$, $f^{(C)}(x)$, and $E_{j}^{(C)}$. In Section 3.1, we study
the set $E_{j}$, namely the set over which $f(x)=f_{j}(x)$. We also study the
effective domain of $F(C)$, over which we look for solutions to Problem 2.7,
and we prove it is convex. In Section 3.2 we seek more or less explicit
formulas for $F_{j}(C)$ and its derivatives. Section 3.3 shows that $F$ is a
_coercive_ function on certain domains, while Section 3.4 shows that it is
continuous everywhere.

###  3.1. Properties of $E_{j}$

The most crucial points of interest will be the intersections between
parabolas, namely where $f_{i}(x)=f_{j}(x)$, or
$C_{i}-(x-x_{i})^{2}=C_{j}-(x-x_{j})^{2}$. We define the unique solution by

(3.1)
$\gamma_{ij}=\displaystyle\frac{C_{i}-C_{j}}{2(x_{j}-x_{i})}+\frac{x_{j}+x_{i}}{2},i\neq
j.$

Observe that if $f_{i}(x)>f_{j}(x)$ for all for $x<\gamma_{ij}$, then
$f_{i}(x)<f_{j}(x)$ for all $x>\gamma_{ij}$, and likewise if
$f_{i}(x)<f_{j}(x)$ for all for $x<\gamma_{ij}$, then $f_{i}(x)>f_{j}(x)$ for
all $x>\gamma_{ij}$. With this in mind, we proceed to the following claim:

######  Lemma 3.1.

$E_{j}$ is either empty, a singleton, or a closed, bounded interval.

###### Proof.

We first note that, since $E_{j}$ is by definition is a set, $E_{j}$ is empty,
a singleton, or contains at least two points. Hence, all we need to prove is
that if $E_{j}$ contains more than one point, then $E_{j}$ is a closed,
bounded interval. So, suppose that $E_{j}$ contains at least two points. Call
two of these points $\alpha,\beta\in E_{j}$. Let $x\in(\alpha,\beta)$. To show
that $E_{j}$ is an interval, we need to show that $x\in E_{j}$, i.e.
$f_{j}(x)\geq f_{k}(x)$ and 0 for all $k\neq j$. Since $f$ is a concave
function and both $f_{j}(\alpha)$ and $f_{j}(\beta)$ are greater than or equal
to 0, $f_{j}(x)\geq 0$. To prove $f_{j}(x)\geq f_{k}(x)$ for all $k\neq j$, it
suffices to show that there is no $\gamma_{jk}\in(\alpha,\beta)$. This can be
shown simply using proof by contradiction: Suppose that there exists a
$\gamma_{jk}\in(\alpha,\beta)$. Then, because of Remark 2.4, either
$f_{j}(\alpha)>f_{k}(\alpha)$ and $f_{j}(\beta)<f_{k}(\beta)$, or
$f_{j}(\alpha)<f_{k}(\alpha)$ and $f_{j}(\beta)>f_{k}(\beta)$. However, since
both $\alpha_{j},\beta_{j}\in E_{j}$, neither of these possibilities can hold.
Thus, $E_{j}$ must contain all the points between $\alpha$ and $\beta$. Next,
to show that $E_{j}$ is bounded, we simply note that
$E_{j}\subseteq{[\gamma_{ij},\gamma_{jk}]}$ for any $i<j<k$, and therefore is
bounded. Lastly, to show $E_{j}$ is closed, we define
$\mathinner{\left\\{x_{k}\right\\}}$ to be a sequence in $E_{j}$ such that
$x_{k}\rightarrow x$. Thus, be definition of $E_{j}$, $f_{j}(x_{k})=f(x_{k})$.
Since both $f_{j}(x)$ and $f(x)$ are continuous, letting
$k\rightarrow\infty,f_{j}(x)=f(x)$. Therefore, if $E_{j}$ contains at least
two points, then $E_{j}$ is a closed, bounded interval. ∎

In light of this, we introduce the following notation:

######  Definition 3.2.

In the case that $E_{j}$ is and interval, we define $\alpha_{j}$ and
$\beta_{j}$ as the left and right end points of $E_{j}$, i.e.
$[\alpha_{j},\beta_{j}]=E_{j}$. We denote the length of $E_{j}$ by
$\mathinner{\\!\left\lVert E_{j}\right\rVert}=\beta_{j}-\alpha_{j}$.

Before proceeding, we observe that Definition 2.6 leaves open the possibility
that $F_{j}(C)=0$. Since we are trying to solve $F(C)=a$ where all the
components of $a$ are positive, the true domain of interest is given as
follows:

######  Definition 3.3.

$\mathbb{G}=\mathinner{\left\\{C_{j}\mathrel{\mathop{\ordinarycolon}}F(C_{j})>0\text{
for }j=1,\dots,n\right\\}}$.

######  Corollary 3.4.

If $C\in\mathbb{G}$, then $E_{j}$ is an interval for all $j=1,\dots,n.$

###### Proof.

First, by definition of $\mathbb{G}$, if $C\in\mathbb{G}$, then $F_{j}(C)>0$
for all $j=1,\dots,n$. Further, by Definition 2.6, it is clear that, in order
for $F_{j}(C)$ to be greater than $0$, then $E_{j}$ cannot be empty or a
singleton. Thus, by 3.1, if $F(C)\in\mathbb{G}$, then $E_{j}$ must be an
interval for all $j=1,\dots,n$. ∎

######  Lemma 3.5.

For $C\in\mathbb{G}$, $\alpha_{j}$ and $\beta_{j}$ are continuous functions of
$C$ given by

(3.2)
$\alpha_{j}=\begin{cases}\displaystyle\frac{C_{j-1}-C_{j}}{2(x_{j}-x_{j-1})}+\frac{x_{j}+x_{j-1}}{2}\text{
if }f_{j}(\alpha_{j})>0,\\\ x_{j}-\sqrt{C_{j}},\text{ if
}f_{j}(\alpha_{j})=0,\end{cases}\beta_{j}=\begin{cases}\displaystyle\frac{C_{j}-C_{j+1}}{2(x_{j+1}-x_{j})}+\frac{x_{j+1}+x_{j}}{2},\text{
if }f_{j}(\beta_{j})>0,\\\ x_{j}+\sqrt{C_{j}},\text{ if
}f_{j}(\beta_{j})=0.\end{cases}$

###### Proof.

As $f=\max{(f_{1},f_{2},\dots,f_{n},0)}$, $\alpha_{j}$ is either
$\gamma_{(j-1)j}$ if $f_{j}(\alpha_{j})>0$ or the leftmost root of $f_{j}(x)$.
Likewise, $\beta_{j}$ is either $\gamma_{j(j+1)}$ or the rightmost root of
$f_{j}(x)$. We already have an equation for $\gamma_{ij}$ in (3.1), so all
that remains is to find when $f_{j}(x)=0$, which we do as follows:

$f_{j}(x)=0\Rightarrow C_{j}-(x-x_{j})^{2}=0\Rightarrow
x=x_{j}\pm\sqrt{C_{j}}.$

The claim follows. ∎

######  Corollary 3.6.

If $C\in{\mathbb{G}}$, then the length of $E_{j}(C)$ is given by

(3.3) $\mathinner{\\!\left\lVert
E_{j}\right\rVert}=k_{j}(C)=\min\\{\gamma_{j,j+1}(C),x_{j}+\sqrt{C_{j}}\\}-\max\\{\gamma_{j-1,j}(C),x_{j}-\sqrt{C_{j}}\\}.$

In fact, $C\in{\mathbb{G}}$ if and only if $k_{j}(C)>0$ for all $j$. In
particular, ${\mathbb{G}}$ is a convex set.

###### Proof.

If $C\in{\mathbb{G}}$, use Lemma 3.5 to see that $\mathinner{\\!\left\lVert
E_{j}\right\rVert}=\beta_{j}-\alpha_{j}=k_{j}(C)>0$. If $C\notin{\mathbb{G}}$,
then at least one $E_{j}$ must be a singleton or empty, i.e. $f_{j}(x)$ is
never the maximum in $f(x)=\max\\{f_{1}(x),\ldots,f_{n}(x),0\\}$. We deduce
that one of the following inequalities must be false:

  1. (1)

$\gamma_{j(j+1)}<\gamma_{(j-1)j}$

  2. (2)

$\gamma_{(j-1)j}>x_{j}+\sqrt{C_{j}}$

  3. (3)

$\gamma_{j(j+1)}<x_{j}-\sqrt{C_{j}}$

hence $k_{j}(C)\leq 0$.

Since $k_{j}(C)$ is a concave function and
${\mathbb{G}}=\\{C\mathrel{\mathop{\ordinarycolon}}k_{j}(C)>0\ \forall j\\}$,
it follows that ${\mathbb{G}}$ is convex. ∎

We make one more simple observation about ${\mathbb{G}}$.

######  Lemma 3.7.

$\mathbb{G}$ is non-empty.

###### Proof.

To show $\mathbb{G}$ is non-empty, by 3.3, we need a $C$ such that
$F_{j}(C)>0$ for all $j=1,\dots,n$. So, let $\epsilon>0$. Choose $\epsilon$
small enough so that, if $C_{j}=\epsilon$ for all $j=1,\dots,n$, then
$f_{j}(\alpha_{j})=0=f_{j}(\beta_{j})$ for all $j=1,\dots,n$, i.e. no $f_{j}$
intersects any $f_{i},i\neq j$ above the $x$-axis. Then
$E_{j}=[x_{j}-\sqrt{\epsilon},x_{j}+\sqrt{\epsilon}]$ for all $j=1,\dots,n$.
Thus, we have

(3.4) $F_{j}(C)=\frac{4}{3}\epsilon^{3/2},$

which is greater than $0$. Thus, $\mathbb{G}$ is non-empty. ∎

Thus, the set $\mathbb{G}$ is both non-empty and convex, which will be useful
in the proof of uniqueness of equilibria.

###  3.2. Computing $F_{j}(C)$ and its Derivatives

In this section, we want to calculate $F_{j}(C)$ and all of its derivatives
for $C\in\mathbb{G}$. So, we have the following:

######  Proposition 3.8.

$F_{j}(C)=C_{j}(\beta_{j}-\alpha_{j})-\frac{1}{3}(\beta_{j}-x_{j})^{3}+\frac{1}{3}(\alpha_{j}-x_{j})^{3}$.

###### Proof.

We can compute this explicitly:

(3.5)
$F_{j}(C)=\int_{\alpha_{j}}^{\beta_{j}}C_{j}-(x-x_{j})^{2}\operatorname{d\\!}x=C_{j}(\beta_{j}-\alpha_{j})-\frac{1}{3}(\beta_{j}-x_{j})^{3}+\frac{1}{3}(\alpha_{j}-x_{j})^{3}.$

In particular, taking into account the different possibilities for
$\alpha_{j}$ and $\beta_{j}$ (see 3.2):

$\begin{split}&F_{j}(C)=\int_{\alpha_{j}}^{\beta_{j}}C_{j}-(x-x_{j})^{2}\operatorname{d\\!}x=C_{j}(\beta_{j}-\alpha_{j})-\frac{1}{3}(\beta_{j}-x_{j})^{3}+\frac{1}{3}(\alpha_{j}-x_{j})^{3}\\\
&\text{If }f_{j}(\alpha_{j}),f_{j}(\beta_{j})>0,\\\
&F_{j}(C)=C_{j}\mathinner{\left(\frac{C_{j}-C_{j+1}}{2(x_{j+1}-x_{j})}+\frac{x_{j+1}+x_{j}}{2}\right)}-\frac{1}{3}\mathinner{\left(\frac{C_{j}-C_{j+1}}{2(x_{j+1}-x_{j})}+\frac{x_{j+1}+x_{j}}{2}-x_{j}\right)}^{3}\\\
&-C_{j}\mathinner{\left(\frac{C_{j-1}-C_{j}}{2(x_{j}-x_{j-1})}+\frac{x_{j}+x_{j-1}}{2}\right)}+\frac{1}{3}\mathinner{\left(\frac{C_{j-1}-C_{j}}{2(x_{j}-x_{j-1})}+\frac{x_{j}+x_{j-1}}{2}-x_{j}\right)}^{3},\\\
&\text{or if }f_{j}(\alpha_{j}),f_{j}(\beta_{j})=0,\\\
&F_{j}(C)=\mathinner{\left(C_{j}(x_{j}+\sqrt{C_{j}})-\frac{1}{3}(x_{j}+\sqrt{C_{j}}-x_{j})^{3}\right)}-\mathinner{\left(C_{j}(x_{j}-\sqrt{C_{j}})-\frac{1}{3}(x_{j}-\sqrt{C_{j}}-x_{j})^{3}\right)}=\frac{4}{3}C_{j}^{3/2}.\end{split}$

The other two cases, when $f_{j}(\alpha_{j})>0,f_{j}(\beta_{j})=0$ and when
$f_{j}(\alpha_{j})=0,f_{j}(\beta_{j})>0$, are straightforward variations of
these two cases. ∎

######  Remark 3.9.

These calculations show that, in particular, $0\leq
F_{j}(C)\leq\frac{4}{3}C_{j}^{3/2}$ for any $C$. This will be useful later.

######  Proposition 3.10.

$\displaystyle\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{j}}}}({C})=\mathinner{\\!\left\lVert
E_{j}\right\rVert}+\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}+\frac{f_{j}(\alpha_{j})}{2(x_{j}-x_{j-1})}$

###### Proof.

With $F(C)$ now calculated, we can easily calculate
$\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{j}}}}$, taking into account
the different possibilities for $\alpha_{j}$ and $\beta_{j}$ as above:

(3.6) $\begin{split}&\text{If }f_{j}(\alpha_{j}),f_{j}(\beta_{j})>0,\\\
\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{j}}}}=&\mathinner{\left(\frac{C_{j}-C_{j+1}}{2(x_{j+1}-x_{j})}+\frac{x_{{j+1}}+x_{j}}{2}\right)}+\mathinner{\left(\frac{C_{j}}{2(x_{j+1}-x_{j})}\right)}-{\frac{1}{2(x_{j+1}-x_{j})}}\mathinner{\left(\frac{C_{j}-C_{j+1}}{2(x_{j+1}-x_{j})}+\frac{x_{{j+1}}-x_{j}}{2}\right)}^{2}\\\
-&\mathinner{\left(\frac{C_{j-1}-C_{j}}{2(x_{j}-x_{j-1})}+\frac{x_{j}+x_{j-1}}{2}\right)}+\mathinner{\left(\frac{C_{i}}{2(x_{j}-x_{j-1})}\right)}-{\frac{1}{2(x_{j}-x_{j-1})}}\mathinner{\left(\frac{C_{j-1}-C_{j}}{2(x_{j}-x_{j-1})}-\frac{x_{j}-x_{j-1}}{2}\right)}^{2}\\\
=&\mathinner{\\!\left\lVert
E_{j}\right\rVert}+\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}+\frac{f_{j}(\alpha_{j})}{2(x_{j}-x_{j-1})},\\\
&\text{or if }f_{j}(\alpha_{j}),f_{j}(\beta_{j})=0,\\\
\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{j}}}}=&\ 2C_{j}.\end{split}$

Again, the other two cases are straightforward adaptations. In any case, we
have:

(3.7)
$\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{j}}}}=\mathinner{\\!\left\lVert
E_{j}\right\rVert}+\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}+\frac{f_{j}(\alpha_{j})}{2(x_{j}-x_{j-1})}.\\\
$

Notice that this holds for every case, as if
$f_{j}(\alpha_{j}),f_{j}(\beta_{j})=0$, then those components of the partial
derivative would be 0, leaving $\mathinner{\\!\left\lVert E_{j}\right\rVert}$,
which is in fact $2C_{j}$ in that case. ∎

######  Theorem 3.11.

$\mathinner{\dfrac{\partial{}F_{i}(C)}{\partial{C_{j}}}}=\displaystyle\frac{-f_{i}(\gamma_{ij})}{2(x_{j}-x_{i})}$
and
$\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{i}}}}=\displaystyle\frac{-f_{j}(\gamma_{ij})}{2(x_{j}-x_{i})}$

###### Proof.

We now seek to calculate all
$\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{i}}}}$ for $i\neq j$.
Without loss of generality, assume $i<j$. So, we only need to calculate
$\displaystyle\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{i}}}}$ and
$\displaystyle\mathinner{\dfrac{\partial{}F_{i}(C)}{\partial{C_{j}}}}$.
Further, if $f_{j}(x)$ does not touch $f_{i}(x)$ on or above the $x$-axis,
then $\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{i}}}}=0$, as slight
changes in $C_{i}$ won’t change $C_{j}$. So we only consider the case if
$f_{j}(\alpha_{j}),f_{j}(\beta_{j})>0$. We simply differentiate:

(3.8)
$\mathinner{\dfrac{\partial{}F_{i}(C)}{\partial{C_{j}}}}=\frac{1}{2(x_{j}-x_{i})}\mathinner{\left(\mathinner{\left(\frac{C_{i}-C_{j}}{2(x_{j}-x_{i})}+\frac{x_{j}-x_{i}}{2}\right)}^{2}-C_{i}\right)}=\frac{-f_{i}(\gamma_{ij})}{2(x_{j}-x_{i})}$
(3.9)
$\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{i}}}}={\frac{1}{2(x_{j}-x_{i})}\mathinner{\left(\mathinner{\left(\frac{C_{i}-C_{j}}{2(x_{j}-x_{i})}-\frac{x_{j}-x_{i}}{2}\right)}^{2}-C_{j}\right)}}=\frac{-f_{j}(\gamma_{ij})}{2(x_{j}-x_{i})}=\frac{-f_{i}(\gamma_{ij})}{2(x_{j}-x_{i})}$

∎

We notice the following:

######  Remark 3.12.

$\mathinner{\dfrac{\partial{}F_{i}(C)}{\partial{C_{j}}}}=\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{i}}}}$

Additionally, we observe that
$\mathinner{\dfrac{\partial{}F_{i}(C)}{\partial{C_{j}}}}$ for any $i$ and $j$
depends on at most three components of $C$, as $\alpha_{j}$ and $\beta_{j}$
each only depend on one $C_{i}$ other than $C_{j}$ itself. With this in mind,
we can say:

######  Lemma 3.13.

$\mathinner{\dfrac{\partial{}F_{i}(C)}{\partial{C_{j}}}}\neq 0$ for at most
three components of $C$: $C_{j-1},C_{j},$ and $C_{j+1}$

###  3.3. Coercivity

We say that the vector field $F(C)$ is coercive provided that
$\displaystyle\lim\limits_{\mathinner{\\!\left\lVert
C\right\rVert}\to\infty}\frac{F(C)\cdot C}{\mathinner{\\!\left\lVert
C\right\rVert}}=\infty$, where $\mathinner{\\!\left\lVert
C\right\rVert}=\mathinner{\left(C_{1}^{2}+\cdots+C_{n}^{2}\right)}^{1/2}$ is
the norm of the vector $C$. This is the usual condition for a vector field to
be surjective, as in the Minty-Browder Theorem [Bro67]. There is no chance for
$F(C)$ to be coercive on all of ${\mathbb{R}}^{n}$, since $F_{j}(C)=0$
whenever $C_{j}\leq 0$. However, this will not matter, is we are always aiming
at restricting to vectors $C$ with non-negative components. Instead, we will
use the following substitute for coercivity to prove existence of solutions.

######  Lemma 3.14.

There exist constants $\delta\in(0,1)$ and $M>0$ such that if
$C_{j}=\max\\{C_{1},...,C_{n}\\}\geq M$, then $F_{j}(C)\geq\delta C_{j}$.

###### Proof.

Choose a $C_{j}$ such that $C_{j}=\\{C_{1},...,C_{n}\\}$. We now seek to
estimate $F_{j}({C})$. Note that $E_{j}$ must be an interval. The smallest
$E_{j}$ could be is if $\beta_{j-1}=\alpha_{j}$ and $\beta_{j}=\alpha_{j+1}$.
In this case, we can compute $\alpha_{j}$ and $\beta_{j}$:

$\alpha_{j}=\frac{C_{j-1}-C_{j}}{2(x_{j}-x_{j-1})}+\frac{x_{j-1}+x_{j}}{2}\text{
, }\beta_{j}=\frac{C_{j}-C_{j+1}}{2(x_{j+1}-x_{j})}+\frac{x_{j}+x_{j+1}}{2},$

and by the maximality of $C_{j}$ we get

$\alpha_{j}\leq\frac{x_{j-1}+x_{j}}{2},\quad\beta_{j}\geq\frac{x_{j}+x_{j+1}}{2}.$

This gives a lower bound on how close $\alpha_{j}$ and $\beta_{j}$ can get to
$x_{j}$. Thus,

$\begin{split}F_{j}({C})\geq\int_{\frac{x_{j-1}+x_{j}}{2}}^{\frac{x_{j}+x_{j+1}}{2}}C_{j}-(x-x_{j})^{2}\operatorname{d\\!}x=C_{j}\frac{x_{j+1}-x_{j-1}}{2}-\frac{1}{3}\mathinner{\left(\mathinner{\left(\frac{x_{j+1}-x_{j}}{2}\right)}^{3}-\mathinner{\left(\frac{x_{j}-x_{j-1}}{2}\right)}^{3}\right)}.\end{split}$

Define

$\bar{\delta}\mathrel{\mathop{\ordinarycolon}}=\frac{\min_{j}\mathinner{\left\\{x_{j+1}-x_{j-1}\right\\}}}{2},\quad\bar{M}\mathrel{\mathop{\ordinarycolon}}=\frac{1}{3}\max_{j}\mathinner{\left(\mathinner{\left(\frac{x_{j+1}-x_{j}}{2}\right)}^{3}-\mathinner{\left(\frac{x_{j}-x_{j-1}}{2}\right)}^{3}\right)}_{+},$

where we recall $x_{+}\mathrel{\mathop{\ordinarycolon}}=\max\\{x,0\\}$. Then

$F_{j}(C)\geq\bar{\delta}C_{j}-\bar{M}.$

Pick any $0<\delta<\min\\{\bar{\delta},1\\}$ and set
$M=(\bar{\delta}-\delta)^{-1}\bar{M}$. If $C_{j}\geq M$, we get
$F_{j}(C)\geq\delta C_{j}$, as desired. ∎

We remark that Lemma 3.14 can be used to prove that $F$ is coercive on
${\mathbb{R}}^{n}_{\geq 0}$, which is defined to be the set of all
$(x_{1},\ldots,x_{n})$ such that $x_{j}\geq 0$ for every $j$. The reason is
that if $C\in{\mathbb{R}}^{n}_{\geq 0}$ and
$C_{j}=\max\\{C_{1},\ldots,C_{n}\\}$, then

$F(C)\cdot C\geq F_{j}(C)C_{j}.$

For the remaining details, let us first define the norm

$\mathinner{\\!\left\lVert
x\right\rVert}_{\infty}=\max\mathinner{\left\\{\mathinner{\\!\left\lvert
x_{1}\right\rvert},\ldots,\mathinner{\\!\left\lvert
x_{n}\right\rvert}\right\\}}$

and note that $\mathinner{\\!\left\lVert
x\right\rVert}_{\infty}=\max\\{x_{1},\ldots,x_{n}\\}$ for
$x\in{\mathbb{R}}^{n}_{\geq 0}$. Recall that

$\mathinner{\\!\left\lVert
x\right\rVert}_{\infty}\leq\mathinner{\\!\left\lVert
x\right\rVert}\leq\sqrt{n}\mathinner{\\!\left\lVert x\right\rVert}_{\infty}.$

Let $M,\delta$ be as in Lemma 3.14. If $C\in{\mathbb{R}}^{n}_{\geq 0}$ and
$\mathinner{\\!\left\lVert C\right\rVert}\geq M\sqrt{n}$, then
$\mathinner{\\!\left\lVert C\right\rVert}_{\infty}\geq M$ and so

$F(C)\cdot C\geq F_{j}(C)C_{j}\geq\delta
C_{j}^{2}=\delta\mathinner{\\!\left\lVert
C\right\rVert}_{\infty}^{2}\geq\frac{\delta}{\sqrt{n}}\mathinner{\\!\left\lVert
C\right\rVert}^{2}.$

Thus $F$ is coercive on ${\mathbb{R}}^{n}_{\geq 0}$. However, we also note
that only Lemma 3.14, and not coercivity per se, will be used directly to
prove the existence of solutions.

###  3.4. Continuity

With coercivity proved, we move to prove that $F(C)$ is continuous. Before
doing so, we need to explore what happens in $E_{j}=\varnothing$.

######  Lemma 3.15.

Let $I\subseteq{1,\dots,n}$ and $\mathcal{C}_{I}$ be the set of all $C$ such
that $E_{j}(C)=\varnothing$ for all $j\in I$. Then $\mathcal{C}_{I}$ is open,
i.e. for every $C\in\mathcal{C}_{I}$, there exists $\epsilon>0$ such that, if
$\mathinner{\\!\left\lVert\tilde{C}-C\right\rVert}<\epsilon$, then
$\tilde{C}\in\mathcal{C}_{I}$.

###### Proof.

Let $C\in{\mathcal{C}}_{I}$. Define
$\tilde{f}^{(C)}(x)=\max\\{f_{j}^{(C)}(x)\mathrel{\mathop{\ordinarycolon}}j\in
I\\}$. Then $f^{(C)}(x)>\tilde{f}^{(C)}(x)$ for all $x$, because all of the
sets $E_{j}(C)$, $j\in I$, are empty. Since all of the $f_{j}^{(C)}(x)$ are
parabolas that go to $-\infty$ as $\mathinner{\\!\left\lvert
x\right\rvert}\to\infty$, there exists $M>0$ such that if
$\mathinner{\\!\left\lvert x\right\rvert}>M$, then
$f^{(C)}(x)-\tilde{f}^{(C)}(x)\geq 1$. Now $f^{(C)}(x)-\tilde{f}^{(C)}(x)$ is
a continuous and positive function, hence it has a minimum $\delta>0$ on the
compact set $[-M,M]$. Without loss of generality, we can assume $\delta\leq
1$; in particular, we have $f^{(C)}(x)-\tilde{f}^{(C)}(x)\geq\delta$ for every
$x$. Now assume
$\mathinner{\\!\left\lvert\tilde{C}_{j}-C_{j}\right\rvert}<\delta/2$ for
$j=1,2,\ldots,n$. Then $\mathinner{\\!\left\lvert
f_{j}^{(\tilde{C})}(x)-f_{j}^{(C)}(x)\right\rvert}<\delta/2$ for every $j$ and
every $x$. It follows that $\mathinner{\\!\left\lvert
f^{(\tilde{C})}(x)-f^{(C)}(x)\right\rvert}<\delta/2$ and
$\mathinner{\\!\left\lvert\tilde{f}^{(\tilde{C})}(x)-\tilde{f}^{(C)}(x)\right\rvert}<\delta/2$,
so

(3.10)
$\begin{split}f^{(\tilde{C})}(x)-\tilde{f}^{(\tilde{C})}(x)&=f^{(\tilde{C})}(x)-f^{(C)}(x)+f^{(C)}(x)-\tilde{f}^{(C)}(x)+\tilde{f}^{(C)}(x)-\tilde{f}^{(\tilde{C})}(x)\\\
&>-\delta/2+\delta-\delta/2=0.\end{split}$

It follows that all the sets $E_{j}(\tilde{C})$, $j\in I$, are empty, and so
$\tilde{C}$ must also be in ${\mathcal{C}}_{I}$. Thus, $\mathcal{C}_{I}$ is
open. ∎

With this in mind, we can now prove continuity:

######  Theorem 3.16.

$F(C)$ is continuous.

###### Proof.

By Lemma 3.1, we know that $E_{j}$ is empty, a singleton, or is a closed,
bounded interval. First, consider the case if $E_{j}=\varnothing$. Then
$C\in\mathcal{C}_{I}$, as defined in Lemma 3.15. Take a sequence of vectors
$C^{(k)}\rightarrow C$. Then, for $k$ sufficiently large,
$C^{(k)}\in\mathcal{C}_{I}$ also. So $F_{j}(C^{(k)})=F_{j}(C)=0$ for every
$j\in I$. Therefore, all the $j\in I$ will not effect the continuity of
$F(C)$. In light of this, we can relabel the sequence
$\mathinner{\left\\{1,\dots,n\right\\}}\setminus I$, calling it
$\mathinner{\left\\{1,\dots,m\right\\}}$. Then $f(x)$ will still be the
$\max\\{f_{1}(x),f_{2}(x),\dots,f_{m}(x),0\\}$, as those indices which may
have been removed will not affect the max. Thus, we can assume for the
remainder of this proof that $E_{j}(C)\neq\varnothing$, i.e. $j=1,\dots,m$.

Next, suppose that $E_{j}(C)$ is an interval. Choose a sequence of vectors
$C^{(k)}\rightarrow C$. By equation (3.3), it is clear that, as $F_{j}(C)$ is
simply a polynomial in terms of $C$, that $F_{j}(C^{(k)})\rightarrow
F_{j}(C)$.

Lastly, suppose $E_{j}$ is a single point. Then, by the contra-positive of
Corollary 3.4, $F_{j}(C)=0$. Again, we choose a sequence of vectors
$C^{(k)}\rightarrow C$. So, we need to show $F_{j}(C^{(k)})\rightarrow 0$. If
$F_{j}(C^{(k)})>0$ then $E_{j}(C^{(k)})$ is an interval whose length is given
by $k_{j}(C^{(k)})$ defined in Corollary 3.6. Further, since $k_{j}(C)$ is a
continuous function of $C$ and $E_{j}(C)$ is a single point, then
$k_{j}(C^{(k)})\rightarrow 0$ as $k\rightarrow\infty$. Therefore, whenever
$E_{j}(C^{(k)})$ is an interval, its length shrinks to zero as
$k\rightarrow\infty$, and since in any other case we have $F_{j}(C^{(k)})=0$,
we deduce that $F_{j}(C^{(k)})\rightarrow 0$ as $k\to\infty$. ∎

We have now sufficiently explored various properties of $f(x)$ and $F(C)$, so
we may now move onto our first goal:

##  4\. Proving Existence of a Solution

In this section, we seek to show that a solution $C$ to the equation $F(C)=a$
exists for all $a$. Before we begin the proof, we state Brouwer Fixed Point
Theorem:

######  Lemma 4.1.

Every continuous function from a nonempty convex compact subset $K$ of a
Euclidean space to $K$ itself has a fixed point.

See e.g. [Flo03]. We now begin the main result:

######  Theorem 4.2.

Let $a\in{\mathbb{R}}^{n}_{\geq 0}$, i.e. let $a=(a_{1},\ldots,a_{n})$ with
$a_{j}\geq 0$ for each $j$. Then there exists a solution $C$ to the equation
$F(C)=a$.

###### Proof.

Let $a$ be given. Define
$L\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{n}_{\geq
0}\rightarrow\mathbb{R}^{n}_{\geq 0}$ as
$L_{i}(C)\mathrel{\mathop{\ordinarycolon}}=\mathinner{\left(C_{i}+a_{i}-F_{i}(C)\right)}_{+}$
(recall $x_{+}\mathrel{\mathop{\ordinarycolon}}=\max\\{x,0\\}$). We will use
the norm

$\mathinner{\\!\left\lVert
C\right\rVert}_{1}\mathrel{\mathop{\ordinarycolon}}=\mathinner{\\!\left\lvert
C_{1}\right\rvert}+\cdots+\mathinner{\\!\left\lvert
C_{n}\right\rvert}=C_{1}+\cdots+C_{n}.$

It is useful to recall that

$\mathinner{\\!\left\lVert
C\right\rVert}_{\infty}\leq\mathinner{\\!\left\lVert C\right\rVert}_{1}\leq
n\mathinner{\\!\left\lVert C\right\rVert}_{\infty}.$

By Lemma 3.14, there exist constants $\delta\in(0,1)$ and $M>0$ such that if
$C_{j}=\mathinner{\\!\left\lVert C\right\rVert}_{\infty}\geq M$, then
$F_{j}(C)\geq\delta C_{j}$. In this case,

$L_{j}(C)\leq\mathinner{\left((1-\delta)C_{j}+a_{j}\right)}_{+}=(1-\delta)C_{j}+a_{j}.$

In general,

$L_{i}(C)\leq C_{i}+a_{i}.$

Summing over all $i$, we deduce that

$\mathinner{\\!\left\lVert L(C)\right\rVert}_{1}\leq\mathinner{\\!\left\lVert
C\right\rVert}_{1}+\mathinner{\\!\left\lVert a\right\rVert}_{1}-\delta C_{j}.$

Recall that $C_{j}=\mathinner{\\!\left\lVert C\right\rVert}_{\infty}\geq
n^{-1}\mathinner{\\!\left\lVert C\right\rVert}_{1}$. Let
$r>\mathinner{\\!\left\lVert a\right\rVert}_{1}n/\delta$. Then we have

(4.1) $\mathinner{\\!\left\lVert
L(C)\right\rVert}_{1}\leq\mathinner{\\!\left\lVert
C\right\rVert}_{1}+\mathinner{\\!\left\lVert a\right\rVert}_{1}-\delta
n^{-1}\mathinner{\\!\left\lVert C\right\rVert}_{1}<\mathinner{\\!\left\lVert
C\right\rVert}_{1}\quad\text{if}~{}\mathinner{\\!\left\lVert
C\right\rVert}_{1}\geq r.$

We now restrict to a compact, convex set. Let $M_{r}$ be the set of all
$x\in\mathbb{R}_{\geq 0}^{n}$ such that $\mathinner{\\!\left\lVert
x\right\rVert}_{1}\leq r$. Then $M_{r}$ is closed and bounded, so $M_{r}$ is
compact, and it is clearly nonempty. Lastly, $M_{r}$ is convex, as it is the
intersection of two convex sets, namely, ${\mathbb{R}}^{\geq 0}$ and the
closed unit ball for the 1-norm. Thus, by Lemma 4.1, any continuous function
$g\mathrel{\mathop{\ordinarycolon}}M_{r}\rightarrow M_{r}$ has a fixed point.

Define a function

(4.2) $P(D)\mathrel{\mathop{\ordinarycolon}}=\begin{cases}\displaystyle
D,\mathinner{\\!\left\lVert D\right\rVert}_{1}\leq r\\\
\frac{rD}{\mathinner{\\!\left\lVert
D\right\rVert}_{1}},\mathinner{\\!\left\lVert
D\right\rVert}_{1}>r.\end{cases}$

Note that $P$ maps ${\mathbb{R}}^{n}_{\geq 0}$ into $M_{r}$. Finally, let
$Q(C)\mathrel{\mathop{\ordinarycolon}}=P(L(C))$ for $C\in M_{r}$. Then $Q$ is
continuous, as it is the composition of continuous functions, and it maps
$M_{r}$ to itself. Thus, by 4.1, $Q$ has a fixed point $C\in M_{r}$. If
$\mathinner{\\!\left\lVert C\right\rVert}_{1}=r$, then

(4.3) $r=\mathinner{\\!\left\lVert
C\right\rVert}_{1}=\mathinner{\\!\left\lVert
Q(C)\right\rVert}_{1}=\mathinner{\\!\left\lVert P(L(C))\right\rVert}_{1}.$

By construction of $P$, this implies

(4.4) $\mathinner{\\!\left\lVert L(C)\right\rVert}_{1}\geq r.$

But this contradicts (4.1). It follows that $\mathinner{\\!\left\lVert
C\right\rVert}_{1}<r$. Thus,

(4.5) $r>\mathinner{\\!\left\lVert
C\right\rVert}_{1}=\mathinner{\\!\left\lVert
Q(C)\right\rVert}_{1}=\mathinner{\\!\left\lVert P(L(C))\right\rVert}_{1}.$

Hence, $P(L(C))=L(C)$, so finally we have

(4.6) $C=Q(C)=P(L(C))=L(C).$

This implies that for each $i$,

$C_{i}=\mathinner{\left(C_{i}+a_{i}-F_{i}(C)\right)}_{+}.$

If $C_{i}>0$, this implies

$C_{i}=C_{i}+a_{i}-F_{i}(C)\Rightarrow F_{i}(C)=a_{i}.$

If $C_{i}=0$, then $F_{i}(C)=0$ and we get

$0=\mathinner{\left(a_{i}\right)}_{+}=a_{i}$

so again $F_{i}(C)=a_{i}$. Hence $F(C)=a$, as desired. ∎

##  5\. Proving Uniqueness of a Solution

By Theorem 4.2, we know that a solution to $F(C)=a$ exists whenever
$a\in{\mathbb{R}}^{n}_{\geq 0}$. It can be observed that, if one of the
$a_{j}=0$, then the solution need not be unique; this occurs, for example, if
one of the parabolas $f_{j}(x)$ lies entirely underneath another parabola
$f_{k}(x)$, in which case adjusting $C_{j}$ slightly will not change the
equality $F_{j}(C)=0$. Notice, however, that in such a case, one can simply
throw out $a_{j}$ and reduce the dimension of the problem. Thus, without any
loss of generality, we can assume $a_{j}>0$ for all $j$. We therefore have a
solution $C\in\mathbb{G}$, where ${\mathbb{G}}$ is defined in Section 3.1.

In this section, we prove that the solution is unique. To do so, we show that
$F(C)$ is strictly monotone in a classical sense. It is interesting to note
that monotonicity is often used to prove uniqueness in mean field games, the
most commonly used notion being the “Lasry-Lions monotonicity condition” (cf.
[LL07]). (For a discussion on alternative monotonicity conditions, see
[GM23].) Here, the coupling is indeed precisely through the density of
players, and therefore the Lasry-Lions monotonicity condition holds. However,
it is interesting to note that the monotonicity of $F(C)$, a function of the
parameters $C_{1},\ldots,C_{n}$, is not a direct consequence the Lasry-Lions
monotonicity condition in any obvious way.

######  Lemma 5.1.

On $\mathbb{G}$, $F(C)$ is strictly monotone, i.e. for any
$C,\tilde{C}\in\mathbb{G},(F(C)-F(\tilde{C})(C-\tilde{C})>0$.

###### Proof.

Recall that $\mathbb{G}$ is convex by Corollary 3.6. To show F is strictly
monotone, then, it suffices to show that

$\sum_{i=1}^{n}\sum_{j=1}^{n}\mathinner{\dfrac{\partial{}F_{i}(C)}{\partial{C_{j}}}}\cdot(v_{i}v_{j})>0\
\forall\ {C},{v}\in\mathbb{R}^{n}\text{ with }v\neq 0.$

Assume for the sake of notation that $v_{0}=0=v_{n+1}$. By Lemma 2.11,

$S\mathrel{\mathop{\ordinarycolon}}=\sum_{i=1}^{n}\sum_{j=1}^{n}\mathinner{\dfrac{\partial{}F_{i}(C)}{\partial{C_{j}}}}v_{i}v_{j}=\displaystyle\sum_{j=1}^{n}\mathinner{\left(\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{j}}}}(v_{j})^{2}+\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{j-1}}}}(v_{{j-1}}v_{j})+\mathinner{\dfrac{\partial{}F_{j}(C)}{\partial{C_{j+1}}}}(v_{j}v_{{j+1}})\right)}.$

Using the derivative formulas from (3.8),

$S=\sum_{j=1}^{n}\mathinner{\left(\mathinner{\\!\left\lVert
E_{j}\right\rVert}+\frac{f_{j}(\alpha_{j})}{2(x_{j}-x_{j-1})}+\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}\right)}\mathinner{\left(v_{j}\right)}^{2}\\\
-\mathinner{\left(\frac{f_{j}(\alpha_{j})}{2(x_{j}-x_{j-1})}\right)}\mathinner{\left(v_{j-1}v_{j}\right)}-\mathinner{\left(\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}\right)}\mathinner{\left(v_{j}v_{j+1}\right)}.$

Since $f_{1}(\alpha_{1})=0=f_{n}(\beta_{n})$, we can rewrite this as

$S=\sum_{j=1}^{n}{\mathinner{\\!\left\lVert
E_{j}\right\rVert}}{v_{j}}^{2}+\sum_{j=2}^{n}\mathinner{\left(\frac{f_{j}(\alpha_{j})}{2(x_{j}-x_{j-1})}\mathinner{\left(v_{j}\right)}^{2}-\mathinner{\left(\frac{f_{j}(\alpha_{j})}{2(x_{j}-x_{j-1})}\right)}\mathinner{\left(v_{j-1}v_{j}\right)}\right)}\\\
+\sum_{j=1}^{n-1}\mathinner{\left(\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}\mathinner{\left(v_{j}\right)}^{2}-\mathinner{\left(\frac{f_{j}(\beta_{j})}{{2(x_{(j+1)}-x_{j})}}\right)}\mathinner{\left(v_{j}v_{j+1}\right)}\right)}.$

Since $f_{j+1}(\alpha_{j+1})$ is either $f_{j}(\beta_{j})$ or $0$, we can
combine to get

$\begin{split}S=\displaystyle&\sum_{j=1}^{n}{\mathinner{\\!\left\lVert
E_{j}\right\rVert}}{v_{j}}^{2}+\sum_{j=1}^{n-1}\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}\mathinner{\left(v_{j+1}\right)}^{2}-\frac{2f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}\mathinner{\left(v_{j}v_{j+1}\right)}+\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}\mathinner{\left(v_{j}\right)}^{2}\\\
=\displaystyle&\sum_{j=1}^{n}{\mathinner{\\!\left\lVert
E_{j}\right\rVert}}{v_{j}}^{2}+\sum_{j=1}^{n-1}\mathinner{\left(\frac{f_{j}(\beta_{j})}{2(x_{j+1}-x_{j})}\right)}\mathinner{\left(v_{j+1}-v_{j}\right)}^{2}>0.\end{split}$

Thus, $F(C)$ is strictly monotone on $\mathbb{G}$. ∎

We now have everything we need to prove uniqueness.

######  Theorem 5.2.

The solution to the equation $F(C)=a$ has a unique solution.

###### Proof.

Suppose, for the sake of contradiction, that $C_{1},C_{2}\in\mathbb{G}$ with
$C_{1}\neq C_{2}$ are both solutions to $F(C)=a$. Then $F(C_{1})=F(C_{2})$, so
$(F(C_{1})-F(C_{2})(C_{1}-C_{2})=0$. But this contradicts 5.1. Thus, the
solution to $F(C)$ is unique. ∎

We have now shown that a unique solution to $F(C)=a$ exists and is unique. We
now move on to numerical approximations.

##  6\. Numerical Simulations

We used Newton’s Method to find approximate solutions to $F(C)=a$. By
algebraically solving $\frac{4}{3}C_{j}^{3/2}=a_{j}$, we obtain an initial
guess $C=C^{(0)}$ such that $F_{j}(C)\leq a_{j}$ for each $j$ (see Remark
3.9). Then we define an approximating sequence $C^{(n)}$ by

(6.1)
$C^{(n+1)}=C^{(n)}-DF\mathinner{\left(C^{(n)}\right)}^{-1}\mathinner{\left(F\mathinner{\left(C^{(n)}\right)}-a\right)},$

where $DF(C)$ is the Jacobian matrix computed in Section 3.2. The residual
$\mathinner{\\!\left\lVert F\mathinner{\left(C^{(n)}\right)}-a\right\rVert}$
can be computed easily; in practice it becomes vanishingly small after only a
few iterations. Below we report on three simulations. In each case, the
initial distribution $m=\sum_{j=1}^{n}a_{j}\delta_{x_{j}}$ is illustrated a
set of $n$ bars located at position $x_{j}$ with height $a_{j}$. Then the
final density $f(x)$ is then graphed over these bars. The function $f(x)$ is
piecewise quadratic with the $j$th parabola centered at $x_{j}$ and
representing the spreading of the population that starts at $x_{j}$; the area
under this parabola above the $x$-axis equals $a_{j}$. This accords with our
geometric intuition about the game: players spread out as much as possible so
as to avoid areas of high population density, and the shape of the final
density matches the quadratic cost $(x-y)^{2}$ of traveling from $x$ to $y$.

Figure 3. Example 1: $x=(1,2.25,3,3.75),a=(0.25,0.25,0.25,0.25)$

In the first simulation (Figure 3), we took an initial population uniformly
distributed over four points:
$m=\frac{1}{4}\mathinner{\left(\delta_{1}+\delta_{2.25}+\delta_{3}+\delta_{3.75}\right)}$.
In this example, the leftmost quarter of the population is far enough away
from rest that it can spread out without colliding with the others, whereas
the other three are close enough to each other that players have less room to
move, resulting in a narrower probability density.

Figure 4. Example 2: $x=(1,2,\dots,10),a=random$

Figure 5. Example 3: $x=random$, $a=random$

In the second simulation (Figure 4), we let $x_{j}-x_{j-1}$ be constant and
let $a$ be random. (Note that we did not normalize the components of $a$ to
sum to 1, but this does not change the shape of the solution.) Notice how in
this sample the probability density has a very interesting contour. For
example, the “bubble” (i.e. piece of parabola) corresponding to $x=9$ barely
appears because it is crowded by neighboring “bubbles.” This demonstrates how
the different weights locally change the shape of the density.

In the last simulation (Figure 5), the vectors $x$ and $a$ are random (and $a$
is normalized to be a probability vector). Although the initial measure is
highly irregular, the final density will still be smooth, though it appears
less so as the number of random points are chosen. Indeed, as the proof of
uniqueness in Section 5 suggests, the problem becomes more and more ill-
conditioned as the points $x_{j}$ get closer together, since the lengths of
the intervals $E_{j}$ necessarily get smaller.

##  7\. Conclusion

In this paper we have studied a mean field game with final cost equal to the
population density of players. We considered the case where the initial
measure is discrete, which seems to be new in the literature. We showed that
the problem reduces to a finite dimensional problem, which can prove is well-
posed by classical methods. Our numerical simulations illustrate the
“smoothing effect” that mean field games with density penalization are
expected to have. For future research, it would be interesting to see if this
approach could be extended to more general mean field games with density-
dependent costs and discrete initial measures.

## References

  * [ACD+20] Yves Achdou, Pierre Cardaliaguet, François Delarue, Alessio Porretta, Filippo Santambrogio, Pierre Cardaliaguet, and Alessio Porretta.  An introduction to mean field game theory.  Mean Field Games: Cetraro, Italy 2019, pages 1–158, 2020.
  * [BFY13] Alain Bensoussan, Jens Frehse, and Phillip Yam.  Mean field games and mean field type control theory.  Springer, 2013.
  * [Bro67] F.E. Browder.  Existence and perturbation theorems for nonlinear maximal monotone operators in Banach spaces.  Bull. Amer. Math. Soc., 73:322–327, 1967.
  * [Car15] Pierre Cardaliaguet.  Weak solutions for first order mean field games with local coupling.  In Analysis and Geometry in Control Theory and its Applications, volume 11 of Springer INdAM Series, pages 111–158. Springer, 2015.
  * [CD17] R Carmona and F Delarue.  Probabilistic Theory of Mean Field Games: vol. I, Mean Field FBSDEs, Control, and Games.  Stochastic Analysis and Applications, Springer Verlag, 2017.
  * [CG15] Pierre Cardaliaguet and P. Jameson Graber.  Mean field games systems of first order.  ESAIM: COCV, 21(3):690–722, 2015.
  * [CGPT15] Pierre Cardaliaguet, P. Jameson Graber, Alessio Porretta, and Daniela Tonon.  Second order mean field games with degenerate diffusion and local coupling.  Nonlinear Differential Equations and Applications NoDEA, 22(5):1287–1317, 2015.
  * [Flo03] Monique Florenzano.  General equilibrium analysis: existence and optimality properties of equilibria.  Springer Science & Business Media, 2003.
  * [G+14] Diogo A Gomes et al.  Mean field games models–a brief survey.  Dynamic Games and Applications, 4(2):110–154, 2014.
  * [GM18] P. Jameson Graber and Alpár R. Mészáros.  Sobolev regularity for first order mean field games.  Annales de l’Institut Henri Poincaré C, Analyse Non Linéaire, 35(6):1557–1576, September 2018.
  * [GM23] P. Jameson Graber and Alpár R. Mészáros.  On monotonicity conditions for mean field games.  Journal of Functional Analysis, 285(9):110095, 2023.
  * [Gra14] P Jameson Graber.  Optimal control of first-order Hamilton–Jacobi equations with linearly bounded Hamiltonian.  Applied Mathematics & Optimization, 70(2):185–224, 2014.
  * [HMC06] Minyi Huang, Roland P Malhamé, and Peter E Caines.  Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the nash certainty equivalence principle.  Communications in Information & Systems, 6(3):221–252, 2006.
  * [LL06a] Jean-Michel Lasry and Pierre-Louis Lions.  Jeux à champ moyen. I–Le cas stationnaire.  Comptes Rendus Mathématique, 343(9):619–625, 2006.
  * [LL06b] Jean-Michel Lasry and Pierre-Louis Lions.  Jeux à champ moyen. II–Horizon fini et contrôle optimal.  Comptes Rendus Mathématique, 343(10):679–684, 2006.
  * [LL07] Jean-Michel Lasry and Pierre-Louis Lions.  Mean field games.  Japanese Journal of Mathematics, 2(1):229–260, 2007.
  * [Mun22] Sebastian Munoz.  Classical and weak solutions to local first-order mean field games through elliptic regularity.  Annales de l’Institut Henri Poincaré C, 39(1):1–39, 2022.
  * [PS17] Adam Prosinski and Filippo Santambrogio.  Global-in-time regularity via duality for congestion-penalized mean field games.  Stochastics, 89(6-7):923–942, 2017.
  * [Zim24] Brady Zimmerman.  A finite dimensional approximation of a density dependent mean field game, 2024.  Bachelor’s Thesis.