11institutetext: KU Leuven, De Nayer Campus, Dept. of Computer Science
J.-P. De Nayerlaan 5, 2860 Sint-Katelijne-Waver, Belgium 22institutetext:
Leuven.AI - KU Leuven Institute for AI, 3000 Leuven, Belgium 33institutetext:
Flanders Make@KU Leuven, 3000 Leuven, Belgium
33email<EMAIL_ADDRESS>
*Authors contributed equally

# EmoCAM: Toward Understanding What Drives CNN-based Emotion Recognition

Youssef Doulfoukar 11 Laurent Mertens 1122 0000-0001-5175-2673 Joost Vennekens
112233 0000-0001-5175-2673

###### Abstract

Convolutional Neural Networks are particularly suited for image analysis
tasks, such as Image Classification, Object Recognition or Image Segmentation.
Like all Artificial Neural Networks, however, they are “black box” models, and
suffer from poor explainability. This work is concerned with the specific
downstream task of Emotion Recognition from images, and proposes a framework
that combines CAM-based techniques with Object Detection on a corpus level to
better understand on which image cues a particular model, in our case EmoNet,
relies to assign a specific emotion to an image. We demonstrate that the model
mostly focuses on human characteristics, but also explore the pronounced
effect of specific image modifications.

###### Keywords:

Artificial Neural Networks Convolutional Neural Networks Explainable AI
Emotion Recognition Computer Vision

##  1 Introduction

Thanks to recent progress, image analysis problems such as Object Detection
using Artifical Neural Networks (ANN) can be more or less considered to be
solved [19, 8]. However, higher-order tasks, such as identifying the emotion
content of an entire image, remain more challenging. Convolutional Neural
Network (CNN) models such as EmoNet [9] present a promising approach in this
area, but its results are not yet completely convincing. This raises the
question to which extent this network is actually picking up meaningful cues
in the images, and to what extent it is learning spurious correlations that
may be present in the private dataset on which it was trained.

ANNs are still considered “black box” models, and the domain that attempts to
untangle how they make the predictions they make, i.e., to improve their
_explainability_ , is a very active one [16, 18, 2]. One of the techniques for
this is Class Activation Maps [20], or CAM, which allows to highlight those
parts of the image that contributed most to a model’s (say, a CNN image
classifier) output. This technique allows to visually inspect individual
images or videos, but does not immediately allow for an automated global
analysis on a corpus level.

To answer the earlier question of what image cues a CNN-based Emotion
Recognition network, in casu EmoNet, most relies on, we propose EmoCAM. Our
framework combines two information streams, namely CAM and Object Detection,
to build a pipeline that allows to determine those object classes that most
contributed to the model’s decision making on a corpus level. Besides better
understanding what object classes the model relies most upon, we also want to
explore the potential of applying minor changes to the input images that steer
the model towards a specific emotion by leveraging the obtained information
from our EmoCAM analysis. Our source code can be found at
https://gitlab.com/EAVISE/lme/emocam.

The remainder of this paper is organized as follows: in Sec. 2 we describe our
proposed framework in detail, followed by Sec. 3 where we look at a concrete
case using the EmoNet network and FindingEmo [11] image dataset; limitations
and roads for future work are explored in Sec. 4 and we present concluding
remarks in Sec. 5.

##  2 Methodology

We start our analysis by applying, for a given CNN model $M$ and corpus $D$,
the following steps to each image $I\in D$, schematically illustrated in Fig.
1.

Figure 1: Schematic illustrating pipeline combining CAM with Object Detection.
Photo by Sander Sammy on Unsplash, https://tinyurl.com/2tc69hfy, Unsplash
license.

First, we process $I$ with the object detection network of our choice, in
casu, YOLOv3 [14] trained on the Open Images dataset.111We use the PyTorch
LightNet [12] implementation, and OpenImages weights available from
https://pjreddie.com/darknet/yolo/. We opted for this particular pretrained
network as other popular Object Detection dataset choices such as PASCAL VOC
(20 classes) and MS-COCO (80 classes) are too restricted in the classes they
propose. By contrast, Open Images, which contains 601 classes, presents a nice
balance between human-related classes (e.g., “human face”, “mouth”, etc.), and
more general classes representing contextual elements (e.g., “car”, “tree”,
etc.). The result of this operation is a list of detected objects and their
corresponding bounding boxes $B$. We filter the YOLOv3 output by keeping only
bounding boxes with an IoU score $>0.005$.

Second, we process $I$ with $M$, and apply a CAM-based technique $C$ to the
last convolutional layer of $M$.222We experimented with other layers, as well
as combination of layers, but the best results were obtained using only the
last layer. This gives us an activation map that we overlay on top of $I$ to
obtain a new image $A$.

Finally, we lay the bounding boxes $B$ on top of $A$, and look for those boxes
$b\in B$ for which the average CAM activation, or _importance_ ,
$C_{Act}>0.3$, with $C_{Act}$ defined as the sum of the CAM activations within
the box divided by the area of the box.333To be clear, for this we take the
activation values from the original grayscale CAM output, not from the
colorized version of the output. The threshold was heuristically determined by
visually inspecting a limited set of images. We refer to these boxes as the
set $B^{*}$, and interpret these as those objects that most contributed to the
model’s decision.

Once we have found $B^{*}$ for every image in our corpus, we then analyse
these data to find associations between object classes and output labels by
constructing an association matrix $M_{A}$, where $a_{ij}\in M_{A}$ represents
the number of images labeled with the $j^{th}$ EmoNet emotion label in which
the $i^{th}$ object class has been detected at least once. By dividing each
column $j$ through by the total number of images labeled with the $j^{th}$
emotion such as to obtain percentages (after doing $\times 100$), we obtain
$M^{\prime}_{A}$ which allows to ignore imbalances in the prediction rates of
the different emotions.

##  3 Results

We tested our proposed approach using the EmoNet model and the FindingEmo
dataset. EmoNet is a model obtained through replacing the last layer of an
AlexNet model pretrained on the ImageNet [4] corpus. This last layer was then
trained on a private dataset of 137,482 images annotated for the emotion they
evoke in the observer with one of 26 custom emotion labels. We use the Python
port by L. Mertens [10] of the original Matlab release.

FindingEmo is an image dataset consisting of 25,869 images annotated for,
a.o., the dominant emotion in the picture, using one of the 24 emotion labels
in Plutchik’s Wheel of Emotions [13]. All images represent multiple people in
various natural settings and with varying degrees of interaction among them.

We first present detailed results for Grad-CAM [17] in Sec. 3.1, and follow
this up with an exploration of the effect of using other CAM-based methods in
Sec. 3.2.444For the CAM analysis, We use the Python packages grad-cam [6] and
captum https://pypi.org/project/captum/. Finally, we briefly explore the
effect on the predicted label of artificially adding certain objects to
images, attempting to answer the question whether the presence of certain
objects can _cause_ a specific label to be predicted.

###  3.1 Results for Grad-CAM

A heatmap depicting $M^{\prime}_{A}$ as obtained using Grad-CAM together with
EmoNet applied on the FindingEmo corpus can be found in Fig. 3. We limit
ourselves to the 25 most prominent Open Images classes (as determined by the
average of the corresponding row in $M^{\prime}_{A}$). A clear conclusion to
be drawn from this graph is that human features do indeed contribute the most
to the decision making, most particularly the human face which, except for
“Clothing”, represents the most important class for each EmoNet label.

Figure 2: Association between Open Images classes and predicted EmoNet label.
Heatmap entries represent the percentage of images labeled with a certain
EmoNet label for which at least one object of the corresponding Open Images
class was detected with high enough importance. “Aest. Appr.” = Aesthetic
Appreciation.

Additionally, some more specific associations do manifest themselves. Clear
examples are the association between “Sports equipment” and “Excitement”, and
“Food” and “Craving”, both of which seem logical. Less clear is, e.g., the
association between “Furniture” and “Interest”, or “Plant” and “Surprise”,
which hint of spurious associations resulting from biases in either or both
the EmoNet training dataset and FindingEmo.

###  3.2 Comparison of CAM Methods

To answer the question to what extent different CAM methods yield different
results, we performed a Representational Similarity Analysis [1] as follows.
For each CAM method $C\in\\{$Grad-CAM, Ablation-CAM[5], LIME[15], LRP[3],
LIFT-CAM[7]$\\}$, we determine the association matrix $M_{A}$ with all Open
Images classes as described in Sec. 2, keeping the same emotion and class
ordering for each $C$555Which ordering is used does not matter, as long as it
is consistently used.. We then flatten each matrix by concatenating all rows,
turning it in to a 1D vector $V_{M_{C}}$. Finally, we construct a matrix $R$
where each entry $R_{C{C^{\prime}}}$ represents the Spearman Correlation rank
between $V_{M_{C}}$ and $V_{M_{C^{\prime}}}$. The resulting matrix is shown in
Fig. 3. All related $p$-values were $<\\!\\!<0.05$, indicating statistical
significance.

Figure 3: RSA analysis of different CAM methods.

The consistently high correlation values between all pairs indicate that
variations in results obtained through different CAM methods can be expected
to be minimal. We did observe both LIFT-CAM and LRP resulting in a notable
association between “Pillow” and “Sexual Desire”. Other than this, the
differences between the methods appear to lie within the relative strengths of
the associations observed, rather than the assocations themselves.

###  3.3 Prediction Stability

To illustrate how the obtained knowledge can be applied to fool the network by
creating an adversarial attack, consider the image shown in Fig. 4. We know
from Sec. 3.1 that there is a high association between the object category
“Sports equipment” and EmoNet label “Excitement”. This inspired us to take an
image labeled with high probabilty as “Joy” (92.9%; Excitement: 1.8%). After
altering this image by pasting a rugby ball on top of the head of one of the
two main subjects, the prediction changes to 66.1% “Excitement” (“Joy”:
30.2%), demonstrating the dramatic effect the presence of a particular object
can have on the model’s output.

Figure 4: Adversarial example. Original image on the left, labeled by EmoNet
as 92.9% “Joy”. Modified part on the right; full modified image labeled as
66.1% “Excitement”. Sources: original photo by istolethetv,
https://tinyurl.com/3r86e3w2, CC 2.0 license; Rugby ball by Peter Griffin,
https://tinyurl.com/dmb77rks, CC0 license.

Note that the position of the pasted object greatly influences the effect it
has. Moving the rugby ball to the immediate right of the subject’s face alters
the predictions to 43.9% “Joy” and 42.1% “Excitement”, while moving it to the
immediate left only alters the predictions by 4% in the same directions
(“Joy”: 34.2%; “Excitement”: 62.1%). Covering the other subject’s head
instead, we obtain 52.0% “Joy” and 30.7% Excitement.

We further investigate this effect by performing the following experiment. For
each $I$ in $D$, we paste a given object $O\in\\{$Rugby ball, Soccer ball,
Lotus flower$\\}$, resized such that its height equals
$0.2~{}\times~{}\mbox{height}(I)$, in $I$ centered at each one of a set of
predefined relative positions $P$ within the image, resulting in $size(P)$
alterations to $I$. The positions and objects considered are illustred in Fig.
5. We then send the altered images through EmoNet, and observe how the
prediction was affected. The results are shown in Fig. 6.

Figure 5: Schematic illustration of “paste object in image” experiment. The
grid on the left illustrates the relative positions within the image the
objects are pasted and centered at, with the considered objects shown on the
right. Sources: Rugby ball, see Fig. 4; Soccer ball by Jean Schecter,
https://tinyurl.com/3xxkkysb, CC 4.0 BY-NC; Lotus flower at
https://pngimg.com/image/69752, CC 4.0 BY-NC.

The results confirm that the model shows high sensitivity to certain objects.
Although differences of up to more than 10% can be observed in Fig. 6,
specifically for the rugby ball, no real tendencies reveal themselves. In
combination with the example in Fig. 4, we hypothesize the differences are not
so much due to the _absolute_ position of the object, but to what it
_occludes_. For the rugby ball, 4 out of 17 positions resulted most often in a
label switch to “Excitement”. For the soccer ball, the number increases to 7.
The Lotus flower clearly results in much less label shifts overall, with not a
single position favoring “Excitement”, confirming the importance of the object
class in effecting a label switch.

Figure 6: Percentage of images in the corpus whose predicted label changed
when a specific object was pasted at a specific position. Refer to Fig. 5 for
the positions and objects. The letters ‘A’ for ‘Admiration’, ‘E’ for
‘Excitement’ and ‘J’ for ‘Joy’ atop each bar indicate the EmoNet label most
often switched to.

##  4 Limitations and Future Work

Although the currently described approach already provides valuable insights,
some limitations are to be noted.

First, the approach is, by definition, heavily dependent on the choice of
Object Detection network and its corresponding classes and performance. The
upside is that, as a plug-and-play component, different Object Detection
networks can be chosen for different tasks, allowing to pick object classes
tailored to the task at hand.

Second, our current implementation does not take into account the size of the
bounding boxes, which can result in suboptimal results. Consider, e.g., the
example shown in Fig. 7. Although the subject’s ear is clearly not the most
important contributing element in the picture, because of the small size of
the “Human ear” bounding box the average CAM activation is nonetheless the
highest, spuriously pushing this object class to the top. Two main paths could
be explored to counter this issue. The most straightforward would be to
develop a scoring function that does take into account the bounding box size,
or the activation distribution within it. Alternatively, segmentation models
could be used instead of bounding box detection models, so as to obtain
clearly delineated zones representing the different objects. Then of course,
our first limitation still applies, i.e., the segmentation classes need to be
relevant to the task at hand.

