# Méthode de réduction de modèle a priori basée sur des formulations idéales en minimum de résidu

Olivier Zahm, Marie Billaud-Friess, Anthony Nouy

► To cite this version:

Olivier Zahm, Marie Billaud-Friess, Anthony Nouy. Méthode de réduction de modèle a priori basée sur des formulations idéales en minimum de résidu. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717098

HAL Id: hal-01717098

https://hal.science/hal-01717098v1

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Méthode de réduction de modèle *a priori* basée sur des formulations idéales en minimum de résidu

Olivier Zahm, Marie Billaud-Friess, Anthony Nouy

LUNAM Université, GeM, UMR CNRS 6183, Ecole Centrale de Nantes, Université de Nantes, France,
(olivier.zahm,marie.billaud-friess,anthony.nouy)@ec-nantes.fr

**Résumé** — La résolution numérique d'équations aux dérivées partielles stochastiques, nécessitant l'approximation de fonctions de nombreux paramètres, se heurte au problème de dimensionnalité et reste difficile voire impossible avec les méthodes de calcul traditionnelles. Les méthodes de réduction de modèle basées sur l'approximation de tenseurs de faible rang permettent de contourner le problème de dimensionnalité. Cependant, les algorithmes disponibles conduisent généralement à des approximations sous-optimales au sens de normes d'intérêt. Nous proposons ici une approche en minimum de résidu avec utilisation de normes dites idéales dans le but d'obtenir la meilleure approximation de la solution dans un ensemble de tenseurs donné au sens de normes choisies.

**Mots clés** — Réduction de modèle, Approximation de tenseurs, Minimum de résidu, Goal-oriented

## 1 Introduction

Ces dernières années, des efforts ont été réalisés pour développer de nouvelles méthodes numériques basées sur une approche fonctionnelle des incertitudes, ces dernières étant intégrées comme dimensions supplémentaires aux modèles. La résolution numérique de ces problèmes à grand nombre de variables reste difficile voire impossible par les méthodes classiques. Dans ce contexte, des méthodes de réduction de modèle, telles que la méthode des bases réduites [8] ou les méthodes PGD [1], ont récemment été introduites pour diminuer la taille des problèmes considérés, et rendre leur résolution accessible.

On s'intéresse ici à un problème linéaire de la forme $Au = b$. Nous proposons de chercher une approximation de la solution dans un ensemble de tenseurs $S$ tel que les ensembles canoniques ou les ensembles de Tucker de faible rang [5]. Le problème de meilleure approximation associé consiste à chercher un élément $u_S \in S$ vérifiant :

$$ \|u - u_S\| = \min_{v \in S} \|u - v\|. \quad (1) $$

où $\| \cdot \|$ désigne une norme d'intérêt. La solution $u$ étant inconnue, on se ramène classiquement à un problème de minimisation de résidu $\|Av - b\|_*$ mesuré avec une certaine norme $\| \cdot \|_*$. Si cette norme est mal adaptée, on observe que l'approximation obtenue est très sous-optimale vis-à-vis de la norme $\| \cdot \|$. Dans le but d'obtenir *a priori* une approximation optimale $u_S$ de la solution par une formulation en minimum de résidu, la norme $\| \cdot \|_*$ devrait vérifier :

$$ \|Av - b\|_* = \|v - u\|. \quad (2) $$

Cependant, le calcul explicite de cette norme de résidu est d'une complexité équivalente au problème de départ, ce qui n'est pas viable numériquement. En s'inspirant de l'approche introduite dans [2, 3], nous proposons de définir l'approximation de la manière suivante :

$$ \min_{v \in S} \|\Lambda^{\delta}(Av - b)\|_* \quad (3) $$

où $\Lambda^{\delta}$ est un opérateur fournissant une approximation $\Lambda^{\delta}(Av - b)$ du résidu $Av - b$ avec une précision relative contrôlée $\delta$. Pour la résolution de ce problème d'approximation, nous proposons un algorithme
itératif dont nous prouvons la convergence vers un voisinage de la meilleure approximation $u_S$. De plus, cette méthode fournit naturellement un estimateur d'erreur *a posteriori* garanti.

