<figure><img src="image_1.png" /></figure>

HAL
open science

# Angular probability density reconstruction by Maximum Entropy

Thi-Lich Nghiem, Viet-Duc Le, Thi-Lan Le, Daniel Delahaye, Pierre Maréchal

► To cite this version:

Thi-Lich Nghiem, Viet-Duc Le, Thi-Lan Le, Daniel Delahaye, Pierre Maréchal. Angular probability density reconstruction by Maximum Entropy. 2024. hal-04646918

HAL Id: hal-04646918

https://hal.science/hal-04646918v1

Preprint submitted on 12 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Angular probability density reconstruction by Maximum Entropy.
## Towards air traffic complexity estimation

Thi-Lich NGHIEM, Viet-Duc LE, Thi-Lan LE, Daniel DELAHAYE and Pierre MARÉCHAL

**Abstract**

This paper presents a method for evaluating local air traffic complexity based on the heading distribution of trajectories within a specific cell and time period. Utilizing Maximum Entropy estimation, we reconstruct the angular distribution by estimating its Fourier coefficients from statistical samples. Fenchel's duality underpins the analysis, leading to a simple estimation algorithm. The study lays the groundwork for creating complexity maps using differential geometry, detailed in a companion paper.

**keywords.** Angular probability density, Maximum Entropy principle, air traffic complexity.

## 1 Introduction

The complexity of air traffic is a critical area of study, particularly due to the increasing intensity of air traffic. This paper proposes a method for evaluating local complexity, partly based on the angular distribution, or heading, of trajectories passing through a given cell during a specified period. The data provides a statistical sample from which we estimate the Fourier coefficients of the desired distribution, ultimately reconstructing this distribution by applying the Maximum Entropy principle. The focus of this paper is on Maximum Entropy, striving to provide a rigorous construction without discretizing the problem, thereby avoiding the introduction of biases inherent in any discretization. Mathematically, Fenchel's duality is the cornerstone of our analysis and enables the development of a simple estimation algorithm. The construction of complexity maps will employ differential geometry tools, which are quite different from the techniques discussed in this paper. Therefore, this aspect will be covered in a companion paper. By leveraging Maximum Entropy estimation, our approach aims to offer a robust and unbiased framework for understanding and assessing the complexities of air traffic patterns. This method not only enhances the accuracy of local complexity evaluations but also provides a solid foundation for further studies in air traffic management.

The rest of our paper is structured as follows: In Section 2, we formulate our Maximum Entropy problem as an optimization problem and introduce our notation. Section 3 provides the necessary mathematical theory for a rigorous treatment of the Maximum Entropy principle in our context. Section 4 presents simulations to validate the method while Section 5 puts in perspective the application to air traffic complexity estimation.

## 2 Setting

For the sake of simplicity, we model airspace as a two-dimensional image, segmented into a Cartesian grid where each intersection point defines the center of a circular zone, which we name a *cell*. Assessing the local complexity is then made through the following steps:

(1) Select the trajectories that pass through the cell and determine, for each one, the *heading*, that is, the angle with respect to some fixed direction, at entry point (or at exit point, or at both entry and exit points). Then estimate, via empirical means, the Fourier coefficients of the underlying angular probability distribution.

(2) Solve a Maximum Entropy problem to obtain a probability distribution that is as *neutral* as possible. From the mathematical viewpoint, this step uses a *Fenchel's duality* approach, that enables to solve an infinite dimensional entropy problem via the resolution of a finite dimensional concave smooth maximization problem.

(3) Use the optimal probabilities to estimate the *local complexity*. This step used tools from differential geometry, and is differed to ulterior publication.

Figure 1 gives an example of sets of trajectories, while Figure 2 shows a cell extracted from the previous map, together with the corresponding angular histogram.
<figure><img src="image_2.png" /><figcaption>Figure 1: Example of trajectories. The red circle indicates the extracted cell displayed in Figure 2.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: A cell and the corresponding angular histogram.</figcaption></figure>
**Computing moments.** All cells are assumed to have the same radius $\delta$. On denoting by $\mathbf{u} = (u_1, u_2)$ the co-ordinate vector of a grid point where local complexity is to be evaluated, we consider all trajectories through the cell

$$S = \{(x, y) | (x - u_1)^2 + (y - u_2)^2 \le \delta^2\}.$$

Let $n$ be the number of trajectories $\mathcal{T}_j$ that intersect the cell, and the set of indices of these trajectories is $J$:

$$n = \text{card} J \quad \text{where} \quad J = \{j \in \mathbb{N} | \mathcal{T}_j \cap S \neq \emptyset\}.$$

Each trajectory that intersects the cell enters and with a heading denoted by $\theta_j$. We now regard the $\theta_j$'s as realizations of a random angle variable $\theta$. We are then interested in estimating the probability density $p(\theta)$ of such a variable. From the angular sampling $\theta_j$, we may build a set of empirical moments. The Fourier coefficients of $p$ are defined as

$$a_l = \frac{1}{\pi} \int_{0}^{2\pi} p(\theta) \cos(l\theta) \, d\theta \quad (l \in \mathbb{N}) \quad \text{and} \quad b_l = \frac{1}{\pi} \int_{0}^{2\pi} p(\theta) \sin(l\theta) \, d\theta \quad (l \in \mathbb{N}^*).$$