Figure 7: An instance were the current EmoCAM approach fails to detect the
most important class. Source: original image by Darius Bashar on Unsplash,
https://tinyurl.com/4eephaxb, Unsplash license.

Third, with regard to our experiment described in Sec. 3.3, a more interesting
approach might be to, instead of, or along with, considering only a fixed set
of positions, paste the object at the center of bounding boxes relating to
specific features such as human heads, thus investigating the effect of
masking specific objects. We also intend to apply EmoCAM to the modified
images to explore if the shift in label is reflected in a shift in focus in
the modified image.

##  5 Conclusion

We propose the novel EmoCAM approach to explaining CNN decisions, specifically
with the downstream task of Emotion Recognition from images in mind. Our
objective is threefold: 1) better understanding what parts of the input image
the model uses to make its decision, 2) allowing to check whether or not the
information used by the model aligns with expectations from a human
perspective, and 3) uncovering potential model biases. We have demonstrated
our approach using the EmoNet model, FindingEmo dataset and multiple CAM
techniques. Using our approach, we found that EmoNet indeed shows a strong
focus on human elements, most notably (parts of) the human face, which is
encouraging as it aligns with our understanding of human emotion recognition
from Psychology. Nevertheless, we also found the model output to be quite
unstable, in that adding specific objects (e.g., a rugby ball) to an image can
dramatically alter its output and steer it towards a specific target emotion
(e.g., “Excitement”).

## References

  * [1] Chapter 27 - representational similarity analyses: A practical guide for functional mri applications. In: Manahan-Vaughan, D. (ed.) Handbook of in Vivo Neural Plasticity Techniques, Handbook of Behavioral Neuroscience, vol. 28, pp. 509–525. Elsevier (2018). https://doi.org/https://doi.org/10.1016/B978-0-12-812028-6.00027-6, https://www.sciencedirect.com/science/article/pii/B9780128120286000276
  * [2] Ali, S., Abuhmed, T., El-Sappagh, S., Muhammad, K., Alonso-Moral, J.M., Confalonieri, R., Guidotti, R., Del Ser, J., Díaz-Rodríguez, N., Herrera, F.: Explainable artificial intelligence (xai): What we know and what is left to attain trustworthy artificial intelligence. Information Fusion 99, 101805 (2023). https://doi.org/https://doi.org/10.1016/j.inffus.2023.101805, https://www.sciencedirect.com/science/article/pii/S1566253523001148
  * [3] Bach, S., Binder, A., Montavon, G., Klauschen, F., Müller, K.R., Samek, W.: On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation. PLoS ONE 10 (2015), https://api.semanticscholar.org/CorpusID:9327892
  * [4] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.: Imagenet: A large-scale hierarchical image database. In: 2009 IEEE conference on computer vision and pattern recognition. pp. 248–255. Ieee (2009)
  * [5] Desai, S., Ramaswamy, H.G.: Ablation-cam: Visual explanations for deep convolutional network via gradient-free localization. In: 2020 IEEE Winter Conference on Applications of Computer Vision (WACV). pp. 972–980 (2020). https://doi.org/10.1109/WACV45572.2020.9093360
  * [6] Gildenblat, J., contributors: Pytorch library for cam methods. https://github.com/jacobgil/pytorch-grad-cam (2021)
  * [7] Jung, H., Oh, Y.: Towards better explanations of class activation mapping. In: 2021 IEEE/CVF International Conference on Computer Vision (ICCV). pp. 1316–1324 (2021). https://doi.org/10.1109/ICCV48922.2021.00137
  * [8] Kaur, R., Singh, S.: A comprehensive review of object detection with deep learning. Digital Signal Processing 132, 103812 (2023). https://doi.org/https://doi.org/10.1016/j.dsp.2022.103812, https://www.sciencedirect.com/science/article/pii/S1051200422004298
  * [9] Kragel, P.A., Reddan, M.C., LaBar, K.S., Wager, T.D.: Emotion schemas are embedded in the human visual system. Science Advances 5(7), eaaw4358 (2019)
  * [10] Mertens, L.: Emonet: A pytorch port. https://gitlab.com/EAVISE/lme/emonet (2022)
  * [11] Mertens, L., Yargholi, E., de Beeck, H.O., den Stock, J.V., Vennekens, J.: Findingemo: An image dataset for emotion recognition in the wild (2024), https://arxiv.org/abs/2402.01355
  * [12] Ophoff, T.: Lightnet: Building blocks to recreate darknet networks in pytorch. https://gitlab.com/EAVISE/lightnet (2018)
  * [13] PLUTCHIK, R.: Chapter 1 - a general psychoevolutionary theory of emotion. In: Plutchik, R., Kellerman, H. (eds.) Theories of Emotion, pp. 3–33. Academic Press (1980). https://doi.org/https://doi.org/10.1016/B978-0-12-558701-3.50007-7, https://www.sciencedirect.com/science/article/pii/B9780125587013500077
  * [14] Redmon, J., Farhadi, A.: Yolov3: An incremental improvement. CoRR abs/1804.02767 (2018), http://arxiv.org/abs/1804.02767
  * [15] Ribeiro, M.T., Singh, S., Guestrin, C.: "why should i trust you?": Explaining the predictions of any classifier. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 1135–1144. KDD ’16, Association for Computing Machinery, New York, NY, USA (2016). https://doi.org/10.1145/2939672.2939778, https://doi.org/10.1145/2939672.2939778
  * [16] Saeed, W., Omlin, C.: Explainable ai (xai): A systematic meta-survey of current challenges and future opportunities. Knowledge-Based Systems 263, 110273 (2023). https://doi.org/https://doi.org/10.1016/j.knosys.2023.110273, https://www.sciencedirect.com/science/article/pii/S0950705123000230
  * [17] Selvaraju, R.R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., Batra, D.: Grad-cam: Visual explanations from deep networks via gradient-based localization. In: 2017 IEEE International Conference on Computer Vision (ICCV). pp. 618–626 (2017). https://doi.org/10.1109/ICCV.2017.74
  * [18] Yang, W., Wei, Y., Wei, H., Chen, Y., Huang, G., Li, X., Li, R., Yao, N., Wang, X., Gu, X., Amin, M.B., Kang, B.: Survey on explainable ai: From approaches, limitations and applications aspects. Human-Centric Intelligent Systems 3, 161 – 188 (2023), https://api.semanticscholar.org/CorpusID:260817864
  * [19] Zhao, Z.Q., Zheng, P., Xu, S.T., Wu, X.: Object detection with deep learning: A review. IEEE Transactions on Neural Networks and Learning Systems 30(11), 3212–3232 (2019). https://doi.org/10.1109/TNNLS.2018.2876865
  * [20] Zhou, B., Khosla, A., Lapedriza, A., Oliva, A., Torralba, A.: Learning deep features for discriminative localization. In: 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). pp. 2921–2929 (2016). https://doi.org/10.1109/CVPR.2016.319


# Filtering Dynamical Systems Using Observations of Statistics

Eviatar Bach<EMAIL_ADDRESS>Department of Environmental Science
and Engineering, California Institute of Technology, Pasadena, California
91125, USA  Department of Computing and Mathematical Sciences, California
Institute of Technology, Pasadena, California 91125, USA  Tim Colonius
Department of Mechanical and Civil Engineering, California Institute of
Technology, Pasadena, California 91125, USA  Isabel Scherl  Department of
Mechanical and Civil Engineering, California Institute of Technology,
Pasadena, California 91125, USA  Andrew Stuart  Department of Computing and
Mathematical Sciences, California Institute of Technology, Pasadena,
California 91125, USA

###### Abstract

We consider the problem of filtering dynamical systems, possibly stochastic,
using observations of statistics. Thus, the computational task is to estimate
a time-evolving density $\rho(v,t)$ given noisy observations of the true
density $\rho^{\dagger}$; this contrasts with the standard filtering problem
based on observations of the state $v$. The task is naturally formulated as an
infinite-dimensional filtering problem in the space of densities $\rho$.
However, for the purposes of tractability, we seek algorithms in state space;
specifically, we introduce a mean-field state-space model, and using
interacting particle system approximations to this model, we propose an
ensemble method. We refer to the resulting methodology as the ensemble
Fokker–Planck filter (EnFPF).

Under certain restrictive assumptions, we show that the EnFPF approximates the
Kalman–Bucy filter for the Fokker–Planck equation, which is the exact solution
to the infinite-dimensional filtering problem. Furthermore, our numerical
experiments show that the methodology is useful beyond this restrictive
setting. Specifically, the experiments show that the EnFPF is able to correct
ensemble statistics, to accelerate convergence to the invariant density for
autonomous systems, and to accelerate convergence to time-dependent invariant
densities for non-autonomous systems. We discuss possible applications of the
EnFPF to climate ensembles and to turbulence modeling.

††preprint: AIP/123-QED

> Data assimilation (DA) is the process of estimating the state of a dynamical
> system using observations. Here, we modify the standard DA setting to allow
> for observations of _statistics_ of a system with respect to its time-
> evolving probability density. We propose a mathematical framework, a
> resulting ensemble method, and present numerical experiments demonstrating
> accelerated convergence of a system to its attractor. We propose further
> applications to problems in climate and turbulence modeling.

##  I Introduction

The goal of this paper is to introduce a filtering methodology that
incorporates statistical information into a (possibly stochastic) dynamical
system. In Secs. I.1–I.3, we present, respectively, a high-level overview of
the problem, discuss the motivation and previous literature, and outline the
paper structure and our contributions.

###  I.1 Assimilating Statistical Observations

We start by presenting a high-level overview of the problem of incorporating
statistical information into a dynamical system; a detailed problem statement
follows in Sec. II.1.

Data assimilation (DA) is overviewed in a number of books, including
Refs.Jazwinski (1970); Kalnay (2002); Law, Stuart, and Zygalakis (2015); Reich
and Cotter (2015). The problem is to estimate the state of a dynamical system
by combining noisy, partial observations with a model for the system. In the
continuous-time DA problem, we have a stochastic differential equation (SDE)

$\displaystyle dv^{\dagger}$
$\displaystyle=f(v^{\dagger},t)\,dt+\sqrt{\Sigma(t)}\,dW,$ (1)
$\displaystyle{v^{\dagger}(0)}$ $\displaystyle{=v^{\dagger}_{0},}$ (2)

with solution $v^{\dagger}\in\mathbb{R}^{d}$, and observations given by

$dz^{\dagger}=h(v^{\dagger}(t),t)\,dt+\sqrt{\Gamma(t)}\,dB,$ (3)

with $z^{\dagger}\in\mathbb{R}^{p}$. The equations for $v^{\dagger}$ and
$z^{\dagger}$ are driven by independent standard Wiener processes $W$ and $B$.
These SDEs, as with all the SDEs in the paper, are to be interpreted in the
Itô sense. Filtering is then the problem of obtaining the best possible
estimate of the posterior density on $v^{\dagger}(t)$ given the past
observations $\\{z^{\dagger}(s)\\}_{s\in[0,t]}$. Throughout the paper, we use
the $\dagger$ superscript to indicate the true quantities, and omit it for
filtered quantities.

Instead of observing a specific trajectory of a dynamical system, as
$\\{z^{\dagger}(t)\\}$ given by Eq. (3) does, one can also consider
_observations of the system’s statistical behavior_ , that is, observations of
functionals of the probability density $\rho^{\dagger}(v,t)$ over
trajectories. This density reflects the randomness from the initial conditions
for $v$ and/or from the Brownian forcing. For a deterministic dynamical system
($\Sigma\equiv 0$), if the initial conditions are random, then
$\rho^{\dagger}(v,t)$ will reflect the changing density over time under the
action of the system’s dynamics, governed by the Liouville equation.111Here,
we use the term Liouville equation for the equation governing evolution of the
density of any ordinary differential equation, not just in the Hamiltonian
setting. If noise is present, the changing density is also affected by the
Brownian noise $W$ and is governed by the Fokker–Planck equation, a
diffusively regularized Liouville equation. In this paper we focus on
observations of $\rho^{\dagger}(v,t)$ defined by replacing Eq. (3) with

$dz^{\dagger}=\left(\int\mathfrak{h}(v,t)\rho^{\dagger}(v,t)\,dv\right)dt+\sqrt{\Gamma(t)}\,dB.$
(4)

Here, $\mathfrak{h}(v,t)$ defines the observed statistics of $v$, $B$ is a
Wiener process, and $z^{\dagger}\in\mathbb{R}^{p}$. The filtering problem is
to estimate a density $\rho(v,t)$ given all the past observations
$\\{z^{\dagger}(s)\\}_{s\in[0,t]}$. As in the observation equation (3), the
observations are finite-dimensional, noisy, and partial. However, since the
observations are now of $\rho^{\dagger}(v,t)$ instead of $v^{\dagger}(t)$, we
must specify the dynamics of $\rho^{\dagger}(v,t)$. This is given by the
Fokker–Planck (FP) or Kolmogorov forward equation,

$\displaystyle\frac{\partial\rho^{\dagger}}{\partial t}$
$\displaystyle=\mathcal{L}^{*}(t)\rho^{\dagger},$ (5a)
$\displaystyle\mathcal{L}^{*}(t)\psi$ $\displaystyle=-\nabla\cdot(\psi
f)+\frac{1}{2}\nabla\cdot\bigl{(}\nabla\cdot(\psi\Sigma)\bigr{)},$ (5b)

