Space-Time-Parameter PCA for Data-Driven Modeling with Application to Bioengineering

Abstract

Principal component analysis is a recognized powerful and practical method in statistics and data science. It can also be used in modeling as a dimensionality reduction tool to achieve low-order models of complex multiphysics or engineering systems. Model-order reduction (MOR) methodologies today are an important topic for engineering design and analysis. Design space exploration or accelerated numerical optimization for example are made easier by the use of reduced-order models. In this chapter, we will talk about the use of higher-order singular value decompositions (HOSVD) applied to spatiotemporal problems that are parameterized by a set of design variables or physical parameters. Here we consider a data-driven reduced order modeling based on a design of computer experiment: from high-dimensional computational results returned by high-fidelity solvers (e.g. finite element ones), the HOSVD allows us to determine spatial, time and parameters principal components. The dynamics of the system can then be retrieved by identifying the low-order discrete dynamical system. As application, we will consider the dynamics of deformable capsules flowing into microchannels. The study of such fluid-structure interaction problems is motivated by the use of microcapsules as innovative drug delivery carriers through blood vessels.

1. Introduction

Manufactured deformable microcapsules are intended to be used as drug carriers within the human vascular network to deliver drugs at specific targets (tumors, etc.). In order to design reliable capsules, one can make help of numerical simulation and high performance computing. The transportation of such capsules into microchannels is a three-dimensional fluid-structure interaction (FSI) problem involving a fluid flow within a confined environment and the deformation of hyperelastic membranes [1, 2]. The behavior of the capsule depends on dimensionless parameters such as the capillary number denoted by $\mathit{Ca}$ and the aspect ratio $a/\mathit{\ell }$ between the capsule radius $a$ and the channel characteristic length $\mathrm{\ell }$. The parameter vector $\mathit{\mu }=\left(\mathit{Ca},a/\mathrm{\ell }\right)$, for which a capsule steady shape exists, lies in a bounded domain $D\subset {\mathbb{R}}^2$. We look for the time evolution of the capsule shape under a Lagrangian description. From an initial shape $\mathit{X}$, we are interested in determining the capsule position $\mathit{x}\left(\mathit{X},t,\mathit{\mu }\right)$ in the microchannel domain $\mathrm{\Omega }\in {\mathbb{R}}^3$ at time $t$ for a parameter vector $\mathit{\mu }\in D$. By denoting $\mathit{u}$ the displacement vector from the initial position, we have

$$\mathit{x}\left(\mathit{X},\mathit{\mu },t\right)=\mathit{X}+\mathit{u}\left(\mathit{X},\mathit{\mu },t\right)$$E1

with $\mathit{u}\left(\mathit{X},\mathit{\mu },0\right)=\mathbf{0}$. The governing equations of the FSI problem include both kinematics and motion equations. At the membrane, we have equilibrium of the mechanical forces (mechanical equilibrium of the membrane and viscous stresses from the fluid). By denoting $\mathit{v}$ the vector field of velocity at the membrane, the system of differential algebraic equations in abstract form reads

$$\dot{\mathit{u}}=\mathit{v},$$E2

$$\mathit{v}={\mathcal{F}}_{\mathit{\mu }}\left(\mathit{u},\mathit{X}\right).$$E3

Practically, there are different candidate computational approaches to discretize this system of equations. First, the initial capsule membrane has to be discretized by using a finite element triangular mesh made of nodes ${\left\{{\mathit{X}}_i\right\}}_{i=1}^I$. Regarding time discretization, in [2], an explicit time scheme is used and the velocity field is determined by the use of a boundary integral method (BIM) coupled with a finite element method (FEM). A numerical stability condition imposes the use of small time steps. For a given parameter $\mathit{\mu }$, the time evolution of capsule dynamics on the time intervals of interest generally requires hours of CPU time. To better understand the membrane behavior with respect to $\mathit{\mu }$, a design of computer experiment (DoCE) is done: from a set of $J$ parameter samples of ${\mathit{\mu }}_j\in D$, $j=1,\dots ,J$, a spatio-temporal solution is computed for each ${\mathit{\mu }}_j$, leading to a database of shape solutions under the form of a third-order tensor

$${\mathcal{T}}_x={\left(\mathit{x}\left({\mathit{X}}_i,{\mathit{\mu }}_j,t^n\right)\right)}_{i=1,\dots ,I,j=1,\dots ,J,n=0,\dots ,N,}\in {\mathbb{R}}^{3I\times J\times \left(N+1\right)}$$E4

using a triangular finite element discretization of the membrane, a time discretization $t^n=n\mathrm{\Delta }t$ (assuming that the time step is constant) and the parameter samples ${\mathit{\mu }}_j$. Typically, for practical computations, $I=O\left(1000\right)$, $N=O\left(10000\right)$ and $J=O\left(100\right)$, so that the tensor database becomes rather huge (about $O\left(10\right)$ gigabytes). Of course, one can only store the solutions at coarser times steps and reduce $N$ to $O\left(100\right)$ but the database remains rather big even in this case.

From this data tensor, one can imagine different use cases leading to different tools:

1. Data exploration and knowledge extraction;

2. Real-time rendering of capsule dynamics for better understanding;

3. Data-driven modeling of capsule dynamics in the whole parameter domain.

First and second items can be achieved by means of data dimensionality reduction. This leads to a lower storage of data in memory as well as a lower numerical complexity of processing and manipulation. In this chapter, we will consider a Higher Order Singular Value Decomposition (HOSVD), which is a generalization of Principal Analysis Component (PCA) to tensors. The third item involves a model-order reduction (MOR) methodology. From computed data, we would like to derive a lightweight dynamical system that reproduces the data and, even more, that is able to accurately estimate solutions for different parameter values $\mathit{\mu }$. Data-driven model-order reduction first makes use of data-dimensionality reduction by a low-order tensor decomposition of the solutions according to some suitable spatial, temporal and parameter reduced bases (see [3]). In our application, this will give the truncated decomposition (where the solutions are here seen as functions):

