# GS-EMA: Integrating Gradient Surgery Exponential Moving Average with
Boundary-Aware Contrastive Learning for Enhanced Domain Generalization in
Aneurysm Segmentation

###### Abstract

The automated segmentation of cerebral aneurysms is pivotal for accurate
diagnosis and treatment planning. Confronted with significant domain shifts
and class imbalance in 3D Rotational Angiography (3DRA) data from various
medical institutions, the task becomes challenging. These shifts include
differences in image appearance, intensity distribution, resolution, and
aneurysm size, all of which complicate the segmentation process. To tackle
these issues, we propose a novel domain generalization strategy that employs
gradient surgery exponential moving average (GS-EMA) optimization technique
coupled with boundary-aware contrastive learning (BACL). Our approach is
distinct in its ability to adapt to new, unseen domains by learning domain-
invariant features, thereby improving the robustness and accuracy of aneurysm
segmentation across diverse clinical datasets. The results demonstrate that
our proposed approach can extract more domain-invariant features, minimizing
over-segmentation and capturing more complete aneurysm structures.

Index Terms— Domain Generalization, Gradient Surgery, Contrastive Learning

##  1 Introduction

The accurate segmentation of cerebral aneurysms is vital for diagnosing and
treating patients effectively. This process is not just about detecting
aneurysms early; it involves precise measurements of their size and shape,
which are critical for formulating treatment plans [1, 2]. However, the
variability in imaging data quality (see Fig.1) from different medical centers
presents a significant challenge, complicating the segmentation process.

This variability necessitates a domain generalization (DG) approach, where a
model trained on data from multiple sources can adapt to new, unseen domains.
The diversity of multi-source data makes DG a daunting challenge in medical
imaging, pushing the need for models that generalize well across different
medical centers and data types.

Fig. 1: Illustration of the variability in imaging data quality from different
medical centers.

Unlike traditional DG approaches such as domain alignment [3], data
augmentation [4], ensemble learning [5], self-supervised learning [6],
disentangled representation learning [7], and others, our method takes a
different approach. We enhance domain generalization by leveraging gradient
surgery exponential moving average (GS-EMA), offering an innovative solution
to address DG challenges.

In deep learning, EMA is a frequently used technique for parameter averaging
in models, aimed at enhancing the generalization performance and stability of
the model. In a teacher-student [8] network setup, the teacher network
undergoes a process of parameter smoothing, driven by the student network.
However, initially, there are no specific conditions set for this transfer.
Consequently, all parameters learned by the student network, whether they are
domain-invariant or domain-specific, are updated into the teacher network at
some rate. This approach poses a challenge as it fails to distinguish between
domain-invariant and domain-specific parameters. To address this issue, we
introduce the concept of gradient surgery.

Deep neural networks are trained using gradient descent, where gradients guide
the optimization process across the landscape defined by the loss function and
training data. The gradient surgery framework [9, 10] aims to resolve
conflicts arising in multi-task learning. The conflicting gradients are
typically averaged to obtain a final gradient for parameter updates. GSMorph
[11] propose alternative methods like normal vector projection to derive the
ultimate gradient for parameter updates. Instead of devising a new projection
method as suggested by others, we approach the problem by analyzing the
relationships between gradients to determine whether EMA parameter updates
should occur.

Additionally, there is a class imbalance problem in 3D data segmentation due
to the small proportion of aneurysms. After multiple downsampling steps, these
small features are more likely to be overlooked in the latent space. To tackle
this, we introduce the concept of boundary-awareness to traditional
contrastive learning [12].

Contributions: Our study introduces innovative techniques that enhance model
adaptability. We integrate gradient surgery with EMA updates, strengthening
the ability of model to learn domain-invariant features. This novel approach
promises to elevate the performance of DG tasks in medical imaging, ensuring
that our model can generalize effectively to new datasets and medical centers.
Additionally, we pioneer the use of boundary-aware contrastive learning,
enabling our model to discern small target features especially for cerebral
aneurysms.

##  2 Methods

Fig. 2 depicts our neural network architecture dedicated for domain
generalization tasks. It initiates with 3D source images, which undergoes
image transformation to produce target images. Once both the source and target
images are obtained, they are separately fed into the encoders of the student
and teacher networks.

After acquiring the latent space features, a boundary-aware contrastive
learning loss is computed. The central notion here is to amalgamate the same
instance subjected to diverse transformations, while distancing different
instances, aiming to grasp instance-aware representations. This contrastive
learning differs from transformation predictions, as it strives to attain
transformation-invariant representations. The latent space features are then
decoded to yield predictions, which are supervised using ground truth.

Within the student network, the green arrow signifies fully supervised
learning for the source images, and the yellow arrow represents the same for
the target images. By analyzing the gradient relationship between these two
losses, a novel GS-EMA strategy is devised to update the parameters of teacher
network. If the gradient angle between the losses is less than 90 degrees, it
indicates that the network has learned domain-invariant features, prompting an
EMA update. Conversely, if the gradient angle exceeds 90 degrees, no EMA
update is performed, as this suggests the network has grasped domain-specific
features, which is not conducive to domain generalization tasks. Ultimately,
after several updates, a teacher network enriched with more domain-invariant
features is achieved, readying it for domain generalization tasks.

Fig. 2: Schematic of the proposed model.

###  2.1 Problem Definition and Data Transformation

Let $\mathcal{X}$ be the input (image) space and $\mathcal{Y}$ be the
segmentation (label) space, a domain is defined as a joint distribution
$P_{XY}$ on $\mathcal{X}\times\mathcal{Y}$. In the context of DG, we have
access to $K$ similar but distinct source domains
$\left\\{\left(x^{k}_{s},y^{k}_{s}\right)\right\\}_{k=1}^{K}$, each associated
with a joint distribution $P_{XY}^{k}$. Note that $P_{XY}^{i}\neq P_{XY}^{j}$
with $i\neq j$ and $i,j\in\\{1,...,K\\}$. The goal of DG is to learn a
predictive model using only source domain data such that the prediction error
on an unseen target domain is minimized.

To enhance the model adaptability to previously unseen data domains, we use
data transformations to simulate the distribution of the target domain data.
The simulated target data is represented as
$\left\\{\left(x_{k}^{t},y_{k}^{t}\right)\right\\}_{k=1}^{K}$. The process of
data transformation encompasses several key steps, including geometric
transformations, intensity alterations, noise injection and smoothing,
histogram shifting, as well as bias field correction. These operations
collectively aim to generate diverse target data, empowering the model with
enhanced generalization capabilities to adapt to various data sources and
target domains.

###  2.2 Gradient Surgery Exponential Moving Average

In a teacher-student network setup, when the student network is tasked with
learning from data originating from different domains, we calculate distinct
losses for each domain in Eq. 2. This allows us to obtain gradient information
specific to each domain. Our fundamental hypothesis is that when the angle
between gradients from different domains is less than 90 degrees, it suggests
that the student network has effectively learned how to extract domain-
invariant features. In such cases, we employ EMA to update the parameters of
student network, subsequently transferring these parameters to the teacher
network. This transfer is performed to better capture universal features.

However, when the angle between gradients from different domains exceeds 90
degrees, it indicates that the student network is primarily focused on
learning domain-specific features. In such scenarios, we abstain from
utilizing EMA for parameter updates and refrain from transmitting these
parameters to the teacher network. This strategic approach ensures that the
student network can efficiently discriminate between features originating from
different domains, enabling it to adapt effectively to the challenges of
multi-task learning.

$\begin{split}\mathcal{L}_{DCE}(p,y)&=\frac{1}{N}\sum_{i=1}^{N}\left(1-\frac{2\left|p\cap
y\right|}{\left|p\right|+\left|y\right|}-y\log p\right)\end{split}$ (1)
$\begin{split}\mathcal{L}_{stu}^{src}=\mathcal{L}_{DCE}(p_{stu}^{s},y^{s}),\
\mathcal{L}_{stu}^{trg}=\mathcal{L}_{DCE}(p_{stu}^{t},y^{t})\end{split}$ (2)

Data: Student network parameters $\theta_{stu}$; Teacher network parameters
$\theta_{tea}$; Loss on source data in student network $\mathcal{L}_{src}$;
Loss on target data in student network $\mathcal{L}_{trg}$; EMA decay
coefficient $\alpha$.

Result: Decide whether updated $\theta_{tea}$ with EMA from $\theta_{stu}$.

1 for  _each mini-batch_ do

2 $\nabla\mathcal{L}_{stu}^{src}\rightarrow g_{src}$;

