# Generalized Empirical Interpolation Method ( GEIM )

## 1. GENERAL FORMULATION OF THE PROBLEM

We consider the following formulation:

$\mathcal{P}: \Omega \times \mathcal{D} \rightarrow \mathbb{K}$

where $\mathcal{P}$ is the problem, $\Omega \subseteq \mathbb{R}^{d}$ the physical domain ($d=2,3$), $\mathcal{D} \subseteq \mathbb{R}^{N_{p}}$ the parametical domain
and $\mathbb{K}$ a field ( $\mathbb{R}$ or $\mathbb{C}$).

In the following, we denote $u \in \mathcal{X}$ a soluttion of the problem $\mathcal{P}$ where $\mathcal{X}$ is an appropriate Banach space and $\mathcal{M}^{bk}$ the set of solutions of the problem $\mathcal{P}$ given by the best available model of $\mathcal{P}$ (the best knowledge model).

 Usually we assume that at least $L^{2}(\Omega)\subseteq\mathcal{X}\subseteq H^{1}(\Omega)$.

## 2. GEIM FORMULATION

GEIM is an algorithm that combines both data assimilation and order reduction:

• Order Reduction

We consider $\mathcal{M}^{Par}=\{u(p_{1}),u(p_{2}),\cdots ,u(p_{M}) \} \subseteq \mathcal{M}^{bk}$ the set $M$ particular solutions of $\mathcal{P}$.

• Data Assimilation

We consider that we have a set of $M$ experimental observation $y^{obs}_{m}(p), 1\leq m\leq M$ assuming that each observation $y^{obs}_{m}(p)$ is such that:

$y^{obs}_{m}(p)=\sigma_{m}(u^{true}(p))$

where $u^{true}(p)$ represents the true physical state of the system and $\sigma_{m}$ are linear functionals representing the sensors.

 $u^{true}(p)$ is not accessible.

We set $\Sigma=\{\sigma_{1},\sigma_{2},\cdots ,\sigma_{M}\}$.

From $\mathcal{M}^{Par}$ and $\Sigma$, we construct a vector subspace generated by the interpolations functions $\mathcal{M}^{GEIM}=\{q_{1},q_{2},\cdots,q_{M}\}$ and then we define an interpolation operator

$\mathcal{I}_{M}(u)=\sum_{j=1}^{M}\alpha_{j}(u)q_{j}$

such as $\sigma_{i}(\mathcal{I}_{M}(u))=\sigma_{i}(u) \forall 1\leq i\leq M$.

## 3. METHODOLOGY

Ideally we want to choose the linear forms $\sigma_{i} \in \Sigma$ and the basic functions $q_{i} \in \mathcal{M}^{bk}$ in an optimal way. To do this, we consider a Glouton algorithm to minimize the interpolation error. The construction of the interpolation functions $q \in \mathcal{X}$ and the selection of the linear forms $\sigma \in \Sigma$ are done recursively.

Given a first form generative function, chosen as the greatest generative functions in norm $\| \cdot \|_{\mathcal{X}}$, we can choose the linear form associated as being that which gives more information on $u$.

$\sigma_{1}=\underset{\sigma \in \Sigma}{argsup}\ |\sigma(u(p_{1}))|$

The basic interpolation function is defined as

$q_{1}=\frac{u(p_{1})}{\sigma_{1}(u(p_{1}))}$

We then define

$\begin{equation} \begin{array}{rl} u(p_{2})&=\underset{u\in \mathcal{M}^{Par}}{argsup}\ \| u-\sigma_{1}(u)q_{1}\|_{\mathcal{X}}\\ \sigma_{2}&=\underset{\sigma \in \Sigma}{argsup}\ |\sigma[u(p_{2})-\sigma_{1}(u(p_{2}))q_{1}]|\\ q_{2}&=\frac{u(p_{2})-\sigma_{1}(u(p_{2}))q_{1}}{\sigma_{2}[u(p_{2})-\sigma_{1}(u(p_{2}))q_{1}]} \end{array} \end{equation}$

and so on by induction.

For $M> 2$, we solve the linear system for the state $u \in \mathcal{M}^{bk}$ to find the coefficients of interpolations associated $(\alpha_{j}^{M-1}(u))_{1\leq j \leq M-1}$.

$\sigma_{i}(u)=\sum_{j=1}^{M-1}\ \alpha_{j}^{M-1}(u)\sigma_{i}(q_{j})\ \forall 1\leq i\leq M-1$

and use the interpolation operator to define the $M^{th}$ generating function $u(p_{M})$, the linear form $\sigma_{M}$ and the interpolating basic function $q_{M}$ in the following way