where $\mathcal{L}^{*}$ is the adjoint of the generator of Eq. (1).222We
define the divergence of a matrix as is standard in continuum mechanics; see
Gurtin (1981)Gurtin (1981) and Gonzalez and Stuart (2008)Gonzalez and Stuart
(2008). The divergence of a matrix $S$ is defined by the identity
$(\nabla\cdot S)\cdot a=\nabla\cdot(S^{T}a)$ holding for any vector $a$.. For
a deterministic system, with $\Sigma\equiv 0$, the Fokker–Planck equation
reduces to the Liouville equation.

An important question is how one would obtain observations of a system’s
statistics for problems of practical relevance. We discuss this in detail in
Sec. I.2.1. For now, we proceed on the assumption that $z^{\dagger}$ solving
Eq. (4) is given.

Now, Eqs. (5) and (4) define a filtering problem for $\rho(v,t)$. This is an
infinite-dimensional filtering problem, in contrast to the finite-dimensional
filtering problem for $v(t)$ defined by Eqs. (1) and (3). We refer to the
filtering problem defined by Eqs. (5) and (4) as the _Fokker–Planck filtering
problem_. Note that both Eqs. (5) and (4) are _linear_ in $\rho^{\dagger}$,
meaning that the solution to the problem can be written using the infinite-
dimensional Kalman–Bucy (KB) filter; see Sec. IV.1 for more details.

Despite the existence of an exact solution to the filtering problem, through
the infinite-dimensional Kalman–Bucy (KB) filter, approximating the Gaussian
conditional density $\rho$ is in most setting computationally intractable
since the mean is a probability density function and the covariance is an
operator. Thus we seek inspiration from the success of ensemble Kalman
filtering Evensen (1994): we work in state space and seek an ensemble that
evolves in time a number of states whose empirical density approximates the
filtered $\rho$. We note that the particle filter similarly substitutes the
problem of evolving a probability density with that of evolving a number of
particles and weights (Crişan, Del Moral, and Lyons, 1999). Furthermore,
derivation of ensemble Kalman methods via a mean-field limit provides a
systematic methodology for the derivation of equal-weight approximate filters
(Calvello, Reich, and Stuart, 2022). We call the resulting method the
_ensemble Fokker–Planck filter_ (EnFPF).

Figure 1: The density of an Ornstein–Uhlenbeck process evolving in time (top
panel). At regular intervals, we make observations of this density and use
them to inform the evolution of an ensemble (bottom panel).

Figure 1 shows a schematic of such an ensemble method. In the top panel is the
true time-varying probability density, in this case of an Ornstein–Uhlenbeck
process. In the bottom panel is an ensemble of states. At regular intervals,
we observe expectations over the density in the top panel. Using these
observations and our model of the system, we evolve the ensemble over the time
interval between the current and next observations.

###  I.2 Motivation and Literature Review

The subject of Kalman filtering and Kalman–Bucy (KB) filtering in infinite-
dimensional spaces is studied in the control theory literature (Curtain and
Pritchard, 1978). We emphasize that although we sketch out the basic
mathematical foundations of the Fokker–Planck filtering problem in Sec. IV,
many interesting mathematical problems in analysis and probability remain open
in this area. To the best of our knowledge, the methodology proposed here is
the first general method for assimilating observations of statistics directly
into a state-space formulation of dynamical systems. Our methodology is built
on the conceptual approach introduced in the feedback particle filter (Yang,
Mehta, and Meyn, 2013; Reich, 2019), and earlier related work (Crisan and
Xiong, 2010), seeking a mean-field model that achieves the goal of filtering
and can be approximated by particle methods (Pathiraja, Reich, and Stannat,
2021); in particular, we seek particle approximations of the mean-field model
inspired by ensemble Kalman methods (Calvello, Reich, and Stuart, 2022).

The problem of recovering a probability density from a finite number of known
moments is called a moment problem. When $\mathfrak{h}$ in Eq. (4) consists of
monomials in $v$, the problem of reconstructing $\rho$ is similar to a moment
problem, with the major difference that $\rho$ evolves in time according to a
dynamical system. Moment problems are typically regularized by a maximum
entropy approach (Abramov, 2010); in the Fokker–Planck filtering problem,
regularization is provided by the system’s dynamics.

Our motivation comes from a number of applications around which we organize
the remainder of our literature review, after first discussing the general
question of how to obtain observations of statistics.

####  I.2.1 Obtaining observations of statistics

In typical applications, one can only observe a single trajectory of a
dynamical system, and thus the statistics of the density will not be directly
available. If we are interested in the statistics of the invariant measure, as
we are for several of the applications discussed below, then for ergodic
systems we have that

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathfrak{h}(v^{\dagger}(t))\,dt=\int\mathfrak{h}(v)\rho^{\dagger}(v)\,dv,$
(6)

where $\rho^{\dagger}$ is the invariant density, and thus an approximation of
the statistics of the invariant measure can be obtained from a long observed
or simulated trajectory of the dynamical system.

For nonautonomous systems, due to lack of ergodicity, observations of the
statistics cannot be made using long time averages. If the nonstationary
forcing is slow enough, however, an adiabatic approximation, in which the fast
scales are considered to be ergodic with an invariant measure parameterized by
the value of the slow forcing, may be justified (Pavliotis and Stuart, 2008;
Drótos, Bódai, and Tél, 2016). If the forcing is periodic, then observations
of the phase-dependent statistics could be obtained by averaging the
observables at a given phase over multiple periods.

For certain systems, invariant statistics may be acquired analytically or by
numerically solving a different set of equations. For example, for the
Navier–Stokes equations, the Reynolds-averaged Navier–Stokes (RANS) equations
can be used to approximate the stationary statistics.

It may be possible to instead formulate a filtering problem using an
observation operator that involves averaging over a finite time window; we
leave this for future work. This problem was considered in Ref. Dirren and
Hakim (2005), but only a heuristic solution was proposed. We note that other
works have made use of observation operators with time-delayed observations
Rey  _et al._ (2014); Bach  _et al._ (2021), albeit for different purposes.

In Secs. I.2.2–I.2.5, we review the possible applications of the ensemble
Fokker–Planck filter.

####  I.2.2 Acceleration of convergence to a (possibly time-dependent)
invariant measure

Acceleration of the time to convergence of dynamical models to their invariant
measure (often referred to as the “spin-up” period or the transient) is of
importance in many fields, including climate (Bryan, 1984; Daron and
Stainforth, 2015; Drótos, Bódai, and Tél, 2017; Deser  _et al._ , 2020) and
other fluid problems (Wang, Xin, and Zhang, 2022), Langevin sampling (Hwang,
Hwang-Ma, and Sheu, 2005; Abdulle, Pavliotis, and Vilmart, 2019), and
turbulence simulation (Nelson and Fringer, 2017).

For a stochastic differential equation with an invariant measure, under
conditions described in Goldys and Maslowski (2005)Goldys and Maslowski
(2005), the convergence to this invariant measure is exponential with an
exponent related to the spectral gap of the corresponding generator.

In this paper, we show that this convergence can be accelerated using the
ensemble Fokker–Planck filter, and this is the primary application we test in
the numerical experiments. In particular, if some statistics of the invariant
measure are known, these statistics can be assimilated into the ensemble,
obtaining an ensemble whose empirical density is closer to the invariant
measure.

To our knowledge, existing methods of accelerating convergence of model
trajectories to the invariant measure have been problem-dependent, as in Bryan
(1984)Bryan (1984). Isik (2013)Isik (2013) and Isik, Takhirov, and Zheng
(2017)Isik, Takhirov, and Zheng (2017) studied a relaxation-based method of
accelerating the convergence to equilibrium of the Navier–Stokes equations,
which bears some resemblance to our approach.

Non-autonomous (also referred to as non-stationary) and random dynamical
systems can have time-dependent attractors, known as pullback attractors, to
which the evolution converges Arnold (1998). A pullback attractor is the set
that the dynamical system approaches when evolved in time from the infinite
past to a fixed time (say time $0$ without loss of generality). We refer to
the probability measure associated with these attractors as time-dependent
invariant measures, following Chekroun, Simonnet, and Ghil (2011)Chekroun,
Simonnet, and Ghil (2011). These objects are of considerable interest for
climate (Ghil, Chekroun, and Simonnet, 2008; Chekroun, Simonnet, and Ghil,
2011; Daron and Stainforth, 2015). The EnFPF can also accelerate convergence
to these invariant measures.

The problem of accelerating convergence to the invariant measure is related to
the problem of controlling the Fokker–Planck equation, where a density is
controlled in order to reach to a specified target distributionAnnunziato and
Borzì (2013), and to statistical control, wherein one aims to return a
perturbed system to its equilibrium statisticsCovington, Qi, and Chen (2023).

Furthermore, the EnFPF could be tested for accelerating the convergence of
sampling algorithms such as Langevin sampling and Markov chain Monte Carlo,
when some statistics of the target density are known _a priori_.

Finally we note that when estimating Koopman or Perron–Frobenius operators, it
is often necessary to have a large number of trajectories from initial
conditions sampled from the invariant measure.

####  I.2.3 Parameter estimation

The EnFPF could be used for jointly updating states and parameters using
statistical observations by adopting a state augmentation approach. Other work
has adapted methods from data assimilation for parameter estimation using
time-averaged statistics, assumed to be close to the statistics on the
invariant measure by ergodicity (Köhl and Willebrand, 2002; Schneider, Stuart,
and Wu, 2021; Mons, Du, and Zaki, 2021; Schneider, Stuart, and Wu, 2022).

####  I.2.4 Correcting for model error

Generally, methods that correct for model error are formulated in terms of
forecast performance at some lead time, e.g., Levine and Stuart (2022) and
Bach and Ghil (2023). If one is instead interested in correcting statistical
properties, one can postulate a parametric form for the model error and use
time-averaged observations to estimate the parameters, as discussed in the
preceding paragraph. Alternatively, the EnFPF could be tested for directly
correcting model error using statistical observations in a similar manner to
the use of classical DA in reducing the impact of model error for forecast
applications Danforth and Kalnay (2008); Chen  _et al._ (2022). The analysis
increments could then be taken to approximate model error corrections, and
training a machine learning model to predict these corrections could be
tested, as has been done for classical DA (Danforth and Kalnay, 2008; Farchi
_et al._ , 2021; Chen  _et al._ , 2022).

Statistical properties have previously been used to learn closure models for
the Navier–Stokes equation using a 3DVar-like scheme Ephrati  _et al._ (2023).

####  I.2.5 Assimilation of time-averaged observations

In paleoclimate, proxy records often represent time averages instead of
instantaneous measurements. Methods have been developed for making use of
time-averaged observations for state estimation in the paleoclimate data
assimilation literature (Dirren and Hakim, 2005). As discussed above, in the
case of slow forcing, time-averaged observations can be used to approximately
track the system’s time-varying statistics, enabling their use in the EnFPF.

###  I.3 Contributions and Paper Outline

The primary contributions of this work are: (i) to establish a framework for
the filtering of stochastic dynamical systems, or dynamical systems with
random initial data, given only observations of statistics; (ii) to introduce
ensemble-based state-space methods for this filtering problem via a mean field
perspective; and (iii) to demonstrate numerically that the proposed methods
are effective at guiding dynamical systems toward observed statistics. (i) is
covered in Sec. II.1 and Sec. IV; (ii) is covered in sections II.2–II.6; and
(iii) is covered in Sec. III.

In Sec. II.1, we outline the Fokker–Planck filtering problem and distinguish
it from the standard filtering problem. In Secs. II.2–II.4, we introduce a
mean-field algorithm and its particle and discrete-time approximations,
culminating in the ensemble Fokker–Planck filter (EnFPF). In Sec. II.6, we
discuss implementation details, including the approximation of the score
function and a square-root ensemble formulation with reduced computational
effort.

In Sec. III, we carry out numerical experiments with several chaotic dynamical
systems, both autonomous and non-autonomous, and based on the Lorenz63,
Lorenz96, and Kuramoto–Sivashinsky models. In particular, we demonstrate that
the EnFPF can accelerate the convergence of these systems to their invariant
densities, using information about the moments of these densities.

In Sec. IV, we provide a justification of our algorithm. We first formulate
the KB filter for densities (Sec. IV.1), which provides a solution to the
Fokker–Planck filtering problem in function space, and analyze some of its
properties in Appendix A. We then propose an ansatz amenable to a mean-field
model (Sec. IV.2) and show its equivalence to the KB filter for densities
under some assumptions (Theorem 1 in Appendix B). We then show how this ansatz
can be approximated by a mean-field model (Sec. IV.3, providing further
details in Appendix C).

Finally, in Sec. V, we give conclusions and outlook for future work.

##  II Problem and Algorithm

In Sec. II.1, we introduce the probabilistic formulation of the standard
filtering problem, and then contrast it with the Fokker–Planck filtering
problem, where data is in the form of statistics. Sec. II.2 demonstrates an
approach to this problem using a mean-field model. In Sec. II.3, we introduce
a particle approximation of the mean-field algorithm, which forms the basis of
the proposed EnFPF.

###  II.1 Problem Statement

####  II.1.1 The Standard Filtering Problem