3 $\nabla\mathcal{L}_{stu}^{trg}\rightarrow g_{trg}$;

4 if  _$\langle g_{src},g_{trg}\rangle\leq 0$_ then

5 $\theta_{tea}^{\prime}=\theta_{tea}\cdot\alpha+(1-\alpha)\cdot\theta_{stu}$;

6

7 else if  _$\langle g_{src},g_{trg}\rangle >0$_ then

8 $\theta_{tea}^{\prime}=\theta_{tea}$;

9

10 Update $\theta_{tea}$ with $\theta_{tea}^{\prime}$;

11 Update $\theta_{stu}$ as needed;

12

Algorithm 1 Gradient Surgery Exponential Moving Average

###  2.3 Boundary-Aware Contrastive Learning

In our study, we tackle the challenge of uneven distribution of classes in the
segmentation of aneurysms by proposing a unique contrastive learning approach
that operates within a teacher-student network configuration. This method
enhances the distinction between matching (positive) and non-matching
(negative) sample pairs by employing a Fourier transformation strategy, which
is particularly adept at isolating high-frequency elements that delineate
boundaries. Transitioning from volume-based to boundary-based analysis ensures
that the presence of small aneurysms is not disproportionately low compared to
larger vessels.

Both the student and teacher branches receive two distinct sets of data: the
original data from the source domain, represented as $x^{s}$, and the
corresponding transformed data $x^{t}$. Consequently, the latent feature
representations from the student network are symbolized as $z_{stu}^{s}$ and
$z_{stu}^{t}$, while those from the teacher network are signified as
$z_{tea}^{s}$ and $z_{tea}^{t}$.

Advancing further, we harness the power of Fourier transformation paired with
a high-frequency filter in Eq. 3 to extract features that are cognizant of the
boundaries within the data. These extracted features from both student and
teacher networks are represented as $b_{stu}^{s}$, $b_{stu}^{t}$,
$b_{tea}^{s}$, and $b_{tea}^{t}$ respectively. Our primary objective within
this feature space is to cultivate instance-specific representations that are
closely aligned when the same instance is encoded differently, while
simultaneously ensuring a clear demarcation between distinct instances,
irrespective of the encoder used.

To clarify the relationships within our contrastive learning framework, we
delineate the instances processed through different encoders as positive pairs
when they originate from the same instance. This includes pairs like
$z_{stu}^{s}$ with $z_{tea}^{s}$, and $z_{stu}^{t}$ with $z_{tea}^{t}$. In
contrast, negative pairs consist of different instances that have been encoded
either by the same or by different encoders, such as $z_{stu}^{s}$ with
$z_{stu}^{t}$, and $z_{tea}^{s}$ with $z_{tea}^{t}$, as well as cross-encoder
pairs like $z_{stu}^{s}$ with $z_{tea}^{t}$, and $z_{tea}^{s}$ with
$z_{stu}^{t}$. These delineations form the basis of our contrastive learning
process.

Moving forward, we apply a Fourier transformation to the volume features to
construct an amplitude map, which is crucial for identifying the salient high-
frequency components that highlight boundaries in Eq.3. A specialized square
mask is then utilized to isolate this high-frequency information in Eq.3.
Here, the Fourier transform and its inverse are denoted by $\mathcal{F}$ and
$\mathcal{F}^{-1}$ respectively. The mask $\mathbb{1}_{mask}$, with value zero
in its center and one at the periphery, has the same shape as $z$.

For boundary features, positive pairs are formed by analogous instances across
the student and teacher networks, such as $b_{stu}^{s}$ with $b_{tea}^{s}$,
and $b_{stu}^{t}$ with $b_{tea}^{s\rightarrow t}$. Conversely, negative pairs
are created by combining features from distinct instances, which may be within
the same network or across both, exemplified by pairs such as $b_{stu}^{s}$
with $b_{stu}^{t}$, and $b_{tea}^{s}$ with $b_{tea}^{t}$, as well as inter-
network pairs like $b_{stu}^{s}$ with $b_{tea}^{t}$, and $b_{tea}^{s}$ with
$b_{stu}^{t}$.

$b=\mathcal{F}^{-1}\left(\mathcal{F}\left(z\right)\cdot\mathbb{1}_{mask}\right)$
(3) $h(u,v)=\frac{u^{T}v}{\|u\|_{2}\|v\|_{2}}$ (4)
$\mathcal{L}_{c}=-\log\frac{\sum_{i=1}^{N_{p}}e^{{\left(h\left(u_{i}^{+},v_{i}^{+}\right)\right)}}}{\sum_{i=1}^{N_{p}}e^{{\left(h\left(u_{i}^{+},v_{i}^{+}\right)\right)}}+\sum_{j=1}^{N_{n}}e^{{\left(h\left(u_{j}^{+},v_{j}^{-}\right)\right)}}}$
(5)

To quantify the similarity of these pairs, we compute the cosine similarity
for each within both the time and frequency domains in Eq.4 and Eq.5. The
similarity for positive pairs is expressed as $h(u_{i}^{+},v_{i}^{+})$ where
$i$ spans all positive pair indices, and the similarity for negative pairs is
articulated as $h(u_{j}^{+},v_{j}^{-})$ where $j$ represents the indices of
all negative pairs. Here, $N_{p}$ stands for the count of positive pairings,
and $N_{n}$ corresponds to the count of negative pairings.

The contrastive learning loss for these high-frequency boundary pairs is then
calculated using the same equation as for the volumetric pairs. By summing up
the volumetric contrastive learning loss with the boundary contrastive
learning loss, we derive a comprehensive boundary-aware contrastive learning
loss. This loss function is designed to finely tune our model to discriminate
between the nuanced features of aneurysms, enhancing its segmentation
performance.

###  2.4 Overall framework and training objective

The loss function consists of two parts. $\mathcal{L}_{DCE}$ fully supervises
the four outputs of the teacher and student networks. BACL includes volume
contrast $\mathcal{L}_{c}^{z}$ and boundary contrast $\mathcal{L}_{c}^{b}$.
The ratio of $\lambda_{1}$ to $\lambda_{2}$ is set at 0.25:0.5.

$\mathcal{L}=\lambda_{1}\cdot(\mathcal{L}_{stu}^{src}+\mathcal{L}_{stu}^{trg}+\mathcal{L}_{tea}^{src}+\mathcal{L}_{tea}^{trg})+\lambda_{2}\cdot(\mathcal{L}_{c}^{z}+\mathcal{L}_{c}^{b})$
(6)

##  3 Experiments

###  3.1 Experimental Setting

Dataset: We tested our method with 3DRA images from 223 patients from the
@neurIST dataset [13]. These images were collected from four distinct medical
institutions, each employing varied scanning equipment and imaging protocols.
Consequently, this dataset exhibits a broad various in both visual
characteristics and resolution. The data diversity can evaluate the robustness
and adaptability of our proposed GS-EMA method.

Implementation details: Our study was conducted on a NVIDIA RTX 3090 GPU. We
utilized the Swin-UNet [14] architecture for both the student and teacher
networks in our framework. The training was set to 100 epochs. To determine
whether to apply EMA updates, we experimented with setting the EMA coefficient
$\alpha$ to either 0.9999 or 0.9. We started with an initial learning rate of
0.001 and adjusted it downwards by multiplying by 0.1 after every ten epochs.
The code will be publicly available soon.

###  3.2 Quantitative Results

| DSC (%) $\uparrow$ | Sen (%) $\uparrow$ | Jac (%) $\uparrow$ | VS (%) $\uparrow$
---|---|---|---|---
nnUNet [15] | 59.61 | 57.51 | 47.38 | 70.91
VASeg [16] | 60.28 | 54.47 | 49.82 | 67.91
FedDG [4] | 64.50 | 64.31 | 54.26 | 74.73
CMDG [17] | 65.01 | 64.10 | 54.11 | 73.38
Ours | 71.89 | 70.88 | 62.36 | 80.00
no EMA | 61.52 | 55.64 | 50.91 | 69.01
EMA | 64.71 | 62.86 | 54.03 | 72.64
GS-EMA | 68.49 | 72.79 | 58.40 | 76.63
BACL-V | 68.49 | 72.79 | 58.40 | 76.63
BACL-B | 70.62 | 75.14 | 60.54 | 78.22
BACL | 71.89 | 70.88 | 62.36 | 80.00

Table 1: Quantitative results including compare with SOTAs and ablation
studies. Critical metrics includes the Dice similarity coefficient (DSC),
Sensitivity (Sen), Jaccard index (Jac) and Volume similarity (VS).

