# Machine Learning in Viscoelastic Fluids via Energy-Based Kernel Embedding

Samuel E. Otto
Department of Mechanical Engineering,
University of Washington, Seattle, EUA
Cassio M. Oishi
Departamento de Matemática e Computação, Faculdade de Ciências e Tecnologia,
São Paulo State University, Presidente Prudente, Brazil
Fabio Amaral
Departamento de Matemática e Computação, Faculdade de Ciências e Tecnologia,
São Paulo State University, Presidente Prudente, Brazil
Steven L. Brunton
Department of Mechanical Engineering,
University of Washington, Seattle, EUA
J. Nathan Kutz
Department of Applied Mathematics,
University of Washington, Seattle, EUA
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>0000-0002-5086-8057 0000-0002-0904-6561 0000-0001-6945-8376
0000-0002-6565-5118 0000-0002-6004-2275

###### Abstract

The ability to measure differences in collected data is of fundamental
importance for quantitative science and machine learning, motivating the
establishment of metrics grounded in physical principles. In this study, we
focus on the development of such metrics for viscoelastic fluid flows governed
by a large class of linear and nonlinear stress models. To do this, we
introduce a kernel function corresponding to a given viscoelastic stress model
that implicitly embeds flowfield snapshots into a Reproducing Kernel Hilbert
Space (RKHS) whose squared norm equals the total mechanical energy. Working
implicitly with lifted representations in the RKHS via the kernel function
provides natural and unambiguous metrics for distances and angles between
flowfields without the need for hyperparameter tuning. Additionally, we
present a solution to the preimage problem for our kernels, enabling accurate
reconstruction of flowfields from their RKHS representations. Through
numerical experiments on an unsteady viscoelastic lid-driven cavity flow, we
demonstrate the utility of our kernels for extracting energetically-dominant
coherent structures in viscoelastic flows across a range of Reynolds and
Weissenberg numbers. Specifically, the features extracted by Kernel Principal
Component Analysis (KPCA) of flowfield snapshots using our kernel functions
yield reconstructions with superior accuracy in terms of mechanical energy
compared to conventional methods such as ordinary Principal Component Analysis
(PCA) with naïvely-defined state vectors or KPCA with ad-hoc choices of kernel
functions. Our findings underscore the importance of principled choices of
metrics in both scientific and machine learning investigations of complex
fluid systems.

Keywords: viscoelastic flow, energy-based inner product, kernel method,
machine learning, reproducing kernel Hilbert space, principal component
analysis

##  1 Introduction

A basic component of quantitative science is the ability to measure
differences in collected data. Quantities such as mass, momentum, and energy
provide unified descriptions of physical systems, and are therefore desirable
quantities to measure when investigating the behavior of a system of interest.
Likewise, in machine learning it is often necessary to endow the space in
which data lies with geometric notions such as distance and angle, or
otherwise to embed the data in a space with these notions. The assumed
geometry, including for instance, the way measurements are normalized, can
drastically affect the outcome of learning, potentially highlighting spurious
features while ignoring important ones that appear insignificant due to a poor
choice of metric. Therefore, in machine learning for physics applications it
is important to ground the ways we compare data in principled physical notions
such as energy.

One of the most important tasks in machine learning involves extracting low-
dimensional features (variables) that describe a system or allow one to
predict quantities of interest. Since its first applications to fluid dynamics
in [1], the Proper Orthogonal Decomposition (POD), also known as principal
component analysis (PCA) or the Karhunen-Loéve (KL) expansion, has been widely
used for encoding high dimensional data in a low dimensional representation
([2, 3]). Other modal analysis techniques such as those reviewed in [4] have
emerged as powerful tools to reduce the dimension of complex flow by splitting
them into simpler components or “modes”.

Many of these modal analysis techniques, including POD, rely on the choice of
an inner product defined on the flow’s state space. The fact that
significantly different results can be obtained using different inner products
motivates the introduction of principled choices based on physics. For
incompressible Newtonian fluid flows the integrated dot product of velocity
fields is a natural inner product whose resulting (squared) norm is the flow’s
kinetic energy. Applying PCA with this inner product extracts the orthogonal
modes and mode coefficients that are optimal for reconstructing flowfields in
an energetic sense. An inner product with analogous properties for isothermal
compressible flows was introduced by Rowley et al. [5] and used to construct
POD-Galerkin reduced-order models. Moreover, this work showed that
compatibility of the inner product with a conserved or dissipated energy
function guarantees that the stability of an equilibrium is preserved by
Galerkin projection. Energy and dimensionally-consistent inner products for
more general classes of compressible flows were later investigated in [6].
Following this same line of work, energy-consistent inner products have been
introduced for magnetohydrodynamics in [7] and for rotating shallow-water
equations in [8].

In the new era of machine learning applications in fluid mechanics [9], there
is growing interest in nonlinear dimensionality reduction techniques rooted in
manifold learning and artificial neural networks [10, 11, 12]. These
methodologies allow for greater dimensionality reduction, especially in
advection-dominated flow problems, by capturing curved manifolds in the state
space that are poorly approximated by low-dimensional subspaces, i.e., by
superpositions of modes [13]. Complementary to learned manifolds,
dimensionality reduction methods condense states along “fibers” determined by
the choice of reduced modeling variables. By allowing these variables to be
nonlinear functions of the state, nonlinear dimensionality reduction methods
also allow for nonlinear fibers that appropriately group states with similar
dynamical behaviors [14].

Kernel Principal Component Analysis (KPCA) was introduced by Schölkopf et al.
[15] and has become one of the most widely used methods for nonlinear
dimensionality reduction. The key insight is that PCA — a linear
dimensionality reduction method — can be applied through the use of a kernel
function in a high or infinite-dimensional Reproducing Kernel Hilbert Space
(RKHS) into which states have been lifted as illustrated in Fig. 1. Recent
works, such as [16, 17, 18], have leveraged the enhanced capabilities of
kernel-based nonlinear dimensionality reduction for modeling complex fluid
flows. However, to the best of our knowledge, all formulations presented to
date have only been used in the context of Newtonian fluid flows. Moreover,
with a wide range of kernel functions to choose from and the ability to
combine kernel functions to produce new ones, it is a matter of great
practical importance to narrow this selection based on principled physical
considerations.

Viscoelastic fluids play a pervasive role across various industrial sectors,
encompassing consumer goods, food, healthcare, and more. These applications
motivate the use of machine learning models and nonlinear dimensionality
reduction to shed light on the complex interplay of inertial and elastic
effects that give rise to unique unsteady nonlinear dynamics in viscoelastic
fluid flows (see [19]). This is especially pertinent in the context of Elastic
Turbulence (ET) ([20, 21]) and Elasto-Inertial Turbulence (EIT) ([22, 23]),
two phenomena posing distinctive challenges in the realm of viscoelastic fluid
dynamics. Roughly speaking, defining the elasticity parameter as
$E=\frac{Wi}{Re}$, where $Re$ represents the Reynolds number and $Wi$ is the
Weissenberg number, ET manifests itself in the turbulent flow of viscoelastic
fluids when $E$ is large, i.e., for inertialess or creeping flows. On the
other hand, EIT characterizes the turbulent behavior of highly viscoelastic
fluids influenced by inertial forces, particularly in inertia-dominated flows
at moderate or high Reynolds numbers. Both ET and EIT, along with the
identification of regime transitions ([24, 25, 26, 27]), present challenges
and opportunities to enhance our understanding of these phenomena through
machine learning models based on data gathered from computational fluid
dynamics simulations. More broadly, machine learning models have the potential
to enhance our ability to optimize and control non-Newtonian fluid systems
arising in a wide range of engineering applications.

${\mathcal{F}}$$\Phi({\mathcal{F}})$RKHS,
${\mathcal{H}}$${\mathcal{E}}({\boldsymbol{q}})=E_{2}$$\|\phi\|_{{\mathcal{H}}}^{2}=E_{2}$${\mathcal{E}}({\boldsymbol{q}})=E_{1}$$\|\phi\|_{{\mathcal{H}}}^{2}=E_{1}$
Figure 1: Lifting states into a Reproducing Kernel Hilbert Space (RKHS) with
squared norm equal to total mechanical energy. We illustrate how constant-
energy surfaces in the state space with energies $E_{1}$ and $E_{2}$ are
lifted onto spheres centered about the origin in the RKHS with squared radii
$E_{1}$ and $E_{2}$. Improperly chosen metrics in the state space can cause
states with different mechanical energies to appear similar. This problem is
rectified by measuring distances between lifted states in the RKHS.

Following the line of work initiated by Rowley et al. [5] and Schölkopf et al.
[15], in this paper we show that the mechanical energy of viscoelastic fluid
flows can be used to formulate kernel functions giving rise to natural metrics
for distances and angles between states that are well-suited for physics-
informed machine learning applications. Specifically, the kernel function
associated with a given viscoelastic stress model implicitly embeds the flow’s
state in a Reproducing Kernel Hilbert Space (RKHS) with squared norm equal to
the mechanical energy. The induced distance function turns the set of states
with finite energy into a complete, separable metric space. Kernel-based
machine learning algorithms (see [28]) rely on the induced geometric relations
between data points lifted into an RKHS for classification, regression, and
dimensionality reduction. As we mentioned above, it is often difficult to make
a principled choice for the kernel function from among a plethora of commonly
used kernels with various hyperparameters and with different ways of combining
kernels to produce new ones. The kernels we introduce have no hyperparameters,
partially resolving this ambiguity for problems involving viscoelastic fluid
flows. We note that the resolution is only partial because there could still
be other kernel functions that are compatible with the energy and turn the
space of states with finite energy into a complete, separable metric space.

Beyond providing principled choices of kernel functions for viscoelastic
flows, we also provide a solution of the preimage problem (see [29, 30]) for
our kernels. That is, the problem of reconstructing the flowfield from its
lifted representation in the RKHS, or its truncated representation in terms of
kernel principal components. We prove that the velocity and matrix square root
of the conformation tensor field can always be reconstructed from the RKHS
representation of a state via a bounded linear operator. As a consequence, we
show that these fields can be linearly reconstructed from truncated kernel
principal components with guaranteed accuracy depending on the truncated
singular values. Other fields can also be linearly reconstructed depending on
the viscoelastic stress model, but the velocity and square root of the
conformation tensor are sufficient to explicitly reconstruct the entire
flowfield.

To illustrate the utility of our kernel functions and the importance of making
principled choices for measuring distances between states, we perform KPCA-
based dimensionality reduction and reconstruction of an unsteady viscoelastic
lid-driven cavity flow. The flow is simulated at both low and moderate
Reynolds numbers and at both low and high Weissenberg numbers using different
viscoelastic stress models. The selection of the specific pairs of Reynolds
and Weissenberg numbers is grounded in a careful consideration of previous
works that have investigated the complexities of viscoelastic flows under
inertial effects ([31, 32, 33]) as well as investigations specifically focused
on problems related to creeping flows ([34, 35]) While ordinary PCA is optimal
for reconstructing the flowfields with respect to the Frobenius norm, we show
that this can lead to poor reconstruction of the flow’s mechanical energy. In
contrast, by matching the kernel function to the viscoelastic stress model,
Reynolds number, and Weissenberg number of the simulation, kernel PCA is able
to extract low-dimensional coordinates that faithfully reconstruct the
flowfield with superior accuracy in an energetic sense.

##  2 Viscoelastic flow models and energy functions

We consider nondimensional equations for viscoelastic flow [19] on a set
$\Omega\subset\mathbb{R}^{d}$ whose state
${\boldsymbol{q}}=({\boldsymbol{u}},{\boldsymbol{c}})$ consists of the
velocity field ${\boldsymbol{u}}:\Omega\to\mathbb{R}^{d}$ and the conformation
tensor field
${\boldsymbol{c}}:\Omega\to\mathbb{S}_{+}^{d}\subset\mathbb{R}^{d\times d}$.
The conformation tensor field takes values in the space $\mathbb{S}_{+}^{d}$
of positive semidefinite real symmetric matrices. We use $\mathbb{S}^{d}$ to
denote the real symmetric matrices and we use $\mathbb{S}_{++}^{d}$ to denote
the real symmetric positive definite matrices. The time evolution of the state
${\boldsymbol{q}}(t)$ is governed by the momentum equation

$\operatorname{\frac{\operatorname{\mathrm{d}}}{\operatorname{\mathrm{d}}t}}{\boldsymbol{u}}+{\boldsymbol{u}}\cdot\operatorname{\nabla}{\boldsymbol{u}}+\operatorname{\nabla}p=\frac{\beta}{\mbox{{Re}}}\operatorname{\Delta}{\boldsymbol{u}}+\frac{1}{\mbox{{Re}}}\operatorname{\nabla\cdot}\boldsymbol{\tau}+{\boldsymbol{f}},$
(2.1)

subject to the incompressibility constraint

$\operatorname{\nabla\cdot}{\boldsymbol{u}}=0,$ (2.2)

and the conformation tensor transport equation

$\overset{\bigtriangledown}{{\boldsymbol{c}}}={\boldsymbol{s}}({\boldsymbol{c}}),$
(2.3)

where $\beta$ is the viscosity ratio satisfying $0<\beta<1$. The upper
convected derivative [36], adopted in Eq. (2.3) is defined as

$\overset{\bigtriangledown}{{\boldsymbol{c}}}=\operatorname{\frac{\operatorname{\mathrm{d}}}{\operatorname{\mathrm{d}}t}}{\boldsymbol{c}}+{\boldsymbol{u}}\cdot\operatorname{\nabla}{\boldsymbol{c}}-{\boldsymbol{c}}\operatorname{\nabla}{\boldsymbol{u}}-(\operatorname{\nabla}{\boldsymbol{u}})^{T}{\boldsymbol{c}}$
(2.4)

These equations are coupled through a stress model

$\boldsymbol{\tau}=-(1-\beta){\boldsymbol{s}}({\boldsymbol{c}}),$ (2.5)

determined by a function
${\boldsymbol{s}}:D({\boldsymbol{s}})\subset\mathbb{S}_{+}^{d}\to\mathbb{S}^{d}$.
There are a wide variety of available stress models taking the general form

${\boldsymbol{s}}({\boldsymbol{c}})=-\frac{h_{0}({\boldsymbol{c}})}{\mbox{{Wi}}}\left[h_{1}({\boldsymbol{c}}){\boldsymbol{c}}-h_{2}({\boldsymbol{c}}){\boldsymbol{I}}+\kappa({\boldsymbol{c}}-{\boldsymbol{I}})^{2}\right],$
(2.6)

where $\kappa$ is a constant and $h_{0}$, $h_{1}$, and $h_{2}$ are scalar-
valued functions of the conformation tensor. Some examples are given in table
1. Using the stress model, the total mechanical energy of the system is given
by

${\mathcal{E}}({\boldsymbol{q}})=\frac{1}{2}\int_{\Omega}\left[\big{|}{\boldsymbol{u}}({\boldsymbol{x}})\big{|}^{2}+\frac{1}{\mbox{{Re}}}\operatorname{tr}(\boldsymbol{\tau})+\frac{(1-\beta)c}{\mbox{{Re}}\mbox{{Wi}}}\right]\operatorname{\mathrm{d}}{\boldsymbol{x}}.$
(2.7)

Here, we are adopting the energy defined in [37] with the addition of a model-
dependent constant $c\geq 0$ selected to ensure that the energy is always non-
negative. More details about the energy estimation for viscoelastic fluids can
be found in [38]. For the models in table 1, this constant may be taken to be
equal to the dimension $d$ of the flow domain. To simplify the notation, we
define

$\theta:=\frac{1-\beta}{\mbox{{Re}}\mbox{{Wi}}}.$