Dans la première partie nous introduisons l'approche fonctionnelle utilisée pour résoudre les équations aux dérivées partielles stochastiques. La projection de Galerkin sur un sous-espace d'approximation nous amène à considérer un système algébrique formulé sur un espace produit tensoriel. La deuxième partie introduit certains outils d'approximation de faible rang de tenseurs. Puis nous définissons le problème en minimum de résidu idéal et nous proposons qu'un algorithme permettant de résoudre ce problème. Enfin, la robustesse de la méthode est illustrée à l'aide de résultats numériques sur une équation de type réaction-advection-diffusion.

Dans ce papier, nous restreignons la présentation de la méthode à un cadre algébrique de dimension finie. Les résultats et algorithmes peuvent cependant s'étendre à un cadre plus général [7].

## 2 Équation aux dérivées partielles stochastique

Nous considérons ici une équation aux dérivées partielles (EDP) stochastique définie sur un domaine $\Omega$. Certains coefficients de l'EDP sont fonction d'un vecteur aléatoire $\xi$ supposé de dimension finie $d \ge 1$. L'espace de probabilité associé à $\xi$ est noté $(\Xi, \mathcal{B}, P_\xi)$, où $\Xi \subset \mathbb{R}^d$ est l'ensemble des réalisations, $\mathcal{B}$ une $\sigma$-algèbre et $P_\xi$ la mesure de probabilité.

Pour une réalisation $\xi \in \Xi$ donnée, la formulation variationnelle associée à l'EDP consiste à trouver $u(\xi) \in V$ vérifiant $a(u(\xi), v; \xi) = b(v; \xi)$ pour tout $v \in V$. Ici, $V$ désigne un espace de Hilbert approprié de fonctions définies sur $\Omega$: on supposera que cet espace ne dépend pas de $\xi$.

Nous considérons la formulation faible associée au problème stochastique. L'espace de variables aléatoires du second ordre à valeurs dans $V$ est noté $L^2(\Xi, P_\xi; V)$. Cet espace possède une structure produit tensoriel, c'est à dire $L^2(\Xi, P_\xi; V) = V \otimes S$ avec $S = L^2(\Xi, P_\xi)$ l'ensemble des variables aléatoires du second ordre à valeurs réelles. La solution faible $u \in V \otimes S$ satisfait :

$$ 
\begin{gathered}
\mathcal{A}(u, v) = \mathcal{B}(v) \quad \forall v \in V \otimes S, \qquad (4) \\
\text{avec } \mathcal{A}(\cdot, \cdot) = \int_{\Xi} a(\cdot, \cdot; \xi) \, dP_{\xi}, \\
\mathcal{B}(\cdot) = \int_{\Xi} b(\cdot; \xi) \, dP_{\xi}.
\end{gathered}
$$

Pour que ce problème soit bien posé, la forme bilinéaire $\mathcal{A}$ sera supposée continue sur $(V \otimes S) \times (V \otimes S)$ et $(V \otimes S)$-coercive (voire faiblement coercive), et la forme linéaire $\mathcal{B}$ continue sur $(V \otimes S)$.

La dicrétisation du problème (4) consiste dans un premier temps à introduire un sous-espace vectoriel de $V \otimes S$. La partie spatiale est classiquement discrétisée avec la méthodes des éléments finis. Le sous espace associé, noté $V^n$, est de dimension $n$. Quant à la discrétisation stochastique plusieurs possibilités existent pour le choix d'approximation : bases de type polynômial ou bien polynômial par morceaux (éléments finis, ondelettes...) [6, 10]. Ce choix dépend de la régularité de la solution. Le sous espace d'approximation de dimension $p$ est noté $S^p$.

Dans un deuxième temps, une approximation de la solution est cherchée dans le sous espace $(V^n \otimes S^p) \subset (V \otimes S)$. Cette approximation $u \in V^n \otimes S^p$ vérifie $\mathcal{A}(u, v) = \mathcal{B}(v)$ pour tout $v \in V^n \otimes S^p$. De manière équivalente ce problème se met sous forme algébrique

$$ Au = b. \qquad (5) $$

L'opérateur $A$ et le second membre $b$ correspondent aux représentants de la restriction de $\mathcal{A}$ et $\mathcal{B}$ au sous espace $(V^n \otimes S^p)$ associé à une certaine base. La solution $u \in \mathbb{R}^n \otimes \mathbb{R}^p$ représente les coordonnées de l'approximation sur la base considérée.