Table.1 includes comparison with state-of-the-art (SOTA) methods and two
ablation studies. Our model outperforms traditional segmentation approaches
like nnUNet [15] and aneurysm-focused VASeg [16], as well as domain-
generalizing methods for medical image segmentation, including CMDG [17] and
FedDG [4]. The ablation study highlights that our GS-EMA algorithm, regulating
EMA updates with gradient relation, surpasses both regular and non-EMA methods
in segmenting aneurysms. It also indicates superior results for BACL when
integrating volume (BACL-V) and boundary (BACL-B) learning, compared to using
either alone.

###  3.3 Visual Inspection

Fig. 3 offers a visual comparison of aneurysm segmentation between our method
and the SOTAs. It is evident from the comparison that our approach is less
prone to over-segmentation while also being able to segment aneurysms more
completely.

Fig. 3: Comparative visualization of SOTAs and ours method on aneurysm
segmentation.

Fig. 4 displays a t-SNE comparison of latent features using EMA and GS-EMA.
The larger overlap achieved by GS-EMA indicates a stronger capability of the
model to extract domain-invariant features.

Fig. 4: The t-SNE visualization of latent features from EMA (left) and GS-EMA
(right).

##  4 Conclusion

In summary, our study introduces an effective GS-EMA algorithm and a boundary-
aware contrastive learning technique for aneurysm segmentation. These methods
outperform existing approaches by minimizing over-segmentations and capturing
more complete aneurysm structures. For future work, we plan to apply our GS-
EMA technique to a wider array of medical imaging datasets for further
validation and enhancement.

## References

  * [1] Qiongyao Liu, Ali Sarrami-Foroushani, Yongxing Wang, Michael MacRaild, Christopher Kelly, Fengming Lin, Yan Xia, Shuang Song, Nishant Ravikumar, Tufail Patankar, Z Taylor, Toni Lassila, and Alejandro F. Frang,  “Hemodynamics of thrombus formation in intracranial aneurysms: an in-silico observational study,”  APL Bioengineering, 2023.
  * [2] Francesco Pappalardo, Giulia Russo, Flora Musuamba Tshinanu, and Marco Viceconti,  “In silico clinical trials: concepts and early adoptions,”  Briefings in bioinformatics, vol. 20, no. 5, pp. 1699–1708, 2019.
  * [3] Wang Lu, Jindong Wang, Haoliang Li, Yiqiang Chen, and Xing Xie,  “Domain-invariant feature exploration for domain generalization,”  Transactions on Machine Learning Research, 2022.
  * [4] Quande Liu, Cheng Chen, Jing Qin, Qi Dou, and Pheng-Ann Heng,  “Feddg: Federated domain generalization on medical image segmentation via episodic learning in continuous frequency space,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 1013–1023.
  * [5] Shujun Wang, Lequan Yu, Kang Li, Xin Yang, Chi-Wing Fu, and Pheng-Ann Heng,  “Dofe: Domain-oriented feature embedding for generalizable fundus image segmentation on unseen datasets,”  IEEE Transactions on Medical Imaging, vol. 39, no. 12, pp. 4237–4248, 2020.
  * [6] Rayan Krishnan, Pranav Rajpurkar, and Eric J Topol,  “Self-supervised learning in medicine and healthcare,”  Nature Biomedical Engineering, vol. 6, no. 12, pp. 1346–1352, 2022.
  * [7] Hanlin Zhang, Yi-Fan Zhang, Weiyang Liu, Adrian Weller, Bernhard Schölkopf, and Eric P Xing,  “Towards principled disentanglement for domain generalization,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2022, pp. 8024–8034.
  * [8] Fengming Lin, Yan Xia, Nishant Ravikumar, Qiongyao Liu, Michael MacRaild, and Alejandro F Frangi,  “Adaptive semi-supervised segmentation of brain vessels with ambiguous labels,”  arXiv preprint arXiv:2308.03613, 2023.
  * [9] Tianhe Yu, Saurabh Kumar, Abhishek Gupta, Sergey Levine, Karol Hausman, and Chelsea Finn,  “Gradient surgery for multi-task learning,”  Advances in Neural Information Processing Systems, vol. 33, pp. 5824–5836, 2020.
  * [10] Lucas Mansilla, Rodrigo Echeveste, Diego H Milone, and Enzo Ferrante,  “Domain generalization via gradient surgery,”  in Proceedings of the IEEE/CVF international conference on computer vision, 2021, pp. 6630–6638.
  * [11] Haoran Dou, Ning Bi, Luyi Han, Yuhao Huang, Ritse Mann, Xin Yang, Dong Ni, Nishant Ravikumar, Alejandro F Frangi, and Yunzhi Huang,  “Gsmorph: Gradient surgery for cine-mri cardiac deformable registration,”  in MICCAI 2023, Cham, 2023, pp. 613–622, Springer Nature Switzerland.
  * [12] Chen Yang, Xiaoqing Guo, Zhen Chen, and Yixuan Yuan,  “Source free domain adaptation for medical image segmentation with fourier style mining,”  Medical Image Analysis, vol. 79, pp. 102457, 2022.
  * [13] Siegfried Benkner, Antonio Arbona, Guntram Berti, Alessandro Chiarini, Robert Dunlop, Gerhard Engelbrecht, Alejandro F Frangi, Christoph M Friedrich, Susanne Hanser, Peer Hasselmeyer, et al.,  “@ neurist: infrastructure for advanced disease management through integration of heterogeneous data, computing, and complex processing services,”  2010, vol. 14, pp. 1365–1377, IEEE.
  * [14] Hu Cao, Yueyue Wang, Joy Chen, Dongsheng Jiang, Xiaopeng Zhang, Qi Tian, and Manning Wang,  “Swin-unet: Unet-like pure transformer for medical image segmentation,”  in European conference on computer vision. Springer, 2022, pp. 205–218.
  * [15] Fabian Isensee, Paul F Jaeger, Simon AA Kohl, Jens Petersen, and Klaus H Maier-Hein,  “nnu-net: a self-configuring method for deep learning-based biomedical image segmentation,”  Nature methods, vol. 18, no. 2, pp. 203–211, 2021.
  * [16] Fengming Lin, Yan Xia, Shuang Song, Nishant Ravikumar, and Alejandro F Frangi,  “High-throughput 3dra segmentation of brain vasculature and aneurysms using deep learning,”  Computer Methods and Programs in Biomedicine, vol. 230, pp. 107355, 2023.
  * [17] Cheng Ouyang, Chen Chen, Surui Li, Zeju Li, Chen Qin, Wenjia Bai, and Daniel Rueckert,  “Causality-inspired single-source domain generalization for medical image segmentation,”  IEEE Transactions on Medical Imaging, vol. 42, no. 4, pp. 1095–1106, 2022.


&\hspace{15ex}+ \sum^n_{i=1} \int_{\Rn} \left( \sum^n_{j=1} \frac{1}{2} \frac{\partial a^{\alpha}_{ij}}{\partial y_j}(t,y) - b_i(t,y) \right) \frac{\partial (p^+ - p^-)}{\partial y_i}(t,y) p(t,y)^- \; \mathrm{d}y \\
&\le -\frac{\alpha}{2} \| \nabla (p^-)(t,\cdot) \|^2_{H} + \bar C(a,b) \sum^n_{i=1} \int_{\Rn} \left| \frac{\partial p^-}{\partial y_i}(t,y) p(t,y)^- \right| \; \mathrm{d}y \\
&\le -\frac{\alpha}{2} \| \nabla (p^-)(t,\cdot) \|^2_{H} + n M \bar C(a,b)^2 \frac{\| p(t,\cdot)^- \|^2_{H}}{2} + \frac{\| \nabla (p^-)(t,\cdot) \|^2_{H}}{2 M} \\
&\le n M \bar C(a,b)^2 \frac{\| p(t,\cdot)^- \|^2_{H}}{2} ,
\end{align*}

and thanks to a routine application of Gronwäll’s inequality, the fact that $\bar p(x) \ge 0$ almost everywhere yields $p(t,x) \ge 0$ for every $t \in [0,T]$ and almost every $x \in \Rn$.

At this step, for every $k \in \mathbb{N}$ choose a cut-off function $\varphi_k \in C^{\infty}_c(\Rn,[0,1])$ such that $\varphi_k(x) = 1$ for every $x \in \overline{B^{\Rn}_k(0)}$, $\textnormal{supp}(\varphi_k) \subseteq B^{\Rn}_{2k}(0)$, and whose first and second derivatives are uniformly bounded. For every $t \in [0,T]$, from the definition of FPE$_f$ and the monotone and dominated convergence theorems (here, we leverage both Theorem <ref> and the fact that $p(t,\cdot) \in L^2(\Rn,\R)$, $t \in [0,T]$) we may compute