Model | $h_{0}({\boldsymbol{c}})$ | $h_{1}({\boldsymbol{c}})$ | $h_{2}({\boldsymbol{c}})$ | $\kappa$ | constraints on $D({\boldsymbol{s}})$ | $c$
---|---|---|---|---|---|---
Oldroyd-B | $1$ | $1$ | $1$ | $0$ | None | $d$
Giesekus | $1$ | $1$ | $1$ | $\alpha$ | None | $d$
FENE-CR | $1$ | $\frac{L^{2}}{L^{2}-\operatorname{tr}({\boldsymbol{c}})}$ | $\frac{L^{2}}{L^{2}-\operatorname{tr}({\boldsymbol{c}})}$ | $0$ | $\operatorname{tr}({\boldsymbol{c}})<L^{2}$ | $d$
FENE-P | $1$ | $\frac{L^{2}}{L^{2}-\operatorname{tr}({\boldsymbol{c}})}$ | $1$ | $0$ | $\operatorname{tr}({\boldsymbol{c}})<L^{2}$ | $d$
Linear PTT | $1+\varepsilon\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})$ | $1$ | $1$ | $0$ | None | $d$
Nonlinear PTT | $\exp{\left[\varepsilon\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})\right]}$ | $1$ | $1$ | $0$ | None | $d$
Table 1: Some common stress models for viscoelastic flow (see [19]) in the
general form of Eq.0(2.6). With the exception of the linear Oldroyd-B model,
the nonlinear models introduce an additional nondimensional parameter, each
defined accordingly: $\alpha$ represents the mobility factor for the Giesekus
model, whereas $L^{2}$ and $\varepsilon$ correspond to the extensibility
parameters for the finitely extensible nonlinear elastic with Peterlin closure
(FENE) and the Phan-Thien–Tanner(PTT) models, respectively.

We denote the space of states with finite total mechanical energies by

${\mathcal{F}}=\left\\{{\boldsymbol{q}}\ :\
{\mathcal{E}}({\boldsymbol{q}})<\infty\right\\}.$

We declare two elements in ${\mathcal{F}}$ to be equal when they agree almost
everywhere on $\Omega$. In this paper we introduce energy-based distance
functions that allow states in ${\mathcal{F}}$ to be compared. The next
section explains how this can be accomplished for a large class of stress
models including those listed in table 1 by embedding ${\mathcal{F}}$ in a
Reproducing Kernel Hilbert Space (RKHS) whose squared norm equals the total
mechanical energy, as illustrated in Fig. 1. This enables a variety of kernel-
based machine learning algorithms (see [28]) to be applied in viscoelastic
flow problems with principled choices for the kernel functions based on flow
physics.

##  3 Energy-based reproducing kernels

First we review the concept of a Reproducing Kernel Hilbert Space (RKHS). An
RKHS ${\mathcal{H}}$ over ${\mathcal{F}}$ consists of functions
$\phi:{\mathcal{F}}\to\mathbb{R}$ where pointwise evaluation
$\phi\mapsto\phi({\boldsymbol{q}})$ is a bounded linear map for every
${\boldsymbol{q}}\in{\mathcal{F}}$. By the Riesz lemma [39], there is a unique
element $K_{{\boldsymbol{q}}}\in{\mathcal{H}}$ so that the value of every
$\phi\in{\mathcal{H}}$ at ${\boldsymbol{q}}\in{\mathcal{F}}$ is given by the
inner product

$\phi({\boldsymbol{q}})=\left\langle K_{{\boldsymbol{q}}},\
\phi\right\rangle_{{\mathcal{H}}}.$ (3.8)

The function $k:{\mathcal{F}}\times{\mathcal{F}}\to\mathbb{R}$ defined by

$k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})=K_{{\boldsymbol{q}}_{2}}({\boldsymbol{q}}_{1})=\left\langle
K_{{\boldsymbol{q}}_{1}},\
K_{{\boldsymbol{q}}_{2}}\right\rangle_{{\mathcal{H}}}$ (3.9)

is called the “reproducing kernel” of ${\mathcal{H}}$ and the map
$\Phi:{\boldsymbol{q}}\mapsto K_{{\boldsymbol{q}}}$ is called the “feature
map”. It is easy to verify that for every finite collection of states
${\boldsymbol{q}}_{1},\ldots,{\boldsymbol{q}}_{m}\in{\mathcal{F}}$ and
coefficients $a_{1},\ldots,a_{m}\in\mathbb{R}$ the kernel satisfies the
positive-definiteness condition

$\sum_{i=1}^{m}\sum_{i=1}^{m}a_{i}a_{j}k({\boldsymbol{q}}_{i},{\boldsymbol{q}}_{j})\geq
0.$ (3.10)

Conversely, any function $k:{\mathcal{F}}\times{\mathcal{F}}\to\mathbb{R}$
satisfying the above positive-definiteness condition is the reproducing kernel
of a unique RKHS thanks to the Moore-Aronszajn theorem (see [40, 41]).

For a large class of stress models, including those in table 1, we show that
there is an RKHS ${\mathcal{H}}$ whose feature map
$\Phi:{\mathcal{F}}\to{\mathcal{H}}$ is injective and respects the total
energy in the sense that

${\mathcal{E}}({\boldsymbol{q}})=\|\Phi({\boldsymbol{q}})\|_{{\mathcal{H}}}^{2}=k({\boldsymbol{q}},{\boldsymbol{q}}).$
(3.11)

This notion of compatibility with the total mechanical energy is illustrated
in Fig. 1. Using the distance between lifted states in the RKHS, the state
space ${\mathcal{F}}$ becomes a metric space with

$d_{{\mathcal{E}}}({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2}):=\left\|\Phi({\boldsymbol{q}}_{1})-\Phi({\boldsymbol{q}}_{2})\right\|_{{\mathcal{H}}}=\sqrt{k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{1})-2k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})+k({\boldsymbol{q}}_{2},{\boldsymbol{q}}_{2})}.$
(3.12)

The injective property of the feature map is required to ensure that
$d_{{\mathcal{E}}}({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})=0$ if and only
if ${\boldsymbol{q}}_{1}={\boldsymbol{q}}_{2}$ in ${\mathcal{F}}$, i.e, almost
everywhere in $\Omega$. We say that an RKHS is an “injective RKHS” with an
“injective kernel function” when the associated feature map is injective. Most
importantly, the metric in Eq. (3.12) can be computed using the kernel
function without doing explicit computations in the abstract, possibly
infinite-dimensional space ${\mathcal{H}}$.

Our main result, stated in the following theorem, provides general conditions
on the stress model ensuring there is an RKHS with injective feature map
compatible with the total energy. Moreover, the theorem provides an explicit
formula for the corresponding reproducing kernel. Before stating the result,
we need some preliminary definitions. Let
$\sigma({\boldsymbol{c}})\subset\mathbb{R}$ denote the set of eigenvalues
(spectrum) of a symmetric matrix ${\boldsymbol{c}}\in\mathbb{S}^{d}$ and let
${\boldsymbol{P}}_{{\boldsymbol{c}}}(\lambda)$ denote the orthogonal
projection onto the eigenspace of ${\boldsymbol{c}}$ corresponding to an
eigenvalue $\lambda\in\sigma({\boldsymbol{c}})$. Then the action of a function
$f:\sigma({\boldsymbol{c}})\to\mathbb{R}$ on the matrix ${\boldsymbol{c}}$ is
defined by

$f({\boldsymbol{c}}):=\sum_{\lambda\in\sigma({\boldsymbol{c}})}f(\lambda){\boldsymbol{P}}_{{\boldsymbol{c}}}(\lambda).$
(3.13)

It is easy to see that $1({\boldsymbol{c}})={\boldsymbol{I}}$,
$\operatorname{Id}_{\mathbb{R}}({\boldsymbol{c}})={\boldsymbol{c}}$,
$(f+g)({\boldsymbol{c}})=f({\boldsymbol{c}})+g({\boldsymbol{c}})$, and
$(fg)({\boldsymbol{c}})=f({\boldsymbol{c}})g({\boldsymbol{c}})$, making
$f\mapsto f({\boldsymbol{c}})$ a commutative algebra homomorphism commonly
referred to as the “functional calculus” of ${\boldsymbol{c}}$. This extends
the usual notions of matrix polynomials, matrix square roots, and matrix
exponentials.

######  Theorem 1.

Suppose that there are (measurable) functions $\\{f_{i}\\}_{i=0}^{\infty}$ on
$\sigma(D({\boldsymbol{s}}))=\bigcup_{{\boldsymbol{c}}\in
D({\boldsymbol{s}})}\sigma({\boldsymbol{c}})$ and nonnegative constants
$\\{c_{i,p}\\}_{i,p=0}^{\infty}$ so that

$h({\boldsymbol{c}}):=-\mbox{{Wi}}\cdot\operatorname{tr}[{\boldsymbol{s}}({\boldsymbol{c}})]+c=\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}c_{i,p}\left[\operatorname{tr}\big{(}f_{i}({\boldsymbol{c}})^{2}\big{)}\right]^{p}$
(3.14)

for every ${\boldsymbol{c}}\in D({\boldsymbol{s}})$. Then the series

$\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2}):=\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}c_{i,p}\left[\operatorname{tr}\big{(}f_{i}({\boldsymbol{c}}_{1})f_{i}({\boldsymbol{c}}_{2})\big{)}\right]^{p}$
(3.15)

converges absolutely for every ${\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2}\in
D({\boldsymbol{s}})$ and the function
$k:{\mathcal{F}}\times{\mathcal{F}}\to\mathbb{R}$ defined by

$k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2}):=\frac{1}{2}\int_{\Omega}\left[{\boldsymbol{u}}_{1}({\boldsymbol{x}})\cdot{\boldsymbol{u}}_{2}({\boldsymbol{x}})+\theta\tilde{k}\big{(}{\boldsymbol{c}}_{1}({\boldsymbol{x}}),{\boldsymbol{c}}_{2}({\boldsymbol{x}})\big{)}\right]\operatorname{\mathrm{d}}{\boldsymbol{x}}$
(3.16)

is a positive-definite kernel satisfying
${\mathcal{E}}({\boldsymbol{q}})=k({\boldsymbol{q}},{\boldsymbol{q}})$.
Moreover, if there is a coefficient $c_{i,p}>0$ with $p\geq 1$, $f_{i}$
injective on $\sigma(D({\boldsymbol{s}}))$, and $p$ odd or $f_{i}$
nonnegative, then the feature map of the corresponding RKHS is injective. We
provide a proof in A.

Furthermore, ${\mathcal{F}}$ is a complete metric space, as we state in the
next theorem. Intuitively, this means that ${\mathcal{F}}$ has no missing
points that can be approached, but never reached. Completeness is important
for understanding limits that appear, for example, when studying the long-time
behavior of the system. The fact that ${\mathcal{F}}$ is complete provides
additional evidence that our kernel-based metric is a natural one for the
states of viscoelastic flows.

######  Theorem 2.

With the same assumptions as Theorem 1, suppose that there is a coefficient
$c_{i,p}>0$ with $p\geq 1$, $f_{i}$ injective on
$\sigma(D({\boldsymbol{s}}))$, and $p$ odd or $f_{i}$ nonnegative. Then
${\mathcal{F}}$ is a complete metric space with metric given by Eq. (3.12) and
$\Phi({\mathcal{F}})$ is a closed subset of ${\mathcal{H}}$. A proof is
provided in B.

To see how these theorems work, we use them to show that the stress models in
table 1 admit injective positive-definite kernel functions turning
${\mathcal{F}}$ into a complete metric space. However, we note that in certain
cases the parameters in the stress models must be constrained. Interestingly,
for the Oldroyd-B, Giesekus, and linear PTT stress models it is also possible
to explicitly construct feature maps $\Psi:{\mathcal{F}}\to L^{2}(\Omega)$
acting point-wise so that

$k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})=\langle\Psi({\boldsymbol{q}}_{1}),\
\Psi({\boldsymbol{q}}_{2})\rangle_{L^{2}(\Omega)}:=\int_{\Omega}\Psi({\boldsymbol{q}}_{1})({\boldsymbol{x}})^{T}\Psi({\boldsymbol{q}}_{2})({\boldsymbol{x}})\operatorname{\mathrm{d}}{\boldsymbol{x}}.$
(3.17)

As a result of the Moore-Aronszajn theorem (see Theorem 4 in A), for these
cases the RKHS ${\mathcal{H}}$ is identified isometrically with the closed
subspace of $L^{2}(\Omega)$ spanned by
$\\{\Psi({\boldsymbol{q}})\\}_{{\boldsymbol{q}}\in{\mathcal{F}}}$. We let
$\operatorname{vec}({\boldsymbol{A}})$ denote any vectorization of a matrix
${\boldsymbol{A}}$ and we recall that
$\operatorname{vec}({\boldsymbol{A}}_{1})^{T}\operatorname{vec}({\boldsymbol{A}}_{2})=\operatorname{tr}({\boldsymbol{A}}_{1}^{T}{\boldsymbol{A}}_{2})$.

######  Example 1.

For the Oldroyd-B model, we have

$h({\boldsymbol{c}})=\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})+d=\operatorname{tr}({\boldsymbol{c}}),$
(3.18)

which can be written in the form of Eq. (3.14) with $f_{0}:x\mapsto\sqrt{x}$
and coefficient $c_{0,1}=1$. Thus, Eq. (3.16) with

$\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2}):=\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)},$
(3.19)

is an injective positive semi-definite kernel function satisfying Eq. (3.11)
for the Oldroyd-B. Moreover, it is easy to see that

$\Psi({\boldsymbol{q}})({\boldsymbol{x}})=\frac{1}{\sqrt{2}}\begin{bmatrix}{\boldsymbol{u}}({\boldsymbol{x}})\\\
\sqrt{\theta}\operatorname{vec}\big{(}\sqrt{{\boldsymbol{c}}({\boldsymbol{x}}})\big{)}\end{bmatrix}\in\mathbb{R}^{d+d^{2}}.$
(3.20)

provides an explicit feature map satisfying Eq. (3.17).

######  Example 2.

For the Giesekus model, we have

$h({\boldsymbol{c}})=\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})+\alpha\operatorname{tr}\left[({\boldsymbol{c}}-{\boldsymbol{I}})^{2}\right]+d=\operatorname{tr}({\boldsymbol{c}})+\alpha\operatorname{tr}\left[({\boldsymbol{c}}-{\boldsymbol{I}})^{2}\right],$
(3.21)

which can be written in the form of Eq. (3.14) with $f_{0}:x\mapsto\sqrt{x}$
and $f_{1}:x\mapsto(x-1)$ with coefficients $c_{0,1}=1$ and $c_{1,1}=\alpha$.
Thus, Eq. (3.16) with

$\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2}):=\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}+\alpha\operatorname{tr}\big{[}({\boldsymbol{c}}_{1}-{\boldsymbol{I}})({\boldsymbol{c}}_{2}-{\boldsymbol{I}})\big{]},$
(3.22)

is an injective positive semi-definite kernel function satisfying Eq. (3.11)
for the Giesekus model with parameter $\varepsilon\geq 0$. Moreover,

$\Psi({\boldsymbol{q}})({\boldsymbol{x}})=\frac{1}{\sqrt{2}}\begin{bmatrix}{\boldsymbol{u}}({\boldsymbol{x}})\\\
\sqrt{\theta}\operatorname{vec}\big{(}\sqrt{{\boldsymbol{c}}({\boldsymbol{x}}})\big{)}\\\
\sqrt{\theta\alpha}\operatorname{vec}\big{(}{\boldsymbol{c}}({\boldsymbol{x}})-{\boldsymbol{I}}\big{)}\end{bmatrix}\in\mathbb{R}^{d+2d^{2}},$
(3.23)

defines an explicit feature map satisfying Eq. (3.17).

######  Example 3.

For the FENE-CR model, we have

$h({\boldsymbol{c}})=\frac{L^{2}\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})}{L^{2}-\operatorname{tr}({\boldsymbol{c}})}+d=\frac{(L^{2}-d)\operatorname{tr}({\boldsymbol{c}})}{L^{2}-\operatorname{tr}({\boldsymbol{c}})}=(L^{2}-d)\sum_{p=1}^{\infty}\frac{1}{L^{2p}}\operatorname{tr}({\boldsymbol{c}})^{p},$
(3.24)

where the geometric series converges when
$|\operatorname{tr}({\boldsymbol{c}})|<L^{2}$, hence on $D({\boldsymbol{s}})$
for this model. This expression takes the form of Eq. (3.14) with
$f_{0}:x\mapsto\sqrt{x}$ and $c_{0,p}=(L^{2}-d)L^{-2p}$ for every $p\geq 1$.
Therefore, if the parameter $L$ of the FENE-CR model is chosen so that
$L^{2}>d$, then Eq. (3.16) with