In the standard filtering problem, we are given state observations
$z^{\dagger}(t)$ of $v^{\dagger}(t)$, defined by Eq. (3), and the dynamics of
$v^{\dagger}(t)$ are given by Eq. (1). The problem is then to find an equation
for the conditional distribution of $v|Z^{\dagger}(t)$, where
$Z^{\dagger}(t)=\\{z^{\dagger}(s)\\}_{s\in[0,t]}$ are the observations
accumulated up to time $t$ under a fixed realization of $B$. The solution to
the filtering problem is given by the Kushner–Stratonovich equation,

$\frac{\partial\rho}{\partial
t}=\mathcal{L}^{*}{(t)}\rho+\Bigl{\langle}h(v,{t})-\mathbb{E}h,\frac{dz^{\dagger}}{dt}-\mathbb{E}h\Bigr{\rangle}_{\Gamma{(t)}}\rho,$
(7)

where $\langle\cdot,\cdot\rangle_{A}\equiv\langle
A^{-1/2}\cdot,A^{-1/2}\cdot\rangle$ is the weighted Euclidean inner product.
Treatments of the standard filtering problem can be found in Jazwinski
(1970)Jazwinski (1970) and Bain and Crisan (2009)Bain and Crisan (2009).

####  II.1.2 The Fokker–Planck Filtering Problem

In this paper, we consider instead noisy observations of
$\rho^{\dagger}(v,t)$: the observation process $z^{\dagger}(\cdot)$ is given
by

$dz^{\dagger}=H(t)\rho^{\dagger}(\cdot,t)\,dt+\sqrt{\Gamma(t)}\;dB.$ (8)

Here, $H(t)$ is a linear operator mapping the space of probability densities
into a finite-dimensional Euclidean space, and the dynamics of
$\rho^{\dagger}$ are given by the Fokker–Planck equation (5). That is, we make
observations of _statistics_ of the dynamical system. We refer to the problem
of finding the conditional density of $v|Z^{\dagger}(t)$, where
$Z^{\dagger}(t)=\\{z^{\dagger}(s)\\}_{s\in[0,t]}$ is given by Eq. (8), as the
_Fokker–Planck filtering problem_. In the following subsection, we propose an
approximation to the solution to this problem in state space.

###  II.2 Mean-Field Equation

Although in Sec. IV.1 we treat the Fokker–Planck filtering problem for more
general $H$, in the rest of what follows we focus on the setting where

$H(t)\rho=\mathbb{E}[\mathfrak{h}(v,t)]=\int\mathfrak{h}(v,t)\rho(v,t)\,dv,$
(9)

for some $\mathfrak{h}$. With this assumption on $H$, Eq. (8) reduces to Eq.
(4). In particular, if $\mathfrak{h}$ is a monomial in $v$, e.g.,
$\mathfrak{h}(v)=v$ or $\mathfrak{h}(v)=\operatorname{vec}(v\otimes v)$, then
$H\rho$ will correspond to moments of $\rho$. We will henceforth use
$\mathbb{E}$ to denote expectation under $\rho$, unless otherwise indicated.

######  Remark 1.

Note that if $\rho^{\dagger}(v,0)=\delta(v-v^{\dagger}_{0})$ for some
$v^{\dagger}_{0}$ and $\Sigma=0$, then the Fokker–Planck filtering problem is
equivalent to the standard filtering problem with
$v^{\dagger}(0)=v^{\dagger}_{0}$, observation operator $\mathfrak{h}$, and
$\Sigma=0$.

Our proposed methodology is to introduce a mean-field model for variable $v$,
depending on its own probability density function $\rho(v,t)$. The mean-field
model is chosen to drive the system toward the observed statistical
information. Algorithms are then based on particle approximation of this
model, leading to ensemble Kalman–type methods. The mean-field model
considered is

$\displaystyle dv$
$\displaystyle=f(v,t)\,dt+\sqrt{\Sigma(t)}\,dW+K(t)\bigl{(}dz^{\dagger}-d\hat{z}\bigr{)},$
(10a) $\displaystyle d\hat{z}$
$\displaystyle=(\mathbb{E}\mathfrak{h})(t)\,dt+\sqrt{\Gamma(t)}dB,$ (10b)
$\displaystyle K(t)$ $\displaystyle=C^{v\mathfrak{h}}(t)\Gamma(t)^{-1},$ (10c)
$\displaystyle C^{v\mathfrak{h}}(t)$
$\displaystyle=\mathbb{E}\bigl{[}\bigl{(}v(t)-\mathbb{E}v(t)\bigr{)}\bigl{(}\mathfrak{h}(v,t)-(\mathbb{E}\mathfrak{h})(t)\bigr{)}^{T}\bigr{]}.$
(10d)

The terms in the mean-field model can be understood intuitively as follows.
The first two terms on the right-hand side of Eq. (10a) are simply the
dynamics of the system (1). The third term resembles the standard nudging
observer term from control theory, with an ensemble Kalman–inspired gain, and
the use of noisy simulated data, as in the stochastic ensemble Kalman filter.

In some problems, we find that it is beneficial to include an additional
score-based term in the model, replacing Eq. (10a) by

$\displaystyle dv$
$\displaystyle=f(v,{t})\,dt+\sqrt{\Sigma{(t)}}\,dW+K{(t)}\Bigl{(}dz^{\dagger}-d\hat{z}\Bigr{)}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad+K{(t)}\Gamma{(t)}K{(t)}^{T}\nabla\log\rho{(v,t)}\,dt.$
(11)

The additional term induces negative diffusion in the equation for the density
of $v$, exactly balancing the diffusion introduced through $z^{\dagger}$ and
$\hat{z}$. We justify equations (10) and (II.2) in detail in Sec. IV by
building on the Fokker–Planck picture in density space.

###  II.3 Particle Approximation of Mean-Field Equation

In order to tractably implement the mean-field equations (10), we use a
particle (or ensemble) approximation. That is, given $J$ particles, we
consider the following interacting particle system for
$\\{v^{(j)}\\}_{j=1}^{\mathsf{J}}$:

$\displaystyle dv^{(j)}$
$\displaystyle=f(v^{(j)},t)\,dt+\sqrt{\Sigma(t)}\,dW^{(j)}+K(t)\bigl{(}dz^{\dagger}-d\hat{z}^{(j)}\bigr{)},$
(12a) $\displaystyle d\hat{z}^{(j)}$
$\displaystyle=(\mathbb{E}^{\mathsf{J}}\mathfrak{h})(t)\,dt+\sqrt{\Gamma(t)}\,dB^{(j)},$
(12b) $\displaystyle K(t)$
$\displaystyle=(C^{v\mathfrak{h}}(t))^{\mathsf{J}}\,\Gamma(t)^{-1}.$ (12c)

Here, $\mathbb{E}^{\mathsf{J}}$ denotes expectation with respect to the
empirical measure formed by equally weighting Dirac measures at the particles
$\\{v^{(j)}\\}_{j=1}^{\mathsf{J}}$; $(C^{v\mathfrak{h}})^{\mathsf{J}}$ denotes
the sample cross-covariance computed using this empirical measure,

$C^{v\mathfrak{h}}(t)=\mathbb{E}^{\mathsf{J}}\bigl{[}\bigl{(}v(t)-\mathbb{E}v(t)\bigr{)}\bigl{(}\mathfrak{h}(v,t)-(\mathbb{E}\mathfrak{h})(t)\bigr{)}^{T}\bigr{]}.$

Note that unlike the ensemble Kalman filter, the predicted observation for
each ensemble member, Eq. (12b), involves the expectation of $\mathfrak{h}$
over the ensemble, instead of the observation operator applied to that
ensemble member.

###  II.4 Discrete-Time Approximation of Mean-Field Equation

A discrete-time analog of Eqs. (12) is given by

$\displaystyle\hat{v}^{(j)}_{i+1}$
$\displaystyle=\Psi_{i}(v_{i}^{(j)})+\xi^{(j)}_{i},$ (13a) $\displaystyle
v^{(j)}_{i+1}$
$\displaystyle=\hat{v}^{(j)}_{i+1}+K_{i+1}(y_{i+1}^{\dagger}-\hat{y}^{(j)}_{i+1}),$
(13b) $\displaystyle\hat{y}^{(j)}_{i+1}$
$\displaystyle=\mathbb{E}^{\mathsf{J}}[\mathfrak{h}_{i+1}(\hat{v}_{i+1})]+\eta^{(j)}_{{i+1}},$
(13c) $\displaystyle K_{i+1}$
$\displaystyle=(\hat{C}_{i+1}^{v\mathfrak{h}})^{\mathsf{J}}((\hat{C}_{i+1}^{\mathfrak{h}\mathfrak{h}})^{\mathsf{J}}+(\Gamma_{d})_{i+1})^{-1},$
(13d)

where $\xi^{(j)}_{i}\sim\mathcal{N}(0,(\Sigma_{d})_{i})$,
$\eta^{(j)}_{i}\sim\mathcal{N}(0,(\Gamma_{d})_{i})$,
$\mathfrak{h}_{i}(v)=\mathfrak{h}(v,t)$, and

$\displaystyle(\hat{C}_{i+1}^{v\mathfrak{h}})^{\mathsf{J}}$
$\displaystyle=\mathbb{E}^{\mathsf{J}}[(\hat{v}_{i+1}-\mathbb{E}^{\mathsf{J}}\hat{v}_{i+1})$
(14)
$\displaystyle\quad\quad\quad\otimes(\mathfrak{h}_{i+1}(\hat{v}_{i+1})-\mathbb{E}^{\mathsf{J}}[\mathfrak{h}_{i+1}(\hat{v}_{i+1})])],$
$\displaystyle(\hat{C}_{i+1}^{\mathfrak{h}\mathfrak{h}})^{\mathsf{J}}$
$\displaystyle=\mathbb{E}^{\mathsf{J}}[(\mathfrak{h}_{i+1}(\hat{v}_{i+1})-\mathbb{E}^{\mathsf{J}}[\mathfrak{h}_{i+1}(\hat{v}_{i+1})])$
(15)
$\displaystyle\quad\quad\quad\otimes(\mathfrak{h}_{i+1}(\hat{v}_{i+1})-\mathbb{E}^{\mathsf{J}}[\mathfrak{h}_{i+1}(\hat{v}_{i+1})])].$

Furthermore, we introduce the following rescalings adopted in Law, Stuart, and
Zygalakis (2015)Law, Stuart, and Zygalakis (2015):

$\displaystyle f(\cdot,t)$
$\displaystyle=(\Psi{{}_{i}}(\cdot)-I\cdot)/\tau,\quad
z_{i+1}^{\dagger}-z_{i}^{\dagger}=\tau y_{i+1}^{\dagger},$ (16)
$\displaystyle\Sigma(i\tau)$
$\displaystyle=(\Sigma_{d})_{i}/\tau\quad\Gamma(i\tau)=\tau(\Gamma_{d})_{i},\quad
i=t/\tau,$

Then, Eqs. (13) can be seen to be a discretization of Eqs. (12) with time step
$\tau$. More justification is given for these rescalings in Salgado,
Middleton, and Goodwin (1988)Salgado, Middleton, and Goodwin (1988) and Simon
(2006)Simon (2006). Note that both
$K{{}_{i+1}}=(\hat{C}_{i+1}^{v\mathfrak{h}})^{\mathsf{J}}\,{(\Gamma_{d})_{i+1}^{-1}}$
and Eq. (13d) are consistent with the continuous-time gain as $\tau\to 0$. We
use the latter, similar to the discrete-time Kalman filter.

###  II.5 Score Function Term

We now discuss further computational issues that arise when Eq. (10a) is
replaced by (II.2). This term involves the score function, defined as
$\nabla\log\rho$, but with an additional preconditioning. If this term is
added to the discrete-time particle version of the filter, Eq. (13b) becomes

$\displaystyle v^{(j)}_{i+1}$
$\displaystyle=\hat{v}^{(j)}_{i+1}+K_{i+1}(y_{i+1}^{\dagger}-\hat{y}^{(j)}_{i+1})$
(17)
$\displaystyle\quad\quad\quad\quad+K_{i+1}(\Gamma_{d})_{i+1}K_{i+1}^{T}(\nabla\log\rho_{i+1})^{\mathsf{J}}(\hat{v}^{(j)}_{i+1}),$

where $(\nabla\log\rho_{i+1})^{\mathsf{J}}$ denotes particle-based
approximation of the score function using
$\\{\hat{v}^{(j)}_{i+1}\\}_{j=1}^{\mathsf{J}}$. If we make the assumption that
the density is Gaussian with mean $\mathbb{E}v$ and covariance $C^{vv}$, the
score function takes on a simple form,

$\nabla\log\rho=-(C^{vv})^{-1}(v-\mathbb{E}v).$ (18)

A natural particle approximation $(\nabla\log\rho)^{\mathsf{J}}$ follows by
replacing the mean and covariance with the corresponding quantities computed
under the empirical measure of the set of particles.

More general kernel-based nonparametric estimators for the score function have
been developed, such as those defined in Zhou, Shi, and Zhu (2020)Zhou, Shi,
and Zhu (2020) and implemented in the kscore package. In the numerical
experiments reported in this paper, we either omit the score term completely,
or use it and employ only the Gaussian approximation.

###  II.6 Implementation

####  II.6.1 Ensemble Square-Root Formulation

In order to make the method scale well to high dimensions, an ensemble square-
root formulation (Tippett  _et al._ , 2003) of Eq. (13) can be used, although
we do not use it in the numerical experiments reported here. The advantage of
this formulation is that the most expensive linear algebra operations are
rewritten in the ensemble space, resulting in favorable computational
complexity when $J$ is much smaller than the state-space dimension $d$ or
observation-space dimension $p$.333Note, however, that in many applications
with a high-dimensional state space, the statistics of interest may be
relatively low-dimensional, such that the regular version of the algorithm
(13) will be feasible.