\begin{align*}
&\int_{\Rn} p(t,y) \mathrm{d}y = \underset{k\to\infty}{\lim} \ \int_{\Rn} p(t,y) \varphi_k(y) \mathrm{d}y \\
&= \underset{k\to\infty}{\lim} \ \bigg( \int_{\Rn} \bar p(y) \varphi_k(y) \; \mathrm{d}y + \sum^n_{i,j=1} \int_{\Rn} \frac{1}{2} a^{\alpha}_{ij}(t,y) \frac{\partial^2 \varphi_k}{\partial y_i \partial y_j}(y) p(t,y) \; \mathrm{d}y \\
&\hspace{45ex} + \sum^n_{i=1} \int_{\Rn} b_i(t,y) \frac{\partial \varphi_k}{\partial y_i}(y) p(t,y) \; \mathrm{d}y \bigg) \\
&= \int_{\Rn} \bar p(y) \; \mathrm{d}y = 1 .
\end{align*}

To conclude, it is clear that we just need to prove that the curve $\mu : [0,T] \to \ProbaSpace(\Rn)$ defined in the statement of the theorem is narrowly continuous. For this, let $t \in [0,T]$ and $(t_k)_{k \in \mathbb{N}} \subseteq [0,T]$ such that $t_k \to t$ for $k \to \infty$. Note that $p \in C(0,T;H)$ implies that $\| p(t_n,\cdot) - p(t,\cdot) \|_{L^2} \to 0$ for $k \to \infty$. In particular, for any function $\varphi \in C_c(\Rn,\R)$, as soon as $k \to \infty$ we infer that

\left| \int_{\Rn} \varphi(y) p(t_n,y) \; \mathrm{d}y - \int_{\Rn} \varphi(y) p(t,y) \; \mathrm{d}y \right| \le \| \varphi \|_{L^2} \| p(t_n,\cdot) - p(t,\cdot) \|_{L^2} \to 0 .

We conclude from the fact that the narrow and weak* topologies coincide in $\ProbaSpace(\Rn)$.

Given that the proof is substantially long, for the sake of clarity we divide it in several step. Below, we adopt the notation we introduced and used in the proof of Theorem <ref>.

1) A second parabolic estimate. In this section, we provide computations by rather considering SFPE$_f$ with $0\neq f \in L^2([0,T]\times\Rn,\R)$. Let $(v_i)_{i \in \mathbb{N}}$ be a countable basis of $H^{2(m + 1)}(\Rn,\R)$ (and in turn of $H^{\ell}(\Rn,\R)$, for $\ell=0,\dots,2m+1$), such that $(v_i)_{i \in \mathbb{N}}$ is orthonormal in $H^1(\Rn,\R)$. Therefore, there is $(\bar x_i)_{i \in \mathbb{N}} \subseteq \R$ such that

\begin{equation} \label{eq:convInitialCond}
\underset{j\to\infty}{\lim} \ \sum^j_{i=1} \bar x_i v_i = \bar p , \quad \textnormal{in} \ H^1(\Rn,\R) .
\end{equation}

Moreover, by revisiting the proof of <cit.>, one easily see that for every $k \in \mathbb{N}$ there exists $x^k \in AC([0,T],\R^k)$ such that the function

p_k(t,\cdot) \triangleq \sum^k_{i=1} x^k_i(t) v_i \in C(0,T;L^2(\Rn,\R)) \cap L^2(0,T;H^1(\Rn,\R))

is the unique solution to the variational problem

\begin{equation} \label{eq:VarProblemDis}
\begin{cases}
\displaystyle \frac{\mathrm{d}}{\mathrm{d}t} ( p_k(t,\cdot) , v_i )_{L^2} + \Bil(t;p_k(t,\cdot),v_i) = (f(t,\cdot),v_i)_H , \quad i=1,\dots,k , \\[10pt]
x^k_i(0)= \bar x_i , \quad i=1,\dots,k ,
\end{cases}
\end{equation}

which in addition satisfies (see also <cit.>)

\begin{equation} \label{eq:CondConv}
\begin{cases}
\displaystyle \frac{\partial p_k}{\partial t}(t,\cdot) = \sum^k_{i=1} \frac{\mathrm{d}}{\mathrm{d}t} x^k_i(t) v_i \in L^2(0,T;H^1(\Rn,\R)) , \\
\displaystyle \frac{\partial p_k}{\partial y_j}(t,\cdot) = \sum^k_{i=1} x^k_i(t) \frac{\partial v_i}{\partial y_j} \in AC(0,T;L^2(\Rn,\R)) , \quad j=1,\dots,n , \\
p_k \to p \quad \textnormal{strongly in $L^2(0,T;H^1(\Rn,\R))$} .
\end{cases}
\end{equation}

Moreover, from the proof of Theorem <ref>, it is clear that each $p_k$ satisfies (<ref>).

At this step, by multiplying each equation in (<ref>) by $\displaystyle \frac{\mathrm{d}}{\mathrm{d}t} x^k_i(t)$ and summing those over $i=1,\dots,k$, thanks to (<ref>) for almost every $t \in [0,T]$ we may compute

\begin{align} \label{eq:DesiredBound}
\bigg\| \frac{\partial p_k}{\partial t}&(t,\cdot) \bigg\|^2_{L^2} = \sum^k_{i=1} \left( \frac{\mathrm{d}}{\mathrm{d}t} x^k_i(t) \right) \frac{\mathrm{d}}{\mathrm{d}s}\bigg|_{s=t} ( p_k(s,\cdot) , v_i)_{L^2} \nonumber \\
&= \sum^k_{i=1} \left( \frac{\mathrm{d}}{\mathrm{d}t} x^k_i(t) \right) \Big( (f(t,\cdot),v_i)_{L^2} - \Bil(t;p_k(t,\cdot),v_i) \Big) \nonumber \\
&\le \| f(t,\cdot) \|_{L^2} \left\| \frac{\partial p_k}{\partial t}(t,\cdot) \right\|_{L^2} - \sum^n_{i,j=1} \int_{\Rn} \frac{1}{2} a^{\alpha}_{ij}(t,y) \frac{\partial p_k}{\partial y_j}(t,y) \frac{\mathrm{d}}{\mathrm{d}t}\left( \frac{\partial p_k}{\partial y_i} \right)(t,y) \; \mathrm{d}y \\
&\quad - \sum^n_{i=1} \int_{\Rn} \left( \sum^n_{j=1} \frac{1}{2} \frac{\partial a^{\alpha}_{ij}}{\partial y_j}(t,y) - b_i(t,y) \right) \frac{\partial}{\partial y_i}\left( \frac{\partial p_k}{\partial t} \right)(t,y) \ p_k(t,y) \; \mathrm{d}y . \nonumber
\end{align}

We are going to bound the last two terms in (<ref>). Specifically, for the first to the last term, thanks to Theorem <ref> and (<ref>) for every $t \in [0,T]$ we obtain that

\begin{align*}
&-\int^t_0 \int_{\Rn} \sum^n_{i,j=1} \frac{1}{2} a^{\alpha}_{ij}(s,y) \frac{\partial p_k}{\partial y_j}(s,y) \frac{\mathrm{d}}{\mathrm{d}s}\left( \frac{\partial p_k}{\partial y_i} \right)(s,y) \; \mathrm{d}y \; \mathrm{d}s = \\
&= -\int^t_0 \int_{\Rn} \sum^n_{i,j=1} \frac{1}{2} \bigg( \frac{\mathrm{d}}{\mathrm{d}s}\left( a^{\alpha}_{ij} \frac{\partial p_k}{\partial y_j} \frac{\partial p_k}{\partial y_i} \right)(s,y) - \dot{a}^{\alpha}_{ij}(s,y) \frac{\partial p_k}{\partial y_j}(s,y) \frac{\partial p_k}{\partial y_i}(s,y) \\
&\hspace{37.5ex} - a^{\alpha}_{ij}(s,y) \frac{\mathrm{d}}{\mathrm{d}s}\left( \frac{\partial p_k}{\partial y_j} \right)(s,y) \ \frac{\partial p_k}{\partial y_i}(s,y) \bigg) \; \mathrm{d}y \; \mathrm{d}s ,
\end{align*}

and since $a^{\alpha}_{ji} = a^{\alpha}_{ij}$ for every $i,j=1,\dots,n$, we may compute