$\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})=\frac{(L^{2}-d)\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}}{L^{2}-\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}},$
(3.25)

defines an injective positive semi-definite kernel function satisfying Eq.
(3.11).

######  Example 4.

For the FENE-P model, we have

$h({\boldsymbol{c}})=\frac{L^{2}\operatorname{tr}({\boldsymbol{c}})}{L^{2}-\operatorname{tr}({\boldsymbol{c}})}=L^{2}\sum_{p=1}^{\infty}\frac{1}{L^{2p}}\operatorname{tr}({\boldsymbol{c}})^{p},$
(3.26)

where the geometric series converges when
$|\operatorname{tr}({\boldsymbol{c}})|<L^{2}$, hence on $D({\boldsymbol{s}})$
for this model. This expression takes the form of Eq. (3.14) with
$f_{0}:x\mapsto\sqrt{x}$ and $c_{0,p}=L^{-2(p-1)}$ for every $p\geq 1$.
Therefore, Eq. (3.16) with

$\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})=\frac{L^{2}\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}}{L^{2}-\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}},$
(3.27)

defines an injective positive semi-definite kernel function satisfying Eq.
(3.11) for the FENE-P model. We note that unlike the FENE-CR model, there are
no constraints on the model parameter $L$.

######  Example 5.

For the linear PTT model, we have

$h({\boldsymbol{c}})=\left[1+\varepsilon\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})\right]\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})+d=\varepsilon
d^{2}+(1-2\varepsilon
d)\operatorname{tr}({\boldsymbol{c}})+\varepsilon\operatorname{tr}({\boldsymbol{c}})^{2},$
(3.28)

which takes the form of Eq. (3.14) with $f_{0}:x\mapsto\sqrt{x}$ and
coefficients $c_{0,0}=\varepsilon d^{2}$, $c_{0,1}=1-2\varepsilon d$, and
$c_{0,2}=\varepsilon$. Therefore, if the parameter $\varepsilon$ of the linear
PTT model satisfies $0\leq\varepsilon\leq(2d)^{-1}$, then Eq. (3.16) with

$\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})=\left[1+\varepsilon\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}-{\boldsymbol{I}}\big{)}\right]\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}-{\boldsymbol{I}}\big{)}+d,$
(3.29)

defines an injective positive semi-definite kernel function satisfying Eq.
(3.11). To provide an explicit feature map for this model, we first observe
that

$\displaystyle\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})$
$\displaystyle=d^{2}+(1-2\varepsilon
d)\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}+\varepsilon\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}^{2}$
(3.30) $\displaystyle=d^{2}+(1-2\varepsilon
d)\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}+\varepsilon\operatorname{tr}\Big{[}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}\otimes\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}\Big{]}$
$\displaystyle=d^{2}+(1-2\varepsilon
d)\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}\big{)}+\varepsilon\operatorname{tr}\Big{[}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\otimes\sqrt{{\boldsymbol{c}}_{1}}\big{)}\big{(}\sqrt{{\boldsymbol{c}}_{2}}\otimes\sqrt{{\boldsymbol{c}}_{2}}\big{)}\Big{]},$

where $\otimes$ denotes the Kronecker product of matrices. Here, we have used
the trace and mixed product properties of the Kronecker product. Thanks to the
above expression,

$\Psi({\boldsymbol{q}})({\boldsymbol{x}})=\frac{1}{\sqrt{2}}\begin{bmatrix}d\sqrt{\theta}\\\
{\boldsymbol{u}}({\boldsymbol{x}})\\\ \sqrt{\theta(1-2\varepsilon
d)}\operatorname{vec}\big{(}\sqrt{{\boldsymbol{c}}({\boldsymbol{x}}})\big{)}\\\
\sqrt{\theta\varepsilon}\operatorname{vec}\big{(}\sqrt{{\boldsymbol{c}}({\boldsymbol{x}})}\otimes\sqrt{{\boldsymbol{c}}({\boldsymbol{x}})}\big{)}\end{bmatrix}\in\mathbb{R}^{1+d+d^{2}+d^{4}},$
(3.31)

defines an explicit feature map satisfying Eq. (3.17) for the linear PTT
model.

######  Example 6.

For the nonlinear PTT model, we have

$\displaystyle h({\boldsymbol{c}})$
$\displaystyle=\exp\left[\varepsilon\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})\right]\operatorname{tr}({\boldsymbol{c}}-{\boldsymbol{I}})+d$
(3.32) $\displaystyle=d+e^{-\varepsilon
d}\left[\operatorname{tr}({\boldsymbol{c}})e^{\varepsilon\operatorname{tr}({\boldsymbol{c}})}-de^{\varepsilon\operatorname{tr}({\boldsymbol{c}})}\right]$
$\displaystyle=d+e^{-\varepsilon
d}\left[\sum_{p=1}^{\infty}\frac{1}{(p-1)!}\varepsilon^{p-1}\operatorname{tr}({\boldsymbol{c}})^{p}-d-d\sum_{p=1}^{\infty}\frac{1}{p!}\varepsilon^{p}\operatorname{tr}({\boldsymbol{c}})^{p}\right]$
$\displaystyle=\underbrace{d\left(1-e^{-\varepsilon
d}\right)}_{c_{0,0}}+\sum_{p=1}^{\infty}\underbrace{\frac{\varepsilon^{p-1}e^{-\varepsilon
d}}{(p-1)!}\left(1-\frac{\varepsilon
d}{p}\right)}_{c_{0,p}}\operatorname{tr}({\boldsymbol{c}})^{p},$

which takes the form of Eq. (3.14) with $f_{0}:x\mapsto\sqrt{x}$. We observe
that if $0\leq\varepsilon\leq d^{-1}$ then all of the coefficients $c_{0,p}$
are non-negative and $c_{0,p}$ is strictly positive for every $p\geq 2$.
Therefore, if the parameter $\varepsilon$ of the nonlinear PTT model satisfies
$0\leq\varepsilon\leq d^{-1}$, then Eq. (3.16) with

$\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})=\exp\left[\varepsilon\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}-{\boldsymbol{I}}\big{)}\right]\operatorname{tr}\big{(}\sqrt{{\boldsymbol{c}}_{1}}\sqrt{{\boldsymbol{c}}_{2}}-{\boldsymbol{I}}\big{)}+d,$
(3.33)

defines an injective positive semi-definite kernel function satisfying Eq.
(3.11).

###  3.1 Comparing states with different parameter values

So far, our kernel functions and the associated distance function on
${\mathcal{F}}$ only allows us to compare states at the same values of the
parameters $\beta$, Re, and Wi. However, this situation is easily rectified by
defining a new kernel function

$\check{k}\left(({\boldsymbol{q}}_{1},\theta_{1}),\
({\boldsymbol{q}}_{2},\theta_{2})\right):=\frac{1}{2}\int_{\Omega}\left[{\boldsymbol{u}}_{1}({\boldsymbol{x}})\cdot{\boldsymbol{u}}_{2}({\boldsymbol{x}})+\sqrt{\theta_{1}\theta_{2}}\tilde{k}\big{(}{\boldsymbol{c}}_{1}({\boldsymbol{x}}),{\boldsymbol{c}}_{2}({\boldsymbol{x}})\big{)}\right]\operatorname{\mathrm{d}}{\boldsymbol{x}},$
(3.34)

which depends on the values of the parameter
$\theta=(1-\beta)/(\mbox{{Re}}\mbox{{Wi}})$. This is a positive-definite
kernel function of a unique RKHS $\check{{\mathcal{H}}}$ over
${\mathcal{F}}\times(0,\infty)$ thanks to the kernel product rule stated in
Lemma 1 in A. Under the same conditions stated in Theorem 1, the associated
feature map
$\check{\Phi}:{\mathcal{F}}\times(0,\infty)\to\check{{\mathcal{H}}}$ is
“conditionally injective”, meaning that

$\Phi_{\theta}:{\boldsymbol{q}}\mapsto\check{\Phi}({\boldsymbol{q}},\theta)$
(3.35)

is injective for every $\theta>0$. To see this, we observe that
$k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})=\check{k}\left(({\boldsymbol{q}}_{1},\theta),\
({\boldsymbol{q}}_{2},\theta)\right)$. Since $k$ is an injective kernel
function, it follows that

$\|\Phi_{\theta}({\boldsymbol{q}}_{1})-\Phi_{\theta}({\boldsymbol{q}}_{2})\|_{\check{{\mathcal{H}}}}=d_{{\mathcal{E}}}({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})>0$

when ${\boldsymbol{q}}_{1}\neq{\boldsymbol{q}}_{2}$ in ${\mathcal{F}}$.

##  4 Reconstructing fields and preimages for kernel functions

Once states are lifted into the RKHS, it is often of practical importance to
be able to reconstruct them. This is referred to as the preimage problem in
kernel-based machine learning [29, 30]. The following theorem provides a
partial solution for the class of viscoelastic kernels presented above. It
describes spatial fields that can be linearly reconstructed from the lifted
representations of states in the RKHS. In particular, Corollary 1 means that
for each of the stress models listed in Table 1, the field
${\boldsymbol{y}}=({\boldsymbol{u}},\sqrt{{\boldsymbol{c}}})$ lies in
$L^{2}(\Omega;\ \mathbb{R}^{d}\times\mathbb{R}^{d\times d})$ and can be
reconstructed via a bounded linear operator acting on $\Phi({\boldsymbol{q}})$
for every ${\boldsymbol{q}}\in{\mathcal{F}}$. For the Oldroyd-B, Giesekus, and
linear PTT models, this property is evident from their explicit feature maps
given above in Eqs. (3.20), (3.23), and (3.31).

######  Theorem 3.

With the same assumptions as Theorem 1, let
${\boldsymbol{a}}_{0}\in\mathbb{R}^{d}$ and
${\boldsymbol{A}}_{i,p}\in\mathbb{R}^{d^{p}\times d^{p}}$ satisfy

$A:=\sqrt{2\|{\boldsymbol{a}}_{0}\|_{2}^{2}+\frac{2}{\theta}\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}\|{\boldsymbol{A}}_{i,p}\|_{F}^{2}}<\infty$
(4.36)

and consider the function $\psi:\mathbb{R}^{d}\times
D({\boldsymbol{s}})\to\mathbb{R}$ defined by

$\psi({\boldsymbol{u}},{\boldsymbol{c}})={\boldsymbol{a}}_{0}^{T}{\boldsymbol{u}}+\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}\sqrt{c_{i,p}}\operatorname{tr}\left[{\boldsymbol{A}}_{i,p}^{T}f_{i}({\boldsymbol{c}})^{\otimes
p}\right].$ (4.37)

Then there is a unique bounded linear operator $R_{\psi}:{\mathcal{H}}\to
L^{2}(\Omega)$ satisfying

$R_{\psi}\Phi({\boldsymbol{q}})=\psi\circ{\boldsymbol{q}}$ (4.38)

for every
${\boldsymbol{q}}=({\boldsymbol{u}},{\boldsymbol{c}})\in{\mathcal{F}}$ and
$\|R_{\psi}\|\leq A$. We provide a proof in C.

######  Corollary 1.

With the same assumptions as Theorem 1, suppose that there is a nonzero
coefficient $c_{i,p}>0$. Then there is a unique bounded linear operator
$R_{i,p}:{\mathcal{H}}\to L^{2}(\Omega;\
\mathbb{R}^{d}\times\mathbb{R}^{d^{p}\times d^{p}})$ satisfying

$R_{i,p}\Phi({\boldsymbol{q}})=\left({\boldsymbol{u}},\
f_{i}({\boldsymbol{c}})^{\otimes p}\right)$ (4.39)

for every
${\boldsymbol{q}}=({\boldsymbol{u}},{\boldsymbol{c}})\in{\mathcal{F}}$. It’s
norm is bounded by $\|R_{i,p}\|\leq\sqrt{2d+\frac{2d^{2p}}{\theta c_{i,p}}}$.

######  Remark 1.

Incidentally, the proof of Theorem 3 shows that the metric space
$({\mathcal{F}},d_{{\mathcal{E}}})$ and the RKHS ${\mathcal{H}}$ are
separable. This is a useful property for studying operators on
${\mathcal{H}}$, such as the covariance operator in kernel principal component
analysis described below.

These results justify linear reconstruction of the field
${\boldsymbol{y}}=({\boldsymbol{u}},\sqrt{{\boldsymbol{c}}})$ using low-
dimensional coordinates obtained via kernel Principal Component Analysis
(KPCA) [15]. Here, there is a probability measure $\mu$ over states in
${\mathcal{F}}$, and we seek a collection of low-dimensional features
describing states from this distribution. For the sake of simplicity, we
consider the uncentered version of KPCA where the covariance operator
$C_{\mu}:{\mathcal{H}}\to{\mathcal{H}}$ is defined by

$C_{\mu}f=\int_{{\mathcal{F}}}\Phi({\boldsymbol{q}})\langle\Phi({\boldsymbol{q}}),f\rangle_{{\mathcal{H}}}\operatorname{\mathrm{d}}\mu({\boldsymbol{q}}).$
(4.40)

If the probability measure has finite average energy, i.e.,
$\int_{{\mathcal{F}}}{\mathcal{E}}({\boldsymbol{q}})\operatorname{\mathrm{d}}\mu({\boldsymbol{q}})<\infty$,
then the covariance operator is self-adjoint, positive semidefinite, and
trace-class by Theorem 4.1 in [42]. Therefore, ${\mathcal{H}}$ admits an
orthonormal basis of eigenvectors $\\{u_{j}\\}_{j=1}^{\infty}$ of $C_{\mu}$
with eigenvalues $\lambda_{j}=\sigma_{j}^{2}$, $\sigma_{j}\geq 0$, arranged in
descending order and satisfying

$\sum_{j=1}^{\infty}\sigma_{j}^{2}<\infty.$

Using KPCA, the principal components

$z_{j}({\boldsymbol{q}}):=\langle u_{j},\
\Phi({\boldsymbol{q}})\rangle_{{\mathcal{H}}}$ (4.41)

can be computed without performing explicit operations in the RKHS. To do
this, KPCA uses a kernel integral operator
$K_{\mu}:L^{2}({\mathcal{F}},\mu)\to L^{2}({\mathcal{F}},\mu)$ defined by

$(K_{\mu}f)({\boldsymbol{q}})=\int_{{\mathcal{F}}}k({\boldsymbol{q}},{\boldsymbol{q}}^{\prime})f({\boldsymbol{q}}^{\prime})\operatorname{\mathrm{d}}\mu({\boldsymbol{q}}).$
(4.42)

This operator is self-adjoint and has the same eigenvalues $\lambda_{j}$ as
$C_{\mu}$, with eigenfunctions $\\{v_{j}\\}_{j=1}^{\infty}$ forming an
orthonormal basis for $L^{2}({\mathcal{F}},\mu)$. Crucially, these
eigenfunctions can be computed using the kernel function, enabling one to
compute the KPCA coordinates

$z_{j}({\boldsymbol{q}})=\sigma_{j}^{-1}\int_{{\mathcal{F}}}v_{j}({\boldsymbol{q}}^{\prime})k({\boldsymbol{q}}^{\prime},{\boldsymbol{q}})\operatorname{\mathrm{d}}\mu({\boldsymbol{q}}^{\prime}).$
(4.43)

Since $C_{\mu}$ is a compact operator on a Hilbert space, it can be
approximated in norm by finite-rank operators. One often uses empirical
covariance operators obtained by drawing finitely many independent samples
$\\{{\boldsymbol{q}}_{i}\\}_{i=1}^{m}$ from $\mu$ and forming the covariance
$C_{\mu_{m}}$ for the empirical measure
$\mu_{m}=m^{-1}\sum_{i=1}^{m}\delta_{{\boldsymbol{q}}_{i}}$ with
$\delta_{{\boldsymbol{q}}_{i}}$ denoting the Dirac measure centered at
${\boldsymbol{q}}_{i}$. With the identification
$L^{2}({\mathcal{F}},\mu_{m})\cong\mathbb{R}^{m}$, which is valid when the
sampled states are distinct, the corresponding kernel integral operator
$K_{\mu_{m}}$ becomes the kernel matrix
${\boldsymbol{K}}\in\mathbb{R}^{m\times m}$ with entries
$[{\boldsymbol{K}}]_{i,j}=m^{-1}k({\boldsymbol{q}}_{i},{\boldsymbol{q}}_{j})$
and eigenvectors ${\boldsymbol{v}}_{j}\in\mathbb{R}^{m}$. The KPCA features
for the empirical distribution are then given by

