POD BASIS INTERPOLATION VIA INVERSE DISTANCE WEIGHTING ON GRASSMANN MANIFOLDS

. An adaptation of the Inverse Distance Weighting (IDW) method to the Grassmann manifold is carried out for interpolation of parametric POD bases. Our approach does not depend on the choice of a reference point on the Grassmann manifold to perform the interpolation, moreover our results are more accurate than those obtained in [7]. In return, our approach is not direct but iterative and its relevance depends on the choice of the weighting functions which are inversely proportional to the distance to the parameter. More judicious choices of such weighting functions can be carried out via kriging technics [23], this is the subject of a work in progress.


1.
Introduction. The present work deals with interpolation on Grassmann manifolds. More precisely, we are interested to the interpolation of Proper Orthogonal Decomposition (POD) bases corresponding to the solutions of parametric evolution partial differential equations. An interesting example is the Navier-Stokes equations where the varying parameter is the Reynolds number. In such problems, the solution depends on a parameter λ, describing a set in R P .
Suppose that the solution, which belongs to a function space H, is known for (λ k ) 1≤k≤N . These solutions are in general computed with high computational costs discretization type methods as Finite Elements, Finite Volumes, Discontinuous Galerkin, etc. Our goal is to use low-dimensional parametric methods to compute an approximation of the solution for a new parameter value λ ∈ R P .
In various problems in fluid mechanics or fluid structure interaction, the POD basis is computed for a given Reynolds number or for other flow parameters. This basis is then used to build a reduced model to predict the flow for other values of these parameters. The question concerning the domain of validity of this reduced model arises then naturally. In Akkari et al. [5], for quasilinear parabolic partial differential equations, this domain of validity was specified according to the solution regularity and the number of modes considered in the reduced model. Similar results was obtained for Navier-Stokes equations [3] and for Burgers equations [4]. Numerical works on parametric sensitivity for fluid structure interaction problems was performed in [18]. In [19], the authors carried out a non-intrusive strategy for building computational vademecums dedicated to real-time simulations of nonlinear thermo-mechanical problems.
These methods are also used in control theory [15,27,12], in medical imaging [10,24], where one often seeks to use the POD basis to represent data although time-interval variation, spatial domain variation or parameter variation may occur.
One classical way to perform such approximations is the use of Reduced Basis (RB) methods. RB-methods consist on the approximation of the manifold solution M, i.e., the set of parametric solutions, by a low-dimensional subspace H m ⊂ H. In this topic, one popular way is the construction of such a subspace via the snapshots (u(λ k )) 1≤k≤N as, for example, H m = span{u(λ k ) : 1 ≤ k ≤ N } and use Galerkin type techniques to compute an approximation u m (λ) ∈ H m of u(λ) ∈ M. We refer the interested reader to [21,11] and the references therein.
Another way to approximate the parametric solution u for a new parameter value λ ∈ R P from given data (u(λ k )) 1≤k≤N via POD techniques is the interpolation on Grassmann manifolds. We describe briefly this method and refer the reader to [7] for more details. This procedure can be summarized as follows: for each k ∈ {1, 2, · · · , N }, consider the rank-m POD basis (ϕ i (λ k )) 1≤i≤m and the corresponding [7], the space H has a large but finite dimension). Now, denote by G m (H) the Grassmann manifold consisting of all m-dimensional subspaces of H. Then the m subspaces H (k) m can be seen as points on this manifold G m (H). For a new parameter value λ, the question is how to compute a good approximation of the rank-m POD basis to the solution u(λ), or in a more relevant way, how to find its corresponding point in the Grassmann manifold G m (H)? Indeed, it is known from the work of Amsallem and Farhat [7] that the appropriate objects to interpolate are not the POD bases but the underlying spanned subspaces. To answer this question, the idea developed in [7] is to choose a reference parameter value λ r ∈ {λ k : 1 ≤ k ≤ N }, to consider the POD subspace H perform an interpolation (see Section 3 for more details). Finally, one returns to the manifold G m (H) via the (reciprocal) exponential mapping from the tangent space T H (r) m G m (H) to get the desired approximation of u(λ) on G m (H). The disadvantage of the method developed in [7] is its dependence of the reference point and the lack of a strategy for choosing this reference point.
In what follows, we will study an interpolation method on Grassmann manifolds based on the Inverse Distance Weighting (IDW).
The present paper is organized as follows: In Section 2, we recall some classical facts on Proper Orthogonal Decomposition (POD) for evolution problems and their based Reduced Order Model. In Section 3, a brief presentation of the Grassmann manifold and the underlying tools used in the paper is carried out. Section 4 is dedicated to the introduction of the Grassmann-IDW algorithm. The convergence of the G-IDW is performed in Section 5, following ideas developed by Karcher [13] and Huiling Le [16]. Finally, numerical simulations on isothermal fluid flows around a circular cylinder are given to compare our approach to that developed by Amsallem et al. [7].
2. Proper Orthogonal Decomposition (POD). The usual mathematical framework for POD is the space L 2 ([0, T ]; H), where [0, T ] and H denote respectively a time-interval and a separable Hilbert space of spatial functions. In this context, the aim of POD is to decompose a given space-time dependent flow field u ∈ L 2 ([0, T ]; H) into an optimal orthonormal system of spatial modes (ϕ k ) k≥1 ⊂ H and their corresponding temporal amplitudes (a k (t)) k≥1 as the following The optimality of bases (ϕ k (x)) k≥1 is to be understood in the sense that any truncated series expansion where · ; · H stands for the inner product in the Hilbert space H, the endowed norm will be denoted by · . It is well known that the operator R is linear, bounded, nonnegative, self-adjoint and compact. Thus, there exists a complete orthonormal basis of H consisting of eigenfunctions of R with corresponding nonnegative eigenvalues. Each nonzero eigenvalue of R has a finite multiplicity and 0 is the only possible accumulation point of the spectrum of R. We denote by (λ k ) k≥1 and (ϕ k ) k≥1 the sequences of positive eigenvalues of R, repeated as many times as their multiplicity and sorted in decreasing order and the associated eigenfunctions, respectively. The importance of the orthonormal basis (ϕ k ) k≥1 lies in the fact that it is optimal in the following sense : for any integer m ≥ 1 and for any rank-m orthonormal family {ψ k , 1 ≤ k ≤ m} in H, it holds : Equivalently, if H m = span {ϕ 1 , ϕ 2 , · · · , ϕ m } and Π V denotes the orthogonal projection operator on a subspace V in H, then for every m−dimensional subspace V m ⊂ H. Furthermore, we have the error of this optimal representation : For sharp POD convergence results, we refer the reader to [8,9], where the authors proved that the POD expansion converges with exponential accuracy for parameterized transient heat equations or Poisson equations with bivariate functions.
On the other hand, the original data u can be completely represented by the eigenfunctions corresponding to positive eigenvalues of the operator R, i.e., There is an alternative characterization of POD basis via correlation operator. Indeed, consider the linear operator then the adjoint operator of A, denoted by A * is defined by It is straightforward to verify that We define the correlation operator K = A * A. A direct computation shows that Moreover, the operator K has the same properties as R and share together the same positive eigenvalues and the same multiplicities. Furthermore, eigenfunctions (ϕ k ) k≥1 of R and eigenfunctions (f k ) k≥1 of K can be found one from the other by : In the finite dimension case, the data u is a real array in R n×Nt In this case, the operator A can be represented by the real matrix On the other hand, the adjoint operator A * can be represented by the matrix A T ∈ M N t ,n (R). It follows that rank-m POD bases are full-rank matrices of n rows and m columns.
2.1. POD based reduced order models. The POD bases provide a suitable framework for Reduced Order Models via the Galerkin's method. Suppose that the given data u(x, t) is the solution of an evolution partial differential equation, as the Navier-Stokes equations for incompressible flows which will be considered in the numerical simulation section. Such equations can be written symbolically as : where F is a differential operator that contains only spatial derivatives and where Ω is the spatial domain under study. To complete the evolution differential problem (6), we need an initial condition and adequate boundary conditions. Notice that for stationary boundary conditions, the POD basis is usually com- homogeneous boundary conditions [25]. These homogeneous boundary conditions are then taken into account in the Hilbert space sequently, every eigenfunction of the POD basis (ϕ k ) 1≤k≤m satisfies the boundary conditions. The reduced order model is then derived by Galerkin projection of the partial differential equations (6) then the vector (a k ) 1≤k≤m satisfies the ordinary differential system Notice that if the POD basis is performed with the L 2 (Ω) inner product, then the m×m matrix M kj = δ kj , the matrix identity. This will be the case in our numerical applications.
3. The Grassmann manifold. Let m, n be two positive integers with m ≤ n.
The set of all R-subspaces of R n of dimension m, is a differentiable manifold of dimension m × (n − m), called the Grassmann manifold and usually denoted by be the set of all matrices of size n × m, whose column vectors are linearly independent; one defines on this set the following equivalence relation where GL m (R) denotes the linear group of degree m, i.e., the set of m×m invertible matrices. The Grassmann manifold can be realized as the set of such equivalence classes. More precisely, The Grassmann manifold is a compact topological space and thanks to the canonical projection, denoted by π : R n×m * → G m (R n ), it can be endowed by a differentiable manifold structure so that π is a submersion ( [20]).
3.1. Metric on the Grassmann manifold. On each point X ∈ R n×m * , one constructs a metric g X which is right-invariant by GL m (R), i.e., and A ∈ GL m (R). This metric induces, by the submersion π, a Riemannian metric on G m (R n ). Moreover, one constructs a Riemannian connection on G m (R n ), which allows the definition of geodesics and the geodesic exponential mapping on G m (R n ) [1].