\begin{align*}
&-\int^t_0 \int_{\Rn} \sum^n_{i,j=1} a^{\alpha}_{ij}(s,y) \frac{\partial p_k}{\partial y_j}(s,y) \frac{\mathrm{d}}{\mathrm{d}s}\left( \frac{\partial p_k}{\partial y_i} \right)(s,y) \; \mathrm{d}y \; \mathrm{d}s \le \\
&\le -\sum^n_{i,j=1} \frac{1}{2} \int_{\Rn} \left( a^{\alpha}_{ij}(t,y) \frac{\partial p_k}{\partial y_j}(t,y) \frac{\partial p_k}{\partial y_i}(t,y) - a^{\alpha}_{ij}(0,y) \frac{\partial p_k}{\partial y_j}(0,y) \frac{\partial p_k}{\partial y_i}(0,y) \right) \; \mathrm{d}y \; \mathrm{d}s \\
&\hspace{30ex} + \frac{\| \dot{a} \|_{L^{\infty}}}{2} \int^t_0 \int_{\Rn} \sum^n_{i,j=1} \left| \frac{\partial p_k}{\partial y_j}(s,y) \right| \left| \frac{\partial p_k}{\partial y_i}(s,y) \right| \; \mathrm{d}y \; \mathrm{d}s .
\end{align*}

Therefore, the convergence (<ref>) and (<ref>) finally yield

\begin{align} \label{eq:IntermediateIneq}
-\int^t_0 \int_{\Rn} \sum^n_{i,j=1} &a^{\alpha}_{ij}(s,y) \frac{\partial p_k}{\partial y_j}(s,y) \frac{\mathrm{d}}{\mathrm{d}s}\left( \frac{\partial p_k}{\partial y_i}(s,y) \right) \; \mathrm{d}y \; \mathrm{d}s \le \nonumber \\
&\le -\frac{\alpha}{2} \| \nabla p_k(t,\cdot) \|^2_{L^2} + C\big( \alpha , \|(a,b)\|_{\Cbasis} \big ) \left( \| \bar p \|^2_{H^1} + \int^t_0 \| f(s,\cdot) \|^2_{L^2} \; \mathrm{d}s \right) ,
\end{align}

where $C\big( \alpha , \|(a,b)\|_{\Cbasis} \big ) > 0$ is a constant which continuously depends on $\alpha$ and $\|(a,b)\|_{\Cbasis}$ uniquely. Below, we will implicitly overload the constant $C\big( \alpha , \|(a,b)\|_{\Cbasis} \big )$.

We now focus on the last term in (<ref>). For this, first thanks to Theorem <ref> one obtains that, for every $i=1,\dots,n$, $\varphi \in C^{\infty}_c(\Rn,\R)$, and almost every $t \in [0,T]$,

\begin{align*}
-\int_{\Rn} \bigg( \sum^n_{j=1} \frac{1}{2} \frac{\partial a^{\alpha}_{ij}}{\partial y_j}(t,y) - b_i(t,y) &\bigg) \frac{\partial \varphi}{\partial y_i}(y) p_k(t,y) \; \mathrm{d}y = \\
&= \int_{\Rn} \frac{\partial}{\partial y_i} \left( \left( \sum^n_{j=1} \frac{1}{2} \frac{\partial a^{\alpha}_{ij}}{\partial y_j} - b_i \right) p_k \right)(t,y) \ \varphi(y) \; \mathrm{d}y .
\end{align*}

Since from (<ref>) we infer the existence of a sequence $(\varphi_{\ell})_{\ell \in \mathbb{N}} \subseteq C^{\infty}_c(\Rn,\R)$ such that $\displaystyle \left\| \frac{\partial p_k}{\partial t}(t,\cdot) - \varphi_{\ell} \right\|_{H^1} \to 0$, for $\ell\to\infty$, again Theorem <ref> and (<ref>) yield

\begin{align*}
-\int^t_0 &\sum^n_{i=1} \int_{\Rn} \left( \sum^n_{j=1} \frac{1}{2} \frac{\partial a^{\alpha}_{ij}}{\partial y_j}(s,y) - b_i(s,y) \right) \frac{\partial}{\partial y_i}\left( \frac{\partial p_k}{\partial s} \right)(s,y) \ p_k(s,y) \; \mathrm{d}y \; \mathrm{d}s = \\
&= \int^t_0 \sum^n_{i=1} \int_{\Rn} \left( \sum^n_{j=1} \frac{1}{2} \frac{\partial^2 a^{\alpha}_{ij}}{\partial y_i \partial y_j}(s,y) - \frac{\partial b_i}{\partial y_i}(s,y) \right) \frac{\partial p_k}{\partial s}(s,y) p_k(s,y) \; \mathrm{d}y \; \mathrm{d}s \\
&\quad+ \int^t_0 \sum^n_{i=1} \int_{\Rn} \left( \sum^n_{j=1} \frac{1}{2} \frac{\partial a^{\alpha}_{ij}}{\partial y_j}(s,y) - b_i(s,y) \right) \frac{\partial p_k}{\partial y_i}(s,y) \frac{\partial p_k}{\partial s}(s,y) \; \mathrm{d}y \; \mathrm{d}s \\
&\le C\big( \alpha , \|(a,b)\|_{\Cbasis} \big )  \int^t_0 \int_{\Rn} \left( |p_k(s,y)| + \sum^n_{i=1} \left| \frac{\partial p_k}{\partial y_i}(s,y) \right| \right) \left| \frac{\partial p_k}{\partial s}(s,y) \right| \; \mathrm{d}y \; \mathrm{d}s \\
&\le \frac{1}{2} \int^t_0 \left\| \frac{\partial p_k}{\partial s}(s,\cdot) \right\|^2_{L^2} \; \mathrm{d}s + C\big( \alpha , \|(a,b)\|_{\Cbasis} \big )  \left( \| \bar p \|^2_{H^1} + \int^t_0 \| f(s,\cdot) \|^2_{L^2} \; \mathrm{d}s \right) ,
\end{align*}

for every $t \in [0,T]$. By summing up this latter inequality with (<ref>) and (<ref>), via a routine Granwäll's inequality argument we infer that, for $k \in \mathbb{N}$ and $t \in [0,T]$,

\begin{align} \label{eq:estimate2}
\int^t_0 \left\| \frac{\partial p_k}{\partial s}(s,\cdot) \right\|^2_{L^2} \; \mathrm{d}s \le C\big( \alpha , \|(a,b)\|_{\C} \big ) \left( \| \bar p \|^2_{H^1} + \int^t_0 \| f(s,\cdot) \|^2_{L^2} \; \mathrm{d}s \right) .
\end{align}

2) First-order-in-time and second-order-in-space regularity. In this section, we provide computations by rather considering SFPE$_f$ with $0\neq f \in L^2([0,T]\times\Rn,\R)$. We first show the inclusion $\displaystyle \frac{\partial p}{\partial t} \in L^2(0,T;L^2(\Rn,\R)) \cong L^2([0,T]\times\Rn,\R)$, this latter identification being true since $L^2(\Rn,\R)$ is a separable Hilbert space. For this, we show the existence of $\xi \in L^2(0,T;L^2(\Rn,\R))$ such that the following holds in $L^2(\Rn,\R)$:

\begin{equation} \label{eq:WeakTimeDer}
\int^T_0 \xi(t,\cdot) \psi(t) \; \mathrm{d}t = -\int^T_0 p(t,\cdot) \frac{\mathrm{d}\psi}{\mathrm{d}t}(t) \; \mathrm{d}t , \quad \textnormal{for all} \ \psi \in C^{\infty}_c([0,T],\R) .
\end{equation}

For this, we first note that, thanks to (<ref>) there exists $\xi \in L^2(0,T;L^2(\Rn,\R))$ such that $\displaystyle\left( \frac{\partial p_k}{\partial t} \right)_{k \in \mathbb{N}} \subseteq L^2(0,T;L^2(\Rn,\R))$ weakly converges to $\xi$, up to a subsequence. Now, for any $\psi \in C^{\infty}_c([0,T],\R)$ and any $\varphi \in L^2(\Rn,\R)$, thanks to the last (convergence) property in (<ref>) we may compute