$z_{j}({\boldsymbol{q}})=m^{-1/2}\sigma_{j}^{-1}\begin{bmatrix}k({\boldsymbol{q}}_{1},{\boldsymbol{q}})&\cdots&k({\boldsymbol{q}}_{m},{\boldsymbol{q}})\end{bmatrix}{\boldsymbol{v}}_{j}.$
(4.44)

The following result says that certain spatial fields such as
$({\boldsymbol{u}},\sqrt{{\boldsymbol{c}}})$ can be linearly reconstructed
with guaranteed accuracy using the leading KPCA features. Here, we arrange the
first $r$ eigenvectors into an operator $U_{r}:\mathbb{R}^{r}\to{\mathcal{H}}$
defined by $(w_{1},\ldots,w_{r})\mapsto w_{1}u_{1}+\cdots w_{r}u_{r}$ and we
denote the vector of leading principal components
${\boldsymbol{z}}_{r}({\boldsymbol{q}})=(z_{1}({\boldsymbol{q}}),\ldots,z_{r}({\boldsymbol{q}}))=U_{r}^{*}\Phi({\boldsymbol{q}})$.

######  Proposition 1.

With the same assumptions as Theorem 1, suppose that there is a nonzero
coefficient $c_{i,p}>0$ and the average energy
$\int_{{\mathcal{F}}}{\mathcal{E}}({\boldsymbol{q}})\operatorname{\mathrm{d}}\mu({\boldsymbol{q}})<\infty$
is finite. For a state vector
${\boldsymbol{q}}=({\boldsymbol{u}},{\boldsymbol{c}})\in{\mathcal{F}}$, the
reconstruction error is bounded by

$\left\|\big{(}{\boldsymbol{u}},\ f_{i}({\boldsymbol{c}})^{\otimes
p}\big{)}-R_{i,p}U_{r}{\boldsymbol{z}}_{r}({\boldsymbol{q}})\right\|_{L^{2}(\Omega)}^{2}\leq
2\left(d+\frac{d^{2p}}{\theta
c_{i,p}}\right)\Big{(}{\mathcal{E}}({\boldsymbol{q}})-\|{\boldsymbol{z}}_{r}({\boldsymbol{q}})\|_{2}^{2}\Big{)}.$
(4.45)

The average reconstruction error for states drawn according to $\mu$ is
bounded by

$\int_{{\mathcal{F}}}\left\|\big{(}{\boldsymbol{u}},\
f_{i}({\boldsymbol{c}})^{\otimes
p}\big{)}-R_{i,p}U_{r}{\boldsymbol{z}}_{r}({\boldsymbol{q}})\right\|_{L^{2}(\Omega)}^{2}\operatorname{\mathrm{d}}\mu({\boldsymbol{q}})\leq
2\left(d+\frac{d^{2p}}{\theta
c_{i,p}}\right)\sum_{j=r+1}^{\infty}\sigma_{j}^{2}.$ (4.46)

We provide a proof in C.

We note that when the centered version of KPCA is used, the reconstruction
will involve an affine, rather than a purely linear map of the leading
principal components.

##  5 A case study: transient lid-driven cavity flow

Let us consider the lid-driven cavity problem defined over a square spatial
domain $[0,1]\times[0,1]$, and driven by the top lid moving to the right.
While the cavity flow problem may not have immediate practical applications,
it holds significant importance as a viscoelastic benchmark (see [43, 44,
45]). This simple geometry can be employed to address challenges associated
with the singularity of the stress field on the corner of the cavity as well
as to study the emergence of elastic turbulence for high Weissenberg number
([46]). In particular, we have illustrated the streamlines for two situations
in Fig. 2 for the lid-driven cavity flow considering the Oldroyd-B model.

To generate the numerical data for our investigation, we employed the same
code designed for solving viscoelastic fluid flows as in [47, 48, 49, 50]. In
this methodology, we employ finite-difference approximations for the
discretization of the momentum and mass equations (2.1) and (2.2), as well as
the constitutive equation (2.3). Managing the interdependence between the
velocity and pressure fields, a key challenge in solving the governing
equations (2.1) and (2.2), is addressed through the application of a
projection scheme. Widely recognized in computational fluid dynamics, this
scheme effectively decouples the velocity and pressure fields, simplifying the
computational process and enhancing the stability and accuracy of our
numerical simulations. Specifically, the momentum equation (2.1) is time-
discretized using a semi-implicit strategy, while an explicit time
discretization is employed for solving the constitutive equation (2.3).
Additionally, for the discretization of the convective terms presented in
equations (2.1) and (2.3), we apply a high-order upwind methodology. In order
to obtain efficient solutions for high Weissenberg number simulations, we
incorporate a stabilization scheme known as log-conformation ([19]). The log-
conformation method proves particularly effective in handling flows with high
Weissenberg numbers, ensuring stability and accuracy in capturing the
intricate rheological behavior of viscoelastic fluids.

In our numerical tests, we use a uniform mesh with a spatial discretization of
$\Delta x=1/80=0.0125$ within the square domain. Time integration was
performed with a fixed time-step of $\Delta t=10^{-5}$, ensuring temporal
stability and capturing the transient dynamics of the viscoelastic fluid flow.
The fixed mesh size and time-step, were selected based on a careful balance
between computational efficiency and the accuracy required for our
investigation.

Since we aim to analyze the dynamical behavior of the flow, we adopt the same
smooth, transient velocity profile for the top lid used in [44]:

$u(x,1,t)=16x^{2}(1-x)^{2}\sin(\pi t).\\\ $ (5.47)

The remaining cavity walls are stationary and we impose no-slip boundary
conditions.

a)
b)

Figure 2: Streamlines (left) and $\operatorname{tr}({\boldsymbol{c}})$ (right)
of the cavity flow using the viscoelastic Oldroyd-B model: a)
$Re=10,Wi=10,\beta=0.9$ and b) $Re=0.001,Wi=5,\beta=0.5$.

###  5.1 Dataset, algorithm and measurement of errors

The data from simulations are organized into a matrix
$\textbf{Y}\in\mathbb{R}^{D\times m}$, where $D$ represents the product of the
number of mesh points and the number of variables in each spatial location.
Each column is a modified state vector

${\boldsymbol{y}}=\begin{pmatrix}u\\\ v\\\ b_{xx}\\\ b_{xy}\\\
b_{yy}\end{pmatrix}\qquad\mbox{where}\qquad{\boldsymbol{b}}=\sqrt{{\boldsymbol{c}}}$
(5.48)

corresponding to one of the $m$ snapshots collected from a simulation. This
definition of a modified state vector is due the fact it can be linearly
reconstructed from the lifted states, as discussed in section 4. We also drop
the $b_{yx}$ entry of the square root conformation tensor since it is equal to
$b_{xy}$ by symmetry. The original state vector ${\boldsymbol{q}}$ can be
obtained from ${\boldsymbol{y}}$ by the relation
${\boldsymbol{c}}={\boldsymbol{b}}^{2}$ [37]. Hence, the matrix Y is
structured as

$\textbf{Y}=\left[\begin{array}[]{cccc}\rule[-4.30554pt]{0.5pt}{10.76385pt}&\rule[-4.30554pt]{0.5pt}{10.76385pt}&&\rule[-4.30554pt]{0.5pt}{10.76385pt}\\\
{\boldsymbol{y}}_{1}&{\boldsymbol{y}}_{2}&\dots&{\boldsymbol{y}}_{m}\\\
\rule[-4.30554pt]{0.5pt}{10.76385pt}&\rule[-4.30554pt]{0.5pt}{10.76385pt}&&\rule[-4.30554pt]{0.5pt}{10.76385pt}\end{array}\right],$
(5.49)

where ${\boldsymbol{y}}_{i}$ is the modified state vector defined by Eq. 5.48
for the $i$th state vector snapshot ${\boldsymbol{q}}_{i}$ gathered from a
simulation.

We use the Kernel Principal Component Analysis (KPCA) algorithm developed in
[15] to extract low-dimensional variables describing the simulation data. In
summary, the kernel matrix ${\boldsymbol{K}}\in\mathbb{R}^{m\times m}$ with
entries
$[{\boldsymbol{K}}]_{i,j}=k({\boldsymbol{q}}_{i},{\boldsymbol{q}}_{j})$ is
centered to form

${\boldsymbol{K}}_{c}=\left({\boldsymbol{I}}-\frac{1}{m}{\boldsymbol{1}}_{m}{\boldsymbol{1}}_{m}^{T}\right){\boldsymbol{K}}\left({\boldsymbol{I}}-\frac{1}{m}{\boldsymbol{1}}_{m}{\boldsymbol{1}}_{m}^{T}\right),$
(5.50)

where ${\boldsymbol{1}}_{m}$ denotes the vector with unit entries in
$\mathbb{R}^{m}$. Computing the symmetric eigenvalue decomposition of the
centered kernel matrix

${\boldsymbol{K}}_{c}={\boldsymbol{V}}{\boldsymbol{\Lambda}}{\boldsymbol{V}}^{T},$
(5.51)

and forming its rank-$r$ truncation
${\boldsymbol{V}}_{r}{\boldsymbol{\Lambda}}_{r}{\boldsymbol{V}}_{r}^{T}$
retaining the largest $r$ eigenvalues yields

${\boldsymbol{Z}}_{r}=\left[\begin{array}[]{cccc}\rule[-4.30554pt]{0.5pt}{10.76385pt}&\rule[-4.30554pt]{0.5pt}{10.76385pt}&&\rule[-4.30554pt]{0.5pt}{10.76385pt}\\\
{\boldsymbol{z}}_{r}({\boldsymbol{q}}_{1})&{\boldsymbol{z}}_{r}({\boldsymbol{q}}_{2})&\dots&{\boldsymbol{z}}_{r}({\boldsymbol{q}}_{m})\\\
\rule[-4.30554pt]{0.5pt}{10.76385pt}&\rule[-4.30554pt]{0.5pt}{10.76385pt}&&\rule[-4.30554pt]{0.5pt}{10.76385pt}\end{array}\right]=\sqrt{{\boldsymbol{\Lambda}}_{r}}{\boldsymbol{V}}_{r}^{T}.$
(5.52)

The columns of this matrix are vectors comprising the leading $r$ kernel
principal components of each snapshot. The kernel matrix is constructed using
Eq. (3.15) for the selected stress model. For example, if the Oldroyd-B model
is chosen, the kernel function Eq. (3.15) is computed using Eq. (3.19). It is
worth noting that the linear PCA reduction is also obtained as a special case
of KPCA using a kernel given by the classical
$L^{2}(\Omega;\mathbb{R}^{d}\times\mathbb{R}^{d\times d})$ inner product of
state vectors.

We reconstruct the modified snapshots via an affine function of the kernel
principal components

${\boldsymbol{\hat{y}}}_{i}={\boldsymbol{r}}_{0}+{\boldsymbol{R}}{\boldsymbol{z}}_{r}({\boldsymbol{q}}_{i}),$
(5.53)

by solving the least-squares problem

$\operatorname*{\arg\\!\min\enskip}_{{\boldsymbol{R}},{\boldsymbol{r}}_{0}}\|\textbf{Y}-{\boldsymbol{r}}_{0}{\boldsymbol{1}}_{m}^{T}-{\boldsymbol{R}}{\boldsymbol{Z}}_{r}\|_{F}^{2},$
(5.54)

where ${\boldsymbol{R}}\in\mathbb{R}^{D\times r}$ and
${\boldsymbol{r}}_{0}\in\mathbb{R}^{D}$. Affine reconstruction of modified
state vectors is justified per the discussion in Section 4. Approximations of
the original state vectors ${\boldsymbol{\hat{q}}}$ are then obtained from
${\boldsymbol{\hat{y}}}$ by the relation
${\boldsymbol{\hat{c}}}={\boldsymbol{\hat{b}}}^{2}$.

In order to study the effect of the kernel function in dimensionality
reduction, we consider several ways of measuring the reconstruction error. The
first way to measure the error is using the ordinary $L^{2}$ norm of the
reconstructed fields

$E_{F}=\frac{\sum_{i=1}^{m}\|{\boldsymbol{y}}_{i}-{\boldsymbol{\hat{y}}}_{i}\|_{L^{2}(\Omega)}^{2}}{\sum_{i=1}^{m}\|{\boldsymbol{y}}_{i}\|_{L^{2}(\Omega)}^{2}}.$
(5.55)

Note that the numerator and denominator can be computed as Frobenius norms of
respective snapshot matrices. The problem with using this notion of error is
that it is not grounded in the mechanical energy of the system, and could
assign improper weight to the kinetic and elastic components of the error. One
approach to rectify this is to measure the error in the total mechanical
energy of the reconstruction by

$E_{E}=\frac{\sum_{i=1}^{m}|{\mathcal{E}}({\boldsymbol{q}}_{i})-{\mathcal{E}}({\boldsymbol{\hat{q}_{i}}})|}{\sum_{i=1}^{m}{\mathcal{E}}({\boldsymbol{q}}_{i})}.$
(5.56)

However, this notion of error can be zero even when the reconstructed and
ground truth snapshots are different. Both issues are resolved by measuring
the error in the kernel Hilbert space associated with the ground-truth stress
model by

$E_{H}=\frac{\sum_{i=1}^{m}\|\Phi({\boldsymbol{q}}_{i})-\Phi({\boldsymbol{\hat{q}}}_{i})\|_{{\mathcal{H}}}^{2}}{\sum_{i=1}^{m}\|\Phi({\boldsymbol{q}}_{i})\|_{{\mathcal{H}}}^{2}}=\frac{\sum_{i=1}^{m}d_{{\mathcal{E}}}({\boldsymbol{q}}_{i},{\boldsymbol{\hat{q}}}_{i})^{2}}{\sum_{i=1}^{m}{\mathcal{E}}({\boldsymbol{q}}_{i})},$
(5.57)

where quantities in the numerator and denominator are computed using Eq.
(3.12), Eq. (3.11), and the kernel function.

In the following two subsections, we study how the choice of kernel affects
the ability of kernel-based dimensionality reduction and reconstruction
methods to capture essential features of the lid-driven cavity flow. To do
this, we simulate the cavity flow using different choices of ground-truth
stress model, Reynolds number, and Weissenberg number. These ground-truth
quantities are used to define the error metrics $E_{E}$ and $E_{H}$ in Eqs.
5.56 and 5.57. Kernel PCA-based dimensionality reduction and reconstruction
are then performed using the ground-truth kernel function and compared to
other ad-hoc choices such as ordinary PCA and kernel functions associated with
incorrect stress models.

###  5.2 Results using the Oldroyd-B model

We begin by studying the case where the linear Oldroyd-B model is used to
define the ground-truth, and we compare against ordinary PCA. This model holds
significant popularity in non-Newtonian fluid mechanics, making it an ideal
starting point for making comparisons. We note that KPCA using the Oldroyd-B
stress model is equivalent to performing PCA using properly weighted modified
state vectors ${\boldsymbol{\tilde{y}}}=\Psi({\boldsymbol{q}})$ defined by the
explicit feature map in Eq. (3.20). To compare, we naïvely perform PCA using
the state vectors defined by Eq. 5.48 without adjusting how the kinetic and
elastic components are weighted. These two methods become nearly identical
when $\theta=1$, with the only difference being that $b_{xy}$ and $b_{yx}$
both appear in the properly weighted modified state vector, while only
$b_{xy}$ appears in the modified state vector used to perform PCA. The results
discussed below show that using the correct weighting, i.e., the kernel
function associated with the Oldroyd-B model is crucial for PCA to properly
capture the most energetic flow structures.