3.2.
Exponential and logarithm mappings, injectivity radius. Let X ∈ G m (R n ) and v ∈ T X G m (R n ). The geodesic of initial point γ(0) = X and initial velocity γ (0) = v, can written in the form The geodesic exponential mapping is defined by exp X (v) = γ(1) and has a unique matricial representation given by: On the other hand, since d 0 X exp X = I, then the inverse function theorem implies that exp X is a local diffeomorphism at 0 X , the null vector of T X G m (R n ).
Definition 3.1. Let X ∈ G m (R n ), the injectivity radius of (G m (R n ), g) at X, is defined by r I X = sup {ρ > 0 : exp X : B(0 X , ρ) → exp X (B(0 X , ρ)) is a diffeomorphism} and the injectivity radius of G m (R n ) is defined by and its value is given by the following result: Let n, m be two positive integers such that min{m, n − m} ≥ 2. Then the injectivity radius of G m (R n ) is equal to π/2. The proof of Theorem 3.2 can be found in [14]. This Theorem gives the radius of the ball where the reciprocal mapping of the geodesic exponential is defined, it is denoted by exp −1 X or log X . More precisely, for 3.3. Sectional curvature and convexity. Definition 3.3. Let ξ, η and ζ be three vector fields on G m (R n ). The curvature tensor on the manifold (G m (R n ), g) is defined by where ∇ is the connection associated to g and [η, ξ] is the Lie bracket of η and ξ.
Using the metric, one can transform this (3, 1)−tensor to an (4, 0)−tensor, denoted again by R: R(u, v, w, z) = g(R(u, v, w), z). When the two vectors u and v are linearly independent, the expression R(u, v, u, v) depends only on span{u, v}, which justifies the definition of sectional curvature (that is a generalization of Gauss curvature).
Definition 3.4. Let X ∈ G m (R n ) and P a plane in T X G m (R n ). The sectional curvature of P is defined by κ(P, X) = R(u, v, u, v), where {u, v} is an orthonormal basis of P.
The following Theorem shows that G m (R n ) has a positive and bounded sectional curvature. The proof of Theorem 3.5 can be found in [6].
Definition 3.6. If is an upper bound of the sectional curvature in G m (R n ), then the convexity radius r c of G m (R n ) is defined by Applying Theorems 3.2 and 3.5, a direct computation gives that the convexity radius of the Grassmann manifold G m (R n ) is r c = π 4 . Definition 3.7. A ball B(X, r) ⊂ G m (R n ) is said to be strongly convex if any pair of points in B(X, r) can be joined by a unique geodesic, achieving the distance between them, and this geodesic is contained in B(X, r).
The next Theorem gives a sufficient condition on the radius of a ball to be strongly convex.
Theorem 3.8. Let X ∈ G m (R n ). If r < r c , then the ball B(X, r) is strongly convex.
The proof of Theorem 3.8 can be found in [22]. Applying Theorem 3.8, if r < π 4 then B(X, r) is strongly convex.