$$\tilde{\mathit{x}}\left(\mathit{X},\mathit{\mu },t\right)=\mathit{X}+\sum _{k=1}^K\sum _{\mathrm{\ell }=1}^L\sum _{m=1}^M{a}_{k\mathrm{\ell }m}{\mathit{\phi }}_k\left(\mathit{X}\right){\psi }_{\mathrm{\ell }}\left(\mathit{\mu }\right){\omega }_m\left(t\right)$$E5

for some expansion coefficients $a_{k\mathrm{\ell }m}$ and some spatial functions ${\mathit{\phi }}_k$, parameter functions ${\psi }_{\mathrm{\ell }}$ and temporal functions ${\omega }_m$ with ${\omega }_m\left(0\right)=0$ ensuring $\mathit{x}\left(\mathit{X},\mathit{\mu },0\right)=\mathit{X}$. The truncation ranks $K$, $L$ and $M$ are expected to be small enough ($K\ll I$, $L\ll J$ and $M\ll N$). Discretized shape solutions returned by the Full-Order computational Model (FOM) are stored into a third-order tensor of data ${\mathcal{T}}_x$. Let ${\tilde{\mathcal{T}}}_x$ denote the truncated tensor expansion related to (5). It reads:

$${\tilde{\mathcal{T}}}_x=A^0\otimes {\mathit{e}}^N+\sum _{k=1}^K\sum _{\mathrm{\ell }=1}^L\sum _{m=1}^M{a}_{k\mathrm{\ell }m}{\mathbf{\Phi }}_k\otimes {\mathbf{\Psi }}_{\mathrm{\ell }}\otimes {\mathbf{w}}_m\approx {\mathcal{T}}_x$$E6

where $A^0\in {\mathbb{R}}^{I\times J}$, ${\left(A^0\right)}_{\mathit{ij}}={\mathit{X}}_i$, ${\mathit{e}}^N={\left(1,\dots ,1\right)}^T\in {\mathbb{R}}^N$, ${\mathbf{\Phi }}_k\in {\mathbb{R}}^I$, ${\mathbf{\Psi }}_{\mathrm{\ell }}\in {\mathbb{R}}^J$ and ${\mathit{w}}_m\in {\mathbb{R}}^N$. In this chapter, we will seamlessly use the functional representation or the tensor one.

From this reduced form, we will then apply a kernel-based Dynamic Mode Decomposition ($k$-DMD, see [4]) to derive a dynamical system able to predict the capsule shape evolution over time for any value of the parameter $\mathit{\mu }$. This will be developed in the next sections.

2. Higher-order singular value decomposition and truncation

2.1 Compact HOSVD

The higher-order singular value decomposition (HOSVD) of a tensor is a specific orthogonal Tucker decomposition. The classical computation of a HOSVD was introduced by L. R. Tucker [5] and further developed by L. De Lathauwer et al. [6]. Robust computations or improvements have been since proposed [6, 7, 8]. For a tensor $\mathcal{T}$ of order $d$, the idea is to compute the singular value decomposition of each factor- $k$ flattening $T_{\left(k\right)}$ of a tensor $\mathcal{T}$, i.e. a “matricisation” of the tensor where the rows of the matrix are related to the $k$-th dimension.

In our case, we consider a third-order tensor and successive spatial, parameter and temporal flattening to find the singular vectors. The spatial tensor flattening $T_{\left(X\right)}\in {\mathbb{R}}^{I\times \left(J.N\right)}$ of the tensor $\left({\mathcal{T}}_x-A^0\otimes {\mathit{e}}^N\right)$ is

$${\left(T_{\left(X\right)}\right)}_{\mathit{ip}}=\mathit{x}\left({\mathit{X}}_i,{\mathit{\mu }}_j,t^n\right)-{\mathit{X}}_i,p=j+\left(n-1\right)J$$E7

(meaning that the $\mathit{\mu }$ and $t$ dimensions are stacked in columns in the matrix). The SVD of $T_{\left(X\right)}$ provides $r_x$ nonzero singular values with $r_x$ corresponding singular orthonormal vectors ${\mathbf{\Phi }}_k$, $k=1,\dots ,r_x$. Similarly the parameter flattening leads to a rank $r_{\mu }$ with $r_{\mu }$ modes ${\mathbf{\Psi }}_{\mathrm{\ell }}$, $\mathrm{\ell }=1,\dots ,r_{\mu }$ and the time flattening a rank $r_t$ with $r_t$ modes ${\mathit{w}}_m$, $m=1,\dots ,r_t$. The tuple $\left(r_X,r_{\mu },r_t\right)$ is the multilinear rank of ${\mathcal{T}}_x$. Then tensor ${\mathcal{T}}_x$ can be written as

$${\mathcal{T}}_x=A^0\otimes {\mathit{e}}^N+\sum _{k=1}^{r_X}\sum _{\mathrm{\ell }=1}^{r_{\mu }}\sum _{m=1}^{r_t}{a}_{k\mathrm{\ell }m}{\mathbf{\Phi }}_k\otimes {\mathbf{\Psi }}_{\mathrm{\ell }}\otimes {\mathit{w}}_m.$$E8

2.2 Approximation

Among the applications, HOSVD can be used to define a low-order approximation of tensors. The so-called truncated HOSVD [9, 10, 11] consists in truncating the expansion (8) at a given multilinear rank $\left(K,L,M\right)$, $K\le r_X$, $L\le r_{\mu }$, $M\le r_t$, leading to (6). Let $\text{mlrank}\left(\mathcal{T}\right)$ denote the multilinear rank of the tensor $\mathcal{T}$. It has been shown that the approximation (6) returns a quasi-optimal solution of the nonlinear non-convex least-square problem

$$\underset{\tilde{\mathcal{T}}}{min}\frac{1}{2}\parallel {\mathcal{T}}_x-\tilde{\mathcal{T}}{\parallel }_F^2$$E9

(here $\parallel .{\parallel }_F$ is the Frobenius norm for tensors) subject to $\text{mlrank}\left(\tilde{\mathcal{T}}\right)=\left(K,L,M\right)$. The truncation ranks can be determined a priori according to the classical relative information content (RIC) of the SVD theory.