The corresponding empirical coefficients are obtained by replacing the probability measure $p(\theta) \, d\theta$ with the empirical measure $(1/n) \sum_{j \in J} \delta_{\theta_j}$. The so obtained coefficients

$$x_l = \frac{1}{\pi n} \sum_{j \in J} \cos(l\theta_j) \quad \text{and} \quad y_l = \frac{1}{\pi n} \sum_{j \in J} \sin(l\theta_j)$$

are regarded as statistical estimators of $a_l$ and $b_l$, respectively. Note in passing that the estimator $x_0$ gives the exact value $1/\pi$ of $a_0$. Clearly, in practice, these estimators will be meaningful for values of $l$ limited to some index $N$, depending on the number of trajectories intersecting the cell $S$. We are then led to consider the *data vector*

$$\mathbf{z} = (x_1, y_1, \dots, x_N, y_N) \in \mathbb{R}^{2N}.$$

**Formulating the Maximum Entropy problem.** If we could dispose of exact Fourier coefficients, we would then maximize the entropy of $p$ under the constraint that the moments of $p$ match these exact values. We would then aim at solving the following constrained optimization problem:

$$\begin{array}{ll} \underset{(\mathcal{P}_\circ)}{\text{Min}} & H(p) \\ \text{s.t.} & 1 = \int_0^{2\pi} p(\theta) \, d\theta, \quad \mathbf{z} = \int_0^{2\pi} p(\theta) \mathbf{T}(\theta) \, d\theta, \end{array}$$

in which $\mathbf{T}$ is the $\mathbb{R}^{2N}$-valued function given by $\mathbf{T}(\theta) := (\cos \theta, \sin \theta, \dots, \cos N\theta, \sin N\theta)$ and $H(p)$ denotes the Boltzmann-Shannon *neg-entropy* of $p$:

$$H(p) = \int_{0}^{2\pi} p(\theta) \ln p(\theta) \, d\theta.$$

However, we only dispose of approximate values of the Fourier coefficients (except for $l=0$, where the value is exact). There are many possibilities for the relaxation of the constraints. Assuming that some information on the covariance of the data vector $\mathbf{z}$ is available, a natural way to relax the second constraint in $(\mathcal{P}_\circ)$ is to penalize the square of the Mahalanobis distance between $\mathbf{z}$ and the expected value of $\mathbf{T}$ under the probability measure $p(\theta) \, d\theta$.

Notice that the exact covariance matrix $\Sigma$ is not known in general: it can only be estimated. Using the Mahalanobis distance requires, in principle, to dispose of a positive definite estimate of $\Sigma$. When the estimated covariance matrix is singular, the inverse does not exist, and then we have to resort to a degenerate version of the Mahalanobis distance, defined via the *degenerate norm*

$$\|\zeta\|_{\Sigma^\dagger}^2 = \begin{cases} \langle \zeta, \Sigma^\dagger \zeta \rangle & \text{if } \zeta \in \text{ran } \Sigma, \\ \infty & \text{otherwise,} \end{cases}$$

in which $\Sigma^\dagger$ denotes the pseudo-inverse of $\Sigma$, $\langle \cdot, \cdot \rangle$ is the standard Euclidean scalar product on $\mathbb{R}^{2N}$ and $\text{ran } \Sigma$ denotes the range of $\Sigma$. The relaxed version of Problem $(\mathcal{P}_\circ)$ reads

$$\begin{array}{ll} \underset{(\mathcal{P})}{\text{Minimize}} & H(p) + \frac{\alpha}{2} \|\mathbf{z} - \mathbf{A}p\|_{\Sigma^\dagger}^2 \\ \text{s.t.} & 1 = \mathbb{I}p, \end{array}$$

in which $\alpha > 0$ balances the intensity of relaxation and $\|\cdot\|_Q$ is the norm associated to $\langle \cdot, \cdot \rangle_Q$. For convenience, we have defined the linear mappings $\mathbb{I}: L^1([0, 2\pi]) \to \mathbb{R}$ and $\mathbf{A}: L^1([0, 2\pi]) \to \mathbb{R}^{2N}$ by

$$\mathbb{I}p = \int_{0}^{2\pi} p(\theta) \, d\theta \quad \text{and} \quad \mathbf{A}p = \int_{0}^{2\pi} p(\theta) \mathbf{T}(\theta) \, d\theta.$$
The infinite value entailed by $\|\cdot\|_{\Sigma^\dagger}$ in Problem (P) corresponds to a sharp constraint. Problem (P) falls under partially finite convex programming [1,2]. A remarkable fact is that, as will appear later on, solving (P) via Fenchel's duality will avoid the computation of $\Sigma^\dagger$. In Section 3 below, we give a detailed account of the dual approach to solving such problems. Fenchel's duality will enable us to solve the infinite dimensional constrained minimization problem (P) via the maximization in finite dimension of some dual function, which turns out to be everywhere finite, concave and smooth.

# 3 Duality in action

Our aim in this subsection is to show how Fenchel's duality can help us solve our optimization problems via the tractable dual problem. We first review a few basic definitions and results. The reader unfamiliar with convex analysis may refer to [3,4].

## 3.1 Review of useful convex analysis

Let $L$ be any real vector space. A function $f: L \to [-\infty, \infty]$ is said to be convex if its epigraph, the set