In Fig. 3, we plot our three error metrics against the number of kernel
principal components (modes) used for reconstruction at different values of
$\theta$. These values were calculated using the triples
$(\mbox{{Re}}=10,\mbox{{Wi}}=10,\beta=0.9)$,
$(\mbox{{Re}}=1,\mbox{{Wi}}=0.5,\beta=0.5)$, and
$(\mbox{{Re}}=0.001,\mbox{{Wi}}=0.5,\beta=0.5)$. Ordinary PCA minimizes the
reconstruction error with respect to the Frobenius norm, so it is no surprise
that its error with respect to $E_{F}$ is lower. However, this norm does not
reflect the energetic content of the flow. When $\theta\neq 1$, reduction and
reconstruction based on the correct Oldroyd-B kernel yields smaller errors
both in the kernel Hilbert space, as measured by $E_{H}$, and in the
mechanical energy of the reconstructed snapshots, as measured by $E_{E}$. As
expected, the two methods have nearly identical performance in all error
metrics when $\theta=1$, as discussed above.

a)
b)
c)

Figure 3: Comparing reconstruction performance using PCA and Kernel PCA (KPCA)
with the Oldroyd-B kernel for simulations performed using the Oldroyd-B stress
model for different values of $\theta$ in rows (a) $\theta=0.001$, (b)
$\theta=1$ and (c) $\theta=1000$. Error is measured using Eq. 5.55 (left), Eq.
5.57 (center), and Eq. 5.56 (right).

Now we select the leading $r=2$ modes and compare how well the transient
dynamics of the flow are captured using the Oldroyd-B kernel and ordinary PCA-
based methods. In Fig. 4 we plot the ground truth mechanical energy and the
mechanical energies of our $2$-mode reconstructions across four values of
$\theta$. In each case, the reduction and reconstruction of the flow using the
correct Oldroyd-B kernel accurately captures the total mechanical energy of
the flow at all times, whereas the naïve PCA-based method fails to captured
the total energy. The error using PCA is most profound for the small $\theta$
cases shown in Fig. 4(a) and Fig. 4(b).

a)
b)
c)
d)

Figure 4: Total mechanical energy ((2.7)) as function of time, reconstructed
using $r=2$ principal components (modes): a) $\theta=0.001$, b) $\theta=0.01$,
c) $\theta=100$ and d) $\theta=1000$.

Finally, figure 5 provides a comparison between simulation and $2$-mode
reconstructions for distribution of mechanical energy in a fixed snapshot at
$t=2.5$ for the $\theta=0.001$ case. Consistent with the previous discussions
concerning the error plots in Fig. 4(a), the reconstructed total mechanical
energy obtained using the Oldroyd-B kernel exhibits an exceptional agreement
with the simulation, even when only two modes are retained. On the other hand,
the naïve PCA-based reconstruction fails to capture the distribution of
mechanical energy in the snapshot.

Figure 5: Spatial distribution of the total mechanical energy ((2.7)) at time
$t=2.5$ for $\theta=0.001$. Results obtained by the simulation (left), PCA
(center) and KPCA (right). The reconstructed values were obtained using $r=2$
modes.

###  5.3 Results using nonlinear stress models

This section explores the use of nonlinear stress models and their associated
kernel functions. We study whether, and to what extent KPCA-based
dimensionality reduction and reconstruction benefits from choosing the kernel
function in accordance with the underlying stress model used to generate the
simulation data. To do this, we fix $\theta=1$ by setting $\mbox{{Re}}=1$,
$\mbox{{Wi}}=0.5$, and $\beta=0.5$, and we simulate the lid-driven cavity flow
using a variety of nonlinear stress models across a range of parameter values
modulating the degree of nonlinearity. As above, we compare the performance of
KPCA-based dimensionality reduction and reconstruction using the kernel
function matching the underlying simulation against kernel functions that do
not. In particular, we are interested in when a nonlinear kernel function is
necessary, or whether there are parameter regimes where we can simply use the
Oldroyd-B kernel, which has the benefit of an explicit bijective feature map
given by Eq. (3.20). The error $E_{H}$ in Eq. 5.57 is always computed using
the kernel associated with the underlying stress model.

In Fig. 6, the simulation was performed using the FENE-P stress model with
parameter $L^{2}=5$, and we compare the corresponding kernel against
Oldroyd-B. While the Oldroyd-B kernel yields better reconstructions in the
$L^{2}(\Omega)$ sense measured by $E_{F}$, the reconstructions using small
numbers of kernel principal components, or “modes”, are worse than those
obtained using the FENE-P kernel function in the energetic sense measured by
$E_{H}$ and $E_{E}$.

a)
b)

Figure 6: Results for the $\theta=1$ case using the nonlinear FENE-P model
with $L^{2}=5$. Row (a) shows reconstruction errors computed using Eq. 5.55
(left), Eq. 5.57 (center) and Eq. 5.56 (right). Row (b) shows the total
mechanical energy ((2.7)) as function of time for fields reconstructed using
$r=2$ kernel principal components (modes).

We repeat this experiment in Fig. 7 using the nonlinear PTT stress model with
parameter $\varepsilon=0.5$ for the simulation. A similar trend is observed;
matching the kernel function to the stress model is beneficial for
reconstruction using small numbers of modes in the energetic sense measured by
$E_{H}$ and $E_{E}$. This improvement is dramatic when $r=2$ modes are used
for reconstruction, as illustrated by the time histories of total mechanical
energy plotted in Fig. 7(b). However, for intermediate numbers of modes, the
simple Oldroyd-B kernel slightly outperforms the nonlinear PTT kernel in all
error metrics.

a)
b)

Figure 7: Analogue of Fig. 6 using the nonlinear PTT model (denoted nl-PTT)
with $\varepsilon=0.5$.

In Fig. 8 we present analogous results for the Giesekus model with $\alpha=1$.
Again, for small numbers of modes, using the corresponding Giesekus kernel
function tends to yield slight improvements over the Oldroyd-B kernel in the
energetic sense measured by $E_{H}$ and $E_{E}$. The performance using the two
kernel functions is nearly the same when larger numbers of modes are used to
reconstruct.

a)
b)

Figure 8: Analogue of Fig. 6 using the Giesekus model with $\alpha=1$.

We study the affect of the parameters in the FENE-P and nonlinear PTT stress
models in Fig. 9. The number of kernel principal components (modes) used for
reconstruction is fixed at $r=2$. These nonlinear stress models become
equivalent to the Oldroyd-B model in the limits $1/L^{2}\to 0$ and
$\varepsilon\to 0$. As expected, the choice of kernel is unimportant at small
values of these parameters and becomes more significant at larger values. In
both cases, it becomes important to use the kernel function matching the
nonlinear stress model of the simulation only when the degree of nonlinearity
is sufficiently high. We do not present analogous results using the Giesekus
model due to numerical challenges associated with simulating the flow at
larger values of the parameter $\alpha$.

a)
b)

Figure 9: The effect of stress model nonlinearity on reconstruction
performance using $r=2$ kernel principal components. Simulations in row (a)
use the FENE-P model with different values of $\frac{1}{L_{2}}$ and
simulations in row (b) use the nonlinear PTT model with different values of
$\varepsilon$. In both cases we have $\theta=1$ and we compare to results
obtained using the Oldroyd-B kernel.

Finally, two-mode reconstructions of the instantaneous total mechanical energy
fields are shown in Fig. 10 at time $t=0.6$ for a simulation performed using
the nonlinear PTT model with parameter $\varepsilon=0.3$. In this case with
$\theta=1$, the dynamic variation is concentrated close to the lid (top of the
spatial domain), which is in contrast to the global variations observed in
Fig.5 at $\theta=0.001$. Here, we see that using the matching kernel yields a
qualitatively more accurate reconstruction of the energy field than using the
Oldroyd-B kernel. We also plot the energy field along horizontal slices at
$y=0.475$ and $y=0.975$, further illustrating that the reconstructions using
the nonlinear PTT kernel function are in closer agreement with the ground
truth than reconstructions obtained using the Oldroyd-B kernel.

a)
b)

Figure 10: Spatial distribution of the total mechanical energy ((2.7)) for a
simulation using the nonlinear PTT model with $\varepsilon=0.3$ and $\theta=1$
at time $t=0.6$. Row (a) compares reconstructed snapshots using $r=2$ kernel
principal components obtained using the Oldroyd-B and nonlinear PTT kernels to
the ground truth. Row (b) shows the mechanical energy distributions across
horizontal slices located at $y=0.475$ (left) and $y=0.975$ (right).

##  6 Discussion

In summary, our theoretical results provide easily verifiable conditions on a
given viscoelastic stress model ensuring that there is a corresponding
positive-definite kernel function compatible with the total mechanical energy
and turning the space of flowfields with finite mechanical energies into a
complete, separable metric space. Moreover, this kernel function is
constructed explicitly from a convergent series representation of the stress
model. Remarkably, these conditions hold for many standard viscoelastic stress
models, yielding principled choices for kernel functions which give rise to
natural measures for distances and angles between flowfield snapshots. These
geometric notions correspond to implicitly-defined lifted representations of
snapshots in a unique reproducing kernel Hilbert space (RKHS) associated with
the kernel function.

The kernel functions and metrics for viscoelastic flows introduced in this
paper enable a variety of machine learning algorithms to be employed to
extract statistical information from flowfield snapshot data. In this paper we
focus primarily on dimensionality reduction using kernel principal component
analysis (KPCA). Here, a major challenge known as the “preimage problem” is to
reconstruct flowfields from their RKHS representations or truncated
representations in terms of the leading kernel principal components. We
provide a solution to the preimage problem for our viscoelastic kernel
functions by showing that the velocity field and square root conformation
tensor field can be linearly reconstructed from RKHS representations. We
provide bounds on the reconstruction quality using truncated representations
given by the leading kernel principal components. More generally, our results
describe which fields can be linearly reconstructed based on the terms in a
series expansion for a given viscoelastic stress model.

We demonstrate the utility of our kernel functions for dimensionality
reduction and reconstruction of snapshots from a lid-driven cavity flow. Here,
the flow is simulated with different underlying stress models and we compare
the reconstructions obtained using the leading kernel principal components
extracted with different choices of kernel function. Our results underscore
the importance of choosing an appropriate metric for reconstruction error.
Ordinary principal component analysis (PCA), which corresponds to a kernel
function given by the $L^{2}$ inner product on the spatial domain, leads to
low reconstruction error in an $L^{2}$ sense, but produces poor
reconstructions as measured by the total mechanical energy. For simulations
performed using the (linear) Oldroyd-B stress model, superior reconstructions
in an energetic sense were obtained using the kernel principal components
extracted using the Oldroyd-B kernel function. We note that KPCA using the
Oldroyd-B kernel is equivalent to ordinary PCA using properly weighted state
vectors based on the square root conformation tensor field.

For simulations performed using nonlinear stress models such as FENE-P and
nonlinear PTT, we compare the reduction and reconstruction performance using
the corresponding kernel functions to the Oldroyd-B kernel. In each case,
using the appropriate nonlinear kernel function was beneficial for
reconstructing flow features in an energetic sense from a small number of
kernel principal components. The benefit of using a nonlinear kernel function
matching the underlying stress model became more pronounced as the parameters
controlling the nonlinearity of the stress model were increased. However, for
our simple lid-driven cavity flow example, this benefit was less significant
than the improvements made by using the Oldroyd-B kernel over naïve PCA. The
performance of the nonlinear kernel function was comparable to using the
simple Oldroyd-B kernel in many cases where the nonlinearity of the stress
model was low or moderate and when more principal components were employed for
reconstruction. This suggests that the Oldroyd-B kernel function (or
equivalently, appropriately modified state vectors) could be a useful default
choice when processing data from viscoelastic flows with only moderately
nonlinear stress models. However, it is possible that accounting for
nonlinearities of the stress model when selecting the kernel function will be
important for capturing the behavior of more complex viscoelastic flows. This
will be a subject for future work.

An exciting avenue for future work involves using our kernel functions to
build low-dimensional data-driven reduced-order models approximating the
dynamics of viscoelastic fluid flows. These low-dimensional dynamical systems
can then be used for a variety of key scientific and engineering tasks
including qualitative analysis of the flow’s dynamics and bifurcations, state
estimation from limited sensor measurements, real-time forecasting, and
feedback control. Promising approaches could combine the dynamics-informed
features extracted using the kernel covariance balancing reduction using
adjoint snapshots (K-CoBRAS) method [18] with modeling techniques such as
Sparse Identification of Nonlinear Dynamics (SINDy) [51, 52] or variants of
Dynamic Mode Decomposition (DMD) [53, 54, 55]. Other approaches could use our
kernel-induced distance metric to formulate loss functions for autoencoder-
based reduced-order modeling methods such as those introduced in [56, 57, 14].
Our kernels may also be of use for classifying flow regimes using support
vector machines (see [58, 28]).

## Acknowledgements

C.M. Oishi and F.V.G. Amaral would like to thank the financial support given
by Sao Paulo Research Foundation (FAPESP) grants #2013/07375-0, #2021/13833-7,
#2021/07034-4 and #2023/06035-2, and the National Council for Scientific and
Technological Development (CNPq), grants #305383/2019-1 and #307228/2023-1.
The authors acknowledge support from the National Science Foundation AI
Institute in Dynamic Systems (grant number 2112085) and from the Army Research
Office (ARO W911NF-19-1-0045).

## References

  * [1] J. L. Lumey, Stochastic Tools in Turbulence, DOVER PUBN INC, 1970.
  * [2] N. Aubry, On the hidden beauty of the proper orthogonal decomposition, Theoret. Comput. Fluid Dynamics 2 (1991) 339––352.
  * [3] L. Sirovich, Analysis of turbulent flows by means of the empirical eigenfunctions, Fluid Dynamics Research 8 (1-4) (1991) 85.  doi:10.1016/0169-5983(91)90033-F.
URL https://dx.doi.org/10.1016/0169-5983(91)90033-F

  * [4] K. Taira, S. L. Brunton, S. T. M. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, L. S. Ukeiley, Modal analysis of fluid flows: An overview, AIAA Journal (2017) 4013–4041doi:doi.org/10.2514/1.J056060.
  * [5] C. W. Rowley, T. Colonius, R. M. Murray, Model reduction for compressible flows using POD and Galerkin projection, Physica D: Nonlinear Phenomena 189 (1-2) (2004) 115–129.
  * [6] E. J. Parish, F. Rizzi, On the impact of dimensionally-consistent and physics-based inner products for pod-galerkin and least-squares model reduction of compressible flows, Journal of Computational Physics 491 (2023) 112387.  doi:https://doi.org/10.1016/j.jcp.2023.112387.
URL https://www.sciencedirect.com/science/article/pii/S0021999123004825

  * [7] A. A. Kaptanoglu, K. D. Morgan, C. J. Hansen, S. L. Brunton, Physics-constrained, low-dimensional models for magnetohydrodynamics: First-principles and data-driven approaches, Phys. Rev. E 104 (2021) 015206.  doi:10.1103/PhysRevE.104.015206.
URL https://link.aps.org/doi/10.1103/PhysRevE.104.015206

  * [8] K. C. Sockwell, Mass conserving Hamiltonian-structure-preserving reduced order modeling for the rotating shallow water equations discretized by a mimetic spatial scheme, Ph.D. thesis, The Florida State University (2019).
  * [9] S. L. Brunton, B. R. Noack, P. Koumoutsakos, Machine learning for fluid mechanics, Annual Review of Fluid Mechanics 52 (2020) 477–508.
  * [10] M. A. Mendez, Linear and nonlinear dimensionality reduction from fluid mechanics to machine learning, Measurement Science and Technology 34 (2022).
URL https://api.semanticscholar.org/CorpusID:251594521

  * [11] P. Dubois, T. Gomez, L. Planckaert, L. Perret, Machine learning for fluid flow reconstruction from limited measurements, Journal of Computational Physics 448 (2022) 110733.  doi:https://doi.org/10.1016/j.jcp.2021.110733.
URL https://www.sciencedirect.com/science/article/pii/S0021999121006288

  * [12] E. Farzamnik, A. Ianiro, S. Discetti, N. Deng, K. Oberleithner, B. Noack, V. Guerrero, From snapshots to manifolds – a tale of shear flows, Journal of Fluid Mechanics 955 (2023) A34.  doi:10.1017/jfm.2022.1039.
  * [13] M. Ohlberger, S. Rave, Reduced basis methods: Success, limitations and future challenges, in: Proceedings of Algoritmy, 2016, pp. 1–12.
  * [14] S. E. Otto, G. R. Macchio, C. W. Rowley, Learning nonlinear projections for reduced-order modeling of dynamical systems using constrained autoencoders, Chaos 33 (11) (2023).
  * [15] B. Schölkopf, A. Smola, K.-R. Müller, Nonlinear component analysis as a kernel eigenvalue problem, Neural Computation 10 (5) (1998) 1299–1319.  doi:10.1162/089976698300017467.
  * [16] M. Salvador, L. Dedè, A. Manzoni, Non intrusive reduced order modeling of parametrized pdes by kernel pod and neural networks, Computers & Mathematics with Applications 104 (2021) 1–13.  doi:https://doi.org/10.1016/j.camwa.2021.11.001.