\begin{align*}
\bigg( \int^T_0 &\xi(t,\cdot) \psi(t) \; \mathrm{d}t , \varphi \bigg)_{L^2} = \int^T_0 \big( \xi(t,\cdot) , \psi(t) \varphi \big)_{L^2} \; \mathrm{d}t = \\
&= \underset{k\to\infty}{\lim} \ \int^T_0 \left( \frac{\partial p_k}{\partial t}(t,\cdot) , \psi(t) \varphi \right)_{L^2} \; \mathrm{d}t = \underset{k\to\infty}{\lim} \ \left( \int^T_0 \frac{\partial p_k}{\partial t}(t,\cdot) \psi(t) \; \mathrm{d}t , \varphi \right)_{L^2} \\
&= -\underset{k\to\infty}{\lim} \ \left( \int^T_0 p_k(t,\cdot) \frac{\mathrm{d}\psi}{\mathrm{d}t}(t) \; \mathrm{d}t , \varphi \right)_{L^2} = -\underset{k\to\infty}{\lim} \ \int^T_0 \left( p_k(t,\cdot) , \frac{\mathrm{d}\psi}{\mathrm{d}t}(t) \varphi \right)_{L^2} \; \mathrm{d}t \\
&= -\int^T_0 \left( p(t,\cdot) , \frac{\mathrm{d}\psi}{\mathrm{d}t}(t) \varphi \right)_{L^2} \; \mathrm{d}t = \left( -\int^T_0 p(t,\cdot) \frac{\mathrm{d}\psi}{\mathrm{d}t}(t) \; \mathrm{d}t , \varphi \right)_{L^2} ,
\end{align*}

and the equivalence (<ref>) readily follows.

Next, we prove that $p \in L^2(0,T;H^2(\Rn,\R))$. For this, we first note that, from Theorem <ref> and the variational problem (<ref>), thanks to our previous computations the following variational equality holds true for almost every $t \in [0,T]$:

\begin{align}
\label{eq:Elliptic1}
&\int_{\Rn} \sum^n_{i,j=1} \frac{1}{2} a^{\alpha}_{ij}(t,y) \frac{\partial p}{\partial y_j}(t,y) \frac{\partial \varphi}{\partial y_i}(y) \; \mathrm{d}y \hspace{20ex} \textnormal{for} \quad \varphi \in H^1(\Rn,\R) \\
&+ \sum^n_{i=1} \int_{\Rn} \left( \sum^n_{j=1} \frac{1}{2} \frac{\partial a^{\alpha}_{ij}}{\partial y_j}(t,y) - b_i(t,y) \right) \frac{\partial \varphi}{\partial y_i}(y) p(t,y) \; \mathrm{d}y = \int_{\Rn} g(t,y) \varphi(y) \; \mathrm{d}y , \nonumber
\end{align}

\begin{equation} \label{eq:Elliptic2}
g \triangleq f - \frac{\partial p}{\partial t} \in L^2([0,T]\times\Rn,\R) .
\end{equation}

Therefore, Theorem <ref> and the classical elliptic regularity theory (see, e.g., <cit.>) imply, for almost every $t \in [0,T]$, both that $p(t,\cdot) \in H^2(\Rn,\R)$ and that

\begin{equation} \label{eq:EllipticBound}
\| p(t,\cdot) \|^2_{H^2} \le C\big( \alpha , \|(a,b)\|_{\Cbasis} \big ) \Big( \| p(t,\cdot) \|^2_{L^2} + \| g(t,\cdot) \|^2_{L^2} \Big) .
\end{equation}

At this step, from our choice for the countable basis $(v_i)_{i \in \mathbb{N}} \subseteq H^{2(m + 1)}(\Rn,\R)$, for almost every $t \in [0,T]$ and every $i \in \mathbb{N}$, in particular there exists $z_i(t) \in \R$ such that

\begin{equation} \label{eq:convInH2}
p(t,\cdot) = \underset{k\to\infty}{\lim} \ \sum^k_{i=1} z_i(t) v_i(\cdot) , \quad \textnormal{in \ $H^2(\Rn,\R)$} ,
\end{equation}

for almost every $t \in [0,T]$. But, up to extracting a subsequence, (<ref>) yields

p(t,\cdot) = \underset{k\to\infty}{\lim} \ \sum^k_{i=1} x^k_i(t) v_i(\cdot) , \quad \textnormal{in \ $H^1(\Rn,\R)$} ,

for almost every $t \in [0,T]$, and from the orthonormality of $(v_i)_{i \in \mathbb{N}}$ in $H^1(\Rn,\R)$ we obtain that $z_i(t) = x^k_i(t)$, for every $k \in \mathbb{N}$, $i=1,\dots,k$ and almost every $t \in [0,T]$. Hence, (<ref>) yields that the mapping $p : [0,T] \to H^2(\Rn,\R)$ is strongly Bochner measurable. In addition, combining the elliptic estimate (<ref>) together with the estimates (<ref>) and (<ref>) finally provides that $p \in L^2(0,T;H^2(\Rn,\R))$ with

\begin{align} \label{eq:estimate3}
\int^t_0 \| p(s,\cdot) \|^2_{H^2} \; \mathrm{d}s &\le \nonumber \\
&\le C\big( \alpha , \|(a,b)\|_{\Cbasis} \big ) \left( \| \bar p \|^2_{H^1} + \int^t_0 \| f(s,\cdot) \|^2_{L^2} \; \mathrm{d}s \right) , \quad t \in [0,T] .
\end{align}

Before moving on with additional regularity properties, we note that, thanks to Theorem <ref> and the fact that $p \in L^2(0,T;H^2(\Rn,\R))$, integrating by parts the variational problem (<ref>) with $f = 0$ readily yields

\int^T_0 \int_{\Rn} \left( \frac{\partial p}{\partial t}(t,y) - (\Kolmo_t)^* p(t,y) \right) \varphi(y) \; \mathrm{d}y \; \mathrm{d}t = 0 , \quad \varphi \in C^{\infty}_c(\Rn,\R) ,

and SFPE$_0$ follows from the previous regularity properties and a density argument.

3) Second-order-in-time and fourth-order-in-space regularity. We now turn to the original setting SFPE$_0$, i.e., SFPE$_f$ with $f = 0$. The first step to further improve the regularity of $p$ consists of formally differentiating (<ref>) with respect to time and studying solutions $q : [0,T]\times\Rn\to\R$ to the new variational problem

\begin{equation} \label{eq:VarProblemReg}
\begin{cases}
\displaystyle \frac{\mathrm{d}}{\mathrm{d}t} ( q(t,\cdot) , \varphi )_{L^2} + \Bil(t;q(t,\cdot),\varphi) + \\[5pt]
\hspace{25ex} + \BilDot(t;p(t,\cdot),\varphi) = 0 , \quad \varphi \in H^1(\Rn,\R) , \\
q(0,\cdot) = (\Kolmo_0)^* \bar p(\cdot) .
\end{cases}
\end{equation}

Problem (<ref>) is well-posed. Indeed, since $p \in L^2(0,T;H^2(\Rn,\R))$, thanks to Theorem <ref> one easily shows that

\BilDot(t;p(t,\cdot),\varphi) = -\int_{\Rn} (\KolmoDot_t)^* p(t,y) \varphi(y) \; \mathrm{d}y ,

for every $\varphi \in H^1(\Rn,\R)$ and almost every $t \in [0,T]$. Therefore, since $(\KolmoDot_t)^* p \in L^2([0,T]\times\Rn,\R)$ and $(\Kolmo_0)^* \bar p \in H^1(\Rn,\R)$, problem (<ref>) fits the setting of Theorem <ref>, thus from our previous computations there exists a unique solution $q \in L^2(0,T;H^2(\Rn,\R))$ to (<ref>), which additionally satisfies $\displaystyle \frac{\partial q}{\partial t} \in L^2([0,T]\times\Rn,\R)$.

At this step, we define the function

\begin{align*}
w : \ &[0,T]\times\Rn\to\R \\
&(t,y) \mapsto \bar p(y) + \left( \int^t_0 q(s,\cdot) \; \mathrm{d}s \right)(y) = \bar p(y) + \int^t_0 q(s,y) \; \mathrm{d}s ,
\end{align*}

where the Lebesgue integral in the last equality stems from the properties of the Bochner integral for functions in $L^2([0,T]\times\Rn,\R)$. Clearly, $w \in C^0(0,T;H^2(\Rn,\R))$, and Fubini theorem yields, for every $\varphi_1 \in L^2(\Rn,\R)$ and $\psi \in C^{\infty}_c([0,T],\R)$,

\int_{\Rn} \left( \int^T_0 \left(\int^t_0 q(s,y) \; \mathrm{d}s \right) \frac{\mathrm{d} \psi}{\mathrm{d}t}(t) \; \mathrm{d}t + \int^T_0 q(t,y) \psi(t) \; \mathrm{d}t \right) \varphi_1(y) \; \mathrm{d}y = 0 ,