$\begin{equation} \begin{array}{rl} \mathcal{I}_{M-1}(u)&=\sum_{j=1}^{M-1}\alpha_{j}^{M-1}(u)q_{j}\\ u(p_{M})&=\underset{u\in \mathcal{M}^{Par}}{argsup\ }\| u-\mathcal{I}_{M-1}(u)\|_{X}\\ \sigma_{M}&=\underset{\sigma \in \Sigma}{argsup\ }|\sigma(u(p_{M})-\mathcal{I}_{M-1}(u(p_{M})) ) |\\ q_{M}&=\frac{u(p_{M})-\mathcal{I}_{M-1}(u(p_{M}))) }{\sigma(u(p_{M})-\mathcal{I}_{M-1}(u(p_{M}) ))} \end{array} \end{equation}$
 For a more stable numerical implementation, the interpolation functions can be orthonormalized by the orthogonalization method of Gram-Schmid.

We define the interpolation matrix $B^{M}=(B_{i,j}^{M})_{ 1\leq i,j \leq M }$ where $B_{i,j}^{M}=\sigma_{i}(q_{j})\ 1\leq i,j \leq M$ $B^{M}$ is a lower triangular matrix with a unitary diagonal and no singular with the other entries $B_{i,j}^{M} \in [-1,1$].

By seting $U^{p}=(U_{i}^{p}=\sigma_{i}(u(p)))_{1\leq i\leq M}$ and $\alpha^{p}=(\alpha_{i}^{p})_{1\leq i\leq M}$ the interpolation coefficients, the resolution of the interpolation problem can be done by solving the linear system

$B^{M} \alpha^{p}= U^{p}$

As $B^{M}$ is an lower triangular matrix with $1$ on the diagonal, we can avoid the inversion of the matrix and solve the system of $M$ equations to find the interpolator $\mathcal{I}_{M}(u)$

$\begin{equation} \left\{ \begin{array}{rl} \alpha_{1}^{p}&=\sigma_{1}(u(p)) \\ \alpha_{2}^{p}&=\sigma_{2}(u(p)-\sigma_{2}(q_{1})\alpha_{1}^{p}\\ \cdots\\ \alpha_{M}^{p}&=\sigma_{M}(u(p))-\sum_{i=1}^{M-1}\sigma_{M}(q_{i})\alpha_{i}^{p} \end{array}\right. \end{equation}$

This gives us the recursive formula for the $M^{th}$ interpolation operator

$\mathcal{I}_{M}(u(p))=\mathcal{I}_{M-1}(u(p))+\frac{\sigma_{M}(u(p)-\mathcal{I}_{M-1}(u(p)))}{\sigma_{M}(u(p_{M})-\mathcal{I}_{M-1}(u(p_{M})))}\times (u(p_{M})-\mathcal{I}_{M-1}(u(p_{M})))$

and the $m^{th}$ interpolation function

$q_{m}=\frac{u(p_{m})-\mathcal{I}_{m-1}(u(p_{m}))}{\sigma_{m}(u(p_{m})-\mathcal{I}_{m-1}(u(p_{m})))}$
 This dependence of the $m^{th}$ interolating basic function of the $m-1^{th}$ interpolation operator is an iterative procedure that could cause a buildup of numerical error. That’s why in practice other more stable algorithms are used.
Algorithm
1. Input: we consider we have already found $u(p_{m})$ and $\sigma_{m}$

2. Set $w=u(p_{m})\in \mathcal{M}^{bk}$

3. for $j=1$ to $m−1$ do:

4. Set :
$\begin{equation} \begin{array}{rl} r_{j,m}=\alpha_{j}(w)&=\sum_{i=1}^{m-1} (B^{M}_{j,i})^{-1} \sigma_{i}(w)\\ &=\sigma_{j}(w)-\sum_{i=1}^{j-1} \sigma_{j}(q_{i})\alpha_{i} \end{array} \end{equation}$

5. Compute $w=w − r_{j,m} q_{j}$

6. end for

7. Set $r_{m,m}=\sigma_{m}(w)$

8. Compute $q_{m}=\frac{w}{r_{m,m}}$

## 4. NUMERICAL IMPLEMENTATION OF GEIM

OFFLINE STAGE
1. Compute particular solutions $\mathcal{M}^{par}$ and linear forms $\Sigma$.

2. Compute interpolation functions $\mathcal{M}^{GEIM}$.

3. Compute $B^{M}$ the matrix interpolation.

ONLINE STAGE
1. Solve interpolation problem: $B^{M} \alpha^{p}= U^{p}$

2. Compute output: $\mathcal{I}_{M}(u)=\sum_{j=1}^{M}\alpha_{j}^{p}(u)q_{j}$