$$ \text{epi } f := \{ (x, \alpha) \in L \times \mathbb{R} \mid f(x) \le \alpha \}, $$

is a convex subset of $L \times \mathbb{R}$. It is said to be proper convex if it never takes the value $-\infty$ and it is not identically equal to $\infty$. A function $g: L \to [-\infty, \infty]$ is said to be concave if $-g$ is convex, and proper concave if $-g$ is proper convex. Notice that $g$ is concave if and only if its hypograph

$$ \text{hypo } g \{(\mathbf{x}, \alpha) | g(\mathbf{x}) \ge \alpha\} $$

is convex. The effective domain of a convex function $f$ is the set

$$ \text{dom } f = \{x \in L \mid f(x) < \infty\}. $$

The effective domain of a concave function $g$ is the set

$$ \text{dom } g = \{x \in L \mid g(x) > \infty\}. $$

The only functions that are both convex and concave are the affine functions and the functions identically equal to $\pm\infty$. Clearly, the domain of each affine function is equal to $L$, both as a convex and concave function. The other remaining functions would give rise to specific notation if they were not so rare.

It is customary to use indicator functions to encode constraints in optimization problems. Recall that the indicator function of a subset $C \subset L$ is the function

$$ \delta_C(x) := \begin{cases} 0 & \text{if } x \in C, \\ \infty & \text{otherwise.} \end{cases} $$

Now, let $L$ and $\Lambda$ be vector spaces paired (in the algebraic sense) by a bilinear mapping

$$ \langle \cdot, \cdot \rangle: L \times \Lambda \to \mathbb{R} \\ (x, \xi) \mapsto \langle x, \xi \rangle. $$

An standard example is $L = \mathbb{R}^d = \Lambda$ with the usual Euclidean scalar product. The convex conjugate of a function $f$ (that is convex or not) is the function

$$ f^*(\xi) = \sup \{ \langle x, \xi \rangle - f(x) \mid x \in X \}, \quad \xi \in \Lambda. $$

The concave conjugate of a function $f$ (that is concave or not) is the function

$$ f_*(\xi) = \inf \{ \langle x, \xi \rangle - f(x) \mid x \in X \}, \quad \xi \in \Lambda. $$

It is remarkable that convex conjugacy acts as an involution on certain classes of functions. For example, if $f: \mathbb{R}^d \to [-\infty, \infty]$ is a lower-semicontinuous proper convex function, then

$$ f^{**} := (f^*)^* = f. $$

Given a convex subset $C \in \mathbb{R}^d$, we call relative interior of $C$ the interior of $C$ with respect to its affine hull aff $C$. Recall that aff $C$ is the smallest affine subspace that contains $C$. The relative interior of $C$ is denoted by ri $C$. For example, if $C$ is a closed segment in $\mathbb{R}^2$, its interior is empty while its relative interior is the segment without its ends. It can be shown that the relative interior of a nonempty convex set is always nonempty.
**Theorem 1 (Fenchel)** Let $f$ and $g$ be functions on $\mathbb{R}^d$ respectively proper convex and proper concave such that

$$ \mathrm{ri} \, \mathrm{dom} f \cap \mathrm{ri} \, \mathrm{dom} g \neq \emptyset $$

Then

$$ \eta := \inf_{\mathbf{x} \in \mathbb{R}^d} \{f(\mathbf{x}) - g(\mathbf{x})\} = \sup_{\xi \in \mathbb{R}^d} \{g_*(\xi) - f^*(\xi)\} $$

and the supremum is attained.

The above theorem asserts equality between the optimal values of two problems, together with attainment in the second one. Is customary to call these underlying optimization problems the *primal problem* and the *dual problem*, respectively. In the above theorem, both the primal and dual are finite dimensional. Moreover, in order to tackle problems such as ($\mathcal{P}$), we must put some linear mapping in the picture. This is the purpose of the next theorem.

**Theorem 2** Let be given:

1. $L$ and $\Lambda$, real vector spaces;

2. $\langle \cdot, \cdot \rangle$, a bilinear form on $L \times \Lambda$;

3. $A: L \to \mathbb{R}^d$, a linear mapping;

4. $F: L \to (-\infty, \infty]$, a proper convex function;

5. $g: \mathbb{R}^d \to [-\infty, \infty)$, a proper concave function.

Assume that $A$ admits a formal adjoint mapping $A^*$, that is, a linear mapping $A^*: \mathbb{R}^d \to \Lambda$ such that $\langle Ax, y \rangle = \langle x, A^*y \rangle$ for every $x \in L$ and every $y \in \mathbb{R}^d$. Then, under the qualification condition

$$ (QC) \quad \mathrm{ri}(\mathrm{dom} F) \cap \mathrm{ri}(\mathrm{dom} g) \neq \emptyset, $$

one has

$$ \eta := \inf_{x \in X} \{F(x) - g(Ax)\} = \max_{\lambda \in \mathbb{R}^d} \{g_*(\lambda) - F^*(A^*\lambda)\}. $$

This theorem is the corner stone of what is referred to as *partially finite convex programming*. Various forms appeared in the literature (see in particular [1, 2]). The selected form is as in [5], where no topological structure on the infinite dimensional side is requested. The optimization problems

$$ \text{Minimize} \quad (F - g \circ A) \quad \text{and} \quad \text{Maximize} \quad (g_* - F^* \circ A^*) $$