URL https://www.sciencedirect.com/science/article/pii/S0898122121003928

  * [17] H. Csala, S. T. M. Dawson, A. Arzani, Comparing different nonlinear dimensionality reduction techniques for data-driven unsteady fluid flow modeling, Physics of Fluids 34 (11) (2022) 117119.  arXiv:https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/5.0127284/16607728/117119\\_1\\_online.pdf, doi:10.1063/5.0127284.
URL https://doi.org/10.1063/5.0127284

  * [18] S. E. Otto, A. Padovan, C. W. Rowley, Model reduction for nonlinear systems by balanced truncation of state and gradient covariance, SIAM J. Sci. Comput. 45 (5) (2023) A2325–A2355.
  * [19] M. A. Alves, P. J. Oliveira, F. T. Pinho, Numerical methods for viscoelastic fluid flows, Annual Review of Fluid Mechanics 53 (2021) 509–541.
  * [20] A. Groisman, V. Steinberg, Elastic turbulence in a polymer solution flow, Nature 405 (2000) 53–55.
  * [21] L. Pan, A. Morozov, C. Wagner, P. E. Arratia, Nonlinear elastic instability in channel flows at low reynolds numbers, Phys. Rev. Lett. 110 (2013) 174502.  doi:10.1103/PhysRevLett.110.174502.
URL https://link.aps.org/doi/10.1103/PhysRevLett.110.174502

  * [22] D. Samanta, Y. Dubief, M. Holzner, C. Schäfer, A. N. Morozov, C. Wagner, B. Hof, Elasto-inertial turbulence, Proceedings of the National Academy of Sciences 110 (26) (2013) 10557–10562.  arXiv:https://www.pnas.org/doi/pdf/10.1073/pnas.1219666110, doi:10.1073/pnas.1219666110.
URL https://www.pnas.org/doi/abs/10.1073/pnas.1219666110

  * [23] Y. Dubief, V. E. Terrapon, B. Hof, Elasto-inertial turbulence, Annual Review of Fluid Mechanics 55 (1) (2023) 675–705.  arXiv:https://doi.org/10.1146/annurev-fluid-032822-025933, doi:10.1146/annurev-fluid-032822-025933.
URL https://doi.org/10.1146/annurev-fluid-032822-025933

  * [24] W. LI, L. XI, M. D. GRAHAM, Nonlinear travelling waves as a framework for understanding turbulent drag reduction, Journal of Fluid Mechanics 565 (2006) 353–362.  doi:10.1017/S0022112006002138.
  * [25] M. D. Graham, Drag reduction and the dynamics of turbulence in simple and complex fluidsa), Physics of Fluids 26 (10) (2014) 101301.  arXiv:https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/1.4895780/13795064/101301\\_1\\_online.pdf, doi:10.1063/1.4895780.
URL https://doi.org/10.1063/1.4895780

  * [26] C. S. Dutcher, S. J. Muller, Effects of moderate elasticity on the stability of co- and counter-rotating Taylor–Couette flows, Journal of Rheology 57 (3) (2013) 791–812.  arXiv:https://pubs.aip.org/sor/jor/article-pdf/57/3/791/16013419/791\\_1\\_online.pdf, doi:10.1122/1.4798549.
URL https://doi.org/10.1122/1.4798549

  * [27] D. G. Thomas, B. Khomami, R. Sureshkumar, Nonlinear dynamics of viscoelastic taylor–couette flow: effect of elasticity on pattern selection, molecular conformation and drag, Journal of Fluid Mechanics 620 (2009) 353–382.  doi:10.1017/S0022112008004710.
  * [28] T. Hofmann, B. Schölkopf, A. J. Smola, kernel methods in machine learning, The Annals of Statistics 36 (3) (2008) 1171–1220.
  * [29] J. T.-Y. Kwok, I. W.-H. Tsang, The pre-image problem in kernel methods, IEEE transactions on neural networks 15 (6) (2004) 1517–1525.
  * [30] P. Honeine, C. Richard, Preimage problem in kernel-based machine learning, IEEE Signal Processing Magazine 28 (2) (2011) 77–88.
  * [31] D. Ritcher, G. Iaccarino, E. S. G. Shaqfeh, Simulations of three-dimensional viscoelastic flows past a circular cylinder at moderate reynolds numbers, Journal of Fluid Mechanics 651 (2010) 415–442.  doi:10.1017/S0022112009994083.
  * [32] U. N. Patel, J. P. Rothstein, Y. Modarres-Sadeghi, Forced oscillations of a cylinder in the flow of viscoelastic fluids, Journal of Fluid Mechanics 975 (2023) A28.  doi:10.1017/jfm.2023.861.
  * [33] F. Hamid, C. Sasmal, Significant influence of fluid viscoelasticity on flow dynamics past an oscillating cylinder, Journal of Fluid Mechanics 975 (2023) A26.  doi:10.1017/jfm.2023.859.
  * [34] V. Dzanic, C. From, E. Sauret, The effect of periodicity in the elastic turbulence regime, Journal of Fluid Mechanics 937 (2022) A31.  doi:10.1017/jfm.2022.103.
  * [35] J. Song, N. Liu, X.-Y. Lu, B. Khomami, Direct numerical simulation of elastic turbulence in the taylor–couette flow: transition pathway and mechanistic insight, Journal of Fluid Mechanics 949 (2022) A49.  doi:10.1017/jfm.2022.801.
  * [36] D. D. Medeiros, H. Notsu, C. M. Oishi, Second-order finite difference approximations of the upper-convected time derivative, SIAM Journal on Numerical Analysis 59 (6) (2021) 2955–2988.  arXiv:https://doi.org/10.1137/20M1364990, doi:10.1137/20M1364990.
URL https://doi.org/10.1137/20M1364990

  * [37] N. Balci, B. Thomases, M. Renardy, C. R. Doering, Symmetric factorization of the conformation tensor in viscoelastic fluid models, Journal of Non-Newtonian Fluid Mechanics 166 (11) (2011) 546–553.
  * [38] A. Lozinski, R. G. Owens, An energy estimate for the oldroyd b model: theory and applications, Journal of Non-Newtonian Fluid Mechanics 112 (2) (2003) 161–176.  doi:https://doi.org/10.1016/S0377-0257(03)00096-X.
URL https://www.sciencedirect.com/science/article/pii/S037702570300096X

  * [39] M. Reed, B. Simon, Functional Analysis, Vol. I of Methods of Modern Mathematical Physics, Academic Press, 1980.
  * [40] A. Berlinet, C. Thomas-Agnan, Reproducing kernel Hilbert spaces in probability and statistics, Springer, 2011.
  * [41] N. Aronszajn, Theory of reproducing kernels, Transaction of the American Mathematical Society 68 (1950) 337–404.
  * [42] H. Q. Minh, V. Murino, Covariances in Computer Vision and Machine Learning, Morgan & Claypool, 2018.
  * [43] R. Sousa, R. Poole, A. Afonso, F. Pinho, P. Oliveira, A. Morozov, M. Alves, Lid-driven cavity flow of viscoelastic liquids, Journal of Non-Newtonian Fluid Mechanics 234 (2016) 129–138.  doi:https://doi.org/10.1016/j.jnnfm.2016.03.001.
URL https://www.sciencedirect.com/science/article/pii/S0377025716300027

  * [44] E. Castillo, R. Codina, First, second and third order fractional step methods for the three-field viscoelastic flow problem, Journal of Computational Physics 296 (2015) 113–137.
  * [45] Y. Chai, J. Ouyang, X. Wang, An efficient stabilized finite element scheme for simulating viscoelastic flows, International Journal for Numerical Methods in Fluids 93 (6) (2021) 1947–1968.
  * [46] R. Comminal, J. Spangenberg, J. H. Hattel, Robust simulations of viscoelastic flows at high weissenberg numbers with the streamfunction/log-conformation formulation, Journal of Non-Newtonian Fluid Mechanics 223 (2015) 37–61.
  * [47] F. Martins, C. Oishi, A. Afonso, M. Alves, A numerical study of the kernel-conformationtransformation for transient viscoelastic fluid flows, Journal of Computational Physics 302 (2015) 653–673.
  * [48] J. Evans, I. Palhares Junior, C. Oishi, Stresses of ptt, giesekus, and oldroyd-b fluids in a newtonian velocity field near the stick-slip singularity, Physics of Fluids 29 (12) (2017).
  * [49] J. Evans, C. Oishi, Transient computations using the natural stress formulation for solving sharp corner flows, Journal of Non-Newtonian Fluid Mechanics 249 (2017) 48–52.
  * [50] J. D. Evans, H. L. França, C. M. Oishi, Application of the natural stress formulation for solving unsteady viscoelastic contraction flows, Journal of Computational Physics 388 (2019) 462–489.
  * [51] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the national academy of sciences 113 (15) (2016) 3932–3937.
  * [52] C. M. Oishi, A. A. Kaptanoglu, J. N. Kutz, S. L. Brunton, Nonlinear parametric models of viscoelastic fluid flows (2023).  arXiv:2308.04405.
  * [53] M. O. Williams, I. G. Kevrekidis, C. W. Rowley, A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition, Journal of Nonlinear Science 25 (6) (2015) 1307–1346.
  * [54] M. O. Williams, C. W. Rowley, I. G. Kevrekidis, A kernel-based method for data-driven Koopman spectral analysis, Journal of Computational Dynamics 2 (2) (2015) 247–265.
  * [55] M. J. Colbrook, The multiverse of dynamic mode decomposition algorithms, arXiv preprint arXiv:2312.00137 (2023).
  * [56] K. Champion, B. Lusch, J. N. Kutz, S. L. Brunton, Data-driven discovery of coordinates and governing equations, Proceedings of the National Academy of Sciences 116 (45) (2019) 22445–22451.
  * [57] K. Lee, K. T. Carlberg, Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders, Journal of Computational Physics 404 (2020).
  * [58] B. E. Boser, I. M. Guyon, V. N. Vapnik, A training algorithm for optimal margin classifiers, in: Proceedings of the fifth annual workshop on Computational learning theory, 1992, pp. 144–152.
  * [59] W. Rudin, Real and Complex Analysis, 3rd Edition, McGraw-Hill Book Co., 1987.
  * [60] R. A. Horn, C. R. Johnson, Topics in Matrix Analysis, Cambridge University Press, 1991.
  * [61] M. Taboga, Properties of the Kronecker product, Published online at https://www.statlect.com/matrix-algebra/Kronecker-product-properties (2021).

##  Appendix A Proof of Theorem 1

The proof of the theorem relies on several preliminary lemmas. First, the
Moore-Aronszajn Theorem is a classical result, stated below for completeness,
allowing us to associate a unique RKHS to any positive-definite kernel
function.

######  Theorem 4 (Moore-Aronszajn [40]).

Let $k:{\mathcal{F}}\times{\mathcal{F}}\to\mathbb{R}$ be a function satisfying
the positive-definiteness condition in Eq. (3.10). Then there is a unique
reproducing kernel Hilbert space ${\mathcal{H}}$ of functions on
${\mathcal{F}}$ whose reproducing kernel is $k$. In particular, the subspace
${\mathcal{H}}_{0}$ consisting of finite linear combinations of elements in
$\\{K_{{\boldsymbol{q}}}\\}_{{\boldsymbol{q}}\in{\mathcal{F}}}$ is dense in
${\mathcal{H}}$ and ${\mathcal{H}}$ is the set of functions that are pointwise
limits of Cauchy sequences in ${\mathcal{H}}_{0}$ with the inner product

$\left\langle\sum_{i=1}^{m}a_{i}K_{{\boldsymbol{q}}_{i}},\
\sum_{j=1}^{n}b_{j}K_{{\boldsymbol{q}}^{\prime}_{j}}\right\rangle_{{\mathcal{H}}_{0}}=\sum_{i=1}^{m}\sum_{j=1}^{n}a_{i}b_{j}k({\boldsymbol{q}}_{i},{\boldsymbol{q}}^{\prime}_{j}).$
(A.58)

The following lemma provides several useful rules for combining positive-
definite kernels to produce new positive-definite kernels.

######  Lemma 1 (Combining kernels [28]).

Let $k_{1},k_{2},\ldots$ be real-valued positive-definite kernel functions on
an arbitrary nonempty set ${\mathcal{X}}$. Then the set of positive-definite
kernels on ${\mathcal{X}}$ is a closed convex cone, that is,

  1. 1.

if $\theta_{1},\theta_{2}\geq 0$, then $\theta_{1}k_{1}+\theta_{2}k_{2}$ is a
positive-definite kernel; and

  2. 2.

if $k(x,x^{\prime}):=\lim_{n\to\infty}k_{n}(x,x^{\prime})$ exists for all
$x,x^{\prime}\in{\mathcal{X}}$, then $k$ is a positive-definite kernel.

The point-wise product $k_{1}k_{2}$ is also a positive-definite kernel.
Consequently, if $a_{0},a_{1},\ldots\geq 0$ are non-negative constants and

$\psi(t)=\sum_{n=0}^{\infty}a_{n}t^{n},$ (A.59)

converges for every $t\in k({\mathcal{X}}\times{\mathcal{X}})$ then $\psi\circ
k$ is a positive-definite kernel.

We also require conditions on the kernel function ensuring that the associated
feature map $\Phi$ is injective. That is, for every pair of distinct states
${\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2}\in{\mathcal{F}}$ we have

$k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{1})-2k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})+k({\boldsymbol{q}}_{2},{\boldsymbol{q}}_{2})=\|\Phi({\boldsymbol{q}}_{1})-\Phi({\boldsymbol{q}}_{2})\|_{{\mathcal{H}}}^{2}>0.$
(A.60)

We call positive-definite kernel functions with the above property “injective
kernels”. In the following lemma, we provide some useful rules for combining
injective kernels.

######  Lemma 2 (Combining injective kernels).

Let $k$ and $k^{\prime}$ be real-valued positive-definite kernel functions on
an arbitrary nonempty set ${\mathcal{X}}$ and suppose that $k$ is injective.
Then the following hold:

  1. 1.

if $\theta_{1}>0$ and $\theta_{2}\geq 0$ then
$\theta_{1}k+\theta_{2}k^{\prime}$ is injective;

  2. 2.

if $p\geq 1$ is odd then $k^{p}$ is injective; and

  3. 3.

if $k$ takes only non-negative values, then $k^{p}$ is injective for every
integer $p\geq 1$.

Consequently, if $k$ is non-negative-valued and
$\psi:k({\mathcal{X}}\times{\mathcal{X}})\to\mathbb{R}$ can be expressed as a
convergent power series

$\psi(t)=\sum_{n=0}^{\infty}a_{n}t^{n},$ (A.61)

with every $a_{n}\geq 0$ and there being some $n\geq 1$ with $a_{n}>0$ then
$\psi\circ k$ is injective.

###### Proof of Lemma 2.

The first statement follows immediately from Eq. (A.60). To prove the second
and third statements, consider an integer $p\geq 1$ and choose two distinct
points $x,y\in{\mathcal{X}}$. Since $t\mapsto t^{p}$ is convex on
$[0,\infty)$, Jensen’s inequality (Theorem 3.3 in [59] with $\mu$ being a sum
of two Dirac measures) yields

$\frac{1}{2}\left[k(x,x)^{p}-2k(x,y)^{p}+k(y,y)^{p}\right]\geq\left(\frac{k(x,x)+k(y,y)}{2}\right)^{p}-k(x,y)^{p}.$
(A.62)

If $k(x,y)\geq 0$ then Eq. (A.60) implies that

$\left(\frac{k(x,x)+k(y,y)}{2}\right)^{p}-k(x,y)^{p}>0$ (A.63)