Eq. (6) can already be used as a compressed representation of the data, allowing for a lower storage complexity and a simpler manipulation, with low information loss if the RIC is high.

3. Reduced-order modeling of capsule dynamics

Eq. (6) provides a summarization of the family of spatio-temporal capsule shape solutions in the time interval $\left[t^0,t^N\right]$. Unfortunately, this algebraic model has no predictability capability for time $t>t^N$. To derive a predictable time-dependent model from the data tensor ${\mathcal{T}}_x$, one has to derive a differential system that approximates the FSI system of Eqs. (2) and (3). The HOSVD reduction can thus be valued in the context of model-order reduction.

Consider of parameter vector of interest $\mathit{\mu }\in D$. The capsule position approximate solution (5) can be rewritten as

$$\tilde{\mathit{x}}\left(\mathit{X},\mathit{\mu },t\right)=\mathit{X}+\sum _{k=1}^K{a}_{\mu ,k}\left(t\right){\mathit{\phi }}_k\left(\mathit{X}\right)$$E10

with

$$a_{\mu ,k}\left(t\right)=\sum _{\mathrm{\ell }=1}^L\sum _{m=1}^M{a}_{k\mathrm{\ell }m}{\psi }_{\mathrm{\ell }}\left(\mathit{\mu }\right){\omega }_m\left(t\right).$$E11

Let ${\mathit{a}}_{\mu }\left(t\right)$ be the vector-valued function

$${\mathit{a}}_{\mu }\left(t\right)={\left(a_{\mu ,1}\left(t\right),\dots ,a_{\mu ,K}\left(t\right)\right)}^T\in {\mathbb{R}}^K.$$E12

We would like to derive a differential system of unknowns ${\mathit{a}}_{\mu }\left(t\right)$. Since $\mathit{x}\left(\mathit{X},\mathit{\mu },0\right)=\mathit{X}$, we have the natural initial condition ${\mathit{a}}_{\mu }\left(0\right)=\mathbf{0}$. Practically, from expansion (8), we can compute coefficients $a_{\mu ,k}\left(t\right)$ at discrete times $t^n$, and have thus access to the list of coefficient vectors

$${\mathit{a}}_{\mu }^n={\mathit{a}}_{\mu }\left(t^n\right),n=0,\dots ,N.$$E13

A dynamical reduced-order model consists in determining (or approximating) a Lipschitz continuous mapping ${\mathit{F}}_{\mu }:{\mathbb{R}}^K\to {\mathbb{R}}^K$ such that

$${\mathit{a}}_{\mu }^0=\mathbf{0},{\mathit{a}}_{\mu }^{n+1}\approx {\mathbf{F}}_{\mu }\left({\mathit{a}}_{\mu }^n\right)\forall n\in 0,\dots ,N-1$$E14

from the data ${\left\{{\mathit{a}}_{\mu }^n\right\}}_{n=0}^N$. We get a low-order discrete dynamical system. Note that, here, we do not search for parameters of a model, but for the equations of the model themselves. Since the problem of finding such a mapping is infinite dimensional, one has to restrict the search to a mapping in a (suitable) finite dimensional functional space.

4. Koopman theory and dynamic mode decomposition

4.1 Koopman operator for discrete dynamical systems

Koopman theory is a powerful mathematical framework that re-expresses a general nonlinear discrete dynamical system as the knowledge of a linear (infinite dimensional) operator, the so-called Koopman operator or compositional operator. Today it is commonly used in machine learning and data-driven model-order reduction methodologies [4, 12]. Let us assume a discrete dynamical system in the form

$${\mathit{a}}^{n+1}=\mathit{F}\left({\mathit{a}}^n\right),n\in \mathbb{N}$$E15

for a Lipschitz continuous mapping $\mathit{F}$ from ${\mathbb{R}}^d$ to ${\mathbb{R}}^d$. Let $g$ be a function of a Banach space $X$, $g:{\mathbb{R}}^d\to \mathbb{R}$. So we have $g\left({\mathit{a}}^{n+1}\right)=\left(g\circ \mathit{F}\right)\left({\mathit{a}}^n\right)$. The Koopman operator related to $\mathit{F}$ is defined as

$$\mathcal{K}g=g\circ \mathit{F}\forall g\in X.$$E16

Then we have

$$g\left({\mathit{a}}^{n+1}\right)=\left(\mathcal{K}g\right)\left({\mathit{a}}^n\right).$$E17

The knowledge of $\mathcal{K}$ includes the knowledge of $\mathit{F}$. Indeed, by taking the particular observables $g_i\left(\mathit{a}\right)=\mathit{a}\cdot {\mathit{e}}_i$, $i=1,\dots ,d$ where ${\mathit{e}}_i$ is the $i$-th vector of the canonical basis of ${\mathbb{R}}^d$, we retrieve ${\left({\mathit{a}}^{n+1}\right)}_i=F_i\left({\mathit{a}}^n\right)$, i.e. the $i$-th equation of (15). Of course, the linear Koopman operator acts on an infinite-dimensional functional space, so it is impossible to determine it exactly. However, one can search for an approximate Koopman operator $\tilde{\mathcal{K}}$ that acts on an approximate finite-dimensional space $\tilde{X}\subset X$.

The concept of (nonlinear) observables is to have a sufficiently large set of independent nonlinear functions of the state vector and measurements of them in order to identify the mapping $\mathit{F}$. A natural question of interest is what are the best observables to choose. There is no absolute answer to this question, and the choice may depend on the underlying Physics. Without any a priori knowledge on the system of equations, one can use basis functions of a universal approximators like polynomials, Fourier or kernel-based functions for example.

4.2 Dynamic mode decomposition

The simplest choice of observables is the linear functions $g_i\left(\mathit{a}\right)=\mathit{a}\cdot {\mathit{e}}_i$, $i=1,\dots ,d$. It leads to the search of a finite-dimensional approximation $A$ of $\mathcal{T}$ from the full state vector data. The matrix $A$ can be searched as the solution of the least square minimization problem