and, for almost every $t \in [0,T]$ and every $\varphi_2 \in C^{\infty}_c(\Rn,\R)$,

\int_{\Rn} w(t,y) \frac{\partial \varphi_2}{\partial y_j}(y) \; \mathrm{d}y = -\int_{\Rn} \left( \frac{\partial \bar p}{\partial y_j}(y) + \int^t_0 \frac{\partial q}{\partial y_j}(s,y) \; \mathrm{d}s \right) \varphi_2(y) \; \mathrm{d}y , \quad j=1,\dots,n .

In particular, from such properties we readily deduce that, for every $j=1,\dots,n$,

\begin{equation} \label{eq:CondDer}
\frac{\partial w}{\partial t} = q \in L^2(0,T;H^2(\Rn,\R)) , \quad \textnormal{with} \quad \frac{\partial}{\partial y_j}\left( \frac{\partial w}{\partial t} \right) = \frac{\partial}{\partial t}\left( \frac{\partial w}{\partial y_j} \right) .
\end{equation}

Thanks to Theorem <ref> and the properties listed in (<ref>),
for every $\varphi \in C^{\infty}_c(\Rn,\R)$ and almost every $t \in [0,T]$, we may compute

\begin{align*}
&\frac{\mathrm{d}}{\mathrm{d}t} ( w(t,\cdot) , \varphi )_{L^2} = ( q(t,\cdot) , \varphi )_{L^2} \\
&= \big( (\Kolmo_0)^* \bar p , \varphi \big)_{L^2} - \int^t_0 \Bil(s;q(s,\cdot),\varphi) \; \mathrm{d}s - \int^t_0 \BilDot(s;p(s,\cdot) \pm w(s,\cdot),\varphi) \; \mathrm{d}s \\
&= \int_{\Rn} (\Kolmo_0)^* \bar p(y) \varphi(y) \; \mathrm{d}y - \\
&\quad -\sum^n_{i,j=1} \int_{\Rn} \frac{1}{2} \left( \int^t_0 \left( a^{\alpha}_{ij}(s,y) \frac{\partial}{\partial s}\left( \frac{\partial w}{\partial y_j} \right)(s,y) + \dot{a}^{\alpha}_{ij}(s,y) \frac{\partial w}{\partial y_j}(s,y) \right) \mathrm{d}s \right) \frac{\partial \varphi}{\partial y_i}(y) \; \mathrm{d}y \\
&\quad - \sum^n_{j=1} \int_{\Rn} \Bigg( \int^t_0 \Bigg( \Bigg( \sum^n_{j=1} \frac{\partial \dot{a}^{\alpha}_{ij}}{\partial y_j}(s,y) - \dot{b}_i(s,y) \Bigg) w(s,y) \\
&\hspace{28.75ex}+ \Bigg( \sum^n_{j=1} \frac{\partial a^{\alpha}_{ij}}{\partial y_j}(s,y) - b_i(s,y) \Bigg) \frac{\partial w}{\partial s}(s,y) \Bigg) \mathrm{d}s \Bigg) \frac{\partial \varphi}{\partial y_i}(y) \; \mathrm{d}y \\
&\quad - \int^t_0 \BilDot(s;(p - w)(s,\cdot),\varphi) \; \mathrm{d}s .
\end{align*}

Therefore, from <cit.> and the fact that

\int_{\Rn} (\Kolmo_0)^* \bar p(y) \varphi(y) \; \mathrm{d}y = -\Bil(0;\bar p,\varphi) ,

we conclude that $w$ solves the following integro-differential variational problem

\begin{equation} \label{eq:VarProblemIntegro}
\begin{cases}
\displaystyle \frac{\mathrm{d}}{\mathrm{d}t} ( w(t,\cdot) , \varphi )_{L^2} + \Bil(t;w(t,\cdot),\varphi) + \\[5pt]
\displaystyle \hspace{10ex} + \int^t_0 \BilDot(s;(p-w)(s,\cdot),\varphi) \; \mathrm{d}s = 0 , \quad \varphi \in H^1(\Rn,\R) , \\[10pt]
w(0,\cdot) = \bar p(\cdot) .
\end{cases}
\end{equation}

Clearly, $p = w$ is solution to (<ref>). In particular, if so we would obtain that $\displaystyle \frac{\partial p}{\partial t} \in L^2(0,T;H^2(\Rn,\R))$ and $\displaystyle \frac{\partial^2 p}{\partial t^2} \in L^2([0,T]\times\Rn,\R)$, and thanks to the choice of the basis $(v_i)_{i \in \mathbb{N}} \subseteq H^{2(m + 1
)}(\Rn,\R)$ and to a higher order elliptic regularity argument applied to (<ref>)–(<ref>) with $f=0$, we would also obtain that $p \in L^2(0,T;H^4(\Rn,\R))$, that is the sought second-order-in-time and fourth-order-in-space regularity.

Hence, to conclude we need to prove that the integro-differential variational problem (<ref>) may have one solution at most. For this, let $z_1 , z_2 \in L^2(0,T;H^2(\Rn,\R))$, such that $\displaystyle \frac{\partial z_1}{\partial t} , \frac{\partial z_2}{\partial t} \in L^2([0,T]\times\Rn,\R)$, satisfy (<ref>). By defining the mapping

\begin{align*}
h : \ &[0,T] \to L^2(\Rn,\R) \\
&t \mapsto \int^t_0 (\KolmoDot_s)^* (z_1 - z_2)(s,\cdot) \; \mathrm{d}s ,
\end{align*}

it is readily checked that $h \in L^2([0,T]\times\Rn,\R)$, in particular satisfying

\begin{align} \label{eq:IntegroBound}
&\int^t_0 \| h(s,\cdot) \|^2_{L^2} \; \mathrm{d}s \le \nonumber \\
&\hspace{5ex}\le C\big( \alpha , \|(a,b)\|_{\Cbasis} \big ) \int^t_0 \left( \int^s_0 \| (z_1 - z_2)(r,\cdot) \|^2_{H^2} \; \mathrm{d}r \right) \; \mathrm{d}s , \quad t \in [0,T] ,
\end{align}

and that

\begin{equation*}
\begin{cases}
\displaystyle \frac{\mathrm{d}}{\mathrm{d}t} ( (z_1 - z_2)(t,\cdot) , \varphi )_{L^2} + \Bil(t;(z_1 - z_2)(t,\cdot),\varphi) = \\[5pt]
\hspace{45ex}= ( h(t,\cdot) , \varphi )_{L^2} , \quad \varphi \in H^1(\Rn,\R) , \\
(z_1 - z_2)(0,\cdot) = 0 .
\end{cases}
\end{equation*}

Therefore, we are in the setting of Theorem <ref> and of our previous computations. In particular, we may combine (<ref>) with (<ref>) to obtain that, for every $t \in [0,T]$,

\begin{align*}
\int^t_0 \| (z_1 - z_2)(s,\cdot) \|^2_{H^2} \; \mathrm{d}s &\le C\big( \alpha , \|(a,b)\|_{\Cbasis} \big ) \int^t_0 \| h(s,\cdot) \|^2_{L^2} \; \mathrm{d}s \\
&\le C\big( \alpha , \|(a,b)\|_{\Cbasis} \big ) \int^t_0 \left( \int^s_0 \| (z_1 - z_2)(r,\cdot) \|^2_{H^2} \; \mathrm{d}r \right) \; \mathrm{d}s ,
\end{align*}

and a routine application of Gronwäll's inequality allows us to conclude that $z_1 = z_2$.

4) Conclusion. Given Theorem <ref>, the regularity of the initial condition $\bar p \in H^{2 m + 1}(\Rn,\R)$, and the choice of the basis $(v_i)_{i \in \mathbb{N}} \subseteq H^{2(m + 1)}(\Rn,\R)$, one may easily iterate the previous computations by induction to define successive time derivatives of $p$ via iteratively considering formal differentiation-in-time of the variational problem (<ref>). In particular, it is clear how to extend properties (<ref>) and define updated integro-differential variational problems (<ref>) by induction, which, thanks to iteratively higher order elliptic regularity arguments applied to (<ref>)–(<ref>) and the elliptic estimate (<ref>), provide that $\displaystyle \frac{\partial^{\ell} p}{\partial t^{\ell}} \in L^2(0,T;H^{2(m+1-\ell)}(\Rn,\R))$, for $\ell=1,\dots,m+1$, and $p \in L^2(0,T;H^{2(m + 1)}(\Rn,\R))$. The conclusion follows.

