%~Mouliné par MaN_auto v.0.34.0 (b0f8651d) 2025-06-26 16:52:10
\documentclass[CRMATH,Unicode,biblatex]{cedram}

\TopicFR{Analyse numérique}
\TopicEN{Numerical analysis}

\addbibresource{CRMATH_Iollo_20241036.bib}

\usepackage{subcaption}
\usepackage{bbm}
\usepackage{mathtools}
\usepackage[np]{numprint}

\newcommand{\fe}{\mathrm{fe}}
\newcommand{\fp}{\mathrm{fp}}
\newcommand{\hf}{\mathrm{hf}}
\newcommand{\KL}{\mathrm{KL}}
\newcommand{\ot}{\mathrm{ot}}
\newcommand{\iid}{\mathrm{iid}}
\newcommand{\reg}{\mathrm{reg}}
\newcommand{\dist}{\mathrm{dist}}
\newcommand{\internal}{\mathrm{int}}

\newcommand*{\dd}{\mathrm{d}}

\DeclareMathOperator*{\Uniform}{Uniform}

%%%---------------------------------------------------------------------------

%%% ARROWS

%%% \to
\renewcommand*{\to}{\mathchoice{\longrightarrow}{\rightarrow}{\rightarrow}{\rightarrow}}

%%% \mapsto
\let\oldmapsto\mapsto
\renewcommand*{\mapsto}{\mathchoice{\longmapsto}{\oldmapsto}{\oldmapsto}{\oldmapsto}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

%%% VERTICAL BAR

%%% ''Evaluated at'' macro
\newcommand\restr[2]{{				
 \left.\kern-\nulldelimiterspace 
 #1 								
 \vphantom{\big\vert} 				
 \right\vert_{#2} 	
}}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}

\title{Point-set registration in bounded domains via the Fokker--Planck equation}
\alttitle{Alignement de nuages de points dans des domaines bornés via l'équation de Fokker--Planck}

\author{\firstname{Angelo} \lastname{Iollo}}
\address{IMB, UMR 5251, Univ. Bordeaux, 33400 Talence, France} 
\address{Centre Inria de l'Université de Bordeaux, Team MEMPHIS, 33400 Talence, France}
\email{angelo.iollo@inria.fr}

\author{\firstname{Tommaso} \lastname{Taddei}}
\address[1]{IMB, UMR 5251, Univ. Bordeaux, 33400 Talence, France}
\address[2]{Centre Inria de l'Université de Bordeaux, Team MEMPHIS, 33400 Talence, France}
\email{tommaso.taddei@inria.fr}

%\begin{abstract}
%We present a point set registration method in bounded domains based on the solution to the Fokker--Planck equation. Our approach leverages (i) density estimation based on Gaussian mixture models; (ii) a stabilized finite element discretization of the Fokker--Planck equation; (iii) a specialized method for the integration of the particles. We review relevant properties of the Fokker--Planck equation that provide the foundations for the numerical method. We discuss two strategies for the integration of the particles and we propose a regularization technique to control the distance of the particles from the boundary of the domain. We perform extensive numerical experiments for two two-dimensional model problems to illustrate the many features of the method.
%\end{abstract}
%
%\begin{altabstract}
%Nous présentons une méthode d’alignement de nuages de points dans des domaines bornés basée sur la solution de l’équation de Fokker--Planck. Notre approche repose sur\tralicstex{}{ }: (i) l’estimation de densité à l’aide de modèles de mélange gaussien\tralicstex{}{ }; (ii) une discrétisation par éléments finis stabilisée de l’équation de Fokker--Planck\tralicstex{}{ }; (iii) une méthode spécialisée pour l’intégration des particules. Nous examinons les propriétés pertinentes de l’équation de Fokker--Planck qui sous-tendent la méthode numérique. Nous discutons deux stratégies pour l’intégration des particules et proposons une technique de régularisation pour contrôler la distance des particules à la frontière du domaine. Nous réalisons de nombreuses expériences numériques pour deux problèmes modèles bidimensionnels afin d’illustrer l'efficacité de la méthode.
%\end{altabstract}

\keywords{\kwd{Point set registration} \kwd{Fokker--Planck} \kwd{model order reduction}}
\altkeywords{\kwd{Alignement de nuages de points} \kwd{équation de Fokker--Planck} \kwd{réduction de modèle}}

\begin{document}

\begin{abstract}
We present a point set registration method in bounded domains based on the solution to the Fokker--Planck equation. Our approach leverages (i) density estimation based on Gaussian mixture models; (ii) a stabilized finite element discretization of the Fokker--Planck equation; (iii) a specialized method for the integration of the particles. We review relevant properties of the Fokker--Planck equation that provide the foundations for the numerical method. We discuss two strategies for the integration of the particles and we propose a regularization technique to control the distance of the particles from the boundary of the domain. We perform extensive numerical experiments for two two-dimensional model problems to illustrate the many features of the method.
\end{abstract}

\begin{altabstract}
Nous présentons une méthode d’alignement de nuages de points dans des domaines bornés basée sur la solution de l’équation de Fokker--Planck. Notre approche repose sur: (i) l’estimation de densité à l’aide de modèles de mélange gaussien; (ii) une discrétisation par éléments finis stabilisée de l’équation de Fokker--Planck; (iii) une méthode spécialisée pour l’intégration des particules. Nous examinons les propriétés pertinentes de l’équation de Fokker--Planck qui sous-tendent la méthode numérique. Nous discutons deux stratégies pour l’intégration des particules et proposons une technique de régularisation pour contrôler la distance des particules à la frontière du domaine. Nous réalisons de nombreuses expériences numériques pour deux problèmes modèles bidimensionnels afin d’illustrer l'efficacité de la méthode.
\end{altabstract}

\maketitle

\section{Introduction}\label{sec:introduction}

Point set registration (PSR) is the process of finding a spatial transformation that aligns two point clouds. Given the domain $\Omega \subset \mathbb{R}^n$ and the reference and the target point clouds $\{X_i^0\}_{i=1}^{N}, \{ X_j^\infty\}_{j=1}^{M} \subset \Omega$, we seek a transformation $\Phi \colon \Omega \to \Omega$ such that the mapped point cloud \smash{$\braces[\big]{\Phi(X_i^0)}_{i=1}^N$} is close --- in the sense of Hausdorff --- to the target point cloud. For several applications, it is also important to identify a (possibly non-autonomous) flow of diffeomorphisms $\Phi \colon \Omega \times \mathbb{R}_+ \to \Omega$ to coherently deform the reference point cloud $\{X_i^0\}_{i=1}^{N}$ into $\{ X_j^\infty\}_{j=1}^{M}$. In this note, we focus on the problem of PSR in moderate-dimensional ($n=2$ or $n=3$) bounded domains; we target applications in model order reduction (MOR) where coordinate transformations are employed to align coherent features of the solution field~\cite{cucchiara2024model,iollo2014advection,iollo2022mapping,taddei2020registration}; we also envision applications to data assimilation tasks that involve the interpolation between point clouds in complex geometries.

We propose a PSR method based on the solution to the Fokker--Planck (FP) equation. First, we model the reference and the target point clouds as independent identically distributed samples from two continuous distributions with densities $\rho_0$ and $\rho_\infty$, respectively. Second, we estimate $\rho_0$ and $\rho_\infty$ using Gaussian mixture models. Third, we solve a suitable FP equation to determine a flow of probability densities $t\mapsto \rho(\,\cdot\,, t)$ that is equal to $\rho_0$ at time $t=0$, and converges to $\rho_\infty$ as $t\to \infty$. Fourth, we exploit the properties of the FP equation to obtain the desired morphing based on the knowledge of $\rho(\,\cdot\,, t)$ at all times.

In this work, we aim to investigate the effectiveness of the FP equation for the transport of point clouds in bounded domains. The FP equation describes the evolution of a probability density $\rho$ under diffusion and drift: it reads as a linear parabolic equation with Neumann boundary conditions on $\partial \Omega$, constant diffusion and time-independent advection. The FP equation has been extensively employed in Information Theory~\cite{gray2011entropy}, stochastic processes in Physics and Chemistry~\cite{van1992stochastic}, and more recently in machine learning for generative modeling~\cite{song2019generative}. The solution to the FP equation is tightly linked to the solution to optimal transportation (OT) problems~\cite{santambrogio2015optimal}: however, while the FP equation is a linear parabolic equation with time-independent coefficients, the evolution equation associated to OT problems is highly nonlinear and much more difficult to handle, despite the many recent contributions to the field~\cite{peyre2019computational}.

We remark that PSR is a well-established field in image processing~\cite{beg2005computing,ma2018nonrigid,myronenko2010point}. Similarly to our approach, several PSR methods (e.g., \cite{ma2018nonrigid,myronenko2010point}) interpret the point clouds as samples from unknown continuous probability distributions and then rely on (Gaussian) mixture models to estimate the densities. However, the vast majority of PSR techniques in computer vision including~\cite{ma2018nonrigid,myronenko2010point} are designed for registration in unbounded domains and are hence not directly applicable to the case of bounded domains. A notable exception is provided by the large deformation diffeomorphic metric mapping (LDDMM, \cite{beg2005computing}) method: similarly to our approach, LDDMM relies on a flow of diffeomorphisms to accomplish the registration task; however, while our method models the point clouds as densities and then resorts to the FP equation to determine the transformation, LDDMM resorts to the solution to a computationally-expensive optimization problem.

The outline of the note is as follows. In Section~\ref{sec:mathematical_background}, we introduce the FP equation and we discuss well-known, elementary properties that are exploited in the remainder of the note. In Section~\ref{sec:methods}, we present the numerical method: first, we introduce the discretization of the FP equation, which is based on the finite element method for space discretization and on the Crank--Nicolson method for time integration; second, we discuss two distinct methods for the transport of the point clouds; third, we present a regularization approach to control the distance of the particles from the boundary; fourth, we briefly review the problem of density estimation through Gaussian mixtures. In Section~\ref{sec:numerics}, we present several numerical results in two dimensions that illustrate the performance of the method. Section~\ref{sec:conclusions} concludes the paper.

\section{Mathematical background}\label{sec:mathematical_background}

We review select properties of the Fokker--Planck (FP) equation that provide the mathematical foundation for the method. We emphasize that the FP equation has been the subject of extensive studies in Analysis and Calculus of Variations; in particular, the material of this section relies heavily on~\cite{gray2011entropy,jordan1998variational,ouhabaz2009analysis,santambrogio2015optimal}. Below, $\Omega$ is a Lipschitz domain in $\mathbb{R}^n$, while $\mathcal{P}_2(\Omega)$ refers to the space of probability density functions (pdfs) with bounded second-order moments; the acronym ``a.s.'' (resp., ``a.e.'') stands for ``almost surely'' (resp., ``almost every'') with respect to the Lebesgue measure in $\Omega$ and the acronym ``i.i.d.'' stands for ``independent identically distributed''.

\subsection{Kullback--Leibler divergence}

Given the densities $\rho,\rho_\infty\in \mathcal{P}_2(\Omega)$, we introduce the Kullback--Leibler (KL) divergence between $\rho$ and $\rho_\infty$ as~\cite{kullback1997information}
\begin{equation}\label{eq:KLdivergence}
D_{\KL}(\rho \parallel \rho_\infty) = \int_{\Omega} \rho(x) \log \parens[\Bigg]{\frac{\rho(x)}{\rho_\infty(x)}} \,\dd x.
\end{equation}
Note that $D_{\KL}(\rho \parallel \rho_\infty)$ is infinite if the support of $\rho_\infty$ does not contain the support of $\rho$. Theorem~\ref{th:KL} summarizes key properties of the KL divergence.

\begin{theo}\label{th:KL}
The following hold.
\begin{enumerate}
\item \label{th:KL_1} Given the densities $\rho,\rho_\infty\in \mathcal{P}_2(\Omega)$, $D_{\KL}(\rho \parallel \rho_\infty) \geq 0$ and $D_{\KL}(\rho \parallel \rho_\infty) =0$ if and only if $\rho=\rho_\infty$ a.s.
\item \label{th:KL_2} The KL divergence~\eqref{eq:KLdivergence} is asymmetric, that is in general $D_{\KL}(\rho \parallel \rho_\infty) \neq D_{\KL}(\rho_\infty \parallel \rho)$.
\item \label{th:KL_3} The KL divergence~\eqref{eq:KLdivergence} is strongly convex, that is, given $\rho_0,\rho_1,\rho_\infty\in \mathcal{P}_2(\Omega)$ and $\lambda\in (0,1)$, we have
\[
D_{\KL}\parens[\big]{\lambda \rho_0 + (1-\lambda) \rho_1 \parallel \rho_\infty}
\leq
\lambda D_{\KL}(\rho_0 \parallel \rho_\infty) + (1-\lambda) D_{\KL}(\rho_1 \parallel \rho_\infty),
\]
where the latter holds with equality if and only if $\rho_0=\rho_1$ a.s.
\item \label{th:KL_4} (Pinsker's inequality) The KL divergence~\eqref{eq:KLdivergence} satisfies the inequality:
\[
\norm{\rho - \rho_\infty}_{L^1(\Omega)}
\leq
\min \parens[\big]{
\sqrt{ 2 D_{\KL}(\rho \parallel \rho_\infty)},
\;
\sqrt{ 2 D_{\KL}(\rho_\infty \parallel \rho)}}.
\]
\end{enumerate}
\end{theo}

We provide the proof of Theorem~\ref{th:KL} in Appendix~\ref{sec:proofs}. We notice that the KL divergence is minimized when $\rho=\rho_\infty$ (cf.\ property~\eqref{th:KL_1}); as discussed below, the strong convexity property and the Pinsker's inequality are the key to prove the existence of paths that link the initial and the target distributions (cf.\ Theorem~\ref{th:FPequation}). We further remark that the KL divergence is not a distance in $\mathcal{P}_2(\Omega)$ since it is not symmetric (cf.\ property~\eqref{th:KL_2}). If we denote by $V=-\log(\rho_\infty)$ the potential associated with the density $\rho_\infty$, we can express~\eqref{eq:KLdivergence} as:
\begin{equation}\label{eq:KLdivergence_revisited}
D_{\KL}(\rho \parallel \rho_\infty) = \int_{\Omega} F[\rho] \,\dd x =
\underbrace{\int_\Omega \rho \log (\rho) \,\dd x}_{\eqqcolon \; \mathrm{(I)}}
\, + \,
\underbrace{\int_\Omega \rho V \,\dd x}_{\eqqcolon \; \mathrm{(II)}},
\quad \text{where} \ \ F[\rho] =
\rho \log (\rho) + \rho V.
\end{equation}
In Bayesian inference, the first term in~\eqref{eq:KLdivergence_revisited} is interpreted as an entropy function, while the second term is interpreted as a potential energy, and the sum of the two terms is interpreted as a free-energy function (see, e.g., \cite{fox2012tutorial}).

\subsection{Fokker--Planck equation}

Let $\rho, \rho_\infty \in \mathcal{P}_2(\Omega)$ and let $V=-\log(\rho_\infty)$; then, we define the Fokker--Planck (FP) equation:
\begin{equation}\label{eq:FPequation}
\begin{cases}
\partial_t \rho + \nabla \cdot \parens{- \rho \nabla V - \nabla \rho} = 0 	& \text{in} \ \Omega \times \mathbb{R}_+, \\
\rho(\,\cdot\,,0) = \rho_0 															& \text{in} \ \Omega, \\
\parens{\rho \nabla V + \nabla \rho} \cdot \mathbf{n} = 0 					& \text{on} \ \partial \Omega \times \mathbb{R}_+, \\
\end{cases}
\end{equation}
where $\mathbf{n}$ denotes the outward normal to $\partial \Omega$. We assume that $\rho_\infty$ is strictly positive in $\Omega$ and is sufficiently smooth to ensure that $\nabla V$ is Lipschitz continuous. We introduce the spaces $\mathcal{H}=L^2(\Omega)$ and $\mathcal{V}=H^1(\Omega)$ and the weak formulation of the FP equation~\eqref{eq:FPequation}: find $\rho \in C(\mathbb{R}_+; \mathcal{H})$ such that $\rho(0) = \rho_0$ and
\begin{equation}\label{eq:weak_form}
\int_\Omega \parens[\Big]{\partial_t \rho(t) \cdot v + \parens[\big]{\nabla \rho(t) + \rho(t) \nabla V} \cdot \nabla v} \,\dd x = 0
\quad
\forall v\in \mathcal{V}, \ \ \text{a.e.} \ t\in \mathbb{R}_+.
\end{equation}

Theorem~\ref{th:FPequation} reviews relevant properties of the FP equation; the proof is postponed to Appendix~\ref{sec:proofs}.

\begin{theo}\label{th:FPequation}
Let\/ $\Omega$ be a Lipschitz domain in $\mathbb{R}^n$; let $\rho_0, \rho_\infty \in \mathcal{P}_2(\Omega)$ satisfy:
\begin{enumerate}\romanenumi
\item $\inf_{x\in \Omega} \rho_\infty (x) > 0$;
\item $\nabla V $ is Lipschitz-continuous in\/ $\Omega$.
\end{enumerate}
The following hold.
\begin{enumerate}
\item \label{th:FPequation_1} There exists a unique weak solution $\rho \in C(\mathbb{R}_+; \mathcal{H})$ to~\eqref{eq:FPequation}.
\item \label{th:FPequation_2} The solution to~\eqref{eq:FPequation} satisfies $\rho(\,\cdot\,,t) \in \mathcal{P}_2(\Omega)$ for all $t>0$.
\item \label{th:FPequation_3} The solution to~\eqref{eq:FPequation} converges to $\rho_\infty$ for $t\to \infty$ in $L^1$, that is, $\lim_{t\to \infty} \norm[\big]{\rho(\,\cdot\,, t) - \rho_\infty}_{L^1(\Omega)} = 0$.
\item \label{th:FPequation_4} Let $u = - \nabla \parens[\big]{\log (\rho) + V}$ and let \smash{$X_1^0,\dotsc,X_N^0 \overset{\iid}{\sim} \rho_0$}. Consider the trajectories
\begin{equation}\label{eq:trajectories_FP}
\begin{cases}
\dot{X}_i(t) = u\parens[\big]{X_i(t), t} 	& t>0, \\
X_i(0) = X_i^0, \\
\end{cases}
\qquad i=1,\dotsc,N.
\end{equation}
Then, for any $t>0$, we have $X_1(t),\dotsc,X_N(t) \overset{\iid}{\sim} \rho(\,\cdot\,, t)$.
\end{enumerate}
\end{theo}

Theorem~\ref{th:FPequation} provides the foundations for the use of the FP equation for point-set registration. The first statement ensures the existence and the uniqueness of the solution to~\eqref{eq:FPequation}; the second statement ensures that $\rho(\,\cdot\,, t)$ is a pdf for all $t>0$; the third statement shows that $\rho(\,\cdot\,, t)$ converges to the target distribution for $t\to \infty$; finally, the fourth statement provides a constructive way to sample from the target distribution and also to generate trajectories that link initial and target pdfs. Notice that the prescribed Neumann boundary conditions in~\eqref{eq:FPequation} ensure that $u\cdot \restr{\mathbf{n}}{\partial \Omega} =0$; this condition ensures that the particles do not exit the domain.

\subsection{Fokker--Planck equation and gradient flows}\label{sec:gradient_flows}

The authors of~\cite{jordan1998variational} showed that the FP equation can be interpreted as the gradient flow of the KL divergence with respect to the Wasserstein metric. This perspective forms the foundation of the ``gradient flow'' method proposed in Section~\ref{sec:methods} to determine the trajectories of particles. We can hence interpret the solution to the time-discrete FP equation as the solution to a sequence of (linearized) optimal transport problems. Below, we provide an informal proof of this statement for $\Omega=\mathbb{R}^n$; we refer to~\cite{jordan1998variational} for the rigorous derivation and for further details.

First, we recall the 2-Wasserstein metric between the densities $\rho_0$ and $\rho_1$:
\begin{equation}\label{eq:wasserstein_metric}
W_2^2(\rho_0,\rho_1) \coloneqq \inf_{X \colon \mathbb{R}^n\to \mathbb{R}^n} \int_{\mathbb{R}^n} \abs[\big]{X(\xi) - \xi}^2 \rho_0(\xi) \,\dd \xi \quad \text{s.t.} \ \ \rho_1\parens[\big]{X(\xi)} \det \nabla X(\xi) = \rho_0(\xi).
\end{equation}
It is possible to show that the optimal transport map $X$ is of the form $X = \texttt{id} + \nabla \psi $ where $\texttt{id}(\xi) = \xi$ is the identity map and $ \psi \colon \mathbb{R}^n \to \mathbb{R}$ (cf.~\cite{brenier1991polar}, see also~\cite[Chapter~1.3.1]{santambrogio2015optimal}). Therefore, we can restate~\eqref{eq:wasserstein_metric} as
\begin{equation}\label{eq:wasserstein_metric_step2}
W_2^2(\rho_0,\rho_1) \coloneqq \inf_{\psi \colon \mathbb{R}^n\to \mathbb{R}} \int_{\mathbb{R}^n} \abs[\big]{\nabla \psi(\xi)}^2 \rho_0(\xi) \,\dd \xi\quad \text{s.t.} \ \
\rho_1\parens[\big]{\xi + \nabla \psi(\xi)} \det \parens[\big]{
\mathbbm{1}+ H[\psi](\xi)} = \rho_0(\xi),
\end{equation}
where $\mathbbm{1}\in \mathbb{R}^{n\times n}$ is the identity matrix and $H[\psi]$ is the Hessian of $\psi$.

We denote by $\braces[\big]{t^{(k)} = k \Delta t}_{k \in \mathbb{N}}$ the time grid, and we denote by $\rho^{(k)}$ the estimate of $\rho$ at time $t^{(k)}$; we assume that $\rho^{(k+1)}$ results from the evaluation of $\rho^{(k)}$ subject to a potential velocity field $\varepsilon \nabla \psi$ with $\varepsilon \ll 1$. We express $\rho^{(k+1)} = \rho^{(k)} + \varepsilon \eta$. Since $\rho^{(k+1)}$ satisfies $\rho^{(k+1)}\parens[\big]{X(\xi)} \det \nabla X(\xi) = \rho^{(k)}(\xi)$ with $X = \texttt{id} + \varepsilon \nabla \psi$, exploiting the Jacobi formula for the derivative of the determinant, we find
\begin{equation}\label{eq:expression_eta}
\eta = - \nabla \cdot (\rho^{(k)} \nabla \psi) + \mathcal{O}(\varepsilon).
\end{equation}

We seek $\psi$ to maximize the first variation of the KL divergence for a fixed Wasserstein distance:
\begin{equation}\label{eq:wasserstein_gradient}
\psi^{\star} \coloneqq \arg\inf_{\psi} \int_{\mathbb{R}^n} F'[\rho^{(k)}] \cdot \eta(\psi) \,\dd x
\quad \text{s.t.} \ \
\begin{dcases}
\int_{\mathbb{R}^n} \abs{\nabla \psi}^2 \rho^{(k)}\,\dd x =c, \\
\eta(\psi) = - \nabla \cdot (\rho^{(k)} \nabla \psi),\\
\end{dcases}
\end{equation}
where $c$ is a given constant and $F'[\rho] = \log(\rho) + 1 + V $. If we integrate by part the objective of~\eqref{eq:wasserstein_gradient}, we find
\[
\int_{\mathbb{R}^n} F'[\rho^{(k)}] \cdot \eta(\psi) \,\dd x = -
\int_{\mathbb{R}^n}
\nabla \cdot \parens[\big]{
\rho^{(k)} F'[\rho^{(k)}]
} \psi
\,\dd x
\]
In conclusion, we can restate~\eqref{eq:wasserstein_gradient} as
\begin{equation}\label{eq:wasserstein_metric_step3}
\psi^{\star} \coloneqq \arg \inf_{\psi} -
\int_{\mathbb{R}^n}
\nabla \cdot \parens[\big]{
\rho^{(k)}
\nabla F'[\rho^{(k)}]
} \psi
\,\dd \xi
\quad \text{s.t.} \ \ \int_{\mathbb{R}^n} \abs{\nabla \psi}^2 \rho^{(k)} \,\dd x =c.
\end{equation}

The optimality conditions for~\eqref{eq:wasserstein_metric_step3} are
\begin{equation}\label{eq:wasserstein_metric_step4}
\begin{dcases}
- \int_{\mathbb{R}^n} \nabla \cdot \parens[\big]{\rho^{(k)} \nabla F'[\rho^{(k)}]} \cdot \delta \psi \,\dd x + \lambda \int_{\mathbb{R}^n}
\rho^{(k)} \nabla \psi \cdot \nabla \delta \psi \,\dd x = 0		&	\forall \delta \psi \in \mathcal{V},
\\
\int_{{\mathbb{R}^n}} \abs{\nabla \psi}^2 \rho^{(k)} \,\dd x =c,
\\
\end{dcases}
\end{equation}
where $\mathcal{V}\coloneqq H^1({\mathbb{R}^n})$. If we integrate by part the second term in the first equation of~\eqref{eq:wasserstein_metric_step4}, we find\footnote{We here exploit the fact that the space $\mathcal{V}$ is dense in $L^2$.}
\[
\nabla \cdot
\parens{
\rho^{(k)} \nabla \psi
} = -\frac{1}{\lambda}
\nabla \cdot
\parens[\big]{
\rho^{(k)} F'[\rho^{(k)}]
},
\]
where the constant $\lambda$ is a function of $c$ that is determined by substituting the latter expression in~\eqref{eq:wasserstein_metric_step4}. We conclude that for an appropriate choice of $c$ in~\eqref{eq:wasserstein_gradient} we have
\[
\rho^{(k+1)} =
\rho^{(k)} -
\Delta t
\nabla \cdot \parens[\big]{
\rho^{(k)} \nabla V +
\nabla \rho^{(k)}
},
\]
which is the explicit Euler discretization of the FP equation.

\subsection{Transport of Gaussian random variables}

We denote by $\mathcal{N}(\,\cdot\,, \mu, \sigma^2)$ the pdf associated with a univariate Gaussian random variable with mean $\mu$ and variance $\sigma^2$. We assume that $\rho_0\coloneqq \mathcal{N}(\,\cdot\,, \mu, \gamma^2)$ and $\rho_\infty=\mathcal{N}(\,\cdot\,, 0, \sigma^2)$. It is possible to show that the solution to the FP equation in $\Omega=\mathbb{R}$ is given by (see, e.g., \cite[Chapter~VIII.4]{van1992stochastic})
\begin{equation}\label{eq:FP_gaussian}
\rho^{\fp}(x,t) =
\mathcal{N}\parens[\big]{x, \mu(t), \Sigma(t)}
\quad \text{where} \ \
\mu(t) = c(t) \, \mu, \ \Sigma(t) = c(t) \gamma^2 +
\parens[\big]{1 - c(t)} \sigma^2,
\end{equation}
with $c(t)\coloneqq \exp \parens[\big]{-\frac{2 t}{\sigma^2}}$. For comparison, we recall that the McCann interpolation between the same two Gaussian random variables, which corresponds to the optimal transport (OT) geodesic (shortest path) in the Wasserstein metric (see, e.g., \cite[Remark~2.31]{peyre2019computational}), satisfies
\begin{equation}\label{eq:OT_gaussian}
\rho^{\ot}(x,t) =
\mathcal{N}\parens[\big]{x, \mu(t), \Sigma(t)},
\quad \text{where} \ \ \mu(t) = (1- t)\mu, \ \Sigma(t) =
\parens[\big]{(1-t)\gamma + t\sigma}^2,
\ \ t\in [0,1].
\end{equation}

We observe that the solution to~\eqref{eq:FPequation} connects $\rho_0$ to $\rho_\infty$ in infinite time at non-constant speed, while the optimal transport map achieves the same goal in finite time and at constant speed. Note, nevertheless, that the convergence of $\rho^{\fp}$ to $\rho_\infty$ is exponentially-fast: this observation motivates the use of adaptive time stepping strategies to integrate~\eqref{eq:weak_form}.

\section{Methodology}\label{sec:methods}

\subsection{High-fidelity discretization}\label{sec:FEM}

We rely on the finite element (FE) method to discretize~\eqref{eq:weak_form} in space, and on the Crank--Nicolson method for time integration. Given the time grid $\braces[\big]{t^{(k)} = k \Delta t}_{k\in \mathbb{N}}$ and the FE space $\widehat{\mathcal{V}} \subset \mathcal{V}$, we define the sequence of approximate solutions $\{ \widehat{\rho}^{(k)} \}_{k\in \mathbb{N}} \subset \widehat{\mathcal{V}}$ such that $\widehat{\rho}^{(0)}= \rho_0$ and
\begin{equation}\label{eq:FEM_weak}
\int_\Omega \parens[\Bigg]{
\frac{\widehat{\rho}^{(k+1)} - \widehat{\rho}^{(k)}}{\Delta t}
}
\cdot v + \frac{1}{2} \parens[\big]{
\nabla \widehat{\rho}^{(k+1)}+ \widehat{\rho}^{(k+1)} \nabla V +
\nabla \widehat{\rho}^{(k)}+ \widehat{\rho}^{(k)} \nabla V
}
\cdot \nabla v \,\dd x = 0
\quad
\forall v\in \widehat{\mathcal{V}},
\end{equation}
for $k=0,1,\dotsc$ We remark that problem~\eqref{eq:FEM_weak} might feature strong advection: for this reason, in the numerical experiments, we resort to a SUPG stabilization of~\eqref{eq:FEM_weak} to avoid spurious oscillations~\cite{brooks1982streamline}; we refer to~\cite[Chapter~12.8.7]{quarteroni2009numerical} for the detailed presentation of the scheme and for the choice of the parameters. We further remark that our method is not guaranteed to preserve positivity; in the numerical experiments, we hence choose the time step and mesh size sufficiently small to ensure positivity through accuracy.

Given the particles $\{ X_j^{(0)} \}_{j=1}^N$, we consider two distinct strategies to generate the approximate trajectories: the \emph{``ODE'' method} and the
\emph{``gradient flow (GF)'' method}. The former is justified by property~\eqref{th:FPequation_4} in Theorem~\ref{th:FPequation}, while the latter is motivated by the discussion in Section~\ref{sec:gradient_flows}.

\begin{enonce*}[remark]{ODE method}
In this approach, we discretize the ODE system
\begin{equation}\label{eq:ODE_refined}
\begin{cases}
\dot{X}_i(\tau) = \widehat{\mathbf{u}}\parens[\big]{{X}_i(\tau), \tau} 	& \tau \in (t^{(k)}, t^{(k+1)}),\\
{X}_i(t^{(k)}) = {X}_i^{(k)}, \\
\end{cases}
\end{equation}
where $\widehat{\mathbf{u}}(\,\cdot\,, \tau) = - \nabla
\parens[\big]{V + \log \widehat{\rho}(\,\cdot\,, \tau)}$ and $\widehat{\rho}(\,\cdot\,, \tau)$ is a suitable approximation of the density in the interval $(t^{(k)}, t^{(k+1)})$. In the numerical experiments, we consider an explicit Euler discretization of~\eqref{eq:ODE_refined} based on a single time step:
\begin{equation}\label{eq:explicit_euler}
X_i^{(k+1)} = X_i^{(k)} +
\Delta t \,
\widehat{\mathbf{u}}^{(k)}(X_i^{(k)}),
\quad \text{where} \ \ \widehat{\mathbf{u}}^{(k)} \coloneqq - \nabla \parens[\big]{V +
\log
\parens{\widehat{\rho}^{(k)}}
},
\ \ i=1,\dotsc,N, \ \ k=0,1,\dotsc
\end{equation}
We also consider an explicit second-order Runge--Kutta (RK2) scheme for~\eqref{eq:ODE_refined} based on the linear interpolation of the density field
\[
\widehat{\rho}(\,\cdot\,, \tau)=
\parens[\Bigg]{1 - \frac{\tau -t^{(k)}}{\Delta t}
}
\widehat{\rho}^{(k)}
\, + \,
\parens[\Bigg]{
\frac{\tau -t^{(k)}}{\Delta t}
}
\widehat{\rho}^{(k+1)},\quad
\tau \in (t^{(k)}, t^{(k+1)}).
\]
To further increase accuracy, we also lower the time step used for~\eqref{eq:ODE_refined}.
\end{enonce*}

\begin{enonce*}[remark]{GF method}
Given the estimates $ \widehat{\rho}^{(k)}, \widehat{\rho}^{(k+1)}$ obtained through~\eqref{eq:FEM_weak}, we compute the potential $\widehat{\psi}^{(k+1/2)} \in \widehat{\mathcal{V}}$ such that
\begin{subequations}\label{eq:GF_method}
\begin{equation}\label{eq:GF_method_a}
\int_\Omega
\parens[\Bigg]{
\frac{\widehat{\rho}^{(k)} +\widehat{\rho}^{(k+1)}}{2}
} \,
\nabla \widehat{\psi}^{(k+1/2)} \cdot \nabla v \,\dd x =
\int_\Omega
\parens{
\widehat{\rho}^{(k+1)} -
\widehat{\rho}^{(k)}} v \,\dd x,
\quad
\forall v\in \widehat{\mathcal{V}},
\end{equation}
and then we estimate the positions of the particles as
\begin{equation}
X_i^{(k+1)} = X_i^{(k)} +
\,
\widehat{\mathbf{u}}^{(k+1/2)}(X_i^{(k)}),
\quad \text{where} \ \
\widehat{\mathbf{u}}^{(k+1/2)} \coloneqq \nabla \psi^{(k+1/2)},
\ \ i=1,\dotsc,N, \ \ k=0,1,\dotsc
\end{equation}
\end{subequations}
Note that~\eqref{eq:GF_method_a} corresponds to the linearization of the Monge--Ampère equation; in future works, we shall exploit this analogy to determine higher-order approximations of the trajectories of the particles. We observe that the linear system associated with~\eqref{eq:GF_method_a} is poorly conditioned when $\inf_{x\in \Omega} \widehat{\rho}^{(k)}$ is small. In the numerical examples, we consider the regularized problem:
\begin{equation}\label{eq:GF_method_regularized}
\int_\Omega
\parens[\big]{
\widehat{\rho}^{(k)}
\nabla \widehat{\psi}^{(k)} \cdot \nabla v \, + \, \epsilon_{\reg} \,
\widehat{\psi}^{(k)} \, v
}
\,\dd x =
\int_\Omega
\parens{
\widehat{\rho}^{(k+1)} -
\widehat{\rho}^{(k)}} v \,\dd x,
\quad
\forall v\in \widehat{\mathcal{V}},
\end{equation}
where $\epsilon_{\reg} >0$ is a regularization parameter that is set equal to $10^{-10}$ in the numerical experiments.
\end{enonce*}

\subsection{Treatment of boundaries}\label{sec:bnd_regularization}

For several applications, we envision that it is necessary to control the distance of the particles from the boundaries of the domain. In our setting, this can be achieved by introducing a regularization in the potential $V$ at the price of introducing a controlled bias. In more detail, we consider the regularized potential
\begin{equation}\label{eq:regularized_potential}
V_\epsilon(x) = - \log \parens[\big]{\rho_\infty(x)}
\; + \;
\frac{\epsilon}{w_\delta(x)},
\end{equation}
where $w_\delta \colon \Omega \to \mathbb{R}_+$ is a (regularized) distance function from the boundary of $\Omega$. The second term in~\eqref{eq:regularized_potential} introduces a repulsive force that is active in the proximity of the boundary and prevents particles from approaching the boundary too closely. The constant $\epsilon>0$ is a regularization parameter whose choice is investigated in the numerical examples, while the constant $\delta>0$ is the smoothing parameter that is introduced below.

In order to define the function $w_\delta$ in~\eqref{eq:regularized_potential}, we first introduce the distance function $w(x) = \max\big(\dist(x,\partial \Omega)$, $\texttt{tol} \big)$, where the user-defined tolerance $\texttt{tol}$ prevents division by zero in~\eqref{eq:regularized_potential}, for the nodes on the obstacle. Since the distance function $w$ is in general only Lipschitz continuous, we apply a smoothing to ensure that $\nabla V$ is sufficiently smooth. We consider a non-conforming approximation of the biharmonic equation for smoothing~\cite{mozolevski2007hp}. We denote by $\texttt{F}_{j}$ the j-th facet of the mesh, by $\mathbf{n}^+$ the positive normal to each facet, by $ \texttt{D}_{k}$ the $k$-th element of the mesh, by $\texttt{I}_{\internal}$ the set of internal facets of the mesh. Given the field $w \colon \Omega \to \mathbb{R}$, we define the positive and the negative limits on the facets of the mesh
\[
w^\pm(x)\coloneqq \lim_{\epsilon\to 0^+} w(x\mp \epsilon \mathbf{n});
\]
then, we introduce the normal jump of the gradient at each facet $\llbracket \nabla w \rrbracket = \mathbf{n}^+ \cdot (\nabla w^+ -
\nabla w^-)$ and the facet average $\{w\} = \frac{1}{2} (w^+ + w^-)$. Exploiting the previous definitions, we define the smoothing problem: find $w_\delta \in \widehat{\mathcal{V}}$ such that $w_{\delta\rvert_{\partial \Omega}}=\texttt{tol}$ and
\begin{subequations}\label{eq:smoothing_problem}
\begin{equation}
a(w_\delta, v)
\,+ \,
\frac{1}{\delta} \int_\Omega (w_\delta - w) \cdot v \,\dd x = 0
\quad
\forall v\in \widehat{\mathcal{V}} \; : \; \restr{v}{\partial \Omega}=0,
\end{equation}
where
\begin{equation}
a(w,v) =
\sum_{k=1}^{N_{\mathrm{e}}}
\int_{ \texttt{D}_{k}} \parens{\Delta w \cdot \Delta v} \,\dd x + \sum_{j \in \texttt{I}_{\internal}} \int_{ \texttt{F}_{j}}
\parens[\Big]{
\beta_j
\llbracket \nabla w \rrbracket \cdot
\llbracket \nabla v \rrbracket +
\tfrac{1}{\beta_j}
\{ \Delta w \} \cdot \{ \Delta v \}} \,\dd s,
\end{equation}
and $\beta_j = \sigma_{\beta} \kappa^2 \abs{\texttt{F}_{j}}^{-1}$ with $ \sigma_{\beta} =10$, and $\kappa$ equal to the polynomial degree of the FE basis. We refer to~\cite{mozolevski2007hp} for a detailed analysis of the formulation.
\end{subequations}

\subsection{Estimate of the target density}\label{sec:target_density}

The application of the FP equation to PSR requires the estimate of the densities $\rho_0$ and $\rho_\infty$. Given the initial and the target point clouds $\{ X_i^{0} \}_{i=1}^N$ and $\{ X_j^{\infty} \}_{j=1}^M$, we fit two Gaussian mixture models (GMMs) to determine the estimates of $\rho_0$ and $\rho_\infty$ that enter in~\eqref{eq:FEM_weak} (see, e.g., \cite[Chapter~8]{hastie2009elements}). We resort to the Matlab function \texttt{fitgmdist} to determine the GMMs associated with the point clouds. We determine the number of mixture models using the Aikake information criterion (AIC, see~\cite[Chapter~7.5]{hastie2009elements}); furthermore, we introduce a regularization of the covariance matrix to robustify the learning procedure\footnote{We consider the regularization value $10^{-2}$; see the Matlab documentation for further details.}.

The use of GMMs for density estimation ensures that $\rho_0$ and $\rho_\infty$ are strictly positive for all $x\in \Omega$, which is crucial to apply the results of Theorem~\ref{th:FPequation}. Furthermore, we recall that for Gaussian distributions $\rho = \mathcal{N}(\,\cdot\,; \mu, \Sigma)$ we have
\begin{equation}\label{eq:log_gaussian}
\log \parens[\big]{
\mathcal{N}(x; \mu, \Sigma)
} = -\frac{1}{2} (x-\mu)^\top \Sigma^{-1} (x-\mu).
\end{equation}
{Since $\Sigma$ is a $n\times n$ symmetric positive definite matrix,} this implies that for GMMs the gradient of $\log(\rho)$ grows linearly with respect to the squared distance from the centers of the mixture models.

\section{Numerical results}\label{sec:numerics}

We present two numerical experiments that involve the transport of point clouds past a bluff body. The examples are designed to illustrate the ability of our method to deform point clouds past two-dimensional objects. As discussed in~\cite{iollo2022mapping}, this tool is central in the development of model reduction techniques based on coordinate transformations. In~\cite{iollo2022mapping}, the point clouds are associated to coherent structures of the solution field (e.g., shock waves) that are displaced due to parameter variations.

\subsection{Transport of Gaussian distributions across a cylinder}\label{sec:numerics_test1}

\subsubsection{Problem setting}

We consider the problem of transporting the Gaussian density $\rho_0 = \mathcal{N}(\,\cdot\,; \mu_0, \Sigma_0)$ with $\mu_0 = [-2,0]$ and $\Sigma_0 = 0.2 \mathbbm{1}$ to $\rho_\infty = \mathcal{N}(\,\cdot\,; \mu_\infty, \Sigma_\infty)$ with $\mu_\infty = [2,0]$ and $\Sigma_\infty = \Sigma_0$ in the domain $\Omega = (-4,4)^2 \setminus \mathcal{B}_{r=0.5}(0)$. Figure~\ref{fig:test1_vis} shows the initial and the target distributions.

\begin{figure}[h!]
\begin{subfigure}{0.49\textwidth}
\centering
\includegraphics[width=.9\textwidth]{imm/test1/rho0.png}
\caption{}\label{fig:test1_vis_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\centering
\includegraphics[width=.9\textwidth]{imm/test1/rho_inf.png}
\caption{}\label{fig:test1_vis_b}
\end{subfigure}
\caption{Transport of Gaussian distributions across a cylinder. (a)~Initial distribution. (b)~Tar\-get distribution.}\label{fig:test1_vis}
\end{figure}

We resort to a P2 FE discretization with $N_{\hf}=\np{6469}$ degrees of freedom; we integrate the equation in the time interval $(0,T=5)$ and we consider the grid
\begin{equation}\label{eq:non_uniform_timegrid}
\braces[\Big]{ t^{(k)} = T \parens[\big]{\tfrac{k}{K}}^{1.5}}_{k=0}^K, \quad K=\np{3000}.
\end{equation}
The non-homogeneous time discretization is motivated by the observation that the dynamics is much faster at the onset of the simulation. The exponent 1.5 was determined empirically through systematic experimentation; in the future, we plan to implement more rigorous adaptive time stepping strategies. The regularized distance function in~\eqref{eq:regularized_potential} is computed using $\delta=10^{-2}$ and $\texttt{tol}=10^{-2}$ and several values of $\epsilon$.

We further remark that $\nabla V$ is computed by differentiating the FE approximation of $V= - \log(\rho_\infty)$: we notice that approximating the density $\rho_\infty$ through its nodal representation $\rho_\infty^{\fe}$ and then computing $\nabla V$ as $\nabla \rho_\infty^{\fe}/\rho_\infty^{\fe}$ leads to significantly less accurate results. Recalling~\eqref{eq:log_gaussian}, for Gaussian distributions, $V= - \log \parens{\rho_\infty}$ is a second-order polynomial whose gradient can be computed exactly. Furthermore, we empirically found that the computation of $\nabla \rho_\infty^{\fe}/\rho_\infty^{\fe}$ is numerically unstable when $\rho_\infty$ is very small.

To validate the PSR procedure, we generate $X_1^0,\dotsc,X_N^0 \overset{\iid}{\sim} \rho_0$ with $N=100$ and we compare the trajectories obtained using our method with $M=100$ i.i.d. samples from the target distribution. Time integration of the particles is performed using the ODE method based on the explicit Euler method~\eqref{eq:explicit_euler} with the same time step that is used for the integration of~\eqref{eq:FEM_weak}. We also considered the GF method and the ODE method with RK2 time integration: since the results of the latter two methods are nearly equivalent to the ones obtained using~\eqref{eq:explicit_euler}, we do not report them here.

\subsubsection{Results}

Figure~\ref{fig:test1_error_analysis} shows the performance of our method. Figure~\ref{fig:test1_error_analysis_a} shows the behavior of the $L^1$ error $\norm[\big]{\rho(\,\cdot\,, t) - \rho_\infty}_{L^1(\Omega)}$ for several values of $\epsilon$: as we increase $\epsilon$ in~\eqref{eq:regularized_potential}, we introduce an increasing bias in the formulation that prevents convergence to the target density. Figure~\ref{fig:test1_error_analysis_b} compares the deformed point cloud at the final time $\{ X_j^{(K)} \}_{j}$ for $\epsilon=0$ with the ``target'' point cloud generated using $\rho_\infty$: we observe that the results are extremely satisfactory for this model problem. Figure~\ref{fig:test1_error_analysis_c} shows the behavior of select trajectories, while Figure~\ref{fig:test1_error_analysis_d} shows the trajectory of one particle for several values of $\epsilon$: we notice that, as we increase $\epsilon$, the particle moves away from the cylinder due to the effect of the repulsive force.

\begin{figure}[h!]
\begin{subfigure}{0.49\textwidth}
\centering
\includegraphics[width=.9\textwidth]{imm/test1/L1error.png}
\caption{}\label{fig:test1_error_analysis_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\centering
\includegraphics[width=.9\textwidth]{imm/test1/initial_final_points.png}
\caption{}\label{fig:test1_error_analysis_b}
\end{subfigure}
\hfill
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\centering
\includegraphics[width=0.9\textwidth]{imm/test1/select_trajectories.png}
\caption{}\label{fig:test1_error_analysis_c}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\centering
\includegraphics[width=.9\textwidth]{imm/test1/detailed_trajectories.png}
\caption{}\label{fig:test1_error_analysis_d}
\end{subfigure}
\caption{Transport of Gaussian distributions across a cylinder. (a)~Behavior of the $L^2$ error $\norm[\big]{\rho(\,\cdot\,, t) - \rho_\infty}_{L^1(\Omega)}$ for several choices of the regularization parameter $\epsilon$. (b)~Deformed point cloud at final time for $\epsilon=0$ and target point cloud. (c)~Behavior of several trajectories for $\epsilon=0$. (d)~Behavior of one trajectory for several values of $\epsilon=0$.}\label{fig:test1_error_analysis}
\end{figure}

\subsection{Point-set registration across a cylinder}\label{sec:numerics_test2}

\subsubsection{Problem setting}

We consider the problem of continuously deforming the point cloud $\{X_i^{0}\}_{i=1}^N$ to match the target point cloud $\{X_j^{\infty}\}_{j=1}^M$. We consider the same domain $\Omega$ of the previous example and we define
\begin{equation}\label{eq:point_clouds}
\begin{dcases}
X_i^0 = \bracks[\Big]{\cos\parens[\big]{\theta_0 + \Delta \theta \tfrac{i-1}{N-1}} - 1, \sin\parens[\big]{\theta_0 + \Delta \theta \tfrac{i-1}{N-1}}} + 0.1 \widetilde{X}_i^0 			& i=1,\dotsc,N;\\
X_j^\infty = \bracks[\Big]{\cos\parens[\big]{\theta_\infty + \Delta \theta \tfrac{j-1}{M-1}}+1, \sin\parens[\big]{\theta_\infty + \Delta \theta \tfrac{j-1}{N-1}}} + 0.1 \widetilde{X}_j^\infty 	& j=1,\dotsc,M;\\
\end{dcases}
\end{equation}
where $N=M=141$, $\theta_0=\frac{\pi}{2}$, $\theta_\infty=\frac{3\pi}{2}$, $\Delta \theta = \pi$ and $\widetilde{X}_i^0,
\widetilde{X}_j^\infty
\overset{\iid}{\sim} \Uniform\parens[\big]{(0,1)^2}$. Figure~\ref{fig:test2_vis_a} shows the reference and the target point clouds. We consider two distinct discretizations of increasing size. The coarse discretization features a P2 FE mesh with $N_{\hf}=\np{23409}$ degrees of freedom and $K=5\cdot 10^3$ time steps; the fine discretization features a P2 FE mesh with $N_{\hf}=\np{93134}$ degrees of freedom and $K=10^4$ time steps. In both cases, we integrate the system in the time interval $(0,T=15)$ and we consider the non-uniform time grid~\eqref{eq:non_uniform_timegrid}.

We rely on the unregularized potential (i.e., $\epsilon=0$ in~\eqref{eq:regularized_potential}). We compare the ODE method based on forward Euler time integration with the GF method. We also run simulations for the ODE method with RK2 time integration with two time steps per time interval: since the results of the latter are comparable with the ones of the ODE method, they are not reported below. We model the point clouds~\eqref{eq:point_clouds} using two GMMs with four components --- we recall that the number of components is selected by maximizing the AIC criterion. Figures~\ref{fig:test2_vis_b} and~\ref{fig:test2_vis_c} show the GMM densities that are obtained by the learning procedure.

\begin{figure}[h!]
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/point_clouds.png}
\caption{}\label{fig:test2_vis_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/rho0.png}
\caption{}\label{fig:test2_vis_b}
\end{subfigure}
\hfill
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/rho_inf.png}
\caption{}\label{fig:test2_vis_c}
\end{subfigure}
\caption{Point-set registration across a cylinder. (a)~Reference (red) and target point (blue) clouds. (b)~Initial condition $\rho_0$. (c)~Target density $\rho_\infty$.}\label{fig:test2_vis}
\end{figure}

\subsubsection{Results}

Figure~\ref{fig:test2_error_analysis_a} shows the behavior of the $L^1$ error $\norm[\big]{\rho(\,\cdot\,, t) - \rho_\infty}_{L^1(\Omega)}$ over time for the two discretizations; Figures~\ref{fig:test2_error_analysis_b} and~\ref{fig:test2_error_analysis_c} show the two FE meshes. The error in the target density decreases by a factor $4$: it is hence inversely proportional to the size of the mesh, as expected for second-order discretizations.

\begin{figure}[h!]
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/L1error.png}
\caption{}\label{fig:test2_error_analysis_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/mesh1.png}
\caption{}\label{fig:test2_error_analysis_b}
\end{subfigure}
\hfill
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/mesh2.png}
\caption{}\label{fig:test2_error_analysis_c}
\end{subfigure}
\caption{Point-set registration across a cylinder. (a)~Behavior of the $L^1$ error $\norm[\big]{\rho(\,\cdot\,, t) - \rho_\infty}_{L^1(\Omega)}$ for both coarse and fine meshes. (b)~Coarse mesh ($N_{\hf}=\np{23409}$). (c)~Fine mesh ($N_{\hf}= \np{93134}$).}\label{fig:test2_error_analysis}
\end{figure}

Figure~\ref{fig:test2_densities} shows the behavior of the solution to~\eqref{eq:FEM_weak} for three time instants. As expected, we notice that the solution is symmetric with respect to the cylinder. We further notice that for small time steps the solution has three peaks: the majority of the mass --- corresponding to the top and the bottom peaks at $t=0.1$ --- rapidly goes through the cylinder, while a small portion of the mass --- corresponding to the weaker peak on the left of the cylinder at $t=0.1$ --- reaches the target density later.

\begin{figure}[h!]
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/density_time1.png}
\caption{$t=0.1$}\label{fig:test2_densities_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/density_time2.png}
\caption{$t=0.5$}\label{fig:test2_densities_b}
\end{subfigure}
\hfill
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/density_time3.png}
\caption{$t=1.5$}\label{fig:test2_densities_c}
\end{subfigure}
\caption{Point-set registration across a cylinder; density profiles for three time instants.}\label{fig:test2_densities}
\end{figure}

Figures~\ref{fig:test2_point_clouds_ODE} and~\ref{fig:test2_point_clouds_GF} show the behaviors of the point clouds for the two particle integration methods discussed in Section~\ref{sec:methods}: we emphasize that both approaches consider the same particles and the same estimates of the initial and target distributions. The results correspond to the fine discretization. Interestingly, the qualitative behavior of the two methods is significantly different, particularly for small values of $t$: in the ODE method the point cloud is nearly split in two parts at $t=0.1$ and all the points proceed towards the target cloud with similar speeds; instead, in the GF method, we can clearly identify three clusters whose elements move towards the target cloud with different speeds. The latter behavior appears to be more consistent with the results for the density in Figure~\ref{fig:test2_densities}.

\begin{figure}[h!]
\begin{subfigure}{.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/ODE/pointclouds_time1.png}
\caption{$t=0.1$}\label{fig:test2_point_clouds_ODE_a}
\end{subfigure}
\hfill
\begin{subfigure}{.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/ODE/pointclouds_time2.png}
\caption{$t=0.5$}\label{fig:test2_point_clouds_ODE_b}
\end{subfigure}
\hfill
\begin{subfigure}{.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/ODE/pointclouds_time3.png}
\caption{$t=1.5$}\label{fig:test2_point_clouds_ODE_c}
\end{subfigure}
\caption{Point-set registration across a cylinder; deformed point clouds $\braces[\big]{X_i(t)}_{i=1}^N$ for three time instants. ODE method.}\label{fig:test2_point_clouds_ODE}
\end{figure}

\begin{figure}[h!]
\begin{subfigure}{.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/GF/pointclouds_time1.png}
\caption{$t=0.1$}\label{fig:test2_point_clouds_GF_a}
\end{subfigure}
\hfill
\begin{subfigure}{.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/GF/pointclouds_time2.png}
\caption{$t=0.5$}\label{fig:test2_point_clouds_GF_b}
\end{subfigure}
\hfill
\begin{subfigure}{.32\textwidth}
\centering
\includegraphics[width=\textwidth]{imm/test2/GF/pointclouds_time3.png}
\caption{$t=1.5$}\label{fig:test2_point_clouds_GF_c}
\end{subfigure}
\caption{Point-set registration across a cylinder; deformed point clouds $\braces[\big]{X_i(t)}_{i=1}^N$ for three time instants. GF method.}\label{fig:test2_point_clouds_GF}
\end{figure}

Figure~\ref{fig:test2_point_clouds_ODE_vs_GF} compares the point clouds obtained with the two methods at the final time. We notice that the ODE method does not provide satisfactory results for this test case, while the GF method is significantly more accurate.

\begin{figure}[h!]
\begin{subfigure}{.49\textwidth}
\centering
\includegraphics[width=.67\textwidth]{imm/test2/ODE/initial_final_points.png}
\caption{ODE}\label{fig:test2_point_clouds_ODE_vs_GF_a}
\end{subfigure}
\hfill
\begin{subfigure}{.49\textwidth}
\centering
\includegraphics[width=.67\textwidth]{imm/test2/GF/initial_final_points.png}
\caption{GF}\label{fig:test2_point_clouds_ODE_vs_GF_b}
\end{subfigure}
\caption{Point-set registration across a cylinder; deformed point clouds $\braces[\big]{X_i(t)}_{i=1}^N$ for $t=T_{\max}$ for the two time integration methods.}\label{fig:test2_point_clouds_ODE_vs_GF}
\end{figure}

\section{Conclusions and perspectives}\label{sec:conclusions}

We developed a PSR method for the interpolation of point clouds in bounded domains based on the Fokker--Planck equation. Our method relies on GMMs for density estimation and on a standard (Eulerian) discretization of the FP equation for the evaluation of the density. We also proposed two distinct strategies for the estimation of the trajectories of the particles, which are both rigorously justified by the analysis of Section~\ref{sec:mathematical_background}.

The numerical experiments of Section~\ref{sec:numerics} showed the potential of the method for PSR. The regularized potential~\eqref{eq:regularized_potential} provides a simple yet effective way to control the distance of the particles from the boundaries at all times; furthermore, in both numerical experiments, the convergence to the target distribution is rapid and monotonic. As regards the motion of particles, the results of Section~\ref{sec:numerics_test2} clearly showed the superiority of the GF approach based on the solution to a linearized OT problem over the ODE approach based on the explicit integration of the characteristics.

Despite the linearity of the equation, the solution to~\eqref{eq:weak_form} using the FE method on a fixed grid is expensive and demands further investigations. We shall devise a specialized solver for the FP equation with emphasis on the accurate approximation of the particles' dynamics. Specifically, we plan to consider (i) adaptive time stepping schemes and (ii) adaptive spatial discretizations to reduce the overall costs: we expect to borrow ideas from deterministic particle methods for the FP equation~\cite{russo1990deterministic,russo1993deterministic,russo1990particle}. In the future, we also plan to investigate variants of~\eqref{eq:regularized_potential} to effectively control the distance of the particles from the boundaries at a reduced impact on the limit distribution. Finally, we plan to assess the robustness of our approach with respect to the choice of the densities $\rho_0$, $\rho_\infty$.

\appendix

\section{Proofs}\label{sec:proofs}

\subsection{Proof of Theorem~\ref{th:KL}}

\begin{proof}[Proof of Theorem~\ref{th:KL}(\ref{th:KL_1})]
We observe that
\begin{equation}\label{eq:silly_identity}
g(t)\coloneqq t - 1 - \log (t) \geq 0 \quad
\forall t \in \mathbb{R}_+.
\end{equation}
To prove~\eqref{eq:silly_identity}, we notice that $t=1$ is a global minimum of $g$ in $\mathbb{R}_+$ and $g(1)=0$. If we define the set $A\coloneqq \braces[\big]{ x\in \Omega \,:\, \rho(x)>0}\subset \Omega$, we find
\[
-D_{\KL}(\rho \parallel \rho_\infty) =
\int_\Omega \rho \, \log \parens[\Bigg]{\frac{\rho_\infty}{\rho}} \,\dd x =
\int_A \rho \, \log \parens[\Bigg]{\frac{\rho_\infty}{\rho}} \,\dd x;
\]
then, exploiting~\eqref{eq:silly_identity} with $t=\frac{\rho_\infty}{\rho}$, we find
\[
-D_{\KL}(\rho \parallel \rho_\infty) =
\int_A \rho \, \log \parens[\Bigg]{\frac{\rho_\infty}{\rho}}
\leq
\int_A \rho \, \parens[\Bigg]{\frac{\rho_\infty}{\rho} - 1} =
\int_A \rho_\infty - 1
\leq
\int_\Omega \rho_\infty - 1 =0,
\]
which proves $D_{\KL}(\rho \parallel \rho_\infty)\geq 0$.

Let $\rho \in \mathcal{P}_2(\Omega)$ satisfy $D_{\KL}(\rho \parallel \rho_\infty)=0$. If we denote by $\phi \colon t \mapsto - \log (t)$ and $f \colon x \mapsto \frac{\rho(x)}{\rho_\infty(x)}$, we find
\[
\phi \parens[\Bigg]{\int_{\Omega} f(x) \rho(x) \,\dd x} = \phi \parens{1} = 0,
\qquad
\int_\Omega \phi \circ f (x) \rho(x) \,\dd x = - D_{\KL}(\rho \parallel \rho_\infty) = 0.
\]
Notice that $\phi \parens[\big]{\int_{\Omega} f(x) \rho(x) \,\dd x} =
\int_\Omega \phi \circ f (x) \rho(x) \,\dd x$: since $f$ is convex, recalling Jensen's inequality, we must have that $f$ is constant, which implies that $\rho= \rho_\infty$ a.s.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{th:KL}(\ref{th:KL_2})]
It suffices to show an example for which $D_{\KL}(\rho \parallel \rho_\infty) \neq D_{\KL}(\rho_\infty \parallel \rho)$. Towards this end, consider $\Omega=(-1,1)$ and
\[
\rho(x) \equiv \frac{1}{2},
\quad
\rho_\infty(x) =
\begin{cases}
0.1 & x\in (-1,0),\\
0.9 & x \in (0,1).
\end{cases}
\]
By straightforward calculations, we find $D_{\KL}(\rho \parallel \rho_\infty) = 0.5108$ and $D_{\KL}(\rho_\infty \parallel \rho) = 0.3681$.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{th:KL}(\ref{th:KL_3})]
We define $\rho_\lambda \coloneqq \lambda \rho_0 + (1-\lambda) \rho_1$. Since the function $t\mapsto t \log(t/c)$ is strictly convex for any $c>0$, we have
\[
\rho_\lambda(x) \log \parens[\Bigg]{\frac{\rho_\lambda(x)}{\rho_\infty(x)}}
\leq
\lambda \rho_0(x) \log \parens[\Bigg]{\frac{\rho_0(x)}{\rho_\infty(x)}} + (1- \lambda) \rho_1(x) \log \parens[\Bigg]{\frac{\rho_1(x)}{\rho_\infty(x)}},
\quad
\forall x\in \Omega, \ \
\lambda\in (0,1);
\]
furthermore, the latter holds with equality if and only if $\rho_0(x)=\rho_1(x)$. Recalling the expression of the KL divergence, we obtain the desired result.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{th:KL}(\ref{th:KL_4})]
We can exploit~\cite[Theorem~3]{gilardoni2010pinsker} to show that $\norm{\rho - \rho_\infty}_{L^1(\Omega)}
\leq
\sqrt{ 2 D_{\KL}(\rho \parallel \rho_\infty)}$; then, the desired result follows from the symmetry of the norm.
\end{proof}