because $t\mapsto t^{p}$ is strictly monotone increasing on $[0,\infty)$. This
proves the third statement. On the other hand, if $k(x,y)<0$ and $p$ is odd
then $k(x,y)^{p}<0$, which immediately implies that

$k(x,x)^{p}-2k(x,y)^{p}+k(y,y)^{p}>0,$ (A.64)

proving the second statement. To prove the final statement, we observe the sum
defining $\psi$ converges absolutely since all of the terms are nonnegative.
Supposing that $a_{m}>0$ for some $m\geq 1$ we have

$\psi\circ k=a_{m}k^{m}+\sum_{n\neq m}a_{n}k^{n},$ (A.65)

where the first term is an injective kernel by the argument above. Since the
sum in the second term converges, it defines a positive-definite kernel by
Lemma 1. Therefore the sum of the two terms is an injective kernel. ∎

The following lemma allows us to convert terms such as
$\operatorname{tr}({\boldsymbol{c}})$ and
$\operatorname{tr}\big{[}({\boldsymbol{c}}-{\boldsymbol{I}})^{2}\big{]}$
appearing in the energy function into injective kernel functions.

######  Lemma 3.

Let $D\subset\mathbb{S}^{d}$ be a set of symmetric matrices and let $f$ be a
real-valued function on $\sigma(D):=\bigcup_{{\boldsymbol{c}}\in
D}\sigma({\boldsymbol{c}})$. With the action of this function on a matrix
${\boldsymbol{c}}\in D$ defined by Eq. (3.13),

$k({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})=\operatorname{tr}\left[f({\boldsymbol{c}}_{1})f({\boldsymbol{c}}_{2})\right]$
(A.66)

is a positive-definite kernel function on $D\times D$. If $f$ is injective on
$\sigma(D)$, then the kernel function is injective.

###### Proof of Lemma 3.

Positive-definiteness is obvious. If $f$ is injective on $\sigma(D)$ then its
inverse $f^{-1}$ can be defined on $f(\sigma(D))$. Since $\sigma(f(D))\subset
f(\sigma(D))$, we have $f^{-1}(f({\boldsymbol{c}}))={\boldsymbol{c}}$ for
every ${\boldsymbol{c}}\in D$, meaning that $f$ is injective on $D$. If
${\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2}$ are distinct elements in $D$, then
$f({\boldsymbol{c}}_{1})\neq f({\boldsymbol{c}}_{2})$ and we have

$k({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{1})-2k({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})+k({\boldsymbol{c}}_{2},{\boldsymbol{c}}_{2})=\|f({\boldsymbol{c}}_{1})-f({\boldsymbol{c}}_{2})\|_{F}^{2}>0.$
(A.67)

This proves that $k$ is an injective kernel function. ∎

Finally, we are ready to prove Theorem 1.

###### Proof of Theorem 1.

The series in Eq. (3.15) converges absolutely because

$\displaystyle\sum_{i,p=0}^{\infty}c_{i,p}\left|\operatorname{tr}\left(f_{i}({\boldsymbol{c}}_{1})f_{i}({\boldsymbol{c}}_{2})\right)\right|^{p}$
$\displaystyle\leq\sum_{i,p=0}^{\infty}\sqrt{c_{i,p}\operatorname{tr}\big{(}f_{i}({\boldsymbol{c}}_{1})f_{i}({\boldsymbol{c}}_{1})\big{)}^{p}}\sqrt{c_{i,p}\operatorname{tr}\big{(}f_{i}({\boldsymbol{c}}_{2})f_{i}({\boldsymbol{c}}_{2})\big{)}^{p}}$
(A.68)
$\displaystyle\leq\sqrt{\sum_{i,p=0}^{\infty}c_{i,p}\operatorname{tr}\big{(}f_{i}({\boldsymbol{c}}_{1})f_{i}({\boldsymbol{c}}_{1})\big{)}^{p}}\sqrt{\sum_{j,q=0}^{\infty}c_{j,q}\operatorname{tr}\big{(}f_{j}({\boldsymbol{c}}_{2})f_{j}({\boldsymbol{c}}_{2})\big{)}^{q}},$

thanks to two applications of the Cauchy-Schwarz inequality. The two terms in
the rightmost product are finite by assumption. By Lemma 1 and Lemma 3 it
follows that $\tilde{k}$ defined by Eq. (3.16) is a positive-definite kernel
function on $D({\boldsymbol{s}})$. It is then easy to see that Eq. (3.16)
satisfies the positive-definiteness condition in Eq. (3.10), and is therefore
a positive-definite kernel on ${\mathcal{F}}$. Since we have

${\mathcal{E}}({\boldsymbol{q}})=\frac{1}{2}\int_{\Omega}\left[\big{|}{\boldsymbol{u}}({\boldsymbol{x}})\big{|}^{2}+\theta
h({\boldsymbol{c}})\right]\operatorname{\mathrm{d}}{\boldsymbol{x}},$ (A.69)

where
$h({\boldsymbol{c}}):=-\mbox{{Wi}}\cdot\operatorname{tr}\big{(}{\boldsymbol{s}}({\boldsymbol{c}})\big{)}+c=\tilde{k}({\boldsymbol{c}},{\boldsymbol{c}})$,
it follows that
${\mathcal{E}}({\boldsymbol{q}})=k({\boldsymbol{q}},{\boldsymbol{q}})$.

Suppose that there is a coefficient $c_{i,p}>0$ with $p\geq 1$, $f_{i}$
injective, and $p$ odd or $f_{i}$ nonnegative. We first show that $\tilde{k}$
is an injective kernel function on $D({\boldsymbol{s}})$. The case when $p$ is
odd follows immediately from Lemma 2(i, ii), and Lemma 3. In the case that
$f_{i}$ is nonnegative, every $f_{i}({\boldsymbol{c}})$ is a positive semi-
definite matrix and it follows that the kernel function

$({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})\mapsto\operatorname{tr}\big{[}f_{i}({\boldsymbol{c}}_{1})f_{i}({\boldsymbol{c}}_{2})\big{]}=\operatorname{tr}\left[\sqrt{f_{i}({\boldsymbol{c}}_{2})}f_{i}({\boldsymbol{c}}_{1})\sqrt{f_{i}({\boldsymbol{c}}_{2})}\right]$
(A.70)

takes only nonnegative values on $\mathbb{S}^{d}\times\mathbb{S}^{d}$. By
Lemma 2(i, iii) and Lemma 3 it follows that $\tilde{k}$ is an injective kernel
function.

Using injectivity of $\tilde{k}$, we prove that $k$ is an injective kernel
function on ${\mathcal{F}}$. We observe that

$|{\boldsymbol{u}}_{1}|^{2}-2{\boldsymbol{u}}_{1}\cdot{\boldsymbol{u}}_{2}+|{\boldsymbol{u}}_{2}|^{2}\geq
0$ (A.71)

by the Cauchy-Schwarz inequality for the dot product on $\mathbb{R}^{d}$ and

$\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{1})-2\tilde{k}({\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2})+\tilde{k}({\boldsymbol{c}}_{2},{\boldsymbol{c}}_{2})\geq
0$ (A.72)

by the Cauchy-Schwarz inequality for the positive-definite kernel $\tilde{k}$.
Thus, if

$k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{1})-2k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})+k({\boldsymbol{q}}_{2},{\boldsymbol{q}}_{2})=0,$
(A.73)

then for almost every ${\boldsymbol{x}}\in\Omega$ we have

$\big{|}{\boldsymbol{u}}_{1}({\boldsymbol{x}})\big{|}^{2}-2{\boldsymbol{u}}_{1}({\boldsymbol{x}})\cdot{\boldsymbol{u}}_{2}({\boldsymbol{x}})+\big{|}{\boldsymbol{u}}_{2}({\boldsymbol{x}})\big{|}^{2}=0,$
(A.74)

which implies that
${\boldsymbol{u}}_{1}({\boldsymbol{x}})={\boldsymbol{u}}_{2}({\boldsymbol{x}})$,
and

$\tilde{k}\big{(}{\boldsymbol{c}}_{1}({\boldsymbol{x}}),{\boldsymbol{c}}_{1}({\boldsymbol{x}})\big{)}-2\tilde{k}\big{(}{\boldsymbol{c}}_{1}({\boldsymbol{x}}),{\boldsymbol{c}}_{2}({\boldsymbol{x}})\big{)}+\tilde{k}\big{(}{\boldsymbol{c}}_{2}({\boldsymbol{x}}),{\boldsymbol{c}}_{2}({\boldsymbol{x}})\big{)}=0,$
(A.75)

which implies that
${\boldsymbol{c}}_{1}({\boldsymbol{x}})={\boldsymbol{c}}_{2}({\boldsymbol{x}})$.
Therefore, ${\boldsymbol{q}}_{1}={\boldsymbol{q}}_{2}$ in ${\mathcal{F}}$,
proving that $k$ is an injective kernel function. ∎

##  Appendix B Proof of Theorem 2

Our proof relies on the following technical lemma.

######  Lemma 4.

Let $x,y$ be vectors in a Hilbert space. If $p\geq 1$ is odd or $\langle
x,y\rangle\geq 0$, then

$\langle x,x\rangle^{p}-2\langle x,y\rangle^{p}+\langle
y,y\rangle^{p}\geq\frac{1}{2^{2(p-1)}}\|x-y\|^{2p}.$ (B.76)

###### Proof of Lemma 4.

When $x=y$, the statement is vacuously true, so we assume that $x\neq y$ and
denote

$\alpha=\frac{2\langle x,y\rangle}{\|x-y\|^{2}}.$

By Jensen’s inequality (Theorem 3.3 in [59] with $\mu$ being a sum of two
Dirac measures) we have

$\langle x,x\rangle^{p}-2\langle x,y\rangle^{p}+\langle y,y\rangle^{p}\geq
2\left[\left(\frac{\langle x,x\rangle+\langle
y,y\rangle}{2}\right)^{p}-\langle
x,y\rangle^{p}\right]=2\left[\left(\frac{1}{2}\|x-y\|^{2}+\langle
x,y\rangle\right)^{p}-\langle x,y\rangle^{p}\right].$ (B.77)

Dividing through by $2(\frac{1}{2}\|x-y\|^{2})^{p}$ gives

$\frac{\langle x,x\rangle^{p}-2\langle x,y\rangle^{p}+\langle
y,y\rangle^{p}}{2^{1-p}\|x-y\|^{2p}}\geq(1+\alpha)^{p}-\alpha^{p}=:f(\alpha),$
(B.78)

and so it remains to lower bound $f(\alpha)$ by a positive constant,
specifically $2^{1-p}$. When $\langle x,y\rangle\geq 0$, we have $\alpha\geq
0$, and it is easy to see that $f(\alpha)\geq 1\geq 2^{1-p}$.

Now we assume that $p\geq 1$ is odd. Differentiating $f(\alpha)$, we find

$f^{\prime}(\alpha)=p(1+\alpha)^{p-1}-p\alpha^{p-1}=p|1+\alpha|^{p-1}-p|\alpha|^{p-1}$
(B.79)

since $p-1$ is even. When $\alpha\geq-\frac{1}{2}$, we evidently have
$f^{\prime}(\alpha)\geq 0$. Likewise, when $\alpha\leq-\frac{1}{2}$ we have
$f^{\prime}(\alpha)\leq 0$. Therefore, for every real $\alpha$ we have

$f(\alpha)\geq f\left(-\tfrac{1}{2}\right)=2(\tfrac{1}{2})^{p}=2^{1-p},$
(B.80)

which completes the proof of the lemma. ∎

With this lemma in hand, we are ready to prove the theorem.

###### Proof of Theorem 2.

Thanks to Theorem 1, ${\mathcal{F}}$ is a metric space with metric
$d_{{\mathcal{E}}}$ defined by Eq. (3.12). Let
$\\{{\boldsymbol{q}}_{n}=({\boldsymbol{u}}_{n},{\boldsymbol{c}}_{n})\\}_{n=1}^{\infty}\in{\mathcal{F}}$
be a Cauchy sequence. Since $L^{2}(\Omega)$ is complete (see Theorem 3.11 in
[59]) and

$d_{{\mathcal{E}}}({\boldsymbol{q}}_{m},{\boldsymbol{q}}_{n})\geq\frac{1}{\sqrt{2}}\|{\boldsymbol{u}}_{m}-{\boldsymbol{u}}_{n}\|_{L^{2}(\Omega)},$
(B.81)

it follows that ${\boldsymbol{u}}_{n}\to{\boldsymbol{u}}$ in $L^{2}(\Omega)$
for a unique velocity field ${\boldsymbol{u}}$.

Considering the term $c_{i,p}>0$ and letting
$\theta=(1-\beta)/(\mbox{{Re}}\mbox{{Wi}})$, we have

$\left(\frac{2}{c_{i,p}\theta}\right)d_{{\mathcal{E}}}({\boldsymbol{q}}_{m},{\boldsymbol{q}}_{n})^{2}\geq\int_{\Omega}\left\\{\left[\operatorname{tr}\left(f_{i}({\boldsymbol{c}}_{m})^{2}\right)\right]^{p}-2\left[\operatorname{tr}\left(f_{i}({\boldsymbol{c}}_{m})f_{i}({\boldsymbol{c}}_{n})\right)\right]^{p}+\left[\operatorname{tr}\left(f_{i}({\boldsymbol{c}}_{n})^{2}\right)\right]^{p}\right\\}\operatorname{\mathrm{d}}{\boldsymbol{x}}\\\
\geq\frac{1}{2^{2(p-1)}}\int_{\Omega}\|f_{i}({\boldsymbol{c}}_{m})-f_{i}({\boldsymbol{c}}_{n})\|_{F}^{2p}\operatorname{\mathrm{d}}{\boldsymbol{x}},$
(B.82)

thanks to Lemma 4. It follows from the completeness theorem for
$L^{2p}(\Omega)$, specifically Theorems 3.11 and 3.12 in [59], that
$f_{i}({\boldsymbol{c}}_{n})$ converges in $L^{2p}(\Omega)$ to a limit
${\boldsymbol{f}}$ and that there is a subsequence
$\\{f_{i}({\boldsymbol{c}}_{n_{k}})\\}_{k=1}^{\infty}$ converging pointwise
almost everywhere to ${\boldsymbol{f}}$ in $\Omega$. Since the function
$f_{i}$ is injective, it follows that the conformation tensor field
${\boldsymbol{c}}=f_{i}^{-1}({\boldsymbol{f}}):\Omega\to\mathbb{S}_{+}^{d}$
satisfies $f_{i}({\boldsymbol{c}}_{n})\to f_{i}({\boldsymbol{c}})$ in
$L^{2p}(\Omega)$ and
${\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}})\to{\boldsymbol{c}}({\boldsymbol{x}})$
for almost every ${\boldsymbol{x}}\in\Omega$.

It remains to show that ${\boldsymbol{q}}=({\boldsymbol{u}},{\boldsymbol{c}})$
satisfies ${\mathcal{E}}({\boldsymbol{q}})<\infty$ and that
$d_{{\mathcal{E}}}({\boldsymbol{q}}_{n},{\boldsymbol{q}})\to 0$. Both of these
are accomplished by means of Lebesgue’s dominated convergence theorem (Theorem
1.34 in [59]). We pass to a further subsequence, still denoted with indices
$n_{k}$, such that

$d_{{\mathcal{E}}}({\boldsymbol{q}}_{n_{k}},{\boldsymbol{q}}_{n_{k+1}})<2^{-k}.$
(B.83)

We let $\tilde{\Phi}:D({\boldsymbol{s}})\to\tilde{{\mathcal{H}}}$ denote the
feature map associated with the reproducing kernel
$\tilde{k}:D({\boldsymbol{s}})\times D({\boldsymbol{s}})\to\mathbb{R}$. We
define a function $G:\Omega\to\mathbb{R}$ by

$G({\boldsymbol{x}})=\big{\|}\tilde{\Phi}({\boldsymbol{c}}_{n_{1}}({\boldsymbol{x}}))\big{\|}_{\tilde{{\mathcal{H}}}}+\sum_{j=2}^{\infty}\big{\|}\tilde{\Phi}({\boldsymbol{c}}_{n_{j}}({\boldsymbol{x}}))-\tilde{\Phi}({\boldsymbol{c}}_{n_{j-1}}({\boldsymbol{x}}))\big{\|}_{\tilde{{\mathcal{H}}}}.$
(B.84)