Notons que le nombre de degrés de liberté de ce système est de $np$, et peut donc facilement atteindre les limites des méthodes de résolution classiques. La structure produit tensoriel de $\mathbb{R}^n \otimes \mathbb{R}^p$ permet d'appliquer des méthodes d'approximation de faible rang de tenseurs.

## 3 Approximation de faible rang

Nous introduisons dans cette partie certains outils d'approximation de faible rang de tenseurs.
## 3.1 Approximation de rang r

Dans $\mathbb{R}^n \otimes \mathbb{R}^p$ l'ensemble des tenseurs de rang $r$, noté $C_r(\mathbb{R}^n \otimes \mathbb{R}^p)$, est une variété classiquement utilisée pour l'approximation de tenseur :

$$ C_r(\mathbb{R}^n \otimes \mathbb{R}^p) = \left\{ v = \sum_{i=1}^{r} v_i^1 \otimes v_i^2 ; \ v_i^1 \in \mathbb{R}^n, v_i^2 \in \mathbb{R}^p \right\} \quad (6) $$

Cet ensemble étant fermé, le problème de meilleure approximation (au sens d'une norme donnée $\| \cdot \|$) d'un tenseur $u \in \mathbb{R}^n \otimes \mathbb{R}^p$ a bien un sens. Nous définissons l'opérateur de "projection" sur $C_r$ par :

$$ \Pi_{C_r}(u) = \arg \min_{v \in C_r} \|u - v\|. \quad (7) $$

L'algorithme de minimisation alternée permet d'approcher numériquement un minimiseur. Le principe est de minimiser alternativement sur les paramètres $v_i^1$ et $v_i^2$ (les détails sont donnés dans [5]). Cet algorithme, inspiré de la méthode des sous-espaces pour les problèmes aux valeurs propres, fournit une approximation optimale.

Cependant, la meilleure approximation d'un certain rang $r$ ne garantit pas d'obtenir une approximation de $u$ à une précision $\delta$ désirée.

## 3.2 Algorithme Greedy, contrôle en précision

Les algorithmes greedy consistent en une succession d'approximations dans un ensemble de faible rang (par exemple de rang 1). Une séquence $\{u_m\}_m$ est construite itérativement : à chaque étape, une correction $w_{m+1} \in \Pi_{C_1}(u - u_m)$ est calculée puis ajoutée à $u_m : u_{m+1} = u_m + w_{m+1}$. La terminologie "greedy" fait référence au fait que la $(m+1)^{\text{ème}}$ itération ne rectifie pas les $m$ corrections précédemment calculées. Une version améliorée de cet algorithme consiste à ajouter à chaque itération une étape de projection sur le sous-espace span($w_i$). On parlera d'algorithme "orthogonal-greedy". La convergence de la suite $u_m$ vers $u$ est démontrée dans [9]. Des stratégies d'update plus sophistiquées sont analysées dans [4].

<caption>**Algorithm 1 Orthogonal greedy**</caption>

1: Initialisation $u_0 = 0$;
2: **while** $e(u_m, u) \le \delta$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Correction $w_{m+1} \in \Pi_{C_1}(u - u_m)$
4: &nbsp;&nbsp;&nbsp;&nbsp;Mise à jour $W_{m+1} = W_m + \text{span}\{w_{m+1}\}$;
5: &nbsp;&nbsp;&nbsp;&nbsp;Projection $u_{m+1} = \arg \min_{v \in W_{m+1}} \|u - v\|$;
6: &nbsp;&nbsp;&nbsp;&nbsp;$m = m + 1$;
7: **end while**

Dans le but d'obtenir une approximation à une précision $\delta$ donnée, l'algorithme est arrêté par la condition

$$ e(u_m, u) = \frac{\|u_m - u\|}{\|u\|} \le \delta. \quad (8) $$

Dans le cas où $u$ n'est pas connu a priori, un estimateur d'erreur de type stagnation peut être utilisé. Pour $p > 1$ fixé, la condition d'arrêt devient :

$$ e(u_m, u_{m+p}) = \frac{\|u_m - u_{m+p}\|}{\|u_{m+p}\|} \le \delta. \quad (9) $$

# 4 Approches en minimum de résidu

Nous appelons $X$ l'espace $(\mathbb{R}^n \otimes \mathbb{R}^p, \| \cdot \|_X)$ où $\| \cdot \|_X$ désigne une norme d'intérêt. La métrique associée à cette norme est noté $R_X$ : elle vérifie $\| \cdot \|_X = \sqrt{\langle R_X \cdot, \cdot \rangle}$, où $\langle \cdot, \cdot \rangle$ est le produit scalaire canonique
sur $\mathbb{R}^n \otimes \mathbb{R}^p$. Nous proposons ici une approche permettant de trouver une meilleure approximation de $u$, solution de (5), au sens de la norme $\| \cdot \|_X$ :

$$ \min_{v \in C_r} \|u - v\|_X \qquad (10) $$

## 4.1 Minimum de résidu idéal et formulation approchée

Comme $u$ n’est pas connu a priori, nous proposons d’utiliser une formulation en minimum de résidu. Les relations suivantes :

$$ \begin{aligned} \|u - v\|_X^2 &= \langle R_X(u - v), (u - v) \rangle \\ &= \langle R_X A^{-1} A(u - v), A^{-1} A(u - v) \rangle \\ &= \langle A^{-*} R_X A^{-1} A(u - v), A(u - v) \rangle \\ &= \| R_Y^{-1} (A v - b) \|_Y^2, \end{aligned} $$

suggèrent d’introduire l’espace $Y = (\mathbb{R}^n \otimes \mathbb{R}^p, \| \cdot \|_Y)$ avec $R_Y = A R_X^{-1} A^*$ la métrique associée à la norme $\| \cdot \|_Y$. Ainsi, le problème (10) est équivalent à :

$$ \min_{v \in C_r} \| R_Y^{-1} (A v - b) \|_Y. \qquad (11) $$

Cette formulation fait intervenir l’inverse de l’opérateur $R_Y$ qu’il n’est pas envisageable de calculer dans la pratique. En s’inspirant de l’approche introduite dans [2, 3], nous proposons de définir l’approximation de la manière suivante :

$$ \min_{v \in C_r} \| \Lambda^{\delta} R_Y^{-1} (A v - b) \|_Y. \qquad (12) $$

Ici, $\Lambda^{\delta}$ est un opérateur d’approximation à $\delta$ près, c’est à dire que pour tout $y \in Y$

$$ \| \Lambda^{\delta}(y) - y \|_{Y} \le \delta \| y \|_{Y} \qquad (13) $$

## 4.2 Algorithme itératif

Pour la résolution de (12) nous proposons un algorithme itératif. Après initialisation de $u^0 = 0$, les suites $\{y^k\}_{k>0}$ et $\{u^k\}_{k>0}$ sont construites comme :

$$ y^k = \Lambda^\delta (R_Y^{-1} (A u^k - b)) \qquad (14) $$

$$ u^{k+1} \in \Pi_{C_r} (u^k - R_X^{-1} A^* y^k). \qquad (15) $$

Les résultats qui suivent sont démontré dans [7] :

- Sous réserve que $2\delta < 1$, que la suite $\{u^k\}_{k>0}$ vérifie

$$ \|u^k - u\|_X \le (2\delta)^k \|u\|_X + \frac{1}{1 - 2\delta} \|u - \Pi_{C_r}(u)\|_X, \qquad (16) $$

ce qui nous assure que $\limsup_k \|u^k - u\|_X \le \frac{1}{1-2\delta} \|u - \Pi_{C_r}(u)\|_X$. En d'autre termes, cet algorithme permet de trouver une approximation quasi-optimale de $u$ dans $C_r$ au sens de la norme $\| \cdot \|_X$.

- La norme de $y^k$ fournit un estimateur d’erreur garanti de $u^k$ via l’encadrement suivant :

$$ \|y^k\|_Y \le \|u^k - u\|_X \le \frac{1}{\sqrt{1 - \delta^2}} \|y^k\|_Y \qquad (17) $$
## 4.3 Remarques computationnelles

La métrique $R_X$ est supposée être de rang 1 (norme $\| \cdot \|_X$ induite), et son inverse calculable facilement (dans la partie suivante 5, nous nous limitons aux opérateurs diagonaux), de sorte que $R_Y = AR_X^{-1}A^*$ soit calculé explicitement. Nous donnons quelques détails sur les étapes (14) et (15) :

- Étape (14) : l'approximation à $\delta$ près de $R_Y^{-1}(Au^k - b)$ est calculée avec l'algorithme 1 en choisisnant la norme $\| \cdot \|_Y$. Comme cette quantité n'est pas connue a priori, l'estimateur d'erreur en stagnation (9) est utilisé pour la condition d'arrêt (avec $p$ suffisamment grand pour assurer une estimation précise).

- Étape (15) : comme $\| \cdot \|_X$ est une norme induite, l'approximation de rang $r$ correspond à une décomposition en valeurs singulières (SVD). Un algorithme greedy fournit dans ce cas une approximation optimale.

# 5 Illustration numérique

Afin d'illustrer la robustesse de notre méthode, nous considérons une équation de réaction-advection-diffusion sur un domaine carré $\Omega = [0, 1]^2$ :

$$ -\Delta u(x, \xi) + c(x, \xi) \cdot \nabla u(x, \xi) + a(\xi)u(x, \xi) = f(x) \quad \text{dans } \Omega, \tag{18} \\ u = 0 \quad \text{sur } \partial\Omega. $$

Le champ advectif $c$ ainsi que le terme de réaction $a$ sont aléatoires et tels que $c(x, \xi) = (x_2 - 1/2, 1/2 - x_1)\xi_1$ et $a = \xi_2$. Les variables aléatoires $\xi_1$ et $\xi_2$ sont indépendantes et ont pour loi de probabilité respectivement la loi uniforme sur $[-350, 350]$ et log-uniforme sur $[0.1, 10]$. Comme décrit sur la figure 1, nous imposons dans deux sous domaines $\Omega_1 = [0.45, 0.55] \times [0.15, 0.25]$ et $\Omega_2 = [0.45, 0.55] \times [0.75, 0.85]$ un terme source $f = 1$ sur $\Omega_1$ et $f = -1$ sur $\Omega_2$.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Géométrie du cas test.</figcaption></figure>

La méthode des éléments finis est utilisée pour la discrétisation spatiale : éléments de type quadrangle à 4 noeuds. La dimensions de $V^n$ s'élève à 1500. La partie stochastique est discrétisé avec une approximation polynômials par moceaux de dimension 72. La dimension totale étant de $10^5$, la solution de référence peut être calculée à partir du système algébrique (5).

## 5.1 Norme adaptée à une quantité d'intérêt

Une particularité de la méthode proposée est de pouvoir choisir la norme $\| \cdot \|_X$ de façon arbitraire. Pour illustrer cet avantage, nous introduisons deux normes. D'abord la norme canonique, notée $\| \cdot \|_2$, dont la métrique associée est $R_X = I_d^n \otimes I_d^p$, où $I_d^n$ (resp. $I_d^p$) est l'opérateur identité de dimension $n$ (resp. $p$).

Supposons que le but de la simulation soit d'estimer avec précision une certaine quantité d'intérêt $Q(u)$ qui ne dépend que des valeurs de la solution $u$ sur un sous domaine $D = [0.15, 0.25] \times [0.45, 0.55] \subset \Omega$. Par exemple, la moyenne spatiale de $u$ sur le domaine $D$ est donnée par :

$$ Q(u) = \frac{1}{|D|} \int_D u(\xi, x) \, dx. \tag{19} $$
Dans ce cas, il semble pertinent d'utiliser une norme qui associe plus d'importance (de poids) au sous domaine $D$. Soit une fonction poids $w : \Omega \to \mathbb{R}$ qui vaut $10^3$ sur $D$ et 1 ailleurs. La norme pondérée $\|\cdot\|_w$ est définie par

$$ \|v\|_w = \sqrt{\sum_{i=1}^{n} \sum_{j=1}^{p} (w(x_i)v_{ij})^2}. \qquad (20) $$

La métrique associée est de la forme $R_X = \text{diag}(w(x_i)^2) \otimes I_d^p$.

Afin de comprendre l'impact du choix de la norme $\|\cdot\|_X$, nous notons $\tilde{u}_w$ (resp. $\tilde{u}_2$) la $\|\cdot\|_w$-meilleure (resp. $\|\cdot\|_2$-meilleure) approximation de $u$ dans $C_r$. D'après la figure 2, il est clair que $\tilde{u}_w$ fournit une estimation bien plus précise de la moyenne et de la variance de $Q(u)$ que $\tilde{u}_2$ ce qui justifie l'utilisation de la norme pondérée. La figure 3 représente les premiers modes déterministes de $\tilde{u}_w$ et $\tilde{u}_2$: elle montre comment le choix de la norme influence la décomposition.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Convergence en fonction du rang $r$ de l'estimation de la moyenne et de la variance de $Q(u)$ par celles de $Q(\tilde{u}_w)$ et $Q(\tilde{u}_2)$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Comparaison des premiers modes déterministes de $\tilde{u}_2$ et $\tilde{u}_w$.</figcaption></figure>

## 5.2 Analyse des algorithmes

Cette partie a pour but d'illustrer les résultats théoriques annoncés au paragraphe 4.2. Nous nous limitons à l'étude du problème de meilleure approximation de rang $r = 10$ de $u$, au sens de la norme pondérée $\|\cdot\|_w$. La figure 4 présente les courbes de convergence de l'algorithme itératif pour différentes valeurs du paramètre $\delta$. On retrouve quantitativement les résultats théoriques, à savoir une première phase de convergence linéaire suivie d'une phase de stagnation dans un voisinage de la meilleure approximation.

Sur la figure 5 la courbe d'erreur $\|u - u^k\|_w$ est clairement encadrée par $\|y^k\|_Y$ et $\frac{1}{\sqrt{1-\delta^2}} \|y^k\|_Y$, ce qui est en accord avec le résultat annoncé.
<figure><img src="image_5.png" /><figcaption>Fig. 4 – Convergence de la suite $u^k$ pour différentes valeurs de $\delta$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Encadrement de l’erreur $\|u - u^k\|_w$ par $\|y^k\|_Y$ et $\frac{1}{\sqrt{1-\delta^2}}\|y^k\|_Y$, avec $\delta = 0.9$.</figcaption></figure>

# 6 Conclusion

Une méthode d’approximation de solution d’EDP stochastique a été proposée. Elle permet d’obtenir une approximation quasi-optimale de faible rang au sens d’une norme choisie. Cette approche ouvre de nombreuses possibilités, comme celle de construire des approximations adaptées à des quantités d’intérêt. La robustesse de la méthode a été montrée pour une équation de type réaction-advection-diffusion à coefficients stochastiques.

# Références

[1] Chinesta F., Ladeveze P. and Cueto E. *A Short Review on Model Order Reduction Based on Proper Generalized Decomposition*, Archives of Computational Methods in Engineering, Volume 18, Number 4, 395-404, 2011.

[2] Cohen A., Dahmen W. and Welper G. *Adaptivity and Variational Stabilization for Convection-Diffusion Equations*, ESAIM : Mathematical Modelling and Numerical Analysis, 2012.

[3] Dahmen W., Huang C. and Schwab C. *Adaptative petrov-galerkin methods for first order transport equations*, IGPM Report 321, RWTH Aachen, 425-467, 2011.

[4] Falcó A. and Nouy A. *Proper generalized decomposition for nonlinear convex problems in tensor Banach spaces*, Numerische Mathematik, Volume 121, Number 3, 503-530, 2012.

[5] Hackbusch W. *Tensor Spaces and Numerical Tensor Calculus*, Springer Series in Computational Mathematics, Vol. 42, 2012.

[6] Nouy A. *Recent Developments in Spectral Stochastic Methods for the Numerical Solution of Stochastic Partial Differential Equations*, Archives of Computational Methods in Engineering, 2009.

[7] Nouy A., Billaud M. and Zahm O. *A tensor approximation method based on ideal minimal residual formula-tions for the solution of high dimensional problems*, in preparation, 2013.

[8] Patera A.T. and Rozza G. Reduced Basis Approximation and A-Posteriori Error Estimation for Parametrized PDEs. MIT-Pappalardo Graduate Monographs in Mechanical Engineering, 2012.

[9] Temlyakov V. *Greedy Approximation*, Cambridge Monographs on Applied and Computational Mathematics, 2011.

[10] Xiu D. *Fast numerical methods for stochastic computations : a review*, Communications in computational physics, Volume 121, Number 3, 503-530, 2009.