To implement this method, we write $(C^{vv})^{\mathsf{J}}=VV^{T}$,
$(C^{v\mathfrak{h}})^{\mathsf{J}}=VY^{T}$, and
$(C^{\mathfrak{h}\mathfrak{h}})^{\mathsf{J}}=YY^{T}$, where the $j$th column
of $V$ and the $j$th column of $Y$ are given by

$\displaystyle V^{(j)}$
$\displaystyle=(v^{(j)}-\mathbb{E}^{\mathsf{J}}v)/\sqrt{J-1},$ (19)
$\displaystyle Y^{(j)}$
$\displaystyle=(\mathfrak{h}(v^{(j)})-\mathbb{E}^{\mathsf{J}}\mathfrak{h})/\sqrt{J-1},$

respectively. Then, $K$ can be written as

$K=VY^{T}W,$ (20)

where
$W=(\Gamma_{d}^{-1}-\Gamma_{d}^{-1}Y(I+Y^{T}\Gamma_{d}^{-1}Y)^{-1}Y^{T}\Gamma_{d}^{-1})$
by the Woodbury identity.

We assume that $\Gamma_{d}^{-1}$ is provided and can be applied cheaply, for
example, if it is diagonal. This is a standard assumption (Tippett  _et al._ ,
2003). With this expression, $K$ can be computed in
$\mathcal{O}(J^{3}+J^{2}p+Jp^{2}+dJp)$.

Note that the Gaussian score function approximation Eq. (18) cannot be applied
in cases when $J<d$ since $(C^{vv})^{\mathsf{J}}$ will be singular. We do not
consider the score function term in the complexity analysis.

The complexity is, thus, a quadratic polynomial in $d$ and $p$, whereas
various ensemble square-root filters can be implemented to be linear in $p$
and $d$. The latter rely on the fact that the in the standard Kalman filter
the updated covariance can be written as $(I-KH)C^{vv}$, where $H$ is the
observation operator. The EnFPF cannot be written in this way. Whether the
EnFPF can be reformulated to be linear in $p$ and $d$ by another approach is a
topic for future research.

####  II.6.2 Code

The open-source Julia code for the EnFPF is available at
https://github.com/eviatarbach/EnFPF. In the numerical experiments that
follow, we compute the Wasserstein distance (explained in Sec. III) using the
Python Optimal Transport library (Flamary  _et al._ , 2021). We used the
parasweep library to facilitate parallel experiments (Bach, 2021).

####  II.6.3 Numerical Methods for the Test Models

In Sec. III, we will present numerical experiments with the Lorenz63,
Lorenz96, and Kuramoto–Sivashinsky models. We integrate the Lorenz63 and
Lorenz96 models using the fourth-order Runge–Kutta method, with a time step of
0.05 for both. We integrate the Kuramoto–Sivashinsky equation in Fourier space
using the exponential time differencing fourth-order Runge–Kutta method
(Kassam and Trefethen, 2005) with 64 Fourier modes and a time step of 0.25.

##  III Numerical Experiments

In this section, we present the results of numerical experiments applying the
discrete-time EnFPF of Sec. II.4 to the Lorenz63, Lorenz96, and
Kuramoto–Sivashinsky systems. The first three subsections are devoted,
respectively, to these three models; a final fourth subsection returns to the
Lorenz63 model, with quasiperiodic forcing.

We found in the experiments that assimilating too often can cause degraded
results for some systems, as opposed to the situation in standard filtering,
where increased assimilation frequency is typically preferred. In the standard
filtering problem, there is a single true trajectory, and under certain
conditions, the filtering distribution will converge to this trajectory in the
limit of zero observational noiseDashti  _et al._ (2013); Law, Stuart, and
Zygalakis (2015). In the non-zero noise case, however, the filtered time-
series (e.g., the maximum a posteriori estimate) will not even generally be a
trajectory of the dynamical system, except in methods such as strong-
constraint 4DVar. Here, we expect that the problem of ensemble members
deviating from being trajectories can be amplified, since the method only aims
to match statistical features of the entire ensemble. Thus, if assimilation is
done too frequently, then ensemble members may be pushed too far from being
trajectories into unphysical or unstable parts of the phase space. In fact, we
found the assimilation frequency to be a key tuning parameter. We refer to a
single forecast–assimilation step (Eqs. (13)) as a _cycle_ , as is common in
the DA literature, and each cycle lasts for $\tau$ time units.

We found, furthermore, that the score term did not consistently improve
filtering performance. In the experiments that follow, we omit the score term
except in the experiments with the Kuramoto–Sivashinsky system in Sec. III.3,
where it leads to clear improvements when used, together with the Gaussian
approximation, in the form of Eq. (18). For both the Lorenz models we found
that the inclusion of the Gaussian approximation of the score degraded
performance and the use of the kernel-based score approximations, based on the
paper of Zhou, Shi, and Zhu (2020)Zhou, Shi, and Zhu (2020), was no better
than simply omitting the term altogether.

In the experiments below, we use a Wasserstein metric to quantify the distance
between the ensemble distribution and the invariant density. We estimate the
invariant density using an ensemble integrated for a sufficiently long time.
We employ the $W_{1}$ Wasserstein metric that allows us to compute distances
between empirical distributions. The code for computing this distance is
readily available (see Sec. II.6.2).

###  III.1 Lorenz63 Model

For the experiments in this subsection, we use the Lorenz (1963)Lorenz (1963)
model

$\displaystyle\frac{dx}{dt}$ $\displaystyle=\sigma(y-x),$ (21)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=x(r-z)-y,$
$\displaystyle\frac{dz}{dt}$ $\displaystyle=xy-\beta z,$

with the standard parameter values $\sigma=10$, $r=28$, and $\beta=8/3$.

####  III.1.1 Assimilating Time-Varying Means and Second Moments

We first verify the ability of the EnFPF to force an ensemble to adopt time-
varying statistics. We do this by applying the EnFPF to a 10-member ensemble,
with noisy statistical observations of the means and uncentered second moments
of the three variables coming from a 100-member ensemble being evolved
concurrently. The difference between the statistics computed over the 10- and
100-member ensembles arise due to both sampling errors and different initial
conditions. The 100-member ensemble (despite having its own sampling error)
better approximates the true statistics of the system, and we view these
100-member ensemble statistics as the truth, based on which we may compute
errors in the statistics of 10-member ensembles. We assimilate observations
every 0.2 time units, with an observation error covariance set to 20% of the
time variability of each statistic computed over the 100-member ensemble.

Figure 2: The impact of filtering on the root-mean-square error (RMSE) in the
mean and second moment in the Lorenz63 model.

Figure 2 shows the resulting error in the means and second moments of the
10-member ensemble, compared with the errors arising from an unfiltered run of
the 10-member ensemble; in both cases the errors are computed by comparison
with the 100-member ensemble. After several cycles, the filter appears to
reach an asymptotic error on the order of the observation error, and this
error is significantly lower than that arising in the unfiltered case.

Means |
---|---
Observation error | Filtered RMSE
10% (0.088) | 0.11
35% (0.31) | 0.40
60% (0.53) | 0.69
85% (0.75) | 0.97
Second moments |
Observation error | Filtered RMSE
10% (2.8) | 20
35% (9.9) | 23
60% (17) | 29
85% (24) | 35
Table 1: The impact of the observation error covariance on filtering
performance. In the first column are the percentages of the standard deviation
of the time variability of each statistic taken to be the observation error,
and in parentheses are the square root of the total variance of the
observation error in the statistic. With no filtering, the RMSE is 2.5 in the
unfiltered means and 73 in the second moments. The RMSE is averaged over 1400
cycles after 100 transient cycles.

Table 1 shows the impact of the observation error covariance magnitude on the
filtering performance. The ‘setup is otherwise the same as that described
above. As expected, the error increases as $\Gamma$ is increased, although
still outperforming the unfiltered ensemble.

####  III.1.2 Accelerating Convergence to the Invariant Density

We now test the ability of the EnFPF to accelerate convergence to the
invariant density. We assimilate observations of fixed statistics of the
invariant density, the means and second moments of the three variables, into a
100-member ensemble. We use the same assimilation frequency and observation
error as in Sec. III.1.1.

Figure 3: The estimated Wasserstein distance to the invariant density in
Lorenz63, in unfiltered and filtered cases. For the filtered case, the first
and second moments are assimilated. Each curve is averaged over 10 different
initializations.

Figure 3 shows the impact of the EnFPF on the convergence to the invariant
density. In this case, we only apply the EnFPF for the first 30 cycles
(indicated by the pink rectangle), and then let the ensemble evolve under the
regular Lorenz63 dynamics. We see that the EnFPF leads to a more rapid
convergence: by the end of the filtering period, the distance is close to the
asymptotic one, while it takes at least 100 cycles for the unfiltered case to
reach the same. Figure 4 visualizes in state space this rapid convergence
toward the invariant density via the EnFPF.

Figure 4: Top panel: an ensemble evolving in time from left to right,
superimposed on the invariant density of Lorenz63 in the $x$–$z$ plane. Orange
corresponds to higher probability density and blue to lower. Bottom panel: the
same but with the EnFPF applied.

####  III.1.3 Impact of Higher-Order Moments

Figure 5: The estimated Wasserstein distance to the invariant density in
Lorenz63, in unfiltered and filtered cases when different moments are
assimilated. The curves are averaged over 25 initial conditions, and the
shaded areas correspond to $\pm$ the standard error over the initializations.
Here, for the filtered cases, the EnFPF is applied at every cycle.

Figure 5 shows the convergence to the invariant measure of Lorenz63 with
different assimilated moments of $x$ and $y$, namely the first, first and
second, and first, second, and third marginal moments. Assimilating the first-
order moments accelerates the convergence to the invariant measure compared to
the unfiltered case. Adding the second- and third-order moments appears to
result in the most rapid initial rate of convergence, and after about 50
cycles assimilating the first two and the first three moments leads to a
similar asymptotic distance to the invariant measure.

###  III.2 Lorenz96 Model

We now test the convergence to the invariant density of the Lorenz
(1996)Lorenz (1996) model

$\frac{\text{d}x_{i}}{\text{d}t}=-x_{i-1}(x_{i-2}+x_{i+1})-x_{i}+F,$ (22)

where the indices $i$ range from 1 to $D$ and are cyclical. We use $F=8$ and
$D=40$ variables. This is a model of an atmospheric latitude circle that is
commonly used in data assimilation experiments.

Figure 6: The estimated Wasserstein distance to the invariant density in
Lorenz96, in unfiltered and filtered cases. For the filtered case, the first
and second moments are assimilated. Here, we show the mean of the Wasserstein
distances corresponding to the marginal density for each variable.

We assimilate the means and second moments of the 40 variables on the
invariant density, with an observation error covariance of 20% of the temporal
variability of the statistics computed over a 100-member ensemble. We
assimilate every 0.05 time units into a 100-member ensemble for 40 cycles.
Figure 6 shows that the convergence toward the invariant density is thereby
significantly accelerated.

###  III.3 Kuramoto–Sivashinsky Model

We now carry out experiments with the Kuramoto–Sivashinsky model, a chaotic
partial differential equation in one spatial dimension,

$u_{t}+u_{xxxx}+u_{xx}+uu_{x}=0,\quad x\in[0,L].$ (23)

We use $L=22$ and periodic boundary conditions, discretized using 64 Fourier
modes (see Sec. II.6.3 for details on the numerical method).

Figure 7: The estimated Wasserstein distance to the invariant density in the
Kuramoto–Sivashinsky equation, in unfiltered and filtered cases. For the
filtered case, the first and second moments are assimilated. Here, we show the
mean of the Wasserstein distances corresponding to the marginal density for
each variable.

We assimilate the means and second moments of the invariant density of the 64
variables in physical space, every 2.0 time units. We assimilate for 30 cycles
using a 100-member ensemble, and again use an observational error covariance
of 20% of the temporal variability. Figure 7 shows the results with and
without the score term included. In both cases, there is an acceleration
compared to the unfiltered case; inclusion of the score term considerably
accelerates convergence.

###  III.4 Time-Dependent Invariant Measures

We now use the Lorenz63 model (Eq. (21)) but with the $r$ parameter subject to
quasiperiodic forcing, as in Daron and Stainforth (2015)Daron and Stainforth
(2015),

$r(t)=28+\sin\left(2\pi
t\right)+\sin\\!\left(\\!\sqrt{3}t\right)+\sin\\!\left(\\!\sqrt{17}t\right).$
(24)

Since this system is non-autonomous, it possesses for each time $t$ a pullback
attractor with a corresponding time-dependent invariant measure, as discussed
in Sec. I.2. The measure at time $t$ can be approximated by the empirical
density at time $t$ of an ensemble initialized sufficiently far back in time,
at $t-T$ for some large $T$. Here, we evolve a 100-member ensemble using
$T=500$ time units to approximate the invariant measures at time $t$. Then, we
evolve the ensemble for the additional time period of $t$ to $t+20$ to obtain
approximations to the invariant measures in this period.

We evolve two separate 100-member ensembles for the same time period $t$ to
$t+20$, but with $T=0$ (no spin-up). We apply the EnFPF to one of these
ensembles and not the other. For the EnFPF, we assimilate every 0.05 time
units with an observation error covariance of 20% of the temporal variability.
We then measure the distance between the empirical densities of these two
ensembles and the one approximating the invariant measure described in the
previous paragraph.