[1]
K. Abraham, Nonparametric Bayesian Posterior Contraction Rates for
Scalar Diffusions with High-frequency Data, Bernoulli, 25 (2019),
pp. 2696–2728.

[2]
C. Aeckerle-Willems and C. Strauch, Sup-norm Adaptive Drift
Estimation for Multivariate Nonreversible Diffusions, The Annals of
Statistics, 50 (2022), pp. 3484–3509.

[3]
N. Aronszajn, Theory of reproducing kernels, Transactions of the
American mathematical society, 68 (1950), pp. 337–404.

[4]
J.-P. Aubin and H. Frankowska, Set-Valued Analysis, Birkhäuser,

[5]
T. Bhudisaksang and A. Cartea, Online Drift Estimation for
Jump-Diffusion Processes, Bernoulli, 27 (2021), pp. 2494–2518.

[6]
R. Bonalli and B. Bonnet-Weill, First-Order Pontryagin Maximum
Principle for Risk-Averse Stochastic Optimal Control Problems, SIAM Journal
on Control and Optimization, 61 (2023), pp. 1881–1909.

[7]
R. Bonalli, T. Lew, and M. Pavone, Sequential convex programming
for non-linear stochastic optimal control, ESAIM: Control, Optimisation and
Calculus of Variations, 28 (2022), p. 64.

[8]
T. Breiten, K. Kunisch, and L. Pfeiffer, Control Strategies for the
Fokker–Planck Equation, ESAIM: Control, Optimisation and Calculus of
Variations, 24 (2018).

[9]
H. Brezis, Functional Analysis, Sobolev Spaces and Partial
Differential Equations, Springer, 2011.

[10]
M. Chipot, Elements of Nonlinear Analysis, Springer, 2000.

[11]
F. Comte and V. Genon-Catalot, Drift Estimation on Non Compact
Support for Diffusion Models, Stochastic Processes and their Applications,
134 (2021), pp. 174–207.

[12]
D. Dacunha-Castelle and D. Florens-Zmirou, Estimation of the
Coefficients of a Diffusion from Discrete Observations, Stochastics: An
International Journal of Probability and Stochastic Processes, 19 (1986),
pp. 263–284.

[13]
L. C. Evans, Partial Differential Equations, American
Mathematical Society, 2010.

[14]
A. Figalli, Existence and uniqueness of martingale solutions for
SDEs with rough or degenerate coefficients, Journal of Functional Analysis,
254 (2008).

[15]
D. Florens-Zmirou, Approximate Discrete-Time Schemes for Statistics
of Diffusion Processes, Statistics, 20 (1989), pp. 547–557.

[16]
A. R. Gallant, The Relative Efficiency of Method of Moments
Estimators, Journal of Econometrics, 92 (1999), pp. 149–172.

[17]
A. R. Gallant, D. Hsieh, and G. E. Tauchen, Estimation of
Stochastic Volatility Models with Diagnostics, Journal of Econometrics, 81
(1997), pp. 159–192.

[18]
H. Garnier and L. Wang, Identification of Continuous-Time Models
from Sampled Data, Springer, 2008.

[19]
V. Genon-Catalot and C. Laredo, Asymptotic Equivalence of
Nonparametric Diffusion and Euler Scheme Experiments, The Annals of
Statistics, 42 (2014), pp. 1145–1165.

[20]
H. J. Kushner and G. G. Yin, Stochastic Approximation and Recursive
Algorithms and Applications, Springer, 2001.

[21]
H. Lavenant et al., Towards a Mathematical Theory of Trajectory
arXiv:2102.09204, 2021.

[22]
J.-F. Le Gall, Brownian Motion, Martingales, and Stochastic
Calculus, Springer, 2016.

[23]
L. Lennart, System Identification – Theory For the User, PTR
Prentice Hall, 1998.

[24]
X. Li et al., Scalable Gradients for Stochastic Differential
Equations, in Int. Conference on Artificial Intelligence and Statistics,
2020, pp. 3870–3882.

[25]
J.-L. Lions, Optimal Control of Systems Governed by Partial
Differential Equations, Springer, 1968.

[26]
J.-L. Lions and E. Magenes, Non-Homogeneous Boundary Value Problems
and Applications: Vol. 2, Springer, 2012.

[27]
R. S. Liptser and A. N. Shiriaev, Statistics of Random Processes:
General Theory, Springer, 1977.

[28]
N. Marie and A. Rosier, Nadaraya Watson Estimator for IID Paths of
Diffusion Processes, Scandinavian Journal of Statistics, 50 (2023),
pp. 589–637.

[29]
U. Marteau-Ferey, F. Bach, and A. Rudi, Non-parametric models for
non-negative functions, Advances in neural information processing systems,
33 (2020), pp. 12816–12826.

[30]
H. Masuda, Convergence of Gaussian Quasi-Likelihood Random Fields
for Ergodic Lévy driven SDE Observed at High Frequency, The Annals of
Statistics, 41 (2013), pp. 1593–1641.

[31]
B. Muzellec, F. Bach, and A. Rudi, Learning psd-valued functions
using kernel sums-of-squares, arXiv preprint arXiv:2111.11306,  (2021).

[32]
S. Nakakita, Parametric Estimation of Stochastic Differential
Equations via Online Gradient Descent.
arXiv:2210.08800, 2022.

[33]
O. Nelles, Nonlinear System Identification, Springer, 2001.

[34]
Y. Nesterov and A. Nemirovskii, Interior-point polynomial algorithms
in convex programming, SIAM, 1994.

[35]
R. Nickl and K. Ray, Nonparametric Statistical Inference for Drift
Vector Fields of Multidimensional Diffusions, The Annals of Statistics, 48
(2020), pp. 1383–1408.

[36]
R. Nickl and J. Söhl, Nonparametric Bayesian Posterior
Contraction Rates for Discretely Observed Scalar Diffusions, The Annals of
Statistics, 45 (2017), pp. 1664–1693.

[37]
J. N. Nielsen, H. Madsen, and P. C. Young, Parameter Estimation in
Stochastic Differential Equations: an Overview, Annual Reviews in Control,
24 (2000), pp. 83–94.

[38]
J. Ridderhof, K. Okamoto, and P. Tsiotras, Nonlinear Uncertainty
Control with Iterative Covariance Steering, in Int. Conference on Decision
and Control, 2019, pp. 3484–3490.

[39]
A. Rudi and C. Ciliberto, Psd representations for effective
probability models, in Advances in Neural Information Processing Systems,

[40]
A. Rudi, U. Marteau-Ferey, and F. Bach, Finding global minima via
kernel approximations, Mathematical Programming,  (2024), pp. 1–82.

[41]
L. Sharrock and N. Kantas, Joint Online Parameter Estimation and
Optimal Sensor Placement for the Partially Observed Stochastic
Advection-Diffusion Equation, SIAM Journal on Uncertainty Quantification,
10 (2022), pp. 55–95.

[42]
I. Steinwart and A. Christmann, Support vector machines, Springer
Science & Business Media, 2008.

[43]
D. W. Stroock and S. Varadhan, Multidimensional Diffusion
Processes, Springer Science & Business Media, 1997.

[44]
S. C. Surace and J.-P. Pfister, Online Maximum-Likelihood
Estimation of the Parameters of Partially Observed Diffusion Processes,
IEEE Transactions on Automatic Control, 64 (2018), pp. 2814–2829.

[45]
M. Uchida and N. Yoshida, Quasi Likelihood Analysis of Volatility
and Nondegeneracy of Statistical Random Field, Stochastic Processes and
their Applications, 123 (2013), pp. 2851–2876.

[46]
F. Van der Meulen and H. Van Zanten, Consistent Nonparametric
Bayesian Inference for Discretely Observed Scalar Diffusions, Bernulli, 19
(2013), pp. 44–63.

[47]
G. Velho, J. Auriol, and R. Bonalli, A Gradient Descent-Ascent
Method for Continuous-Time Risk-Averse Optimal Control.
arxiv:2306.12878, 2024.

[48]
H. Wendland, Scattered Data Approximation, Cambridge University
Press, 2004.

[49]
W. Xu et al., Infinitely Deep Bayesian Neural Networks with
Stochastic Differential Equations.
arXiv:2102.06559, 2022.

[50]
J. Yong and X. Y. Zhou, Stochastic Controls: Hamiltonian Systems
and HJB equations, Springer Science and Business Media, 1999.

[51]
N. Yoshida, Polynomial Type Large Deviation Inequalities and
Quasi-Likelihood Analysis for Stochastic Differential Equations, Annals of
the Institute of Statistical Mathematics, 63 (2011), pp. 431–479.