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}\)