\subsection{Proof of Theorem~\ref{th:FPequation}}

\begin{enumerate}
\item[\eqref{th:FPequation_1}] The existence and uniqueness of the solution to~\eqref{eq:weak_form} can be studied using the theory of abstract parabolic problems. We refer to~\cite[Chapter~10, Theorem~10.11]{salsa2015partial} for a detailed proof.
\item[\eqref{th:FPequation_2}] Exploiting~\cite[Theorem~10.11]{salsa2015partial}, we find that $\rho(\,\cdot\,, t) \in L^2(\Omega)$ for all $t>0$; it hence suffices to show that $\int_\Omega \rho(x,t) \,\dd x = 1$ for all $t>0$ and that $\rho(x,t) \geq 0$ for a.e.\ $x\in \Omega$ and $t>0$. If we consider $v=1$ in~\eqref{eq:weak_form}, we find
\[
0=
\int_\Omega \partial_t \rho(x,t) \,\dd x =
\frac{\dd}{\dd t}
\int_\Omega \rho(x,t) \,\dd x \quad \Longrightarrow \quad
\int_\Omega \rho(x,t) \,\dd x =
\int_\Omega \rho_0(x) \,\dd x = 1,
\]
which proves conservation of mass. The proof of positivity is significantly more technical and is provided in~\cite[Corollary~4.3]{ouhabaz2009analysis}.
\item[\eqref{th:FPequation_3}] Recalling notation introduced in~\eqref{eq:KLdivergence_revisited}, since $\rho \geq 0$ a.e., we find that
\[
\frac{\dd}{\dd t}D_{\KL}\parens[\big]{\rho(\,\cdot\,,t) \parallel \rho_\infty} = - \int_\Omega \rho(x,t) \abs[\big]{\nabla F'\bracks[\big]{\rho(x,t)}}^2 \,\dd x.
\]
Provided that $c\coloneqq \inf_{x\in \Omega} \rho_\infty(x) $ is strictly positive and $\rho_0 \in L^2(\Omega)$, we have that $D_{\KL}(\rho_0 \parallel \rho_\infty)<\infty$: we have indeed
\[
\begin{split}
D_{\KL}(\rho_0 \parallel \rho_\infty)
	& \leq \int_\Omega \rho_0 \log \parens[\bigg]{\frac{\rho_0}{c}}
\\	& = \int_\Omega \rho_0 \log \rho_0 \,\dd x - \log(c)
\\	& \le \int_\Omega (\rho_0 ^2 + \rho_0) \,\dd x - \log(c)
\\	& = \norm{\rho_0}_{L^2(\Omega)}^2 +1 - \log(c),
\end{split}
\]
where in the second inequality we used the identity $t\log (t) \leq t^2+t$ that is valid for all $t\in \mathbb{R}_+$. Therefore, since $\rho(\,\cdot\,, t)\in \mathcal{P}_2(\Omega)$ for all $t>0$, we find that the KL divergence is monotonically-decreasing and uniformly bounded. Exploiting the Pinsker's inequality, we also find that $t\mapsto \norm[\big]{\rho(\,\cdot\,, t)}_{L^1(\Omega)}$ is uniformly bounded. Then, we can apply~\cite[Theorem~2.3]{alves2024strong} to conclude that $\rho(\,\cdot\,,t) \to \rho_\infty$ in $L^1$ for $t\to \infty$.
\item[\eqref{th:FPequation_4}] The result exploits the Lagrangian interpretation of the FP equation, which is extensively discussed in~\cite[Chapter~4]{santambrogio2015optimal} --- the precise result is provided in~\cite[Theorem~4.4]{santambrogio2015optimal}.
\end{enumerate}

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography

\end{document}