Grassmannian IDW (IDW-G).
Let {y(λ 1 ), · · · , y(λ N )} be N given points on the Grassmann manifold G m (R n ) corresponding to the set of parameters Λ = {λ 1 , · · · , λ N } where λ i ∈ R P , for every i ∈ {1, · · · , N } and P is a positive integer. In what follows, we will set y(λ i ) = y i . For every λ ∈ R P , we define the mapping associated to the points y 1 , · · · , y N by: where α 1 (λ), · · · , α N (λ) are nonnegative real numbers satisfying N i=1 α i (λ) = 1 and d is the Riemannian distance on G m (R n ). The weights (α i (λ)) 1≤i≤N may depend on the parameters (λ i ) 1≤i≤N . Now, consider the minimization problem For every λ ∈ R P , find y(λ) ∈ G m (R n ) such that : Since the Grassmann manifold G m (R n ) is compact and F λ is a continue mapping on G m (R n ) then Problem (P λ ) has at least one solution. From now on, we suppose that the integers m and n satisfy the following assumption: Theorem 4.1. Suppose that the points y 1 , · · · , y N are contained in a certain ball B(y * , r), where y * ∈ G m (R n ) and 0 < r < π 4 √ 2 . Then, for every λ ∈ R P , Problem (P λ ) has a unique solution y(λ) in B(y * , r) characterized by ∇F λ ( y(λ)) = 0, where Proof. We carry out an adaptation of the proof given by Karcher in the framework of the Grassmann manifold, [13] . The key point of the proof is the fact that the radius of the ball B(y * , r) where Problem (P λ ) has its unique solution satisfies 0 < r < π 4 √ 2 . Indeed, (i) The convexity radius of G m (R n ) is r c = π 4 . As mentioned before, the injectivity radius of G m (R n ) is r I = π 2 and its sectional curvature belongs to [0, 2], then the upper bound of the sectional curvature is = 2. Therefore, it follows that Moreover, Theorem 3.8 asserts that for 0 < r < r c , the ball B(y * , r) is strongly convex. (ii) For 0 < r < π 4 √ 2 , the mapping y → d 2 (y, z) is strictly convex on B(y * , r), for any z dans B(y * , r), [13,2]. Whence, (P λ ) has a unique solution in a ball of radius equal to min{ π 4 √ 2 , r c } = π 4 √ 2 , which achieves the proof.
is not other than the wighted barycenter of y 1 · · · , y N in the space V . For example, with the following choice of weights where S(λ) = N i=1 1 d(λ, λ i ) p and p > 0, we find the standard Inverse Distance Weighting (IDW) algorithm in vector spaces. In the Grassmann manifold G m (R n ), we will denote this algorithm by IDW-G.
In the Grassmann manifold G m (R n ), the equation (10) can not be solved in an explicit way, unlike in vector spaces. For this, we suggest the following iterative algorithm: For sake of simplicity, we denote Then, Lemma 5.1. If {y 1 , · · · , y N } ⊂ B(y * , r) then for any y ∈ B(y * , r), there exists 0 < t y ≤ 1 2 such that, if ∇F λ (y) = 0, it holds: 1. F λ • γ y is strictly decreasing on ]0, t y ], where γ y (t) = exp y (−t∇F λ (y)).
Proof. Let us define: Applying Lemma 4 in [16], we get that . On the other hand, since Therefore, d dt (F λ • γ y )(t) < 0, since ∇F λ (y) = 0 (11) and consequently, the function t −→ (F λ • γ y )(t) is strictly decreasing on ]0, t y ]. A direct integration of (11) yields to At this stage, we can state and show the following convergence result.
Proof. If there is * ∈ N such that ∇F λ ( y * ) = 0, then by Theorem 4.1, y * is the solution of (P λ ) and by construction y = y * for any ≥ * , which ends the claim. Now, if ∇F λ ( y ) = 0 pour tout ∈ N, then by Lemma 5.1, we get Hence, the sequence (F λ ( y )) is convergent. On one hand, let y ∈ B(y * , r) be an adherence value of the sequence ( y ) , then (F λ ( y )) converges to F λ ( y), by continuity. On the other hand, since then, by the continuity of the mappings y −→ t y and y −→ ∇F λ y on B(y * , r), we get 0 ≥ ∇F λ ( y) and thus ∇F λ ( y) = 0. Hence, the sequence ( y ) has y as unique adherence value. Which achieves the proof.
The IDW-G algorithm can be summarized as follows: ii. The new parameter: λ ∈ R P . iii. A tolerance number: ε.
if Res ≤ ε then Return y = R else Compute again R and Res end Algorithm 1: IDW-G algorithm.
In Figures 4 and 5, the interpolation algorithm developed in [7], our Grassmann IDW algorithm and the full solution computed by the finite elements method are referenced respectively by "Grass", "IDW-G" and "Full". length-height space variables will be denoted by (x, y) and the planar velocity field will be denoted by (u, v). At every point (x, y) on the boundary of the domain, n(x, y) will denote the unit normal outward vector to the boundary. Finally, the euclidean inner product in the plane will be denoted by "·". The boundary conditions are as follows : • On the vertical boundaries : a horizontal inflow (u(0, y), v(0, y)) = (U, 0) is applied on the left boundary. For the right boundary, a homogeneous output boundary condition is imposed. More precisely, we have on the right boundary : σ(u, v) · n = (0, 0), where σ denotes the stress tensor. We assume that the fluid dynamics is governed by the classical adimensional incompressible Navier-Stokes equations. The Reynolds number is defined by Re = U D ν , where ν is the fluid viscosity.
To built the POD bases, the used snapshots are computed with the free Finite Elements software Fenics using Taylor-Hood P 2 /P 1 elements. The mesh including 17290 elements is non uniform and refined near the cylinder. To ensure homogeneous boundary conditions, the velocity field and the pressure are splitted into the mean and fluctuating parts as follows: u(x, t, Re) = u(x) + u (x, t, Re) and p(x, t, Re) = p(x) + p (x, t, Re), where u(x), p(x) are the mean parts of the velocity field and the pressure respectively, and POD is then applied to the fluctuating velocity u and the fluctuating pressure p . By keeping the only most POD energetic vectors : the temporal coefficient a i (i = 1, . . . , N u ) and b i (i = 1, . . . , N p ) are computed by solving a set of differential equations of small size N u + N p . This reduced order model, that is very fast to solve, is obtained by applying Galerkin projection of the momentum equations onto each velocity and pressure POD vectors Φ u i , ∇Φ p i respectively (for more details, we refer the interested reader to [17] [26]). For this application, we fix N u = N p = 8. A bunch of velocity and pressure POD bases corresponding to five Reynolds numbers (Re bunch = [70, 80, 90, 150, 180]) are computed. This set of POD bases is interpolated on the Grassmann manifold by the method developed by Amsallem "Grass" and by our approach "IDW-G", to obtain a valid basis for different Reynolds numbers (Re int = 140, 160,170,190). For each Reynolds number of the bunch, 150 snapshots uniformly distributed on the time of two oscillations period of the wake are chosen to build the spatial POD basis. For these Reynolds numbers, the flow oscillates with Von Karman streets and a vertex shedding is observed (see Fig 2). The availability of the two interpolation approaches to compute the drag coefficient C D and the lift coefficient C L are tested. Theses coefficients are expressed such as: where F D , F L are the drag and the lift forces exerted by the fluid on the cylinder respectively and ρ is the density of the fluid.
In the IDW-G interpolation method, one interesting parameter is the power p appearing in the weight coefficient given in (9). The influence of this power p is drawn on Figure (Fig. 3). It is observed that the value p = 1 does not give satisfactory results. For higher values of the power p, the results are more accurate. Hereafter, we will choose p = 3. Figures (Fig. 4) and (Fig. 5)    of the ROM built with the two interpolations approaches and predictions of the full model. This figures show that our IDW-G proposed method gives more accurate results than the the Grassmann one, moreover the drag and lift coefficients are properly recovered. 7. Conclusion. We have adapted the Inverse Distance Weighting method on the Grassmann manifold to perform interpolation of parametric POD bases. The advantages of our approach is to overcome the choice of a reference point on the Grassmann manifold and offers a better quality of results comparing to those obtained by Amsallem et al. [7]. In return, our approach is not direct but iterative, moreover its relevance depends on the choice of the weighting function: The numerical experiments have shown that the results are more accurate for p ≥ 2.
A more judicious choice of such weighting functions can be carried out via kriging technics [23], this is the subject of a work in progress. Moreover, we notice  that the proposed method can be applied not only for interpolation but also for extrapolation, as it can be observed for the value Re = 190, which is outside of the interpolation set. This observation will be studied more carefully in a future work.