$$\underset{A\in {\mathcal{M}}_d\left(\mathbb{R}\right)}{min}\frac{1}{2}\parallel Y-\mathit{AX}{\parallel }_F^2$$E18

where $X=\left[{\mathit{a}}^0,{\mathit{a}}^1,\dots ,{\mathit{a}}^{N-1}\right]$ and $Y=\left[{\mathit{a}}^1,{\mathit{a}}^2,\dots ,{\mathit{a}}^N\right]$. Assuming $N\ge d$, the solution is given by $A={\mathit{YX}}^T{\left({\mathit{XX}}^T\right)}^{-1}$. The least square problem (18) can be eventually regularized for better conditioning by a Tykhonov regularization term [13, 14]. This practical approach of Koopman operator approximation is referred to as the dynamic mode decomposition (DMD) [4, 12]. This provides a linear dynamical model

$${\mathit{a}}^{n+1}=A{\mathit{a}}^n$$E19

starting from a given initial condition ${\mathit{a}}^0$. The solution of (19) which is ${\mathit{a}}^n=A^n{\mathit{a}}^0$ is bounded for any initial condition ${\mathit{a}}^0$ as soon as $\rho \left(A\right)\le 1$.

5. Kernel-based identification of dynamical systems

In the case of a strongly linear dynamical system, the linear model (19) can be not accurate enough. We have to include suitable nonlinear observables in the data and the model. In this section, observables are selected from kernel-based approximations [15]. Then we use the variant kernel-DMD ($k$-DMD, [4]) approach to identify $\mathit{F}$.

A real-valued function $k$ on ${\mathbb{R}}^d\times {\mathbb{R}}^d$ is called a positive definite kernel function if it is symmetric and if the following property holds:

$$\forall m\in {\mathbb{N}}^{\star },\forall {\left\{{\mathit{z}}_i\right\}}_{i=1}^m\in {\left({\mathbb{R}}^d\right)}^m,\forall {\left\{{\alpha }_i\right\}}_{i=1}^m\in {\mathbb{R}}^m,\sum _{i=1}^m\sum _{j=1}^m{\alpha }_i{\alpha }_{j}k\left({\mathit{z}}_i,{\mathit{z}}_j\right)\ge 0.$$E20

In other words, the square matrix $\mathit{K}={\left(k\left({\mathit{z}}_i,{\mathit{z}}_j\right)\right)}_{i=1,\dots ,m,j=1,\dots ,m}$ is positive semi-definite. A standard kernel function is for example the Gaussian one

$$k\left(\mathit{z},{\mathit{z}}^{\prime }\right)=exp\left(-\frac{1}{2}\frac{\parallel \mathit{z}-{\mathit{z}}^{\prime }{\parallel }^2}{{\sigma }^2}\right)$$E21

for a given parameter $\sigma >0$.

5.1 Kernel-based interpolation

Kernel functions can be used for interpolation in spaces of arbitrary dimension. Let $g:{\mathbb{R}}^d\to \mathbb{R}$ be a continuous function and assume that we know the values of $g\left({\mathit{z}}_i\right)$ at particular points ${\mathit{z}}_i$, $i=1,\dots ,m$. Then one can define an interpolator $\mathcal{I}g$ of $g$ defined as

$$\mathcal{I}g\left(\mathit{z}\right)=\sum _{j=1}^m{\alpha }_{j}k\left(\mathit{z},{\mathit{z}}_j\right)$$E22

where the coefficient vector $\mathit{\alpha }={\left({\alpha }_1,\dots ,{\alpha }_m\right)}^T$ is determined such that the interpolation property

$$\mathcal{I}g\left({\mathit{z}}_i\right)=g\left({\mathit{z}}_i\right)\forall i\in \left\{1,\dots ,m\right\}.$$E23

holds. The interpolation conditions clearly lead to the solution of the symmetric linear square system of size $m$

$$\mathit{K\alpha }=\mathit{b}$$E24

where $\mathit{b}={\left(g\left({\mathit{z}}_1\right),\dots ,g\left({\mathit{z}}_m\right)\right)}^T)$. Assuming that $\mathit{K}$ is positive definite, then Eq. (24) has a unique solution. Let

$$k_i\left(\mathit{z}\right)=k\left(\mathit{z},{\mathit{z}}_i\right),i=1,\dots ,m$$E25

and $V=\text{span}\left(k_1,\dots ,k_m\right)$. Considering any function $\tilde{g}\in V$, it is easy to check that $\tilde{g}=\mathcal{I}\tilde{g}$. One can derive the interpolation error

$$\begin{array}{l}\parallel g-\mathcal{I}g{\parallel }_{L^{\infty }}=\parallel g-\tilde{g}+\tilde{g}-\mathcal{I}g{\parallel }_{L^{\infty }}\\ =\parallel g-\tilde{g}+\mathcal{I}\tilde{g}-\mathcal{I}g{\parallel }_{L^{\infty }}\\ \le \parallel g-\tilde{g}{\parallel }_{L^{\infty }}+\mid \Vert \mathcal{I}\mid \Vert \parallel g-\tilde{g}{\parallel }_{L^{\infty }}\end{array}$$E26

so that

$$\parallel g-\mathcal{I}g{\parallel }_{L^{\infty }}\le \left(1+|\Vert \mathcal{I}|\Vert \right)\underset{\tilde{g}\in V}{inf}\parallel g-\tilde{g}{\parallel }_{L^{\infty }}.$$E27

The interpolation error is controlled by the best approximation error multiplied by a stability constant depending on the norm of the interpolation operator.

5.2 Use of kernel features and $k$-DMD

Let us go back to the parameterized dynamical system of interest (14) and consider a point cloud ${\left\{{\mathit{a}}_{\left(j\right)}\right\}}_{j=1}^m$ of sample states in the admissible reduced state space $\mathcal{X}\subset {\mathbb{R}}^K$. The functions