Figure 8: The estimated Wasserstein distance to the invariant density in the
non-autonomous Lorenz63 model, in unfiltered and filtered cases. For the
filtered case, $\mathbb{E}[x^{i}]$, $\mathbb{E}[y^{i}]$, and
$\mathbb{E}[z^{i}]$ for $i=1,2,3$ are assimilated.

Figure 8 shows that the convergence to the time-dependent invariant measures
is indeed accelerated by the EnFPF, reaching a comparable asymptotic distance
to the invariant measure in less than half the time.

##  IV Justification of Algorithm

###  IV.1 Kalman–Bucy (KB) Filter for Densities

Since both the Fokker–Planck equation (5) and the observation equation (8) are
linear, and since all noise is additive Gaussian, the conditional probability
measure over densities, $\rho|Z^{\dagger}(t)$, is a Gaussian. This filtering
problem can be solved using a Kalman–Bucy filter in Hilbert space, posing
significant challenges because it involves finding a sequence of probability
measures on an infinite-dimensional space of functions (densities).

We start by defining the Hilbert space
$\mathcal{H}=L^{2}(\mathbb{R}^{d},\mathbb{R})$ with inner product

$\langle a,b\rangle_{\mathcal{H}}\equiv\int ab\,dv.$ (25)

We consider density functions $\rho\in\mathcal{H}$, and we require that
$\rho(v,t)\to 0$ as $v\to\infty$. Note that we will sometimes use this inner
product in situations where one of the arguments is only locally square
integrable; in particular we will need to use the constant function
$\mathbbm{1}(v)=1$. To distinguish them from the Hilbert space inner product,
we denote the standard Euclidean inner product in $\mathbb{R}^{p}$ as
$\langle\cdot,\cdot\rangle_{\mathbb{R}^{p}}$ and the weighted Euclidean inner
product, defined for any strictly positive-definite and symmetric
$A\in\mathbb{R}^{p\times p}$, as $\langle\cdot,\cdot\rangle_{A}\equiv\langle
A^{-1/2}\cdot,A^{-1/2}\cdot\rangle_{\mathbb{R}^{p}}$.

Recall definition Eq. ((5)b) of the adjoint of the generator $\mathcal{L}$. We
are given the dynamics and observation equations (5) and (8):

$\displaystyle d\rho^{\dagger}(v,t)$
$\displaystyle=\mathcal{L}^{*}(t)\rho^{\dagger}(v,t)\,dt,$ (26) $\displaystyle
dz^{\dagger}(t)$
$\displaystyle=H(t)\rho^{\dagger}(v,t)\,dt+\sqrt{\Gamma(t)}dB.$ (27)

Then, given all observations up to time $t$,
$Z^{\dagger}(t)=\\{z^{\dagger}(s)\\}_{s\in[0,t]}$, the filtering distribution
is given by

$\rho(\cdot,t)|Z^{\dagger}(t)\sim\mu(t)\equiv\mathcal{N}\bigl{(}m(t),C(t)\bigr{)},$
(28)

where $\mathcal{N}$ is a Gaussian measure on $\mathcal{H}$ with mean $m(t)$
and covariance operator $C(t)$. For notational simplicity, we have dropped the
explicit dependence of $m(t)$, $C(t)$, and $\rho(t)$ on $v$. Here $C\in
L(\mathcal{H},\mathcal{H})$ is necessarily self-adjoint and trace class
Bogachev (1998); that is, $\operatorname{tr}(C)<\infty$. In what follows, the
expectation $\mathbb{E}_{\mu}$ is defined with respect to the measure $\mu$ on
the space of $L^{2}$ densities $\rho$.

Using Theorem 7.10 in Falb (1967)Falb (1967), the KB filter for this system
can be written as

$\displaystyle dm(t)$ $\displaystyle=\mathcal{L}^{*}(t)m(t)\,dt$ (29a)
$\displaystyle\quad\quad\quad+C(t)H^{*}(t)\Gamma(t)^{-1}\bigl{(}dz^{\dagger}(t)-H(t)m(t)\bigr{)}\,dt,$
$\displaystyle dC(t)$
$\displaystyle=\mathcal{L}^{*}(t)C(t)\,dt+C(t)\mathcal{L}(t)\,dt$ (29b)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad-C(t)H^{*}(t)\Gamma(t)^{-1}H(t)C(t)\,dt,$
$\displaystyle m(0)$ $\displaystyle=m_{0},C(0)=C_{0},$ (29c)

where

$C(t)=\operatorname{cov}({\rho(t)}-m(t),{\rho(t)}-m(t)),$ (30)

and

$\operatorname{cov}(x,y)\equiv\mathbb{E}_{\mu}[x\otimes
y]-\mathbb{E}_{\mu}[x]\otimes\mathbb{E}_{\mu}[y].$ (31)

The outer-product $x_{1}\otimes y_{1}$ is defined by the identity

$(x_{1}\otimes y_{1})x=x_{1}\langle y_{1},x\rangle_{\mathcal{H}}$ (32)

holding for all $x\in\mathcal{H}$. Note that Falb (1967)Falb (1967) requires
boundedness of $\mathcal{L}^{*}$, but the results have been extended to
unbounded operators (Curtain, 1975). However, we still require boundedness of
$H$. For the rest of the paper, we will assume that $H$ takes the form in Eq.
(9).

The adjoint operator $H^{*}$ is then given by

$H^{*}(t)u=\langle\mathfrak{h}(v,t),u\rangle_{\mathbb{R}^{p}},$ (33)

for $u\in\mathbb{R}^{p}$. Note that, formally, $H^{*}(t)u$ is to be viewed as
a function of $v$, in the space $\mathcal{H}.$

In general, the solution of Eq. (29a), $m(t)$, will not be normalized.
However, in Appendix A we show that normalization is preserved under certain
conditions on the initializations $m_{0}$ and $C_{0}$ from Eq. (29c). However,
$m(t)$ is not guaranteed to be non-negative for all $v$ and $t$ and, thus,
cannot be a proper probability density. Nonetheless, we can still consider
integrals against it.

###  IV.2 Ansatz and Relation to KB Filter for Densities

Solving the KB filter equations directly is intractable. We therefore seek an
equation that is amenable to a mean-field model, which, in turn, can be
approximated by a particle system. We propose the following ansatz for the
density of $v|Z^{\dagger}(t)$:

$\frac{\partial\rho}{\partial
t}=\mathcal{L}^{*}{(t)}\rho+\Bigl{\langle}\mathfrak{h}(v,{t})-H{(t)}\rho,\frac{dz^{\dagger}}{dt}-H{(t)}\rho\Bigr{\rangle}_{\Gamma{(t)}}\rho.$
(34)

Note the similarity to the Kushner–Stratonovich (KS) equation (7). Although
the solutions to this equation do not match the KB filter for densities in
general, we show in Theorem 1 that they coincide in observation space for
linear $f$ and $\mathfrak{h}$, under additional assumptions detailed there.
The proof sketch is provided in Appendix B.

###  IV.3 Mean-Field Approximation

We would now like to find a mean-field model which has, as its FP equation,
Eq. (34). We postulate the following form:

$\displaystyle dv$
$\displaystyle=f(v,{t})\,dt+\sqrt{\Sigma{(t)}}\,dW+a(v,\rho,{t})\,dt$ (35)
$\displaystyle\qquad+K(v,\rho,{t})\Bigl{(}dz^{\dagger}-H{(t)}\rho{(v,t)}\,dt-\sqrt{\Gamma{(t)}}dB\Bigr{)}.$

Specifically, we aim to choose the pair of functions $(a,K)$ so that the
Fokker–Planck equation for $v$ governed by this mean-field model coincides
with Eq. (34). In Appendix C, we detail the choices which achieve this, and
after making a further approximation of $K$, we obtain equations (10) with
(10a) replaced by (II.2). However, as explained there, in many cases use of
Eq. (10a), which corresponds to setting $a\equiv 0$ and using a simple
approximation of $K$, leads to algorithms that empirically perform well.

##  V Conclusions

In this paper we introduce the Fokker–Planck filtering problem, which consists
of estimating the evolving probability density of a (possibly stochastic)
dynamical system given noisy observations of expectations evaluated with
respect to it. We provide a solution to this problem using the KB filter in
Hilbert space, and introduce an ensemble algorithm, the ensemble Fokker–Planck
filter (EnFPF), that approximates it under conditions on the dynamics and
observables. We also show, through numerical experiments, that this method can
be used to accelerate convergence to the invariant measure of dynamical
systems, and that this acceleration phenomenon applies beyond the conditions
on the dynamics and observables required to provably link the KB filter and
the mean-field model underlying our proposed ensemble method.

Future work will test this method on higher-dimensional models, such as
turbulent channel flows and ocean models. Other future directions, as
described in Sec. I.2, include (i) the testing of this method as an approach
to counteract model error, (ii) use in parameter estimation, and (iii) use in
the acceleration of sampling methods such as Langevin dynamics and Markov
chain Monte Carlo when some statistics of the target density are known.
Furthermore, many of the numerical results require deeper understanding; these
include the impact of the assimilation frequency, the score term, and the
incorporation of higher-order moments, or other observables, on the filter
performance. Finally, on the theoretical side, there is a considerable need
for deeper analysis.

## Acknowledgments

EB is supported by the the Foster and Coco Stanback Postdoctoral Fellowship.
AS is supported by the Office of Naval Research (ONR) through grant
N00014-17-1-2079. TC and AS acknowledge recent support through ONR grant
N00014-23-1-2654. EB and AS are also grateful for support from the Department
of Defense Vannevar Bush Faculty Fellowship held by AS. We thank Tapio
Schneider, Dimitris Giannakis, and two anonymous referees for helpful
comments.