We observe that $G\in L^{2}(\Omega)$ because

$\sqrt{\frac{\theta}{2}}\|G\|_{L^{2}(\Omega)}\leq\sqrt{{\mathcal{E}}({\boldsymbol{q}}_{n_{1}})}+\sum_{j=2}^{\infty}d_{{\mathcal{E}}}({\boldsymbol{q}}_{n_{j}},{\boldsymbol{q}}_{n_{j-1}})<\infty.$
(B.85)

Moreover, by construction we have

$\tilde{k}\big{(}{\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}}),{\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}})\big{)}=\big{\|}\tilde{\Phi}({\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}}))\big{\|}_{\tilde{{\mathcal{H}}}}^{2}\leq
G({\boldsymbol{x}})^{2}$ (B.86)

for every $k\geq 1$. Therefore, the dominated convergence theorem (Theorem
1.34 in [59]) yields

$\lim_{k\to\infty}\int_{\Omega}\tilde{k}\big{(}{\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}}),{\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}})\big{)}\operatorname{\mathrm{d}}{\boldsymbol{x}}=\int_{\Omega}\tilde{k}\big{(}{\boldsymbol{c}}({\boldsymbol{x}}),{\boldsymbol{c}}({\boldsymbol{x}})\big{)}\operatorname{\mathrm{d}}{\boldsymbol{x}}\leq\|G\|_{L^{2}(\Omega)}^{2}<\infty.$
(B.87)

Combined with the $L^{2}(\Omega)$ convergence of ${\boldsymbol{u}}_{n_{k}}$ to
${\boldsymbol{u}}$, this gives

$\lim_{k\to\infty}{\mathcal{E}}({\boldsymbol{q}}_{n_{k}})={\mathcal{E}}({\boldsymbol{q}})<\infty,$
(B.88)

meaning that ${\boldsymbol{q}}\in{\mathcal{F}}$. Next we observe that for
almost every ${\boldsymbol{x}}\in\Omega$,

$\left|\tilde{k}({\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}}),{\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}}))-2\tilde{k}({\boldsymbol{c}}_{n_{k}}({\boldsymbol{x}}),{\boldsymbol{c}}({\boldsymbol{x}}))+\tilde{k}({\boldsymbol{c}}({\boldsymbol{x}}),{\boldsymbol{c}}({\boldsymbol{x}}))\right|\\\
\leq 4G({\boldsymbol{x}})^{2}$ (B.89)

by the Cauchy-Schwarz inequality in $\tilde{{\mathcal{H}}}$. Therefore, by
another application of the dominated convergence theorem we obtain

$d_{{\mathcal{E}}}({\boldsymbol{q}}_{n_{k}},{\boldsymbol{q}})\to 0.$ (B.90)

Since ${\boldsymbol{q}}_{n}$ is Cauchy, we obtain
$d_{{\mathcal{E}}}({\boldsymbol{q}}_{n},{\boldsymbol{q}})\to 0$, proving that
${\mathcal{F}}$ is a complete metric space.

Finally, if $\Phi({\boldsymbol{q}}_{n})\to\phi$ in ${\mathcal{H}}$ then
$\\{{\boldsymbol{q}}_{n}\\}_{n=1}^{\infty}$ is a Cauchy sequence in
${\mathcal{F}}$. Since ${\mathcal{F}}$ is a complete metric space, there
exists ${\boldsymbol{q}}\in{\mathcal{F}}$ satisfying

$d_{{\mathcal{E}}}({\boldsymbol{q}}_{n},{\boldsymbol{q}})=\|\Phi({\boldsymbol{q}}_{n})-\Phi({\boldsymbol{q}})\|_{{\mathcal{H}}}\to
0.$ (B.91)

Since limits in a Hilbert space are unique, we must have
$\phi=\Phi({\boldsymbol{q}})$, proving that $\Phi({\mathcal{F}})$ is closed in
${\mathcal{H}}$. ∎

##  Appendix C Proof of Theorem 3 and Proposition 1

First, we establish a lemma relating powers of the trace to the trace of
iterated Kronecker products.

######  Lemma 5.

For every ${\boldsymbol{A}},{\boldsymbol{B}}\in\mathbb{R}^{d\times d}$ and
integer $p\geq 0$ we have

$\left[\operatorname{tr}({\boldsymbol{A}}{\boldsymbol{B}})\right]^{p}=\operatorname{tr}\left[{\boldsymbol{A}}^{\otimes
p}{\boldsymbol{B}}^{\otimes p}\right],$ (C.92)

with the convention that ${\boldsymbol{A}}^{\otimes 0}=1$.

###### Proof of Lemma 5.

The cases $p=0$ and $p=1$ are trivial. Suppose that the result holds for a
given $p>1$. Then we have

$\displaystyle\left[\operatorname{tr}({\boldsymbol{A}}{\boldsymbol{B}})\right]^{p+1}$
$\displaystyle=\operatorname{tr}({\boldsymbol{A}}{\boldsymbol{B}})\operatorname{tr}\left[{\boldsymbol{A}}^{\otimes
p}{\boldsymbol{B}}^{\otimes p}\right]$ (C.93)
$\displaystyle=\operatorname{tr}\left[({\boldsymbol{A}}{\boldsymbol{B}})\otimes\big{(}{\boldsymbol{A}}^{\otimes
p}{\boldsymbol{B}}^{\otimes p}\big{)}\right]$
$\displaystyle=\operatorname{tr}\left[\big{(}{\boldsymbol{A}}\otimes{\boldsymbol{A}}^{\otimes
p}\big{)}\big{(}{\boldsymbol{B}}\otimes{\boldsymbol{B}}^{\otimes
p}\big{)}\right]=\operatorname{tr}\left[{\boldsymbol{A}}^{\otimes(p+1)}{\boldsymbol{B}}^{\otimes(p+1)}\right].$

The first line holds by the induction hypothesis, the second line is due to
the trace property of the Kronecker product, and the third line follows from
the mixed product property of the Kronecker product [60, 61]. Therefore the
stated result holds for all integers $p\geq 0$ by induction on $p$. ∎

###### Proof of Theorem 3.

With a countable set of indices

$\mathbb{I}=\\{(i,p,j,k)\in\mathbb{Z}^{4}\ :\ i\geq 0,\ p\geq 0,\ 1\leq
j,k\leq d^{p}\\}$ (C.94)

we consider the space $L^{2}(\Omega\times\mathbb{I})$ with inner product

$\langle f,\
g\rangle_{L^{2}(\Omega\times\mathbb{I})}:=\sum_{(i,p,j,k)\in\mathbb{I}}\int_{\Omega}f({\boldsymbol{x}},i,p,j,k)g({\boldsymbol{x}},i,p,j,k)\operatorname{\mathrm{d}}{\boldsymbol{x}}.$
(C.95)

Note that the sum and the integral can be exchanged thanks to Fubini’s theorem
(Theorem 8.8 in [59]). Defining $\Psi:{\mathcal{F}}\to H:=L^{2}(\Omega;\
\mathbb{R}^{3})\times L^{2}(\Omega\times\mathbb{I})$ by

$\Psi({\boldsymbol{q}})=\frac{1}{\sqrt{2}}({\boldsymbol{u}},\
\sqrt{\theta}\tilde{\Psi}({\boldsymbol{q}})),\qquad\tilde{\Psi}({\boldsymbol{q}})({\boldsymbol{x}},i,p,j,k)=\sqrt{c_{i,p}}\left[f_{i}({\boldsymbol{c}}(x))^{\otimes
p}\right]_{j,k}$ (C.96)

and applying Lemma 5, we obtain

$k({\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2})=\langle\Psi({\boldsymbol{q}}_{1}),\
\Psi({\boldsymbol{q}}_{2})\rangle_{H}=\frac{1}{2}\int_{\Omega}\left\\{{\boldsymbol{u}}_{1}({\boldsymbol{x}})^{T}{\boldsymbol{u}}_{2}({\boldsymbol{x}})+\theta\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}c_{i,p}\operatorname{tr}\left[f_{i}({\boldsymbol{c}}_{1}({\boldsymbol{x}}))^{\otimes
p}f_{i}({\boldsymbol{c}}_{2}({\boldsymbol{x}}))^{\otimes
p}\right]\right\\}\operatorname{\mathrm{d}}{\boldsymbol{x}}$ (C.97)

for every ${\boldsymbol{q}}_{1},{\boldsymbol{q}}_{2}\in{\mathcal{F}}$.
Therefore, by the Moore-Aronszajn theorem (Theorem 4), there is a unique to a
linear isometry $U:{\mathcal{H}}\to H$ satisfying
$U\Phi({\boldsymbol{q}})=\Psi({\boldsymbol{q}})$ for every
${\boldsymbol{q}}\in{\mathcal{F}}$.

Incidentially, the above shows that the metric space
$({\mathcal{F}},d_{{\mathcal{E}}})$ is separable. This is because
${\mathcal{F}}$ is isometric (via the map $U\Phi)$ to a subset of $H$, which
is a separable Hilbert space.

For every $f\in L^{2}(\Omega\times\mathbb{I})$, let
${\boldsymbol{f}}_{i,p}:\Omega\to\mathbb{R}^{d^{p}\times d^{p}}$ be defined by
$[{\boldsymbol{f}}_{i,p}({\boldsymbol{x}})]_{j,k}=f({\boldsymbol{x}},i,p,j,k)$.
The linear map $T:H\to L^{2}(\Omega)$ defined by

$T({\boldsymbol{u}},f)=\sqrt{2}{\boldsymbol{a}}_{0}^{T}{\boldsymbol{u}}+\sqrt{\frac{2}{\theta}}\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}\operatorname{tr}\left({\boldsymbol{A}}_{i,p}^{T}{\boldsymbol{f}}_{i,p}\right)$
(C.98)

clearly satisfies

$TU\Phi({\boldsymbol{q}})=T\Psi({\boldsymbol{q}})={\boldsymbol{a}}_{0}^{T}{\boldsymbol{u}}+\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}\sqrt{c_{i,p}}\operatorname{tr}\left[{\boldsymbol{A}}_{i,p}^{T}f_{i}({\boldsymbol{c}})^{\otimes
p}\right]=\psi\circ{\boldsymbol{q}}$ (C.99)

for every
${\boldsymbol{q}}=({\boldsymbol{u}},{\boldsymbol{c}})\in{\mathcal{F}}$. The
operator $T$ is bounded because

$\displaystyle\|T({\boldsymbol{u}},f)\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\int_{\Omega}\left|\sqrt{2}{\boldsymbol{a}}_{0}^{T}{\boldsymbol{u}}({\boldsymbol{x}})+\sqrt{\frac{2}{\theta}}\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}\operatorname{tr}\left[{\boldsymbol{A}}_{i,p}^{T}{\boldsymbol{f}}_{i,p}({\boldsymbol{x}})\right]\right|^{2}\operatorname{\mathrm{d}}{\boldsymbol{x}}$
(C.100)
$\displaystyle\leq\int_{\Omega}\left|\sqrt{2}\|{\boldsymbol{a}}_{0}\|_{2}\|{\boldsymbol{u}}({\boldsymbol{x}})\|_{2}+\sqrt{\frac{2}{\theta}}\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}\|{\boldsymbol{A}}_{i,p}\|_{F}\|{\boldsymbol{f}}_{i,p}({\boldsymbol{x}})\|_{F}\right|^{2}\operatorname{\mathrm{d}}{\boldsymbol{x}}$
$\displaystyle\leq\underbrace{\left(2\|{\boldsymbol{a}}_{0}\|_{2}^{2}+\frac{2}{\theta}\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}\|{\boldsymbol{A}}_{i,p}\|_{F}^{2}\right)}_{A^{2}<\infty\mbox{
by
assumption}}\underbrace{\int_{\Omega}\left(\|{\boldsymbol{u}}({\boldsymbol{x}})\|_{2}^{2}+\sum_{i=0}^{\infty}\sum_{p=0}^{\infty}\|{\boldsymbol{f}}_{i,p}({\boldsymbol{x}})\|_{F}^{2}\right)\operatorname{\mathrm{d}}{\boldsymbol{x}}}_{\|({\boldsymbol{u}},f)\|_{H}^{2}}$

thanks to two applications of the Cauchy-Schwarz inequality. It follows that
the linear operator $R_{\psi}:=TU$ is bounded, with operator norm
$\|R_{\psi}\|\leq A$. Since the span of vectors
$\\{\Phi({\boldsymbol{q}})\\}_{{\boldsymbol{q}}\in{\mathcal{F}}}$ is dense in
${\mathcal{H}}$ by the by the Moore-Aronszajn theorem (Theorem 4), the bounded
linear operator $R_{\psi}$ is uniquely defined by the relation
$R_{\psi}\Phi({\boldsymbol{q}})=\psi\circ{\boldsymbol{q}}$ for all
${\boldsymbol{q}}\in{\mathcal{F}}$. ∎

###### Proof of Proposition 1.

Let $P_{r}=U_{r}U_{r}^{*}$ denote the orthogonal projection onto the span of
$u_{1},\ldots,u_{r}$ in ${\mathcal{H}}$. Thanks to Corollary 1, we have

$\displaystyle\left\|\big{(}{\boldsymbol{u}},\
f_{i}({\boldsymbol{c}})^{\otimes
p}\big{)}-R_{i,p}U_{r}{\boldsymbol{z}}_{r}({\boldsymbol{q}})\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\|R_{i,p}\Phi({\boldsymbol{q}})-R_{i,p}P_{r}\Phi({\boldsymbol{q}})\|_{L^{2}(\Omega)}^{2}$
(C.101)
$\displaystyle\leq\|R_{i,p}\|^{2}\|(I-P_{r})\Phi({\boldsymbol{q}})\|_{{\mathcal{H}}}^{2}$
$\displaystyle=\|R_{i,p}\|^{2}\left(\|\Phi({\boldsymbol{q}})\|_{{\mathcal{H}}}^{2}-\|P_{r}\Phi({\boldsymbol{q}})\|_{{\mathcal{H}}}^{2}\right).$

Since
$\|\Phi({\boldsymbol{q}})\|_{{\mathcal{H}}}^{2}={\mathcal{E}}({\boldsymbol{q}})$
and
$\|P_{r}\Phi({\boldsymbol{q}})\|_{{\mathcal{H}}}^{2}=\|U_{r}^{*}\Phi({\boldsymbol{q}})\|_{2}^{2}=\|{\boldsymbol{z}}_{r}({\boldsymbol{q}})\|_{2}^{2},$
we obtain Eq. (4.45) using the bound on $\|R_{i,p}\|$ stated in Corollary 1.
Integrating, applying Parseval’s theorem (Theorem II.6 in [39]), and using
Fubini’s theorem (Theorem 8.8 in [59]) to exchange summation and integration
yields

$\displaystyle\int_{{\mathcal{F}}}\|(I-P_{r})\Phi({\boldsymbol{q}})\|_{{\mathcal{H}}}^{2}\operatorname{\mathrm{d}}\mu({\boldsymbol{q}})$
$\displaystyle=\int_{{\mathcal{F}}}\sum_{j=r+1}^{\infty}\left|\langle u_{j},\
\Phi({\boldsymbol{q}})\rangle_{{\mathcal{H}}}\right|^{2}\operatorname{\mathrm{d}}\mu({\boldsymbol{q}})$
(C.102) $\displaystyle=\sum_{j=r+1}^{\infty}\int_{{\mathcal{F}}}\langle
u_{j},\
\Phi({\boldsymbol{q}})\rangle_{{\mathcal{H}}}\langle\Phi({\boldsymbol{q}}),u_{j}\rangle_{{\mathcal{H}}}\operatorname{\mathrm{d}}\mu({\boldsymbol{q}})$
$\displaystyle=\sum_{j=r+1}^{\infty}\langle u_{j},\
C_{\mu}u_{j}\rangle_{{\mathcal{H}}}=\sum_{j=r+1}^{\infty}\sigma_{j}^{2}.$

Combining this with Eq. (C.101) and the bound on $\|R_{i,p}\|$ in Corollary 1
yeilds Eq. (4.46), completing the proof. ∎