$$k_j\left(\mathit{a}\right)=k\left(\mathit{a},{\mathit{a}}_{\left(j\right)}\right)$$E28

can be seen as features and thus be used as suitable nonlinear observables to approximate the Koopman operator. From any known full state vector ${\mathit{a}}_{\mu }^n$ at time $t^n$, we build the vector of observables

$$\mathit{\kappa }\left({\mathit{a}}_{\mu }^n\right)={\left(k_1\left({\mathit{a}}_{\mu }^n\right),k_2\left({\mathit{a}}_{\mu }^n\right),\dots ,k_m\left({\mathit{a}}_{\mu }^n\right)\right)}^T.$$E29

By definition of the Koopman operator, we have $k_i\left({\mathit{a}}_{\mu }^{n+1}\right)=\left(k_i\circ \mathit{F}\right)\left({\mathit{a}}_{\mu }^n\right)=\left(\mathcal{K}{k}_i\right)\left({\mathit{a}}_{\mu }^n\right)$. Then we look for a finite-dimensional approximation $A_{\mu }$ of the Koopman operator $\mathcal{K}$ in the sense

$$\mathit{\kappa }\left({\mathit{a}}_{\mu }^{n+1}\right)\approx A_{\mu }\mathit{\kappa }\left({\mathit{a}}_{\mu }^n\right)\forall n\in \left\{0,\dots ,N-1\right\}$$E30

The matrix $A_{\mu }$ is searched as the minimum of the least square problem

$$\underset{A\in {\mathcal{M}\mathcal{M}}_m\left(\mathbb{R}\right)}{min}\frac{1}{2}\parallel Y_{\mu }-A{X}_{\mu }{\parallel }_F^2$$E31

where the entry and output data matrices are now $X_{\mu }=\left[\mathit{\kappa }\left({\mathit{a}}_{\mu }^0\right),\mathit{\kappa }\left({\mathit{a}}_{\mu }^1\right),\dots ,\mathit{\kappa }\left({\mathit{a}}_{\mu }^{N-1}\right)\right]$ and $Y_{\mu }=\left[\mathit{\kappa }\left({\mathit{a}}_{\mu }^1\right),\mathit{\kappa }\left({\mathit{a}}_{\mu }^2\right),\dots ,\mathit{\kappa }\left({\mathit{a}}_{\mu }^N\right)\right]$. We get the dynamical system $\mathit{\kappa }\left({\mathit{a}}_{\mu }^{n+1}\right)=A_{\mu }\mathit{\kappa }\left({\mathit{a}}_{\mu }^n\right)$ with a specified initial condition $\mathit{\kappa }\left({\mathit{a}}_{\mu }^0\right)$.

Let us emphasize that the computational variables are now the $\mathit{\kappa }\left({\mathit{a}}_{\mu }^n\right)$. But we still need the full state variables ${\mathit{a}}^n$ to determine the displacements or the capsule shapes. The full state can be retrieved for example by interpolation: taking $g_i\left(\mathit{a}\right)=\mathit{a}\cdot {\mathit{e}}_i$, we get

$${\left({\mathit{a}}_{\mu }\right)}_i^{n+1}\approx \mathcal{I}{g}_i\left({\mathit{a}}_{\mu }^{n+1}\right)=\sum _{j=1}^m{\left({\mathit{\alpha }}_j\right)}_i{k}_i\left({\mathit{a}}_{\mu }^{n+1}\right)$$E32

with interpolation coefficients vectors ${\left\{{\mathit{\alpha }}_j\right\}}_{j=1}^m$ such that

$${\mathit{a}}_{\left(i\right)}=\sum _{j=1}^m{\mathit{\alpha }}_j{k}_j\left({\mathit{a}}_{\left(i\right)}\right)\forall i\in \left\{1,\dots ,m\right\}$$E33

(linear system of dimension $m\times K$). The coefficient vectors ${\mathit{\alpha }}_j$ can be computed once for all.

5.3 Including the full state vector and the constants into the features

The $k$-DMD approach can be improved by including the full-state vector ${\mathit{a}}_{\mu }$ itself into the input feature vector. Moreover, if the kernel functions are not able to perfectly reproduce the constant ones, for accuracy reasons it is justified to add the constants into the features. Considering the augmented vector $\mathit{\eta }\left({\mathit{a}}_{\mu }\right)$ of observables

$$\mathit{\eta }\left({\mathit{a}}_{\mu }\right)=\left[\begin{array}{c}{\mathit{a}}_{\mu }\\ \mathit{\kappa }\left({\mathit{a}}_{\mu }\right)\\ 1\end{array}\right]\in {\mathbb{R}}^{K+m+1},$$E34

we look for a dynamical system in the form

$${\mathit{a}}_{\mu }^{n+1}=A_{\mu }\mathit{\eta }\left({\mathit{a}}_{\mu }^n\right)\approx {\mathit{F}}_{\mu }\left({\mathit{a}}_{\mu }^n\right)$$E35

with a constant rectangular matrix $A_{\mu }$ of size $K\times \left(K+m+1\right)$ to identify. The first advantage is that the output vector is the full state itself. The second one is the number of elements of $A_{\mu }$ which is less than $m^2$ as soon as $K\left(K+m+1\right)\le m^2$. In this case, there are less coefficients to identify. The input and output matrices now are

$$X_{\mu }=\left[\mathit{\eta }\left({\mathit{a}}_{\mu }^0\right),\mathit{\eta }\left({\mathit{a}}_{\mu }^1\right),\dots ,\mathit{\eta }\left({\mathit{a}}_{\mu }^{N-1}\right)\right],Y_{\mu }=\left[{\mathit{a}}_{\mu }^1,{\mathit{a}}_{\mu }^2,\dots ,{\mathit{a}}_{\mu }^N\right].$$E36

Assuming $N\ge \left(K+m\right)$, the rectangular matrix $A_{\mu }$ solution of the least square problem

$$\underset{A\in {\mathcal{M}}_{K,K+m+1}\left(\mathbb{R}\right)}{min}\frac{1}{2}\parallel Y_{\mu }-A{X}_{\mu }{\parallel }_F^2$$E37