## References

  * Jazwinski (1970) A. H. Jazwinski,  _Stochastic Processes and Filtering Theory_ , Mathematics in Science and Engineering No. 64 (Academic Press, Inc., New York, 1970).
  * Kalnay (2002) E. Kalnay,  _Atmospheric Modeling, Data Assimilation and Predictability_ (Cambridge University Press, New York, 2002).
  * Law, Stuart, and Zygalakis (2015) K. Law, A. Stuart, and K. Zygalakis,  _Data Assimilation: A Mathematical Introduction_, Texts in Applied Mathematics No. 62 (Springer International Publishing, 2015).
  * Reich and Cotter (2015) S. Reich and C. Cotter,  _Probabilistic Forecasting and Bayesian Data Assimilation_ (Cambridge University Press, Cambridge, 2015).
  * Gurtin (1981) M. E. Gurtin,  _An Introduction to Continuum Mechanics_ , Mathematics in Science and Engineering, Vol. 158 (Academic Press, 1981).
  * Gonzalez and Stuart (2008) O. Gonzalez and A. M. Stuart,  _A First Course in Continuum Mechanics_, Cambridge Texts in Applied Mathematics (Cambridge University Press, Cambridge, 2008).
  * Evensen (1994) G. Evensen, “Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics,” Journal of Geophysical Research: Oceans 99, 10143–10162 (1994).
  * Crişan, Del Moral, and Lyons (1999) D. Crişan, P. Del Moral, and T. J. Lyons, “Interacting Particle Systems Approximations of the Kushner-Stratonovitch Equation,” Advances in Applied Probability 31, 819–838 (1999), 1428398 .
  * Calvello, Reich, and Stuart (2022) E. Calvello, S. Reich, and A. M. Stuart, “Ensemble Kalman Methods: A Mean Field Perspective,” (2022), arxiv:2209.11371 [cs, math] .
  * Curtain and Pritchard (1978) R. F. Curtain and A. J. Pritchard,  _Infinite Dimensional Linear Systems Theory_, Lecture Notes in Control and Information Sciences, Vol. 8 (Springer-Verlag, Berlin/Heidelberg, 1978).
  * Yang, Mehta, and Meyn (2013) T. Yang, P. G. Mehta, and S. P. Meyn, “Feedback Particle Filter,” IEEE Transactions on Automatic Control 58, 2465–2480 (2013).
  * Reich (2019) S. Reich, “Data assimilation: The Schrödinger perspective,” Acta Numerica 28, 635–711 (2019).
  * Crisan and Xiong (2010) D. Crisan and J. Xiong, “Approximate McKean–Vlasov representations for a class of SPDEs,” Stochastics 82, 53–68 (2010).
  * Pathiraja, Reich, and Stannat (2021) S. Pathiraja, S. Reich, and W. Stannat, “McKean–Vlasov SDEs in Nonlinear Filtering,” SIAM Journal on Control and Optimization 59, 4188–4215 (2021).
  * Abramov (2010) R. V. Abramov, “The multidimensional maximum entropy moment problem: A review of numerical methods,” Communications in Mathematical Sciences 8, 377–392 (2010).
  * Pavliotis and Stuart (2008) G. A. Pavliotis and A. M. Stuart,  _Multiscale Methods: Averaging and Homogenization_, Texts in Applied Mathematics, Vol. 53 (Springer, New York, NY, 2008).
  * Drótos, Bódai, and Tél (2016) G. Drótos, T. Bódai, and T. Tél, “Quantifying nonergodicity in nonautonomous dissipative dynamical systems: An application to climate change,” Physical Review E 94, 022214 (2016).
  * Dirren and Hakim (2005) S. Dirren and G. J. Hakim, “Toward the assimilation of time-averaged observations,” Geophysical Research Letters 32 (2005), 10.1029/2004GL021444.
  * Rey  _et al._ (2014) D. Rey, M. Eldridge, U. Morone, H. D. I. Abarbanel, U. Parlitz, and J. Schumann-Bischoff, “Using waveform information in nonlinear data assimilation,” Physical Review E 90, 062916 (2014).
  * Bach  _et al._ (2021) E. Bach, S. Mote, V. Krishnamurthy, A. S. Sharma, M. Ghil, and E. Kalnay, “Ensemble Oscillation Correction (EnOC): Leveraging Oscillatory Modes to Improve Forecasts of Chaotic Systems,” Journal of Climate 34, 5673–5686 (2021).
  * Bryan (1984) K. Bryan, “Accelerating the Convergence to Equilibrium of Ocean-Climate Models,” Journal of Physical Oceanography 14, 666–673 (1984).
  * Daron and Stainforth (2015) J. D. Daron and D. A. Stainforth, “On quantifying the climate of the nonautonomous Lorenz-63 model,” Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 043103 (2015).
  * Drótos, Bódai, and Tél (2017) G. Drótos, T. Bódai, and T. Tél, “On the importance of the convergence to climate attractors,” The European Physical Journal Special Topics 226, 2031–2038 (2017).
  * Deser  _et al._ (2020) C. Deser, F. Lehner, K. B. Rodgers, T. Ault, T. L. Delworth, P. N. DiNezio, A. Fiore, C. Frankignoul, J. C. Fyfe, D. E. Horton, J. E. Kay, R. Knutti, N. S. Lovenduski, J. Marotzke, K. A. McKinnon, S. Minobe, J. Randerson, J. A. Screen, I. R. Simpson, and M. Ting, “Insights from Earth system model initial-condition large ensembles and future prospects,” Nature Climate Change 10, 277–286 (2020).
  * Wang, Xin, and Zhang (2022) Z. Wang, J. Xin, and Z. Zhang, “DeepParticle: Learning invariant measure by a deep neural network minimizing Wasserstein distance on data generated from an interacting particle method,” Journal of Computational Physics 464, 111309 (2022).
  * Hwang, Hwang-Ma, and Sheu (2005) C.-R. Hwang, S.-Y. Hwang-Ma, and S.-J. Sheu, “Accelerating diffusions,” The Annals of Applied Probability 15, 1433–1444 (2005).
  * Abdulle, Pavliotis, and Vilmart (2019) A. Abdulle, G. A. Pavliotis, and G. Vilmart, “Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations,” Comptes Rendus Mathematique 357, 349–354 (2019).
  * Nelson and Fringer (2017) K. S. Nelson and O. B. Fringer, “Reducing spin-up time for simulations of turbulent channel flow,” Physics of Fluids 29, 105101 (2017).
  * Goldys and Maslowski (2005) B. Goldys and B. Maslowski, “Exponential ergodicity for stochastic reaction–diffusion equations,” in  _Stochastic Partial Differential Equations and Applications - VII_ , edited by G. D. Prato and L. Tubaro (Chapman and Hall/CRC, 2005).
  * Isik (2013) O. R. Isik, “Spin up problem and accelerating convergence to steady state,” Applied Mathematical Modelling 37, 3242–3253 (2013).
  * Isik, Takhirov, and Zheng (2017) O. R. Isik, A. Takhirov, and H. Zheng, “Second order time relaxation model for accelerating convergence to steady-state equilibrium for Navier–Stokes equations,” Applied Numerical Mathematics 119, 67–78 (2017).
  * Arnold (1998) L. Arnold,  _Random Dynamical Systems_, Springer Monographs in Mathematics (Springer, Berlin, Heidelberg, 1998).
  * Chekroun, Simonnet, and Ghil (2011) M. D. Chekroun, E. Simonnet, and M. Ghil, “Stochastic climate dynamics: Random attractors and time-dependent invariant measures,” Physica D: Nonlinear Phenomena 240, 1685–1700 (2011).
  * Ghil, Chekroun, and Simonnet (2008) M. Ghil, M. D. Chekroun, and E. Simonnet, “Climate dynamics and fluid mechanics: Natural variability and related uncertainties,” Physica D: Nonlinear Phenomena Euler Equations: 250 Years On, 237, 2111–2126 (2008).
  * Annunziato and Borzì (2013) M. Annunziato and A. Borzì, “A Fokker–Planck control framework for multidimensional stochastic processes,” Journal of Computational and Applied Mathematics 237, 487–507 (2013).
  * Covington, Qi, and Chen (2023) J. Covington, D. Qi, and N. Chen, “Effective Statistical Control Strategies for Complex Turbulent Dynamical Systems,”  (2023), arxiv:2307.15637 [math-ph] .
  * Köhl and Willebrand (2002) A. Köhl and J. Willebrand, “An adjoint method for the assimilation of statistical characteristics into eddy-resolving ocean models,” Tellus A: Dynamic Meteorology and Oceanography 54, 406–425 (2002).
  * Schneider, Stuart, and Wu (2021) T. Schneider, A. M. Stuart, and J.-L. Wu, “Learning stochastic closures using ensemble Kalman inversion,” Transactions of Mathematics and Its Applications 5, tnab003 (2021).
  * Mons, Du, and Zaki (2021) V. Mons, Y. Du, and T. A. Zaki, “Ensemble-variational assimilation of statistical data in large-eddy simulation,” Physical Review Fluids 6, 104607 (2021).
  * Schneider, Stuart, and Wu (2022) T. Schneider, A. M. Stuart, and J.-L. Wu, “Ensemble Kalman inversion for sparse learning of dynamical systems from time-averaged data,” Journal of Computational Physics 470, 111559 (2022).
  * Levine and Stuart (2022) M. Levine and A. Stuart, “A framework for machine learning of model error in dynamical systems,” Communications of the American Mathematical Society 2, 283–344 (2022).
  * Bach and Ghil (2023) E. Bach and M. Ghil, “A Multi-Model Ensemble Kalman Filter for Data Assimilation and Forecasting,” Journal of Advances in Modeling Earth Systems 15, e2022MS003123 (2023).
  * Danforth and Kalnay (2008) C. M. Danforth and E. Kalnay, “Using Singular Value Decomposition to Parameterize State-Dependent Model Errors,” Journal of the Atmospheric Sciences 65, 1467–1478 (2008).
  * Chen  _et al._ (2022) T.-C. Chen, S. G. Penny, J. S. Whitaker, S. Frolov, R. Pincus, and S. Tulich, “Correcting Systematic and State-Dependent Errors in the NOAA FV3-GFS Using Neural Networks,” Journal of Advances in Modeling Earth Systems 14, e2022MS003309 (2022).
  * Farchi  _et al._ (2021) A. Farchi, P. Laloyaux, M. Bonavita, and M. Bocquet, “Using machine learning to correct model error in data assimilation and forecast applications,” Quarterly Journal of the Royal Meteorological Society 147, 3067–3084 (2021).
  * Ephrati  _et al._ (2023) S. Ephrati, A. Franken, E. Luesink, and B. Geurts, “Data-assimilation closure for large-eddy simulation of quasi-geostrophic flow on the sphere,”  (2023), arxiv:2312.12858 [physics] .
  * Bain and Crisan (2009) A. Bain and D. Crisan,  _Fundamentals of Stochastic Filtering_, Stochastic Modelling and Applied Probability (Springer, New York, NY, 2009).
  * Salgado, Middleton, and Goodwin (1988) M. Salgado, R. Middleton, and G. C. Goodwin, “Connection between continuous and discrete Riccati equations with applications to Kalman filtering,” IEE Proceedings D (Control Theory and Applications) 135, 28–34 (1988).
  * Simon (2006) D. Simon,  _Optimal State Estimation: Kalman, H ∞, and Nonlinear Approaches_ (John Wiley & Sons, Ltd, 2006).
  * Zhou, Shi, and Zhu (2020) Y. Zhou, J. Shi, and J. Zhu, “Nonparametric Score Estimators,” in  _Proceedings of the 37th International Conference on Machine Learning_ (PMLR, 2020) pp. 11513–11522.
  * Tippett  _et al._ (2003) M. K. Tippett, J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker, “Ensemble Square Root Filters,” Monthly Weather Review 131, 1485–1490 (2003).
  * Flamary  _et al._ (2021) R. Flamary, N. Courty, A. Gramfort, M. Z. Alaya, A. Boisbunon, S. Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, L. Gautheron, N. T. H. Gayraud, H. Janati, A. Rakotomamonjy, I. Redko, A. Rolet, A. Schutz, V. Seguy, D. J. Sutherland, R. Tavenard, A. Tong, and T. Vayer, “POT: Python Optimal Transport,” Journal of Machine Learning Research 22, 1–8 (2021).
  * Bach (2021) E. Bach, “parasweep: A template-based utility for generating, dispatching, and post-processing of parameter sweeps,” SoftwareX 13, 100631 (2021).
  * Kassam and Trefethen (2005) A.-K. Kassam and L. N. Trefethen, “Fourth-Order Time-Stepping for Stiff PDEs,” SIAM Journal on Scientific Computing 26, 1214–1233 (2005).
  * Dashti  _et al._ (2013) M. Dashti, K. J. H. Law, A. M. Stuart, and J. Voss, “MAP estimators and their consistency in Bayesian nonparametric inverse problems,” Inverse Problems 29, 095017 (2013).
  * Lorenz (1963) E. N. Lorenz, “Deterministic Nonperiodic Flow,” Journal of the Atmospheric Sciences 20, 130–141 (1963).
  * Lorenz (1996) E. N. Lorenz, “Predictability: A problem partly solved,” in  _Proceedings of a Seminar Held at ECMWF on Predictability, 4-8 September 1995_ , Vol. 1 (ECMWF, Shinfield Park, Reading, 1996) pp. 1–18.
  * Bogachev (1998) V. I. Bogachev,  _Gaussian Measures_ , Mathematical Surveys and Monographs No. 62 (American Mathematical Society, 1998).
  * Falb (1967) P. L. Falb, “Infinite-dimensional filtering: The Kalman–Bucy filter in Hilbert space,” Information and Control 11, 102–137 (1967).
  * Curtain (1975) R. Curtain, “A Survey of Infinite-Dimensional Filtering,” SIAM Review 17, 395–411 (1975).

##  Appendix A Properties of the KB Filter for Densities

Lemma 1 and Remark 2 below give the conditions under which $m(t)$ and
$\rho(t)\sim\mathcal{N}(m(t),C(t))$ will be normalized. The function
$\mathbbm{1}$ is defined as $\mathbbm{1}(v)\equiv 1$ for all $v$.

######  Lemma 1.

Assume that $\rho(0)\sim\mu(0)=\mathcal{N}(m_{0},C_{0})$ with

$\begin{cases}\langle m_{0},\mathbbm{1}\rangle_{\mathcal{H}}=1,\\\
C_{0}\mathbbm{1}=0.\end{cases}$ (36)

Then, for $m(t)$ and $C(t)$ satisfying equations (29a)–(29c),

  1. (a)

$C(t)\mathbbm{1}=0$ for all $t\geq 0$, and

  2. (b)

$\langle m(t),\mathbbm{1}\rangle_{\mathcal{H}}=1$ for all $t\geq 0$.

###### Proof.

(Sketch)

  * (a)

Since $\mathcal{L}\mathbbm{1}=0$, we have

$\frac{d}{dt}(C\mathbbm{1})=\mathcal{L}^{*}C\mathbbm{1}-CH^{*}\Gamma^{-1}HC\mathbbm{1}.$
(37)

Assuming uniqueness of the solution to Eq. (29b) for the evolution of $C(t)$,
we deduce that $C(t)\mathbbm{1}=0$ solves Eq. (37).

  * (b)

Applying Itô’s lemma to $\langle m,\mathbbm{1}\rangle_{\mathcal{H}}$ (the Itô
correction does not appear due to linearity of the inner product),

$\displaystyle\frac{d}{dt}\langle m,\mathbbm{1}\rangle_{\mathcal{H}}$
$\displaystyle=\langle\mathcal{L}^{*}m,\mathbbm{1}\rangle_{\mathcal{H}}+\left\langle
CH^{*}\Gamma^{-1}(dz^{\dagger}-Hm),\mathbbm{1}\right\rangle_{\mathcal{H}},$
(38) $\displaystyle=\langle
m,\mathcal{L}\mathbbm{1}\rangle_{\mathcal{H}}+\left\langle
H^{*}\Gamma^{-1}(dz^{\dagger}-Hm),C\mathbbm{1}\right\rangle_{\mathcal{H}},$
$\displaystyle=0,$

since $\mathcal{L}\mathbbm{1}=0$, $C$ is self-adjoint by construction, and
$C\mathbbm{1}=0$ by (a). Now assuming uniqueness of the equation (29a) for
$m(t)$ we find that $\langle m(t),\mathbbm{1}\rangle_{\mathcal{H}}=1$ solves
Eq. (38).

∎

######  Remark 2.

If the conditions in Eq. (36) hold then $\langle\rho(t),\mathbbm{1}\rangle=1$
for $t\geq 0$ almost surely, where $\rho(t)\sim\mu(t)=\mathcal{N}(m(t),C(t))$.
This is because $\mathbbm{1}$ is in the null-space of both the symmetric
operator square-root of $C(t)$, $\sqrt{C(t)}$, and

$\rho(t)=m(t)+\sqrt{C(t)}\xi,$ (39)