are respectively referred to as the *primal* and *dual* problems. The function $D := g_* - F^* \circ A^*$ appearing in the dual problem is called the *dual function*. The theorem asserts the equality between the optimal values of the primal and dual problems, together with *dual attainment*. The next result will provide conditions that will guarantee *primal attainment* as well.

**Theorem 3 (Primal attainment)** With the notation and assumptions of the previous theorem, assume in addition that

$$ (QC^*) \quad \mathrm{ri} \, \mathrm{dom} g_* \cap \mathrm{ri} \, \mathrm{dom}(F^* \circ A^*) \neq \emptyset. $$

Suppose further that

(a) $F^{**} = F$ and $g^{**} = g$;

(b) there exists $\bar{\lambda}$ dual optimal and $\bar{x} \in \partial F^*(A^*\bar{\lambda})$ such that $F^* \circ A^*$ has gradient $A^*\bar{x}$ at $\bar{\lambda}$.

Then $\bar{x}$ is *primal optimal*.

The latter result provides not only a condition for primal attainment, but it also makes appear as a watermark the possibility of a link between primal and dual solutions. The bi-conjugate relationships in Assumption (a) are central in the theorem, and the difficulty in our problem is to prove that the entropy satisfies this property. It turns out that, in our context, it is possible to compute the conjugate of the Shannon entropy by *conjugating through the integral sign*.

Generally speaking, an integral functional is a functional of the form

$$ H(u) = \int_{\Omega} h(u(\omega), \omega) d\mu(\omega), \quad u \in L. $$

Here, $\Omega$ is assumed to be endowed with a $\sigma$-algebra and with a measure denoted by $\mu$; the function $h$ is called the integrand, and the argument $u$ is assumed to pertain to some space of measurable functions $L$. In our context, it is enough to consider such functional on the familiar space $L = L^1([0, 2\pi])$, implicitly endowed with the Borel
σ-algebra and the Lebesgue measure. Moreover, dependence of h in its second argument is not vital here, and we are only interested, in this paper, in the case where $h(u(\omega), \omega) = h_\circ(u(\omega))$, with

$$h_{\circ}(t) = \begin{cases} t \ln t & \text{if } t > 0, \\ 0 & \text{if } t = 0, \\ \infty & \text{if } t < 0. \end{cases}$$

Clearly, $h_\circ$ is a lower semi-continuous convex proper function and, as such, it satisfies $h_\circ^{**} = h_\circ$. Conjugating $H$ is elegantly performed by conjugating the integrand, as we shall see now. Following Rockafellar, we say that 'a space $L$ of measurable functions is decomposable if it is stable under bounded alteration on sets of finite measure.' Otherwise expressed, $L$ is decomposable if and only if it contains all functions of the form

$$1_T u_\circ + 1_{T^c} u,$$

in which $T$ has finite measure, $u_\circ$ is a measurable function such that the set $u_\circ(T)$ is bounded, and $u$ is any member of $L$. Here, $1_T$ denotes the characteristic function of $T$: $1_T(x)$ equals 1 if $x \in T$ and $1_T(x)$ equals zero otherwise; and $T^c$ denotes the complement of $T$. It is easy to see that $L^p$-spaces are decomposable, which includes our workspace $L^1([0, 2\pi])$.

**Theorem 4 (Rockafellar)** Let $L$ and $\Lambda$ be spaces of measurable functions on $\Omega$ paired by means of the standard integral bilinear form

$$\langle f, \varphi \rangle = \int_{\Omega} f(\omega) \varphi(\omega) d\omega.$$

Let $H$ be the integral functional of integrand $h_\circ$, that is,

$$H(u) = \int_{\Omega} h_{\circ}(u(\omega)) d\omega,$$

with $h_\circ$ proper convex and lower semi-continuous. Assume that $L$ is decomposable and that $H$ has nonempty effective domain. Then

$$H^*(\varphi) = \int h^*(\varphi(\omega)) d\omega$$

for every $\varphi \in \Lambda$, and $H^*$ is convex on $\Lambda$.

Applying the latter theorem with $L = L^1([0, 2\pi])$, $\Lambda = L^\infty([0, 2\pi])$ and $h_\circ$ the integrand of the Shannon entropy defined above, we see that, in our case,

$$H^*(\varphi) = \int_0^{2\pi} h_\circ^*(\varphi(\omega)) d\omega.$$

As can be seen from the appendix, the function $h_\circ^*$ is given on $\mathbb{R}$ by $h_\circ^*(\tau) = \exp(\tau - 1)$. Finally, since $\Lambda = L^\infty([0, 2\pi])$ is also decomposable, we obtain that

$$H^{**}(u) = \int_{0}^{2\pi} h_{\circ}^{**}(u(\omega)) d\omega = \int_{0}^{2\pi} h_{\circ}(u(\omega)) d\omega = H(u).$$

We conclude that our entropy satisfies the bi-conjugacy relationship requested in Theorem 3.

Before returning to our specific problem, let us state one more theorem, in which an explicit relationship between primal and dual solutions is obtained.