is computed as $A_{\mu }=Y_{\mu }{X}_{\mu }^T{\left(X_{\mu }{X}_{\mu }^T\right)}^{-1}$.

5.4 Summary and whole algorithm

We give a summary of the model-order reduction algorithm:

1. Input data: third-order tensor ${\mathcal{T}}_x$ (4) made of capsule shape solutions of size $3I\times J\times \left(N+1\right)$.

2. HOSVD + truncated approximation: compute (8) and get the truncated approximation with the truncated multilinear ranks $\left(K,L,M\right)$:

   $${\tilde{\mathcal{T}}}_x=A^0\otimes {\mathit{e}}^N+\sum _{k=1}^K\sum _{\mathrm{\ell }=1}^L\sum _{m=1}^M{a}_{k\mathrm{\ell }m}{\mathbf{\Phi }}_k\otimes {\mathbf{\Psi }}_{\mathrm{\ell }}\otimes {\mathit{w}}_m$$E38

3. Online stage: choose a parameter vector $\mathit{\mu }$. From (13), compute the data coefficients ${\mathit{a}}_{\mu }\left(t^n\right)={\mathit{a}}_{\mu }^n$ (see Eq. (13)) from (11). Choose a kernel function $k\left(.,.\right)$, choose $m$ and the centers ${\mathit{a}}_{\left(j\right)}$, $j=1,\dots ,m$. Assemble the observables $\mathit{\eta }\left({\mathit{a}}_{\mu }^n\right)$ and assemble the matrices $X_{\mu }$ and $Y_{\mu }$ (Eq. (36)). Compute the $k$-DMD matrix $A_{\mu }=Y_{\mu }{X}_{\mu }^T{\left(X_{\mu }{X}_{\mu }^T\right)}^{-1}$. Get the reduced-order dynamical system (35) with ${\mathit{a}}_{\mu }^0=\mathbf{0}$ as initial value.

6. Numerical results for capsule dynamics

The algorithm is applied to a problem of a deformable capsule flowing into a square-based microchannel (typical for microfluidic channels created by soft lithography) with a mean inflow velocity $V$ [1, 2]. Related works and variant ROM approaches on this topic can be found in [14, 16]. The capsule dynamics is simulated by the full-order fluid–structure interaction solver Cite [1]. The fluid solver is based on the solution of the Stokes equations and a nonlinear Neo-Hookean law is used for the membrane. The initial capsule is spherical, corresponding to the shape at rest. The sphere is discretized with a finite element mesh made of $I=2562$ vertices.

6.1 Study for a particular parameter couple $\left(\mathit{Ca},a/\mathrm{\ell }\right)$

The objective of this subsection being to illustrate the methods on an example and show how to apply them, we only consider the snapshot FOM solutions for $\left(\mathit{Ca},a/\mathrm{\ell }\right)=\left(\mathrm{0.1,0.9}\right)$ for the sake of simplicity and brevity. Figure 1 shows both membrane shapes and positions in the channel at different instants.

Figure 1.

Example of a microcapsule dynamics within a square-base channel for ($\mathit{Ca},a/\mathrm{\ell }$) = (0.1, 0.90). Shapes and locations of the initially spherical capsule are shown at $t=\mathrm{0}$ (in transparency), 0.4, 2.8, 5.2, 7.6.

At each time $t^n=n\mathrm{\Delta }t$, where the time step is equal to $\mathrm{\Delta }t=0.04$, a snapshot is saved and stored in the database. Note that all time quantities are non-dimensionalized by the factor $\mathrm{\ell }/V$. The resulting generated data matrix is used as the entry matrix for the ROM learning process. Then a truncated SVD is applied to get the spatial POD modes. In Figure 2, the four first eigenmodes ${\mathrm{\Phi }}^k$ are plotted (more precisely this is a superposition of each mode onto the original spherical shape for a better visualization and understanding of their influence). Based on, the graphics $k\mapsto 1-\mathit{RIC}\left(k\right)$ plotted in log scale in Figure 3, we decide to use a truncation rank $K$ equal to $K=10$, returning a relative information content of about $1-3.5\times {10}^{-5}$.

Figure 2.

SVD: Four first spatial principal components computed by the HOSVD. Each mode has been added on the initial spherical shape and amplified by a factor 2 for better visualization. Higher-order modes show oscillations at the rear of the capsule.

Figure 3.

SVD: Plot of $k\mapsto 1-\mathrm{RIC}\left(k\right)$, where $\mathrm{RIC}\left(k\right)$ represents the relative information content at truncation rank $k$.

Then a reduced-order dynamical system for the capsule time evolution is searched. In this example, we compare two models: the first one is the affine approximation (denoted by ROM-A)

$${\mathit{a}}_{\mu }^{n+1}=A_{\mu }{\mathit{a}}_{\mu }^n+{\mathit{b}}_{\mu }$$E39

with the matrix $A$ and the vector $\mathit{b}$ to identify. It is equivalent to consider the vector of features $\mathit{\eta }\left(\mathit{a}\right)={\left(\mathit{a},1\right)}^T$. The second nonlinear model is built from the observable vector $\mathit{\eta }\left(\mathit{a}\right)={\left(\mathit{a},\mathit{\kappa }\left(\mathit{a}\right),1\right)}^T$ with the recurrent time scheme

$${\mathit{a}}_{\mu }^{n+1}=A_{\mu }\mathit{\eta }\left({\mathit{a}}_{\mu }^n\right)$$E40

(ROM-B). For ROM-B, the Gaussian kernel function (21) is used. The standard deviation parameter $\sigma$ in (21) is chosen as $\sigma ={max}_n\parallel {\mathit{a}}_{\mu }^n\parallel =\parallel {\mathit{a}}_{\mu }^{N+1}{\parallel }_2$. In both cases ROM-A and ROM-B, the determination of the matrix $A_{\mu }$ by minimization of the least square problem (37) leads to a very small residual. In Figure 4, the logarithm of the relative time residual