where $\xi\sim\mathcal{N}(0,\mathbb{I})$, with $\mathbb{I}$ being the
identity. Thus

$\displaystyle\langle\rho(t),\mathbbm{1}\rangle_{\mathcal{H}}$
$\displaystyle=\langle
m(t),\mathbbm{1}\rangle_{\mathcal{H}}+\langle\sqrt{C(t)}\xi,\mathbbm{1}\rangle_{\mathcal{H}},$
(40) $\displaystyle=1+\langle\xi,\sqrt{C(t)}\mathbbm{1}\rangle_{\mathcal{H}},$
$\displaystyle=1.$

This explains the importance of the conditions in Eq. (36): they ensure that
$\rho(t)$ is normalized.

##  Appendix B Theorem 1

######  Theorem 1.

Assume that:

  1. 1.

The system dynamics $f$ and $\mathfrak{h}$ are linear in state space:
$f(v,{t})=\mathsf{L}^{T}v$ and $\mathfrak{h}(v,{t})=\mathsf{H}v$, with
injective $\mathsf{H}$.

  2. 2.

$\Sigma=0$.

  3. 3.

$\rho(0)$ is chosen such that its mean $\mathsf{m}(0)$ and covariance
$\mathsf{C}(0)$ satisfy

$\begin{cases}\mathsf{H}\mathsf{m}(0)=Hm_{0},\\\
\mathsf{H}\mathsf{C}(0)\mathsf{H}^{T}=HC_{0}H^{*}.\end{cases}$ (41)

  4. 4.

$m(t)$ stays in the subspace

$\mathcal{S}\equiv\left\\{u\in\mathcal{H}\,\Bigl{|}\int|u(v)|v_{i}v_{j}dv<\infty\;\forall
i,j\in\\{1,\ldots,d\\}\right\\}$

and $C(t)$ stays in $L(\mathcal{S},\mathcal{S})$, the space of bounded linear
operators from $\mathcal{S}$ into itself.

Then, under the same noise realization for $Z^{\dagger}$,
$\mathsf{H}\mathsf{m}(t)=Hm(t)$ and
$\mathsf{H}\mathsf{C}(t)\mathsf{H}^{T}=HC(t)H^{*}$ will hold for $t\geq 0$,
where $\mathsf{m}(t)$ and $\mathsf{C}(t)$ are the mean and covariance of
$\rho(t)$ obtained from Eq. (34), and $m(t)$ and $C(t)$ are given by the KB
filter for densities (29a)–(29c).

###### Proof.

(Sketch)

We give here the outlines of a proof, but a rigorous proof, as well as
analysis of whether the equivalence holds in any setting more general than the
above restrictive conditions, will require considerably more work.

We consider the evolution of the mean and covariance of the KB filter for
densities (Eqs. (29a) and (29b)) projected into observation space,

$\displaystyle d(Hm)$
$\displaystyle=H\mathcal{L}^{*}m\,dt+HCH^{*}\Gamma^{-1}(dz^{\dagger}-Hm\,dt),$
(42a) $\displaystyle d(HCH^{*})$
$\displaystyle=H\mathcal{L}^{*}CH^{*}\,dt+HC\mathcal{L}H^{*}\,dt$
$\displaystyle\quad\quad\quad\quad\quad\quad-HCH^{*}\Gamma^{-1}HCH^{*}\,dt,$
(42b)

where $H(t)=H$ is not time-dependent because
$\mathfrak{h}(v,t)=\mathfrak{h}(v)=\mathsf{H}v$. These equations now describe
the time evolution of the finite-dimensional quantities $Hm$ and $HCH^{*}$.

Now, imposing $f(v,{t})=\mathsf{L}^{T}v$ and $\mathfrak{h}(v,{t})=\mathsf{H}v$
on the ansatz (Eq. (34)), the time evolution of $\rho$ can be entirely
characterized by its mean and covariance, and we obtain the following
equations for them:

$\displaystyle d\mathsf{m}$
$\displaystyle=\mathsf{L}^{T}\mathsf{m}\,dt+\mathsf{C}\mathsf{H}^{T}\Gamma^{-1}(dz^{\dagger}-\mathsf{H}\mathsf{m}\,dt),$
(43a) $\displaystyle d\mathsf{C}$
$\displaystyle=\mathsf{L}^{T}\mathsf{C}\,dt+\mathsf{C}\mathsf{L}\,dt-\mathsf{C}\mathsf{H}^{T}\Gamma^{-1}\mathsf{H}\mathsf{C}\,dt,$
(43b)

where $\mathsf{m}\equiv\mathbb{E}[v]$ and
$\mathsf{C}\equiv\mathbb{E}[(v-\mathsf{m})(v-\mathsf{m})^{T}]$. A similar
calculation is made in, e.g., section 7.4 of Jazwinski (1970)Jazwinski (1970).
In observation space, we have that

$\displaystyle d(\mathsf{H}\mathsf{m})$
$\displaystyle=\mathsf{H}\mathsf{L}^{T}\mathsf{m}\,dt+\mathsf{H}\mathsf{C}\mathsf{H}^{T}\Gamma^{-1}(dz^{\dagger}-\mathsf{H}\mathsf{m}\,dt),$
(44a) $\displaystyle d(\mathsf{H}\mathsf{C}\mathsf{H}^{T})$
$\displaystyle=\mathsf{H}\mathsf{L}^{T}\mathsf{C}\mathsf{H}^{T}\,dt+\mathsf{H}\mathsf{C}\mathsf{L}\mathsf{H}^{T}\,dt$
$\displaystyle\quad\quad\quad\quad\quad\quad-\mathsf{H}\mathsf{C}\mathsf{H}^{T}\Gamma^{-1}\mathsf{H}\mathsf{C}\mathsf{H}^{T}\,dt.$
(44b)

We would now like to show that $\mathsf{H}\mathsf{m}(t)=Hm(t)$ and
$\mathsf{H}\mathsf{C}(t)\mathsf{H}^{T}=HC(t)H^{*}$ for all $t\geq 0$. We do
this by showing that the RHS of Eqs. (42a) and (42b) are equal to the RHS of
Eqs. (44a) and (44b) at time $t^{*}$ if
$\mathsf{H}\mathsf{m}(t^{*})=Hm(t^{*})$ and
$\mathsf{H}\mathsf{C}(t^{*})\mathsf{H}^{T}=HC(t^{*})H^{*}$. Together with the
initial conditions (41) and uniqueness, this proves the theorem.

It follows immediately that

$\displaystyle
HC(t^{*})H^{*}\Gamma^{-1}\\!\left[\frac{dz^{\dagger}}{dt}-Hm(t^{*})\right]$
(45)
$\displaystyle\qquad\qquad=\mathsf{H}\mathsf{C}(t^{*})\mathsf{H}^{T}\Gamma^{-1}\\!\left[\frac{dz^{\dagger}}{dt}-\mathsf{H}\mathsf{m}(t^{*})\right],$

and that

$HC(t^{*})H^{*}\Gamma^{-1}HC(t^{*})H^{*}=\mathsf{H}\mathsf{C}(t^{*})\mathsf{H}^{T}\Gamma^{-1}\mathsf{H}\mathsf{C}(t^{*})\mathsf{H}^{T}.$
(46)

Note that

$\displaystyle\mathsf{H}\mathsf{m}(t^{*})$ $\displaystyle=Hm(t^{*}),$ (47)
$\displaystyle=\mathsf{H}\int vm(t^{*})\,dv,$

which implies that

$\mathsf{m}(t^{*})=\int vm(t^{*})\,dv,$ (48)

because $\mathsf{H}$ was assumed to be injective.

We proceed with the rest of the terms. For the first term of the RHS of Eq.
(42a),

$\displaystyle H\mathcal{L}^{*}m$ $\displaystyle=\mathsf{H}\int
v\mathcal{L}^{*}m\,dv,$ (49) $\displaystyle=-\mathsf{H}\int
v\nabla\cdot(mf)\,dv,$ $\displaystyle=-\mathsf{H}\mathsf{L}^{T}\int
v\nabla\cdot(mv)\,dv,$ $\displaystyle=\mathsf{H}\mathsf{L}^{T}\int vm\,dv,$
$\displaystyle=\mathsf{H}\mathsf{L}^{T}\mathsf{m},$

where the fourth line follows from integration by parts and the last from Eq.
(48). Note that the boundary term in the integration by parts vanishes from
assumption 4. Thus,

$H\mathcal{L}^{*}m=\mathsf{H}\mathsf{L}^{T}\mathsf{m}.$ (50)

It remains to show that
$H\mathcal{L}^{*}C(t^{*})H^{*}=\mathsf{H}\mathsf{L}^{T}\mathsf{C}(t^{*})\mathsf{H}^{T}$.
We have that for any $u$,

$HC(t^{*})H^{*}u=\mathsf{H}\int
vC(t^{*})H^{*}u\,dv=\mathsf{H}\mathsf{C}(t^{*})\mathsf{H}^{T}u.$ (51)

Since $\mathsf{H}$ was assumed to be injective,

$\int vC(t^{*})H^{*}u\,dv=\mathsf{C}(t^{*})\mathsf{H}^{T}u.$ (52)

Then, for any $w$,

$\displaystyle H\mathcal{L}^{*}C(t^{*})H^{*}w$ $\displaystyle=\mathsf{H}\int
v\mathcal{L}^{*}C(t^{*})H^{*}w\,dv$ (53) $\displaystyle=-\mathsf{H}\int
v\nabla\cdot(C(t^{*})H^{*}w\mathsf{L}^{T}v)\,dv$
$\displaystyle=\mathsf{H}\mathsf{L}^{T}\int vC(t^{*})H^{*}w\,dv$
$\displaystyle=\mathsf{H}\mathsf{L}^{T}\mathsf{C}(t^{*})\mathsf{H}^{T}w$

where the third line follows from integration by parts (with the boundary term
vanishing by the same argument as above) and the last line from Eq. (52).
Taking the adjoint demonstrates that
$HC(v,t^{*})\mathcal{L}H^{*}=\mathsf{H}\mathsf{C}(t^{*})\mathsf{L}\mathsf{H}^{T}$,
completing the proof. ∎

##  Appendix C Mean-Field Approximation

We omit the function arguments until the end of the subsection, for brevity.
Using Eq. 3.30 from Calvello, Reich, and Stuart (2022)Calvello, Reich, and
Stuart (2022), we know that the FP equation of Eq. (35) when $f(v,{t})=0$ and
$\Sigma=0$ is

$\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle=-\nabla\cdot(\rho(a-KH\rho))-\left\langle\nabla\cdot(\rho
K^{T}),\frac{dz^{\dagger}}{dt}\right\rangle$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+\nabla\cdot\bigl{(}\nabla\cdot(\rho
K\Gamma K^{T})\bigr{)}.$ (54)

We now match the terms of Eqs. (54) and (34) to make them equal. By matching
the terms involving $dz^{\dagger}/dt$, we obtain that

$\Gamma^{-1}({\mathfrak{h}}-H\rho)\rho=-\nabla\cdot(\rho K^{T}),$ (55)

and matching the rest of the terms,

$-\rho\langle{\mathfrak{h}}-H\rho,H\rho\rangle_{\Gamma}=-\nabla\cdot(\rho(a-KH\rho))+\nabla\cdot(\nabla\cdot(\rho
K\Gamma K^{T})).$ (56)

Substituting Eq. (55) into Eq. (56), we obtain

$\displaystyle\langle\nabla\cdot(\rho K^{T}),H\rho\rangle$
$\displaystyle=\nabla\cdot(\rho KH\rho)$ (57)
$\displaystyle=-\nabla\cdot(\rho(a-KH\rho))$
$\displaystyle\quad\quad\quad+\nabla\cdot(\nabla\cdot(\rho K\Gamma K^{T})).$

Setting the term in the divergence to 0, we obtain

$a=K\Gamma K^{T}\nabla\log\rho.$ (58)

This is the origin of the score function term discussed in subsection II.5.

We propose a test function $\psi(v)=v-\mathbb{E}v$, take the outer product of
it with both sides of Eq. (55), and integrate by parts, obtaining the identity

$\mathbb{E}K=\mathbb{E}[\psi({\mathfrak{h}}-H\rho)^{T}]\Gamma^{-1}=C^{v\mathfrak{h}}\Gamma^{-1},$
(59)

where
$C^{v\mathfrak{h}}{(t)}\equiv\mathbb{E}[(\mathfrak{h}(v,{t})-H\rho)(\mathfrak{h}(v,{t})-H\rho)^{T}]$.

Fixing the value of the gain $K$ to its expectation (the constant gain
approximation discussed in Calvello, Reich, and Stuart (2022)Calvello, Reich,
and Stuart (2022)), we then obtain

$K{(t)}=C^{v\mathfrak{h}}{(t)}\Gamma{(t)}^{-1}.$ (60)

Thus, the mean-field model is

$\displaystyle dv$
$\displaystyle=f(v,{t})\,dt+\sqrt{\Sigma{(t)}}\,dW+K{(t)}\Bigl{(}dz^{\dagger}-d\hat{z}\Bigr{)}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad+K{(t)}\Gamma{(t)}K{(t)}^{T}\nabla\log\rho{(v,t)}\,dt,$
$\displaystyle d\hat{z}$
$\displaystyle={(\mathbb{E}\mathfrak{h})(t)}\,dt+\sqrt{\Gamma{(t)}}dB,$

which gives Eqs. (10), with (10a) replaced by (II.2).