**Theorem 5 (Primal-dual relationship)** With the notation and assumptions of Theorem 2 assume in addition that $dom D$ has nonempty interior, that $H$ is an integral functional of integrand $h$ such that conjugacy through the integral sign is permitted. Assume that, as in Theorem 3, $H^{**} = H$ and $g^{**} = g$. Assume finally that the conjugate integrand $h^*$ is differentiable throughout $\mathbb{R}$, and that there exists some dual-optimal vector $\tilde{\lambda}$ in $int dom D$. If $\tilde{u}(\omega) := h^{*'}([\Lambda^*\tilde{\lambda}](\omega), \omega)$ belongs to $L$, then $\tilde{u}$ is a primal solution.

We are now ready to get back to our specific entropy problem.
## 3.2 Computing maximum entropy densities

We first note that Problem (P) can be written as

Minimize $H(p) - g_\circ(A_\circ p)$

over the space $L^1([0, 2\pi])$, in which

$$A_\circ p = (\mathbb{I}p; A_\circ p) \in \mathbb{R} \times \mathbb{R}^{2N} = \mathbb{R}^{1+2N}$$

and $g_\circ$ is the function on $\mathbb{R}^{1+2N}$ given by

$$g_{\circ}(\eta_{\circ}, \eta) = -\frac{\alpha}{2} \|\mathbf{z} - \eta\|_{\Sigma^{\dagger}}^{2} - \delta_{\{1\}}(\eta_{\circ}).$$

Straightforward computations show that the adjoint mapping $A_\circ^*: \mathbb{R}^{1+2N} \to L^\infty([0, 2\pi))$ is given by

$$A_\circ^*(\lambda_\circ, \lambda)(\theta) = \lambda_\circ + \langle \lambda, \mathbf{T}(\theta) \rangle.$$

Moreover,

$$(g_\circ)_*(\lambda_\circ, \lambda) = \lambda_\circ + \langle \mathbf{z}, \lambda \rangle - \frac{1}{2\alpha} \|\lambda\|_{\Sigma}^2. \quad (1)$$

Accounting for the fact that, as we have seen above, the entropy can be conjugated through the integral sign, the dual problem reads:

$$\text{Maximize} \quad D(\lambda_\circ, \lambda) := \lambda_\circ + \langle \mathbf{z}, \lambda \rangle - \frac{1}{2\alpha} \|\lambda\|_{\Sigma}^2 - \exp(\lambda_\circ - 1) \int_0^{2\pi} \exp \langle \lambda, \mathbf{T}(\theta) \rangle d\theta.$$

The function $D$ is concave, finite and differentiable throughout $\mathbb{R}^{1+2N}$. Its stationary points must satisfy the system

$$\left\{ \begin{array}{l} 0 = 1 - \exp(\bar{\lambda}_\circ - 1) \int_0^{2\pi} \exp \langle \bar{\lambda}, \mathbf{T}(\theta) \rangle d\theta, \\ \mathbf{z} = \mathbf{z} - \frac{1}{\alpha} \Sigma \bar{\lambda} - \exp(\bar{\lambda}_\circ - 1) \int_0^{2\pi} \mathbf{T}(\theta) \exp \langle \bar{\lambda}, \mathbf{T}(\theta) \rangle d\theta, \end{array} \right.$$

which reduces to

$$\mathbf{z} = \mathbf{z} - \frac{1}{\alpha} \Sigma \bar{\lambda} - \frac{\int_{0}^{2\pi} \mathbf{T}(\theta) \exp \langle \bar{\lambda}, \mathbf{T}(\theta) \rangle d\theta}{\int_{0}^{2\pi} \exp \langle \bar{\lambda}, \mathbf{T}(\theta) \rangle d\theta}.$$

It should be noticed that the above system is also the optimality system of the problem

$$(\tilde{D}) \quad \begin{array}{l} \text{Maximize} \quad \langle \mathbf{z}, \mathbf{z} \rangle - \frac{1}{2\alpha} \|\lambda\|_{\Sigma}^2 - \ln \int_0^{2\pi} \exp \langle \lambda, \mathbf{T}(\theta) \rangle d\theta \\ \text{s.t.} \quad \lambda \in \mathbb{R}^{2N}. \end{array}$$

**Proposition 6** The function

$$\tilde{D}(\lambda) := \langle \mathbf{z}, \mathbf{z} \rangle - \frac{1}{2\alpha} \|\lambda\|_{\Sigma}^2 - \ln \int_{0}^{2\pi} \exp \langle \lambda, \mathbf{T}(\theta) \rangle d\theta$$

to be maximized in Problem $(\tilde{D})$ is concave, smooth and everywhere finite. Its gradient is given by

$$\nabla \tilde{D}(\lambda) = \mathbf{z} - \frac{1}{\alpha} \Sigma \lambda - \frac{\int_{0}^{2\pi} \mathbf{T}(\theta) \exp \langle \lambda, \mathbf{T}(\theta) \rangle d\theta}{\int_{0}^{2\pi} \exp \langle \lambda, \mathbf{T}(\theta) \rangle d\theta}.$$

The function $h_\circ^*(\tau) = \exp(\tau - 1)$ obviously meets the requirements of Theorem 5. Provided we can obtain a dual solution $(\tilde{\lambda}_\circ, \tilde{\lambda})$, the optimal density is then given by

$$\bar{p}(\theta) = \exp[\tilde{\lambda}_\circ - 1 + \langle \tilde{\lambda}, \mathbf{T}(\theta) \rangle] = \frac{\exp\langle \tilde{\lambda}, \mathbf{T}(\theta) \rangle}{\int_{0}^{2\pi} \exp\langle \tilde{\lambda}, \mathbf{T}(\theta) \rangle d\theta}, \quad (2)$$

where $\tilde{\lambda}$ maximizes the function $\tilde{D}$.
**Remark 7** The primal dual relationship obtained above shows that the solution to our maximum entropy problem exists and belongs to the exponential family. The particular form in Equation (2) is well known. In the unrelaxed form of Problem (P), it can be obtained using classical inequalities (see e.g. [6], cited in [7]). However, we wish to stress some advantages one may have to resort to Partially Finite Convex Programming. First, we notice that the latter formalism can deal with the problems with relaxed constraints in a very elegant way. More importantly, by using our approach, we learn that the appropriate value of the Lagrange multiplier maximizes the well behaved dual function, which offers a powerful way to derive a computationally efficient algorithm: solving a nonlinear system is difficult in general, unless this nonlinear system is known to derive from Fermat's principle applied to the maximization of an everywhere finite smooth concave function. Finally, we see that the proposed framework may be used with a variety of entropy functional (such as Burg [8], Renyi [9], Tsallis [10], etc.). Although this goes beyond the scope of this paper, it certainly deserves to be mentioned.

## 3.3 Parameter selection

In the maximum entropy problem (P), the attach term

$$ \frac{\alpha}{2} \|z - A_p p\|_{\Sigma}^2 $$

depends of the choice of the highest frequency N to be accounted for, on the estimation of the covariance matrix Σ and on the choice of the regularization parameter α. We now describe our strategy for their empirical determination.

For a given value of N, the estimation of Σ can be made by means of the standard formula:

$$ \hat{\Sigma} = \frac{1}{n} \sum_j (Z_j - \bar{Z}) (Z_j - \bar{Z})^\top, $$

in which

$$ Z_j = \left( \frac{1}{\pi} \cos(\theta_j), \frac{1}{\pi} \sin(\theta_j), \dots, \frac{1}{\pi} \cos(N\theta_j), \frac{1}{\pi} \sin(N\theta_j) \right) $$

and

$$ \bar{Z} = \left( \frac{1}{n\pi} \sum_{j \in J} \cos(\theta_j), \frac{1}{n\pi} \sum_{j \in J} \sin(\theta_j), \dots, \frac{1}{n\pi} \sum_{j \in J} \cos(N\theta_j), \frac{1}{n\pi} \sum_{j \in J} \sin(N\theta_j) \right). $$

The choice of α and N is a bit more intricate due to their mutual dependency.

On the one hand, in the literature on inverse problems, one can find a variety of parameter selection rules. Among these rules, the Morozov discrepancy principle states that the regularization parameter α should be such that the corresponding solution $\bar{p}$ would give a residual $\|z - A\bar{p}\|$ equal to a number strictly greater than, but close to, the estimated size ρ of the error on the data. Therefore, in order to implement the Morozov principle, an estimate of the error size should be known, a reasonable assumption that we make in the context of this paper.

On the other hand, it is intuitively clear that the information provided by higher frequencies tend to be irrelevant: in the spirit of the maximum entropy strategy, we are looking for the most unbiased information on the underlying random variable (the heading θ) so that fine (high resolution) details are unwanted; moreover, the empirical Fourier coefficients are expected to become poorer approximations of the true ones as the frequency l increases, in particular when the number of angular samples is low. It is therefore natural to put a limit N on the highest frequency to be accounted for. The dimension of the dual variable will then be 2N. In the context of simulated random angles, we observe that the Kullback-Leibler entropy of the true density p₀ relative to the reconstructed density p decreases as N increases, and stabilizes beyond some value of N. Otherwise expressed, beyond a certain value of N, the gain in information provided by additional Fourier coefficients is negligible or non-existent. Figure 3 illustrates that there is no gain beyond N = 50 in our simulation, indicating that higher frequencies do not substantially improve the approximation of the true distribution. This finding supports the principle of using a limited number of Fourier coefficients to avoid overfitting and unnecessary computational complexity in the estimation process.

We may then alternatively estimate α and, using the Morozov discrepancy principle for choosing α and our information theoretic criterion for choosing N. The computation of the Morozov parameter α is performed using a dichotomic approach.

In practice, our alternating strategy proves to be easy to implement, and it provides reasonable estimates after very few iterations. The simulations presented in the following section demonstrate the ease of implementation and validate our approach.
<figure><img src="image_4.png" /><figcaption>Figure 3: Kullback–Leibler relative entropy versus highest frequency N.</figcaption></figure>

## 4 Simulations

In order to validate our estimation method, we now propose a few simulations. We first build a prior distribution, from which a sample {$\theta_1, \dots, \theta_n$} is drawn. We then compute the corresponding empirical Fourier coefficients. We finally apply the strategy described above for estimating the searched angular probability density.

**Building original distributions.** We choose prior densities as sums of quasi-Dirac distributions. To this end, we recall the definition of the bump function:

$$ \psi_{\circ}(\theta) = \begin{cases} \exp\left(-\frac{1}{1-\theta^2}\right) & \text{if } \theta \in (-1, 1), \\ 0 & \text{otherwise.} \end{cases} $$

It is well known that $\psi_{\circ}$ is a $C^{\infty}$-function, whose support is the interval $[-1, 1]$. We can emulate the behaviour of a Dirac distribution by letting

$$ \psi_{\beta}(\theta) := \frac{1}{\beta c} \psi_{\circ}\left(\frac{\theta}{\beta}\right) \quad \text{in which} \quad c := \int \psi_{\circ}(\theta) d\theta. $$

Here, $\beta$ is a strictly positive parameter. The closer $\beta$ is to zero, the more $\psi_{\beta}$ concentrates about the origin. It is readily seen that the support of $\psi_{\beta}$ is the interval $[-\beta, \beta]$. We may now define, for the purpose of our simulations, densities of the form

$$ p_{\circ}(\theta) = \sum_{k=1}^{p} c_k \psi_{\beta_k}(\theta - \theta_k), \quad (3) $$

in which the $c_k$'s are positive numbers that sum up to 1. The translation parameters $\theta_k$ are chosen in such a way that the support of $p_{\circ}$ is contained in $[0, 2\pi]$. We then simulate angular samples that follow such a distribution.

Figure 4 shows two examples of densities in the family defined in Equation (3), with $p = 2$ and $p = 5$ with all width parameters $\beta_k$ equal to 0.1, together with angular samples drawn randomly according to the density and Maximum Entropy estimations.

In our reconstructions, we exploit the fact that the dual function $\tilde{D}$ to be maximized is everywhere finite, concave and smooth. The BFGS algorithm is, of course well adapted for solving Problem (D) in our context. Figure 4 shows the reconstruction of the above two-peak density, in which the parameters have been adjusted according to the empirical strategy described in the previous section. After evaluating the highest frequency as $N = 40$, the Morozov principle was applied via a simple dichotomy. Figure 4 shows the resulting reconstruction.

In all situations that were tested, the optimal density closely mirrors the shape of the original density, successfully capturing the positions and general form of the peaks. As expected, the maximum entropy method produces reconstructions that tend to smooth out the peaks of the original distribution. Our strategy thus avoids underestimating the complexity (or disorder) of local traffic.
<figure><img src="image_5.png" /><figcaption>Figure 4: Prior densities (in red) with drawn samples (in green) and Maximum Entropy reconstructions, with 2-peak (first line) or 5-peak (second line) scenarios.</figcaption></figure>

# 5 Towards Air Traffic Complexity Estimation

The application of the maximum entropy method offers a novel approach for modeling and visualizing the complexities inherent in air traffic management.

Various measures of complexity can be derived from Maximum Entropy probability densities, not only from entrance densities but also from exit densities and especially from the discrepancy between both in a given cell. In fact, a large discrepancy would indicate that aircraft tend to change heading in this cell, reflecting a high level of complexity due to the unpredictability of the trajectories. Another criterion would involve counting the number of peaks in the entrance density, indicating potential conflicting routes. Ideally, a single peak suggests that all aircraft tend to follow the same direction, resulting in a low level of conflict. Naturally, the complexity must also be weighted by the intensity of the traffic, measured here by the number n of trajectories in the cell under consideration.

In the context of aircraft trajectories, complexity measures are used to understand and quantify the level of predictability and potential conflicts. The goal is to use these measures to analyze the patterns and behaviors of aircraft trajectories. Higher complexity, indicated by convergence, multiple peaks, or intense traffic, could imply a need for more sophisticated air traffic management strategies to ensure safety and efficiency. These measures help in assessing and optimizing air traffic control systems.

The potential applications of our methodology in real air traffic management are promising. By integrating actual traffic data into this framework, future research could refine and validate the models to enhance their practical utility. Additionally, expanding this approach to incorporate other variables such as altitude and speed could provide a more comprehensive tool for complexity assessment, supporting more advanced air traffic control strategies.

# A Computing conjugates

In this appendix, we give computation details for the conjugate functions involved in our maximum entropy problem. We start with the conjugate of $h_○$. This is a classical computation that we include for completeness. By definition,

$$h_○^*(\tau) = \sup_{t \in \mathbb{R}} \{\tau t - h_○(t)\}.$$

Suppose that the supremum is attained at $\bar{t} \in \mathbb{R}_+^*$. Then

$$\tau - h'_○(\bar{t}) = \tau - \ln \bar{t} - 1 = 0,$$

so that $\bar{t} = e^{\tau^{-1}}$. We see a posteriori that $\bar{t}$ maximizes the function $x \mapsto \tau t - h_○(t)$, since the maximum can't be attained on $\mathbb{R}_-$. Therefore,

$$h_○^*(\tau) = \tau \bar{t} - h_○(\bar{t}) = \tau e^{\tau^{-1}} - e^{\tau^{-1}}(\tau - 1) = e^{\tau^{-1}}.$$
Next, we deal with the computation of the concave conjugate of the function $g_\circ$. Clearly,

$$g_{\circ}(\eta_{\circ}, \eta) = -\delta_{\{1\}}(\eta_{\circ}) + g(\eta) \quad \text{with} \quad g(\eta) := -\frac{\alpha}{2} \|\mathbf{z} - \eta\|_{\Sigma^{\dagger}}^{2},$$

so that $(g_\circ)_*(\lambda_\circ, \lambda) = (-\delta_{\{1\}})_*(\lambda_\circ) + g_*(\lambda)$. We readily see that $(-\delta_{\{1\}})_*(\lambda_\circ) = \lambda_\circ$. We will show that

$$g_*(\lambda) = \phi(\lambda) := \langle \lambda, \mathbf{z} \rangle - \frac{1}{2\alpha} \langle \lambda, \Sigma\lambda \rangle, \quad (4)$$

from which Equation (1) will follow. Clearly, $g(\eta)$ is a closed proper concave function. Therefore, in order to obtain Equation (4), it suffices to show that $\phi_* = g$. By definition,

$$\phi_*(\eta) = \inf_{\lambda} \left\{ \langle \lambda, \eta \rangle - \langle \lambda, \mathbf{z} \rangle + \frac{1}{2\alpha} \langle \lambda, \Sigma\lambda \rangle \right\}.$$

The function between braces is convex, everywhere finite and differentiable. Its infimum will be attained whenever its gradient can vanish, which corresponds to the condition

$$\eta - \mathbf{z} + \frac{1}{\alpha} \Sigma \lambda = 0,$$

that is, $\Sigma\lambda = \alpha(\mathbf{z} - \eta)$. If $\mathbf{z} - \eta \in \text{ran}\Sigma$ (which always happens when $\Sigma$ is non-singular), then the gradient vanishes at any $\lambda$ such that $\Sigma\lambda = \alpha(\mathbf{z} - \eta)$, e.g. at $\tilde{\lambda} = \alpha\Sigma^{\dagger}(\mathbf{z} - \eta)$. As a matter of fact,

$$\Sigma\tilde{\lambda} = \alpha\Sigma\Sigma^{\dagger}(\mathbf{z} - \eta) = \alpha\text{Proj}_{\text{ran}\Sigma}(\mathbf{z} - \eta) = \alpha(\mathbf{z} - \eta),$$

in which $\text{Proj}_{\text{ran}\Sigma}(\cdot)$ denotes the orthogonal projection onto $\text{ran}\Sigma$. Therefore,

$$\begin{align*}
\phi_*(\eta) &= \langle \alpha\Sigma^{\dagger}(\mathbf{z} - \eta), \mathbf{z} - \eta \rangle + \frac{1}{2\alpha} \langle \alpha\Sigma^{\dagger}(\mathbf{z} - \eta), \alpha\Sigma\Sigma^{\dagger}(\mathbf{z} - \eta) \rangle \\
&= -\alpha \|\eta - \mathbf{z}\|_{\Sigma^{\dagger}}^{2} + \frac{\alpha}{2} \|\eta - \mathbf{z}\|_{\Sigma^{\dagger}}^{2} \\
&= -\frac{\alpha}{2} \|\eta - \mathbf{z}\|_{\Sigma^{\dagger}}^{2}.
\end{align*}$$

If $\mathbf{z} - \eta \notin \text{ran}\Sigma$ (which can only happen if $\Sigma$ is singular), then pick $\lambda_\circ \in \text{ker}\Sigma \setminus \{\mathbf{0}\}$. For every $t \in \mathbb{R}$,

$$\langle t\lambda_\circ, \eta - \mathbf{z} \rangle + \frac{1}{2\alpha} \langle t\lambda_\circ, \Sigma(t\lambda_\circ) \rangle = t \langle \lambda_\circ, \eta - \mathbf{z} \rangle.$$

But $\langle \lambda_\circ, \eta - \mathbf{z} \rangle$ cannot vanish. Otherwise, we would have:

$$\eta - \mathbf{z} \in (\text{ker}\Sigma)^{\perp} = \text{ran}\Sigma^{\top} = \text{ran}\Sigma,$$

in contradiction with the working condition $\mathbf{z} - \eta \notin \text{ran}\Sigma$. Consequently, in this case,

$$\phi_*(\eta) = -\infty = -\frac{\alpha}{2} \|\eta - \mathbf{z}\|_{\Sigma^{\dagger}}^{2}.$$

Putting things together, we obtain the formula in Equation (1).

## References

[1] J. M. Borwein and A. S. Lewis, “Partially finite convex programming, part i: Quasi relative interiors and duality theory,” *Mathematical Programming*, vol. 57, no. 1, pp. 15–48, 1992.

[2] ——, “Partially finite convex programming, part ii: Explicit lattice models,” *Mathematical Programming*, vol. 57, no. 1, pp. 49–83, 1992.

[3] R. T. Rockafellar, *Convex analysis*. Princeton university press, 1970, vol. 18.

[4] C. Zalinescu, *Convex analysis in general vector spaces*. World scientific, 2002.

[5] P. Maréchal, “On the principle of maximum entropy on the mean as a methodology for the regularization of inverse problems,” *Probability Theory and Mathematical Statistics*, pp. 481–492, 1998.

[6] A. Kagan, Y. Linnik, and C. Radhakrishna Rao, *Characterization problems in mathematical statistics*. Wiley, 1973.
[7] R. Gatto and S. R. Jammalamadaka, "The generalized von Mises distribution," *Statistical Methodology*, vol. 4, no. 3, pp. 341-353, 2007.

[8] J. P. Burg, "The relationship between maximum entropy spectra and maximum likelihood spectra," *Geophysics*, vol. 37, no. 2, pp. 375-376, 1972.

[9] A. Rényi, *Probability theory*. Courier Corporation, 2007.

[10] C. Tsallis, R. Mendes, and A. R. Plastino, "The role of constraints within generalized nonextensive statistics," *Physica A: Statistical Mechanics and its Applications*, vol. 261, no. 3-4, pp. 534-554, 1998.