Figure 4.

Matrix identification. Log of the normalized residual $n\mapsto {log}_{10}\left(\parallel {\mathit{a}}_{\mu }^{n+1}-A_{\mu }\mathit{\eta }\left({\mathit{a}}_{\mu }^n\right){\parallel }^2/\parallel {\mathit{a}}_{\mu }^{n+1}{\parallel }^2\right)$ for both ROM-A and ROM-B.

$$n\mapsto {log}_{10}\left(\parallel {\mathit{a}}_{\mu }^{n+1}-A_{\mu }\mathit{\eta }\left({\mathit{a}}_{\mu }^n\right){\parallel }^2/\parallel {\mathit{a}}_{\mu }^{n+1}{\parallel }^2\right)$$E41

is plotted for each ROM model. One can observe values between ${10}^{-14}$ and ${10}^{-8}$. The residual for ROM-B appears to be slightly smaller than that of ROM-A thanks to the added nonlinear terms. A surprising result is that the affine ROM-A model returns rather accurate results whereas the fluid–structure interaction problem is intrinsically nonlinear. In order to study the stability of the model ROM-A, in Figure 5 we plot the complex eigenvalues of the square matrix $A_{\mu }$. We observe that all the eigenvalues have a modulus less or equal to one. One of the eigenvalues is equal to $1$ exactly up to Double Precision roundoff errors, meaning that there is a physical invariant in the system. It is known that the capsule volume is kept constant during time, because of the incompressibility of the fluid. For ROM-A, since ${\mathit{a}}_{\mu }^{n+1}=A_{\mu }{\mathit{a}}_{\mu }^n+\mathit{b}$ and ${\mathit{a}}_{\mu }^0=\mathbf{0}$, we have

Figure 5.

Matrix identification. Eigenvalues of the computed matrix $A_{\mu }$ plotted in the complex plane for the ROM-A model. One of the eigenvalue is $1$ up to round-off error.

$${\mathit{a}}_{\mu }^n=A_{\mu }^n{\mathit{a}}_{\mu }^0+\sum _{k=0}^{n-1}{A}_{\mu }^k\mathit{b}=\sum _{k=0}^{n-1}{A}_{\mu }^k\mathit{b}.$$E42

The matrix $A_{\mu }$ is observed to be diagonalizable in $\mathrm{ℂ}$. There is an invertible matrix $P_{\mu }$ such that $A_{\mu }=P_{\mu }{\mathrm{\Lambda }}_{\mu }{P}_{\mu }^{-1}$ where ${\mathrm{\Lambda }}_{\mu }$ is the diagonal matrix of the eigenvalues. Since it is observed that $\rho \left({\mathrm{\Lambda }}_{\mu }\right)=1$, we have

$$\parallel {\mathit{a}}_{\mu }^n\parallel \le \text{Cond}\left(P_{\mu }\right)n\parallel \mathit{b}\parallel ,\forall n\in \mathbb{N},$$E43

showing that the coefficients in the PCA space grow at most linearly in time.

In Figure 6, we compare the computed capsule shapes and positions in the channel for the computed FOM capsules obtained at different times: $t=0$, $0.4$, $2.8$, $5.2$ and $7.6$. What can be observed is that the ROM-B model returns very satisfactory results where the shape solutions fully overlap the FOM ones’at the eye norm’. Finely, we plot in Figure 7 the time evolution of the errors in the capsule 3D shape of the ROM solutions as compared to the FOM solutions. The difference between the shapes are quantified by ${\epsilon }_{\text{Shape}}\left(t\right)$, the ratio between the modified Hausdorff distance (MHD) computed between the FOM shape ${\mathcal{S}}_{\mathrm{FOM}}\left(t\right)$ and the ROM shape ${\mathcal{S}}_{\mathrm{ROM}}\left(t\right)$ and the channel characteristic length $\mathrm{\ell }$:

Figure 6.

Sequence of cross-section capsule shapes and positions in the microchannel from the initial spherical shape shown in light green at the beginning of the channel: Comparison of the FOM solutions (gray dots) and of the solutions computed from the dynamical k-DMD reduced-order model (dark green solid line) at the same instants as in Figure 1: $t=\mathrm{0,\; 0.4,2.8,5.2,7.6}$.

Figure 7.

(a) Comparison of the time evolution of the shape error between the affine DMD model (ROM-A with $K=10$) and the kernel-based one (ROM-B with $K=10$, $M=5$). One can observe a maximum error less than 0.1% in both cases. (b) Sensitivity analysis of the parameter $\sigma$ ($\parallel {\mathit{a}}^{N+1}{\parallel }_2=667.14$.)

$${\epsilon }_{\text{Shape}}\left(t\right)=\frac{\mathrm{MHD}\left({\mathcal{S}}_{\mathrm{FOM}}\left(t\right),{\mathcal{S}}_{\mathrm{ROM}}\left(t\right)\right)}{\mathrm{\ell }}.$$E44

The modified Hausdorff distance measures the maximum value of the mean distance between the two shapes to compare [17]. The ROM-A and ROM-B models return very accurate solutions with maximum 0.1% error. It is also observed that the ROM-B models is slightly more accurate than the affine approximation.

6.2 HOSVD on the whole data tensor and error measurements on the whole database

Now the consider the whole database made of 55 samples in the parameter domain. In Figure 8 we plot the location of the 55 chosen samples in the plane $\left(\mathit{Ca},a/\mathrm{\ell }\right)$. The design zone of interest corresponds to capsule shapes that can reach a steady state after a certain time.

Figure 8.

Samples of the design of experiment in the parameter space $\left(\mathit{Ca},a/\mathrm{\ell }\right)$. The zone of interest corresponds to capsule shapes that reach a steady state after a certain time.

A SVD is first performed on the $\mu$-flattening of the data tensor ${\mathcal{T}}_x$. In Figure 9, we plot the four first parameter (normalized) eigenmodes in the parameter domain. These modes give an idea on the dependency of the capsule shapes with respect to the parameters. To complete the analysis, we show in Figure 10 the spectrum of the singular values for the $\mu$-flattening matrix. One can observe a rather fast decay of the singular values especially for the ten first modes.

Figure 9.

HOSVD: The first four parameter eigenmodes in the parameter domain, computed from the $\mu$-flattening of the data cube.

Figure 10.

Spectrum of the singular values of the $\mu$-flattening matrix.

Next, we perform the SVD of the time $t$-flattening matrix of ${\mathcal{T}}_x$. The SVD provides us temporal eigenmodes. In Figure 11, the four first temporal eigenmodes ${\omega }_m$, $m=1,\dots ,4$ are plotted. The first one appears to be the linear function while the others add details especially in the transient zone of the capsule dynamics. The spectrum of the singular values again shows a fast decay especially for the six first modes.

Figure 11.

HOSVD. Left: First four temporal eigenmodes computed from the SVD decomposition of the $t$-flattening data tensor. Right: Spectrum of the singular values.

To conclude this section, we have tested the accuracy of both ROM-A and ROM-B on the whole database. For each sample, we have derived a ROM model, i.e. a low-order dynamical system formulated in the PCA space. Then we have compared the ROM solution to the FOM solution by calculating ${\epsilon }_{\text{Shape}}$ between the two capsule shapes. In Figure 12, the heat maps of ${\epsilon }_{\text{Shape}}$ are plotted for ROM-A and ROM-B. One can observe a uniform accuracy over the whole parameter domain, with errors less than 0.1%, thus showing the efficiency of the approach. Reported computational speedups are between 5000 and 10,000 using ROM models. A computer with two Intel Xeon GOLD 6130 CPU (2.1 Ghz) has been used for the numerical tests.

Figure 12.

Heat maps of the modified Hausdorff distance between the FOM solutions and the ROM ones at dimensionless time $\mathit{Vt}/\mathrm{\ell }=10$. Left: ROM-A, right: ROM-B. Errors are less than 0.1%.

Figure 13.

Interpretation of the method as a recurrent neural network (RNN) in the PCA space.

7. Discussion

Having shown how to apply the affine DMD and $k$-DMD models and the very high precision in prediction that they offer, we now would like to provide some further comments and remarks on model order reduction and DMD-type approaches.

7.1 Kolmogorov $n$-width

The method is efficient if the spectrum of singular values decays rapidly, leading to a small truncation rank $K$. If the spectrum decays slowly, there are two possible reasons for that: either the entropy (variety) of information in the data is high, or the solutions do not live in a linear space but rather on a nonlinear manifold. To fix the problem, one can proceed by performing a preliminary clustering of the data, scanning the parameter dimension. One can either use standard clustering techniques such as $K$-means, or a multidimensional scaling (MDS) approach. Then for each cluster, one can consider again a HOSVD and reduced-order approach suitable for each element of the cluster.

7.2 Selection of the kernel functions and kernel interpolation points

As already mentioned, the choice of the kernel function depends on the applications, on the behavior of solutions and/or on the underlying Physics. Without any a priori information, one can use universal approximation kernels like the Gaussian one. The accuracy of the results will also strongly depends on the choice of the kernel interpolation points ${\mathit{a}}_{\left(j\right)}$. The sampling ${\left\{{\mathit{a}}_{\left(j\right)}\right\}}_{j=1}^m$ has to correctly fill in the admissible space, or at least the state-space trajectory of interest. There are different possible strategies. A first candidate is the use of a clustering approach applied to the state-space data. The points ${\mathit{a}}_{\left(j\right)}$ are then the centroids of each cluster. But one can consider more sophisticated approaches like a greedy iterative approach that controls the interpolation error on the data. At each iteration an interpolation point ${\mathit{a}}_{\left(j\right)}$ is added at the location of worst interpolation error, considering all the sample solutions.

7.3 Interpretation in terms of a recurrent neural network

Let us remark that the approach can be reinterpreted as a (supervized) two-layer recurrent artificial neural network (RNN) (Figure 13) [18]. The first layer consists in generating the features $k_i\left({\mathit{a}}_{\mu }\right)$. The second layer is a linear matrix–vector multiplication using the matrix $A_{\mu }$.

8. Conclusions

In this chapter, the higher-order singular value decomposition has been proved to be a flexible and valuable tool in the data-driven reduced-order modeling of solutions of space–time-parameter problems, which are today at the heart of many industrial applications. The methodology has been tested on a problem of fluid–structure interaction of a deformable microcapsule flowing into a microchannel. Stokes equations have been used in the fluid whereas a nonlinear hypereleastic law has been used for the membrane. Different shape solutions computed by the full-order model have been stored into a third-order tensor. First, HOSVD allows us to compute spatial, temporal and parameter principal components and at the same time to compress the data. We get a low-order representation of the solutions with a shared spatial reduced basis. Spatial principal components are observed to provide suitable details in the shape solutions. The modes are arranged in decreasing order of importance according to the relative information content criterion. Next, additional ingredients such as kernel approximation and kernel-based dynamic mode decomposition are used to determine a reduced-order dynamical system for any parameter vector in the admissible parameter domain. The resulting low-dynamical system can be seen as an encoded recurrent neural network set into the PCA space. The approach allows us to explore the different shape solutions and visualize their evolution in the channel in real time.

Acknowledgments

The authors would like to thank Prof. Pierre Villon (UTC) for fruitful discussions on this topic. The project received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (Grant agreement No. ERC-2017-COG - MultiphysMicroCaps).

Conflict of interest

The authors declare no conflict of interest.

Nomenclature

SVD	Singular value decomposition
HOSVD	Higher-order singular value decomposition
PCA	Principal component analysis
MOR	Model order reduction
FOM	Full order model
ROM	Reduced order model
FSI	Fluid–structure interaction
DoCE	Design of computer experiment
POD	Proper orthogonal decomposition
RIC	Relative information content
DMD	Dynamic mode decomposition
k-DMD	Kernel-based dynamic mode decomposition
MHD	Modified Hausdorff distance
RNN	Recurrent neural network