Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models

This paper revisits a powerful discretization technique, the Proper Generalized Decomposition -PGD-, illustrating its ability for solving highly multidimensional models. This technique operates by constructing a separated representation of the so-lution, such that the solution complexity scales linearly with the dimension of the space in which the model is deﬁned, instead the exponentially-growing complexity characteristic of mesh based discretization strategies. The PGD makes possible the eﬃcient solution of models deﬁned in multidimensional spaces, as the ones encountered in quantum chemistry, kinetic theory description of complex ﬂuids, genetics (chemical master equation), ﬁnancial mathematics, ... but also those, classically deﬁned in the standard space and time, to which we can add new extra-coordinates (parametric models, ...) opening numerous possibilities (optimization, inverse identiﬁcation, real time simulations, ...).


Introduction
Many models encountered in science and engineering are defined in multidimensional spaces, as the ones involved in quantum chemistry, kinetic theory descriptions of materi-als (including complex fluids), the chemical master equation governing many biological processes (e.g.cell signaling), models of financial mathematics (e.g.option pricing), among many others.These models exhibit the redoubtable curse of dimensionality when usual mesh-based discretization techniques are applied.Other times, standard models can become multidimensional if some of the parameters that they involve are considered as new coordinates.This possibility is specially attractive when these coefficients are not well known, they have a stochastic nature, or when one is interested in optimization or inverse identification.
The difficulty related to the solution of multidimensional models is quite obvious and it needs the proposal of new appropriate strategies able to circumvent the curse of dimensionality.One possibility lies in the use of sparse grids [12].However, as argued in [1], the use of sparse grid is restricted to models with moderate multidimensionality (up to 20).Another technique able to circumvent, or at least alleviate, the curse of dimensionality consists of using a separated representation of the unknown field (see [44] [10] for some numerical elements on this topic).Basically, the separated representation of a generic function u(x 1 , • • • , x D ) (also known as finite sum decomposition) writes: Remark 1 The coordinates x i , i = 1, • • • , D, are defined in spaces of moderate dimension, i.e. x i ∈ Ω i ⊂ R di , d i ≤ 3. Thus, the dimension of the model results P i=D i=1 d i .Eventually, one of these coordinates could be the time t ∈ I ⊂ R + . 2 This kind of representation is not new, it was widely employed in the last decades in the framework of quantum chemistry.In particular the Hartree-Fock (that involves a single product of functions) and post-Hartree-Fock approaches (as the MCSCF that involves a finite number of sums) made use of a separated representation of the wavefunction [13] [19].
We proposed recently a technique able to construct, in a way completely transparent for the user, the separated representation of the unknown field involved in a multidimensional partial differential equation.This technique, originally described and applied to multi-bead-spring FENE models of polymeric liquids in [3], was extended to transient models of such complex fluids in [4].More complex models (involving different couplings and non-linearities) based on the reptation theory of polymeric liquids were analyzed in [34].
In the current life when one is buying a car there are many choices concerning the color, the power, ... In some cases the number of possible variants reaches astronomical values: 10 50 in the case of a highly appreciated luxury car.Thus, the construction of a high dimensional matrix representing all the possible variants is simply impossible (we must recall that the presumed number of elementary particles in our universe is of around 10 80 !). Obviously, the only possibility of following the sales of such a car is using a sparse storage structure.This is well known and widely used in practice.However, our present interest is not only in storing information, but in solving models defined in multidimensional spaces.
When the solution of highly multidimensional models explores the whole space there is no possibility to follow or represent such a system.This is a well known phenomenon in quantum systems.However, one could track the evolution of such a system in a time interval of moderate length by using a "sparse" representation.In our knowledge two possibilities exist, the first one consists in using Monte Carlo simulations within a stochastic framework.The second one, within a purely deterministic framework, lies in using separated representations like the one expressed by Eq. (1) where the number of sums increases as the solution explores the more and more larger domains.On the other hand, if the solution only explores a small enough subdomain of the whole domain the use of sparse grids or separated representations are some appealing candidates for discretizing such models.
Coming back to models defined in spaces of moderate dimension (d×D, d = 1, 2, 3) but whose solutions evolve in large time intervals, if one uses standard incremental time-discretizations, in the general case (models involving time-dependent parameters, non-linear models, ...), one must solve at least a linear system at each time step.When the time step becomes too small as a consequence of stability requirements, and the simulation time interval is large enough, standard incremental simulation becomes inefficient.To illustrate this scenario, one could imagine the simple reaction-diffusion model that describes the degradation of plastic materials, where the characteristic time of the chemical reaction involved in the material degradation is in the order of some microseconds and the one related to the diffusion of chemical substances (that also represents the material degradation characteristic time itself) is of the order of years.In this case standard incremental techniques must be replaced by other more efficient strategies.
One possibility consists again in performing a separated representation of the unknown field, that in the present case reduces to: that allows, as we describe later, to non-incremental time integration strategies, which can reduce spectacularly the CPU time.This space-time separated representation is not a new proposal.In fact such decompositions were proposed many year ago by Pierre Ladeveze as an ingredient of the powerful non-linear-non-incremental LATIN solver that he proposed in the 80s.During the last twenty years many works were successfully accomplished by the Ladeveze's group.The interested reader can refer to [29] [30] [32] and the references therein.In the radial approximation approach (the name given in the pioneer works of Ladeveze) functions depending on space and the ones depending on time were a priori unknown, and they were computed by an appropriate minimization technique.
In what follows we are reporting the recent advances in the solution of multidimensional models by applying the Proper Generalized Decomposition (PGD), that is, a separated representation of the unknown fields.

Motivating the use of separated representations
From a historical point of view, separated representations were extensively used.The analytical solution of PDEs (elliptic, parabolic and hyperbolic) that the reader can find in any book devoted to the solution of PDEs starts assuming a separated representation: the method of separation of variables.After some manipulations, simple but sometimes quite technical, the final solution could be some times found in the form expressed by Eq. (1).
Another well established and widely employed technique allowing to define a separated representation of a given space-time function is based on the application of a proper orthogonal decomposition.We are illustrating the main ideas related to this technique.
Let u(x, t) be the solution of a certain transient model (in what follows x ∈ Ω ⊂ R d , d = 1, 2, 3, and t ∈ I ⊂ R + ).We are also assuming that this field is known in a discrete manner, that is, at some points x i (the nodes of a mesh or a grid) and at certain times Now, we introduce the notation u p i ≡ u(x i , tp) and construct the matrix Q that contains the snapshots: The proper orthogonal decomposition (POD) of this discrete field consists in solving the eigenvalue problem: that results in Nn couples eigenvalue-eigenvector When the field evolves smoothly, the magnitude of the eigenvalues decreases very fast, fact that reveals that the evolution of the field can be approximated from a reduced number of modes (eigenvectors).Thus, if we define a cutoff value ǫ (ǫ = 10 −8 × λ 1 in practice, λ 1 being the highest eigenvalue) only a reduced number of modes are retained.Let R (R << Nn) be the number of modes retained, i.e.
the eigenvalues are assumed ordered).Thus, one could write: where for the sake of clarity the space modes φ i (x) will be, from now on, denoted as X i (x).Eq. ( 5) represents a natural separated representation (also known as finite sums decomposition).These modes could be now used to solve other "similar" problems, that is, models involving slight changes in the boundary conditions, model parameters, ... [39] [33] [46].Other possibility is to compute the reduced basis from the standard transient solution within a short time interval (with respect to the whole time interval in which the model is defined) and then solve the remaining part of the time interval by employing the reduced basis.Obviously, both strategies induce the introduction of an error whose evaluation, control and reduction is a challenging issue.
One possibility to construct an adaptive reduced approximation basis, that should be the best reduced approximation basis for the treated problem, consists in alternating a reduction step (based on the application of the proper orthogonal decomposition) and an enrichment stage to improve the quality of the reduced approximation basis in order to capture all the solution features.We proposed recently an enrichment technique based on the use of some Krylov's subspaces generated by the equation residual.This technique known as "a priori" model reduction was originally proposed in [47], widely described in [48] and successfully applied for solving complex fluid flows within the kinetic theory framework [2] [17] and for speeding up thermomechanical simulations [18].However, some difficulties were noticed in the application of this strategy: (i) the enrichment based on the use of the Krylov's subspaces is far to be optimal in a variety of models (e.g. the wave equation); (ii) the incremental nature of the algorithm; ...
From the previous analysis we can conclude: (i) the transient solution of numerous models can be expressed using a very reduced number of products each one involving a function of time and a function of space; and (ii) the functions involved in these functional products should be determined simultaneously by applying an appropriate algorithm to guarantee robustness and optimality.
In what follows we are illustrating the simplest strategy able to compute these separated functional couples.

Illustrating the Proper Generalized Decomposition
In this section we are illustrating the discretization of partial differential equations using a separated representation (radial approximation in the Ladeveze's terminology) of the unknown field.
Let us consider the advection-diffusion equation with the following initial and boundary conditions where a is the diffusion coefficient and v the velocity field, The aim of the separated representation method is to compute N couples of functions {(X i , T i )} i=1,...,N such that {X i } i=1,...,N and {T i } i=1,...,N are defined respectively in Ω and [0, Tmax] and the solution u of this problem can be written in the separate form The weak form of problem (6) yields: Find u(x, t) verifying the boundary conditions (7) such that for all the functions u ⋆ (x, t) in an appropriate functional space.
We compute now the functions involved in the sum (8).We suppose that the set of functional couples {(X i , T i )} i=1,...,n with 0 ≤ n < N are already known (they have been previously computed) and that at the present iteration we search the enrichment couple (R(t), S(x)) by applying an alternating directions fixed point algorithm that after convergence will constitute the next functional couple (X n+1 , T n+1 ).Hence, at the present iteration, n, we assume the separated representation The weighting function u ⋆ is then assumed as Introducing (10) and ( 11) into (9) it results We apply an alternating directions fixed point algorithm to compute the couple of functions (R, S): -Computing the function S(x).
First, we suppose that R is known, implying that R ⋆ vanishes in (11).Thus, Eq. ( 12) writes where The weak formulation (13) is satisfied for all S ⋆ , therefore we could come back to the associated strong formulation that one could solve by using any appropriate discretization technique.
-Computing the function R(t).¿From the function S(x) just computed, we search R(t).In this case S ⋆ vanishes in (11) and Eq. ( 12) reduces to where all the spatial functions can be integrated in Ω.Thus, by using the following notations equation ( 16) reads As Eq. ( 18) holds for all S ⋆ , we could come back to the strong formulation which is a first order ordinary differential equation that can be solved easily (even for extremely small time steps) from its initial condition.
These two steps must be repeated until convergence, that is, until verifying that both functions reach a fixed point.If we denote by R (q) (t) and R (q−1) (t) the computed functions R(t) at the present and previous iteration respectively, and the same for the space functions: S (q) (x) and S (q−1) (x), the stoping criterion used in this work writes: where 10 −8 represents the square root of the machine precision.We denote by Q n+1 the number of iterations for solving this non-linear problem to determine the enrichment couple of functions X n+1 (x) and T n+1 (t).After reaching convergence we write X n+1 (x) = S(x) and T n+1 (t) = R(t).The enrichment procedure must continue until reaching the convergence of the enrichment global procedure at iteration N , when the separated representation of the unknown field writes: The more usual global stopping criteria are: -For models whose exact solution u ref is known: -For models whose exact solution is not known: with ǫ a small enough parameter (ǫ = 10 −8 in our simulations and the L 2 -norm applies in the whole space-time domain).
Discussion on the space-time separated representations The just proposed strategy needs for the solution of about N × Q space and time problems (with and N the number of functional couples needed to approximate, up to the desired precision, the searched solution).Thus, one must compute N × Q d-D problems, d = 1, 2, 3, whose complexity depends on the spatial mesh considered, and also N ×Q 1D problems (defined in the time interval I) that only need the solution of an ordinary differential equation from its initial condition.Obviously, even for extremely small time steps, the solution of these transient 1D problems does not introduce any difficulty.
If instead the separated representation just discussed, one performs a standard incremental solution, P dD models, d = 1, 2, 3, must be solved (P being the number of time steps, i.e.P = Tmax/∆t, where the time step ∆t must be chosen for verifying the stability conditions).
In all the analyzed cases N and Q are of the order of tens that implies the solution of about hundred three-dimensional problems defined in Ω, instead the thousands (or even millions) needed for solving those models using standard incremental solvers.
A first comparison between both kind of approaches (the one based on the separated representation and the one based on standard incremental strategies) was presented in [6].
3 Recent advances 3.1 Some advances in kinetic theory, quantum chemistry and cell signaling

Solving the Fokker-Planck equation
The first works focusing in the solution of multidimensional models by applying the PGD described in the previous section concerned the modeling of polymeric liquids.These fluids are composed of an unimaginable number of molecules, that are modeled as a series of dumbbells (on which the different forces arising from the flow drag, interactions, solvent molecules bombardment, ... apply) connected by linear or nonlinear springs that represent the stiffness between each two consecutive dumbbells.Let q 1 , • • • , q D be the vectors defining the orientation and extension of each one of these springs.A kinetic theory approach of the molecular nature of such fluids consists of introducing a distribution function Ψ (x, t, q 1 , • • • , q D ) that represents the fraction of molecules that at position x and time t have a conformation given by vectors the distribution function is then defined in a space of dimension 3 + 1 + 3 × D, i.e.
Remark 2 In the case of linear springs Ωq = R 3 , for non-linear springs of finite extension (let b be the maximum spring extension) Ωq reduces to the ball of radius b, B(0, b).Finally, in the case of rigid rods of unit length, Ωq reduces to the surface of the unit sphere S(0, 1). 2 The evolution of the molecular configuration qi , ∀i, is due to both advective and diffusive terms, the first ones arising from the presence of a prescribed flow with nonzero gradient of velocities.The system should verify the momentum balance at any time, i.e. the equilibrium of forces in absence of inertia terms.For the particular expressions of qi the interested reader can refer to [11] that constitutes a huge catalogue of physical systems.
The balance equation for the distribution function (also known as the Fokker-Planck equation) writes: where D Dt denotes the material derivative.The solution of this multidimensional equation by applying the PGD assumes a separated representation of the distribution function The case of steady-state multi-bead-finite-extension-spring models of polymeric liquids in homogeneous flows was addressed in [3].Later, we considered the transient case, in which the time was introduced as an additional coordinate being careful on its upwinding discretization [4].The case of short fiber suspensions in which the beads connector is inextensible was addressed in [40].All these models are linear even when the spring behavior is strongly non-linear (the resulting Fokker-Planck equations (24) are linear with respect to the unknown function Ψ ).
In more complex systems like liquid crystalline polymers or entangled polymeric systems involving double reptation and convective constraint release, as the one proposed in [23], the resulting models are no more linear and a linearization procedure is compulsory.In [34] we considered the simplest one, a fixed point strategy for solving some models associated with entangled polymeric systems.
Finally the case of non-homogeneous flows was deeply analyzed in [40] and [35].The main issue in treating non-homogenous flows using the fully separated representation (25) including physical space, time and conformation spaces lies in the convective stabilization issues, because the stabilization in the conformation space depends on the gradient of velocities that is defined in the physical space.Thus, the adequate separated representation of the stabilization term remains an open issue.
In [17] we addressed the difficulty related to the reduced modeling of systems composed of moving particles within the framework of Brownian dynamics simulations.We provend that classical model reduction techniques based on the use of the proper orthogonal decomposition fails, and that POD or PGD based discretizations perfectly run as soon as the models are rewritten as Brownian configuration fields.Thus, the issue related to updated Lagrangian simulations (the reference configuration is changing in time) -of capital importance in material forming simulations, finite transformations thermomechanic or molecular dynamics -moved to the first plane.We come back to this challenging issue later.

Solving the Schrödinger equation
After treating kinetic theory models we were interested in models arising from quantum chemistry.These models are redoubtable because the equations that govern the electronic distribution, the Schrödinger equation or its fully relativistic counterpart -the Dirac's equation-are defined in spaces whose dimensionality scales with the number of elementary particles involved in the quantum system.The first one, within the Born-Oppenheimer hypothesis (the wavefunction depends perimetrically on the nuclei positions), writes: where i = √ −1, is the Planck's constant divided by 2π, me is the electron mass, De and Dn are the number of electrons and nuclei in the system respectively, is the wavefunction, being xe the space in which the electron e lives (xe ∈ R 3 ), and Xn the position of nuclei n.The differential operator ∇e refers to the gradient with respect to the coordinates xe.Finally V ee ′ and Ven represent the electron-electron and electron-nuclei Coulomb potentials respectively.
The Schrödinger and Dirac equations are defined in a space that scales with the number of elementary particles involved in the system, each one living in R 3 .Thus, when one considers few particles D, the time-independent Schrödinger equation is defined in a space of dimension 3 × D. However, we have proven in [19] that in many systems the real issue is not the model multi-dimensionality, that could be circumvented by using the PGD, but the Pauli's exclusion principle that enforces the anti-symmetry of the wavefunction.For enforcing this constraint, the use of Slater's determinants is the simplest alternative, but the number of terms that these determinants involve explodes as the number of concerned particles grows.The interested reader can refer to [19] [5] for more details concerning the treatment of quantum systems by using the PGD.

Solving the chemical master equation
We also considered recently in [7] the treatment of the chemical master equation, modeling, among many other physics, the cell signaling.In chemical systems involving some (or several) species each one composed by a few number of individuals, the state of the system is given by a probability (the deterministic approach fails when the number of individuals is too reduced for introducing the concept of chemical concentration): and #s i refers to the number of individuals of each species s i .Obviously #s i is changing in time because the different chemical reactions that are producing or destroying individuals in the population s i .The chemical master equation writes: where a j are the propensities (the chemical reaction j operates growing the population z from z − v j : z − v j → z, and simultaneously population z decreases originating individuals of other populations).
Despite the fact that this equation does not imply partial derivatives other that the temporal one, we assumed a separated representation and we proceed in a standard manner as described in [7].

Further considerations on non-linear models
For illustrating the treatment of non linear models we are treating the following simple nonlinear parabolic problem, that we considered in [6] 8 > > < > > : where Ω ⊂ R d , d ≥ 1, Tmax > 0 and a > 0 is the diffusion coefficient.To build-up the approximated solution of ( 29) by using a separated representation, we propose two alternatives -An incremental linearization -A Newton linearization which we describe in sections below.

Incremental linearization
We look for writing the solution of problem (29) in the separated form We suppose that at iteration n, with n < N , the n first modes (X i , T i ), i = 1, • • • n, are already known and that at present iteration we search the new enrichment functional product R(t) • S(x) such that the updated approximation writes The weak form of problem (29) writes The alternating directions scheme proceed by calculating S(x) from the temporal function R(t) just computed, and then, updating R(t) from the just computed S(x), as we described in section 2. The iteration procedure should continue until reaching convergence.Here, the novelty is the treatment of the non-linear term u 2 .The simplest possibility consists in computing this term at the previous iteration, that is, assuming at the present iteration the following approximation of the non-linear term where the notation introduced in section 2 is used again, and where Φ(x) and Γ (x) are given by The associated strong form writes ¿From this solution S(x), we can update the temporal function R(t), by solving its associated strong form where

Newton linearization
¿From now on we denote by u n the solution computed at iteration n, i.e.
Now, after linearization, the solution at the next iteration can be written as u where e u is the solution of the problem whose weak formulation writes: To compute both functions R(t) and S(x) we apply again the alternating directions method deeply described in the previous sections.

Discussion
Both procedures converge but no significant differences in the number of required iterations were noticed.The convergence rate and the computing time were similar.
In the case of linear models, if we solve the problem and then apply the POD (for a given precision) we obtain When the numerical solution is computed using the PGD described in section 2, the solution obtained using N sums u P GD,N is very close to u P OD,N even if the space functions X P OD i and X i , and the time functions T P OD i and T i are very different.In the case of non-linear models the situation is radically different.Even when the exact solution can be represented by a single functional product, i.e.
the non linear solver produces a solution composed of many sums with N > 1.The main reason is that the number of sums is in this case subsidiary of the convergence rate of the non-linear solver.
In [42] we analyzed many other linearization schemes.When we considered the improved fixed point, in which the non-linear term is approximated at iteration q of the enrichment step (the one for computing the couple R(t) andS(x)) according to: then we proved, in the case described above, that the solver converges after computing the first functional couple.In that sense the solver is optimal (only one couple was required for representing the solution that contains a single functional couple) but the computing time is similar to the one required by using the standard fixed point or the Newton strategy previously described.

PGD tensor form
The procedure described in section 2 can be generalized by using a tensor notation [8] [42].
We assume that the discrete problem, that we writes formally as U * T AU = U * T B, can be written in a separated form: where A i , B i and u i involve only the coordinate x i .The separated representation of A and B comes directly from the differential operators involved in the PDE weak form.
At iteration n, vectors u i j , ∀i ≤ n and ∀j ≤ D are assumed known.Now we are looking for an enrichment: where R i , i = 1, • • • , D, are the unknown enrichment vectors.We assume the following form of the test field: Introducing the enriched approximation into the weak form, the following discrete form results: For alleviating the notation we define: where This sum only contains known fields.Thus Eq. ( 49) can be written as: This problem is strongly non linear.To solve it, a method of alternated directions can be applied.The idea is, starting with the trial vectors R (0) update them using an appropriate strategy.The simplest alternatives consist of: The last strategy converges faster but the advantage of the first one is the possibility of updating each vector simultaneously making use of a parallel computing platform.The fixed point of this iteration algorithm allows defining the enrichment vectors When we look for vector R k assuming known all the others R i , i = k, the test field reduces to: The resulting discrete weak form writes: Making use of the arbitrariness of R * K the following linear system can be easily obtained: which can be easily solved.

Decomposition optimality: residual minimization
The alternating direction strategy that we used in section 2 as well as in section 3.3 converges very fast in the case of symmetric differential operators.However, when we include strong non-linearities, advection terms, parameters as additional coordinates, ... the alternating directions strategy fails (the computed functional couples do not improve noticeably the solution).As Pierre Ladeveze proposed many years ago in the framework of the LATIN method, in that case a more efficient strategy consists of minimizing the residual [42]: We denote by ., .the scalar product and by .its associated norm.Using this notation the residual norm writes: The minimization problem with respect to R k reads: The stoping criterion is defined from the residual norm: Despite the fact that this strategy is much more efficient than the alternating directions strategy for non-symmetric operators, it is far to be optimal.Thus, in some particular cases, even if the exact solution accepts a separated representation, the PGD coupled with the residual minimization produces a solution that contains much more functional products that the exact one.
We are at present exploring other strategies, based on concepts of information theory (maximum entropy, ...).Other possibility consists in looking at each iteration for some functional products (instead of a single one at each iteration).This strategy improves the convergence rate reducing the number of sums in the computed decomposition, but at present no strategy allows for a general optimal decomposition.
The reason of this lack of optimality was understood by performing a numerical analysis (see [21] in the present issue).However, when we tried to enforce the enrichment maximizing the residual diminution the numerical procedure exhibits locking.
Previously we identified two major issues of the PGD, we cited the convective stabilization and the treatment of models defined in moving domains, now we should add a third issue that concerns the decomposition optimality.

Non homogeneous essential boundary conditions and complex domains
In [25] we addressed another different topic, the one related to the enforcement of non-homogeneous boundary conditions and the treatment of complex geometries.

Non homogeneous essential boundary conditions
The enforcement of homogeneous boundary conditions is quite simple, it suffices to enforce at each enrichment stage that the first and last components of vectors R i , i = 1, • • • , D, vanish.However, consider the 3D problem: where x = (x, y, z) ∈ Ω = (0, L) 3 , with In [25] we proposed to determine a function ψ, regular enough (continuous in Ω and ∆ψ ∈ L 2 (Ω)) such that ψ verifies the essential boundary condition, i.e. ψ(x) = ug, x ∈ Γ .This function can be constructed by using the transfinite interpolation method [45].Now, by applying the change of variable u = ψ + Υ , the problem (60) writes: where Υ (x) = 0, x ∈ Γ .Thus, the solution of problem (62) reduces to the ones that we considered previously subjected to homogenous essential boundary conditions, where the unknown field is searched in the separated form: The solution procedure simplifies if function ψ is given in a separated form, i.e.
This decomposition can be efficiently performed by invoking the SVD or its multidimensional counterpart.

Complex domains
In [25] we also addressed the question of complex domains.Until now, the whole domain Ω was considered as the direct product of the domains in which each coordinate is defined, Due to this difficulty when we considered in section 3.1.1models defined in the physical (x) and conformation spaces (q 1 , • • • q D ) we considered the decomposition [3] where x ∈ Ω ⊂ R 3 .Thus, this decomposition is not affected at all by the complexity of Ω.However, a full decomposition of the physical space would result in: that can be used without major difficulties if Ω is a parallelepiped, i.e.Ω = Ωx × Ωy × Ωz.However, when it is not the case, a special treatment is required.We considered again the problem (60) now with homogeneous boundary conditions on ∂Ω ≡ Γ , with Ω quite complex in the sense discussed above.
Thus, we must solve the problem with the unknown field u extended to the domain ω by using u = Φ • Υ .Function Φ is regular enough, defined in the extended domain ω and it vanishes on Γ ≡ ∂Ω in order to enforce the boundary conditions.In [25] we proposed the use of the R-functions [45] for computing function Φ than should be rewritten again in a separated form by invoking the SVD.Function Υ is searched again in the separated form

Localization, interfaces and FE coupling
In this section we are modifying the technique just presented in order to capture localized behaviors such as high gradients or weak and strong discontinuities.This kind of behaviors induces too many functional products within a standard separated representation.However, it is well known that they can be easily taken into account within the finite element framework, by using adaptive local mesh refinement in the first case or any of the techniques able to represent weak and strong discontinuities (meshes compatible with the interfaces or enriched approximations when the meshes are not compatible with the interfaces).
Thus, the main goal is how combining separated representations, able to approximate the solution in the most part of the domain where the model is defined, reducing significantly the computational cost, and a localized finite element description, that even if it is expensive in nature, it only applies in a small region of the whole domain where the solution is expected to exhibit localized behavior.To this end, we proposed in [9] a multi-scale approach that somewhat resembles the s-version of the finite element method by J. Fish [24] or the multi-scale FEM proposed by Rank [43], but in this case in a global (Ritz) and separated basis function setting.
For the sake of simplicity we consider a generic 2D problem where K and L are two differential operators.We are assuming that the first one only involves derivatives with respect to the x-coordinate, the second one involving the derivatives with respect to the other coordinate.This model is defined in Ω = Ωx × Ωy and we imagine that the solution u(x, y) is smooth enough everywhere except in a small region Ω l ⊂ Ω.Now, we could imagine an approximation combining both, the separated representation and a finite element approximation, where the last one only applies in Ω l .We assume a mesh on Ω l composed of N l nodes.We also assume that the finite element shape functions related to the finite element approximation cannot be expressed from the tensor product of the one dimensional bases employed to build-up the separated representation to avoid rank deficiency.
Now, the approximation in the whole domain can be written as: where u SR (x, y) is defined in the whole domain Ω = Ωx ×Ωy whereas the finite element enrichment u F E (x, y) is only defined within Ω l .In order to ensure the continuity of the resulting approximation we must enforce the nullity of the enrichment u F E (x, y) on the boundary of Ω l , ∂Ω l .The resulting approximation writes: where N j (x, y) are the standard finite element shape functions and u j the associated weights.Because of the contribution of u SR (x, y) in Ω l , u j does not correspond to the values of the unknown field u(x, y) at the nodal positions (x j , y j ), .
Remark 3 We are assuming that the enrichment is performed by using standard finite elements, but in fact any kind of compact support approximation could be used.In order to represent interfaces involving weak discontinuities one could proceed within the standard finite element method by ensuring that the interface coincides with the elements edges.If the interface passes across the elements an enriched version of the finite element method should be considered (e.g. the extended finite element [49]). 2 In [9] we illustrated the application of this enriched separated approximation when model solutions exhibit localization or for gluing subdomains in which different separated representations are defined, ...

On the coupling of local and global models: the local problem globalization
We are considering again the model given by Eq. ( 6) in absence of advection, i.e. v = 0, and in a one-dimensional physical space Ω: with the following initial and boundary conditions We are assuming that the source term depends on the local value of r fields where the time evolution of the r fields C i (x, t) is governed by r simultaneous ordinary differential equations (the so-called kinetic model).For the sake of simplicity we consider the linear case, the non-linear one reduces to a sequence of linear problems by applying an appropriate linearization strategy [6].The system of linear ODEs writes at each point x ∈ Ω: We are assuming that the kinetic coefficients α ij evolve smoothly in Ω, because in practical applications these coefficients depend on the solution of the diffusion problem, u(x, t).For the sake of simplicity and without loss of generality, from now on we assume those coefficients constant (they were assumed evolving linearly in the physical space in [20]).Now, we are describing three possible coupled solution of Eqs. ( 72) and (75).
1.The simplest strategy consists in using a separated representation of the global problem solution (72) whereas the local problems are integrated in the whole time interval at each nodal position (or integration point).Obviously, this strategy implies the solution of r ordinary differential equations at each node (or integration point).Moreover, the resulting fields C i (x, t), r = 1, • • • , r, don't have a separated structure, and by this reason before of injecting these fields into the global problem (72) we should separate them by invoking, for example, the singular value decomposition (SVD) leading to: As soon as the source term has a separated structure, the procedure illustrated in previous sections can be applied again for computing the new trial solution of the global problem that writes: Thus, this coupling strategy requires the solution of many local problems (for all species) and at all nodal positions (or integration points).Moreover, after these solutions (that we recall that could be performed in parallel) a singular value decomposition must be applied in order to separate these solutions prior to inject them in the PGD solver of the global problem (72).2. The second coupling strategy lies in globalizing the solution of the local problems.
Thus, we assume that the field related to each species can be written in a separated form: and now, we apply the procedure described in section 2 to build-up the reduced separated approximation, i.e. for constructing all the functions involved in (78).Thus, instead of solving the r ODEs in Eq. (76) (that define r one-dimensional problems) at each nodal position (or integration point), we should solve only r higher dimensional coupled models defined in the physical space and time.Obviously, if the number of nodes (or integration points) is important (mainly when 3D physical spaces are considered) the present coupling strategy could offer significant CPU time savings.This strategy allows computing directly a separated representation, and then, with respect to the previous one, the application of the SVD is avoided.However, if the number of species is high, the computational efforts can become important, because the space-time separated solver must be apply to each species.3. The third alternative, that in our opinion is the more appealing one for solving models involving many species, as large as one wants, implies the definition of a new variable C(x, t, c), that as we can notice contains an extra coordinate c, with discrete nature, and that takes integer values: . Thus, we have increased the dimensionality of the problem, but now, only a single problem should be solved, instead of one for each species as was the case when using the previous strategy.This increase of the model dimensionality is not dramatic because as argued in the first section of this work, the separated representation allows circumventing the curse of dimensionality, allowing for fast and accurate solutions of highly multidimensional models.Now, the issue is the derivation of the governing equation for this new variable C(x, t, c) and the separated representation constructor able to define the approximation: As this strategy was retained in our simulations in [20] we are focusing in its associated computational aspects in the next section.

Fully globalized local models
The third strategy just referred implies the solution of a single multidimensional model involving the field C(x, t, c).This original introduction deserves some additional comments.The first one concerns the discrete nature of the kinetic equations Now, by introducing C(x, t, c), such that C(x, t, i) ≡ C i (x, t), the kinetic equations could be written as: where Lc is an operator in the c-coordinate.
If for one second we try to discretize Eq. ( 81) by finite differences, we could write at each node (x k , tp, i): where represents the discrete form of the c-operator at point i.Now, we come back to the separated representation constructor for defining the approximation: For defining such approximation one should repeat the procedure deeply illustrated in section 2. As the operator here involved is less standard we are summarizing the main steps.
We assume that the first n iterations allowed computing the first n sums of Eq. ( 84) and now, we look for the enrichment R(x) Obviously, due to the discrete character of the third coordinate, an integration quadrature consisting of r points, c 1 = 1, • • • , cr = r will be considered later.Now, for computing the three enrichment functions we are considering again (as in section 2) an alternating directions strategy, that proceeds in three steps (that are repeated until reaching convergence): 1. Assuming functions S(t) and W (c) known, the trial function C * (x, t, c) writes Thus the weak form (87) reads: where S ′ = dS dt and (T C q ) ′ = dT C q dt .Now, time integrals and the ones involving the c-coordinate can be performed.The ones involving the coordinate c write: where as just mentioned c i = i, ∀i, and similar expressions can be derived for the integrals involved in the right hand member.Thus, finally it results: where the coefficient ξ x contains all the integrals in the time and c-coordinates related to the left hand member of Eq. ( 88) and F x (x) all the integrals appearing in the right hand member.The strong form related to Eq. (91) writes: whose algebraic nature describes the fact that the kinetic model is local and then it does not involve space derivatives.2. Assuming functions R(x) and W (c) known, the trial function C * (x, t, c) writes Thus the weak form (87) reads: Now, integrals defined in the physical space Ω must be computed, but this task does not involve additional difficulties.Finally it results: where coefficients ξ t and υ t contain all the integrals in the space and the ccoordinate related to the left hand member of Eq. ( 93) and F t (t) the associated integrals appearing in the right hand member.The strong form related to Eq. ( 94) writes: whose first order differential nature results from the first order time derivatives involved by the kinetic model.3. Assuming functions R(x) and S(t) known, the trial function C * (x, t, c) writes Thus the weak form (87) reads: After performing integration in space and time, it results: where coefficients ξ c and υ c contain all the integrals in the space and time related to the left hand member of Eq. ( 96) and F c (c) the associated integrals appearing in the right hand member.The strong form related to Eq. (97) writes: that results in the algebraic system: Again, its algebraic nature comes from the nature of the kinetic model.

Remark 4
For the sake of simplicity we illustrated the solution procedure within the alternating directions framework, but the solutions reported in [20] were computed by enforcing the residual minimization described in section 3.4.2

Error estimation
In this section we are summarizing an error estimator procedure proposed in [8] and based on the use of primal and dual formulations.
For this purpose we assume a generic multidimensional model whose weak form writes a(u, u * ) = b(u * ) (100) x d not necessarily one-dimensional.From now on, the form (100) will be referred as primal form.
The discrete counterpart of (100) reads where and The discrete separated representation of u at iteration nu writes: where u j i is the discrete (nodal) form of u j i (x i ) Now, we are interested in a certain function of u, o(u), of physical interest (the model output).In what follows we assume that the operator defining the output is linear.Thus, we could write If this operator accepts a separated representation, that is: then, the output can be evaluated from: Now, the error on the output, can be evaluated by solving the so-called dual problem (see [8] and the references therein): whose discrete form writes: where O was already defined and where A T is given by A good error estimation needs for an accurate solution of the dual problem.Within the separated representation framework the solution of the dual problem (109) (assumed accurate enough) can be written as: Now, the error in the output can be evaluated from: whose discrete counterpart writes: that is easily computed from: 3.9 Parametric models The case of parametric models was addressed in [42].In this section we revisit the main ideas in the treatment of such models.For this purpose we consider again the heat equation but in which the thermal diffusivity is unknown.The most usual strategy in that case lies in solving that equation for many values of the thermal diffusivity (Monte Carlo procedure).If the statistical distribution of the realizations of the thermal diffusivity represents accurately the real thermal diffusivity distribution, then from the computed temperature fields we can infer different outputs (in a statistical sense).However, this problem can also be solved in a fully deterministic way if the unknown parameter is considered as a new coordinate (as the physical ones -space and time -).Thus, the dimensionality of the problem increases, but this fact is not a serious handicap for the PGD.
We consider again the one-dimensional heat equation where the thermal diffusivity is assumed depending on the temperature field, i.e. k(u): where coefficients a and b are assumed unknown or badly known.
If we introduce the diffusivity expression into the heat equation ( 115) it results: The goal in the solution of that equation is the calculation of the temperature at each point and time, and for any value of the parameters a and b within theirs domains of variability, i.e. u(t, x, a, b): In [42] we considered the approximations of the different functions T i (t), X i (x), A i (a) and B i (b) performed by using standard one dimensional piecewise linear finite element shape functions on a uniform mesh consisting of 500 nodes in each 1D-domain Ω t , Ωx, Ωa and Ω b .If this problem is solved using a mesh-based strategy in the whole domain the complexity scales with 500 4 .However, the Proper Generalized Decomposition needed less than one minute for solving this problem using a personal laptop.
There is no restriction on the number of parameters that can be transformed into additional model coordinates.In [42] we considered the thermal diffusivity, the source term and the initial condition as additional coordinates.

Remark 5
The introduction of model parameters as extra-coordinates has not an impact in the computational efforts, because as the departure model (the heat transfer equation in the present case) does not involve derivatives with respect to those parameters, when we computes the functions related to these coordinates (associated to the model parameters) only algebraic problems must be solved.Moreover, the use of the residual minimization described in section 3.4 increases significantly the convergence rate of the PGD. 2

Solvers coupling and new integration procedures
As the PGD operates by solving the different functions involved in the separated representation independently, one could imagine that different solvers could be applied for solving the problems in the different coordinates.
A direct consequence of this fact was the definition of new integration procedures.In [28] we considered the integration of transient models by applying the Boundary Element Method (BEM).In what follows we revisit the main ideas on this non-incremental boundary element strategy.

Non-incremental boundary element discretizations
The Boundary Element Method allows for an efficient solution of partial differential equations whose kernel functions are known.The heat equation is one of these candidates.When the model involves large physical domains and time simulation intervals the amount of information that must be stored increases significantly.
We proposed in [28] an alternative radically different that leads to a separated solution of the space and time problems within a non-incremental integration strategy.The technique is based on the use of a space-time separated representation of the unknown field that introduced in the residual weighting formulation allows to define a separated solution of the resulting weak form.The spatial step can be then treated by invoking the standard BEM for solving the resulting steady state problem defined in the physical space.Then, the time problem that results in an ordinary first order differential equation is solved using any standard appropriate integration technique (e.g.backward finite differences).
In principle this procedure opens new possibilities for integrating transient models (linear or non-linear) whose space-time fundamental solution is not known but whose steady state fundamental solution is known.
We are illustrating the strategy for constructing these functional products in the case of an academic transient problem, the transient linear heat equation: with the following initial and boundary conditions where a is the diffusion coefficient.The weak formulation yields: Find u(x, t) verifying the boundary conditions (7) such that for all the functions u ⋆ (x, t) in an appropriate functional space.
We compute now the functions involved in the sum (8).We suppose that the set of functional couples {(X i , T i )} i=1,...,n with 0 ≤ n < N are already known (they have been previously computed) and that at the present iteration we search the enrichment couple (R(t), S(x)) by applying an alternating directions fixed point algorithm that after convergence will constitute the next functional couple (X n+1 , T n+1 ).Hence, at the present iteration, n, we assume the separated representation The weighting function u ⋆ is then assumed as Introducing ( 10) and ( 11) into ( 9) it results We apply an alternating directions fixed point algorithm to compute the couple of functions (R, S).Following the procedure described in section 2, functions S(x) and R(t) satisfy: and We can notice that Eq. ( 125) defines a steady-state elliptic equation with constant coefficients.Being the Green solution associated to that equation known, one could apply the BEM for solving it.The only issue in applying the BEM for solving Eq. ( 125) is the integration of the right hand member that implies volumetric integrations that need the definition of appropriate approximation everywhere in the domain.In [28] the approximation was performed by using a moving least square technique.The main advantage in using this strategy is the possibility of solving models for which the Green solution of the space model is known even when it is not the case for the space-time model.
The interested reader can refer to [28] for some numerical experiments and convergence analysis on this non-incremental boundary element technique.

Miscellaneous: High resolution solvers, multi-physics and mixed formulations
Separated representations involving ordinary models (transient or steady state; 2D or 3D) allow in several cases impressive computing time savings.In the transient case is evident as discussed previously (see the discussion at the end of section 2).In the case of rectangular domains is also possible to separate the different space coordinates x, y and z (in the case of more complex domains the separated representation can also be applied by using appropriate strategies, as the one proposed in [25] and here summarized in section 3.5).
If a fully separated representation can be applied, only 1D problems have to be solved.These boundary-value problems ca be solved efficiently even for extremely large number of nodes distributed in the 1D interval, if the boundary value problem is transformed into a initial value problem by using the technique proposed in [15] [14] [16] that avoids the solution of any linear system.The combination of the PGD and these advanced 1D boundary value solvers allows for impressive CPU time savings of several orders of magnitude.
In [18] we considered linear computational homogenization of heterogeneous materials.A simple thermal model must be solved in a quite simple domain (a cube) but in order to capture all the microstructure details a extremely fine mesh is needed (high resolution).The use of a fully separated representation allows considering millions of nodes in each space direction, solving models, at present, beyond the finite elements expectations.
The same kind of techniques were applied for solving multi-physics models arising in the composites manufacturing processes, where curing kinetics are coupled with the non-linear thermal and thermo-mechanical behaviors [41].The simplest coupling lies in the use of a fixed point strategy, but the LATIN framework or the monolithic approaches could be more efficient alternatives.
Recently we have also solved in a fully separated space description the incompressible infinitesimal strain elasticity equations (that are the same the ones that govern the flow of a Newtonian incompressible fluid -the Stokes equations-).The treatment of these mixed formulations needs to satisfy the LBB conditions.At present there are no mathematical results on this topic, but in our numerical experiments we noticed that the one-dimensional meshes used for approximating the displacements field must be richer than the ones used for approximating the pressure field (Lagrange multiplier associated with the incompressibility constraint) and moreover, we must enrich the pressure separated representation one time for, at least, two displacement enrichments.The understanding of these observation represents a work in progress.

Future developments
The separated representation methodology allows in many cases circumventing the redoubtable curse of dimensionality, or at least alleviating its impact.Thus, the PGD seems to be an appealing and powerful strategy for solving models defined in high dimensional spaces including physical (space and time) and configurational (also known as conformational) coordinates.Thus, the solution of many physical models suffering the curse of dimensionality illness is now possible, opening new possibilities in quantum chemistry, the kinetic theory description of materials within the statistical mechanics framework, the solution of the chemical master equation that is present in many branches of biology (genetics, cell-signaling, ...), several models encountered in financial mathematics, ...However, there are numerous other possibilities based in transforming usual models (defined in moderate dimensions, generally space and time) into higher dimensional models that include all the desired model parameters as new model coordinates defined in theirs respective intervals.Thus, if the increase in the model dimensionality is no more a serious drawback (and it is the case if the PGD is applied) one could solve a parametric equation only one time and then particularize the solution for different choices of the parameters.This procedure only implies post-processing, the computation is performed off-line and only one time !It is easy to infer the impact that such procedure could constitute in the context of optimization or inverse identification.There are many other possibilities, it suffices to change our mind, our way of modeling the physics, without trying to enforce the models of physics to be three-dimensional: the models of a physics taking place in the space and time could be multi-dimensional (the mathematics accepts it!).
In this section we are addressing some comments in different potential and exciting application of the Proper Generalized Decomposition.

-Optimization and inverse identification
Imagine that we are interested in optimizing a thermal process by choosing the optimal thermal diffusivity k (that in what follows we consider constant), or that we are interested in identifying this thermal diffusivity from the experimental data recorded by a thermocouple placed at a certain position.In both cases, when one uses any standard strategy, the solution of many direct problems is needed, as well as a minimization strategy able to search the optimum value with respect to a certain cost function.
Our approach consists in assuming the thermal diffusivity as a new coordinate of the model, as described in section 3.9.By applying the PGD we can compute the general solution that represents the value at each position x, time t and for each value of the thermal diffusivity k.
As soon as a trial thermal diffusivity k trial is generated by the minimization algorithm, the temperature field becomes defined at each point x ∈ Ωx and time that allows an easy and fast calculation of the cost function and its gradient.If we denote the cost function C(u), its gradient reads: that can be easily determined from the cost function expression and from the derivative of Eq. ( 128): Moreover, for simple cost functions one could expect to find directly the extremum points, that is, the points at which the derivative with respect to the design variables (in the present case the thermal diffusivity) vanish: Let K be the set of such extremum points (K = {k 1 , • • • , k M }).Now, it suffices calculating the cost functions for each value of the thermal diffusivity k ∈ K: C (that could be identified to the domain at the initial time, i.e.Ω ref x = Ωx(t = 0)).Now, if the kinematics is known, the function given the position at time t of a point that occupied the position X at time t = 0 can be easily derived: x = x(X, t).Now, the partial differential equation defined in Ωx(t) × Ω t can be transformed in its counterpart now defined in Ω ref x × Ω t .This procedure is quite standard in computational mechanics when one proceeds in the Lagrangian framework.Now, with the problem defined in (X, t) ∈ Ω ref x × Ω t , the natural separated representation writes: Because the global weak space-time formulation is first written, and then the functions of space and time searched after integration in the time and space domains respectively (we don't enforce the balance equation at any particular time step as is the case when using incremental strategies).

-Towards real time simulations
Many problems in engineering needs for very fast solutions, sometimes in the real time range.As the models in engineering are quite complex (non-linear, involving thousands or even millions degrees of freedom, ...) the real time constraint is an intractable issue in numerical simulation by using the standard procedures, despite the impressive recent progresses in the computational resources.Parallel platforms, high performance computing ... do not suffice, at least at present, for reaching the real time simulation requirements.
One area in which we are specially interested is the one that concerns surgery simulators.In these applications, real time is not a caprice, real time is needed to be a valuable surgery tool.Living tissues are non-linear, anisotropic, geometrically complex, and the applied loads move on their surfaces ... all them making difficult fast simulations.
In [36] [37] and [38] we explored the use of model reduction techniques based on the use of Proper Orthogonal Decompositions (POD) combined with an advanced non-linear solver (the asymptotic numerical method) that avoids the necessity of recomputing the tangent matrix at each load updating.Despite these advanced computational ingredients, real time requirements were reached with some difficulties.
At present we are considering the application of the PGD.We could solve the non-linear elasticity problem (involving material and geometrical non-linearities) including the moving load applying in the organ surface as a new model coordinate.Let p the vector representing the applied force on the organ surface.The point at which that force applies is designated by y.As we are considering quasi-static behaviors we don't need to specify the dependence y(t).Now, the non-linear elasticity problem is solved assuming the separated representation of the displacement field given by: that represents the displacement u at point x when a load p applies at position y.Now, as soon as the force and its location are given, p g and y g , the displacement field is calculated The construction of this general solution is made once, off line, and then used in-line as a simple post-processing allowing to fulfill real time constraints.

-Homogenization
Many processes and/or materials exhibit many scales.Microstructure in materials induces different space scales.Processes can also induce different time scales.Thus, for example when one considers processes like ultrasonic welding [31] the scale related to the load evolution is of some microseconds, whereas the characteristic time related to the welding process itself is of some seconds.When these scales are well separated (turbulence is a nice couner-example!) we could consider separated (and independent) coordinates related to the different scales.For the sake of simplicity we consider two time scales.The fastest one implies a time τ being t the one related to the slowest one.The characteristic times of both scales define the scale factor ǫ: t = ǫτ .Thus, the time dependence of a generic field u(t) can be rewritten as u( t, τ ) and then that taking into account that t = t reduced to: and so on for the higher order derivatives.Thus, the partial differential equations involved in the thermo-mechanical model of the process can be rewritten making use of the two new time coordinates.In general, for solving the resulting model an asymptotic expansion of the different fields is performed, and then by identifying the equation at the different orders the model is solved and the effects of the microscopic scale appears naturally in the macroscopic one.However, in some situations this procedure becomes quite technical.
The use of the PGD avoids the asymptotic expansion because the multidimensional problem that results by introducing the different independent time coordinates can be efficiently solved without major difficulties by applying the PGD.Thus, the solution is searched under the separated form: In general this procedure can be extended to models involving more than two time scales, for the ones involving many space scales (rich microstructures with many different characteristic lengths) and for the ones combining many space and time scales.

-Convective stabilization
It is well-known that standard finite element (Galerkin) methods do not work well for convection-diffusion or convection-diffusion-reaction equations, since they lead to unstable, oscillating, solutions [22].The first stabilization methods including upwinding in the convective term in order to produce an artificial diffusion that eventually lead to stable solutions were published in [26].Among the very numerous methods that have been proposed for the stabilization of convection-diffusion equations, the streamline-upwind/Petrov-Galerkin (SUPG) method [27] is one of the most extended.The major drawback of this method is that it introduces some cross-wind artificial numerical diffusion when applied to two-or three-dimensional equations.The use of PGD could circumvent, or at least alleviate this drawback.It has been noticed that the use of such an approximation leads to a sequence of different onedimensional problems, for which "exact" SUPG stabilizations exist.In this way, we can apply a standard SUPG method to a sequence of one-dimensional problems to obtain a properly stabilized solution to a two-or three-dimensional convectiondiffusion(-reaction) problem.
For simplicity, we consider the steady-state 2D convection-diffusion(-reaction) equation.This equation is given by with n = 2, 3 and with boundary conditions where u is the scalar quantity to be transported and also the unknown field of the problem, v is the advective velocity, a > 0 the diffusivity, assumed constant, σ the reaction term and s(x) a volumetric source term.The function u D denotes the prescribed value of u on the Dirichlet portion of the boundary given by Γ D and t denotes the value of the normal diffusive flux on the Neumann boundary In what follows, for the sake of simplicity, we shall omit the reaction term of the equation, since it does not imply any special difficulty for the technique here developed.
The weak form of the problem given by Eqs.(138)-( 139) is: which is very often expressed compactly with the help of the following notation: giving rise to the following compact form of the equation: The general form of the consistent stabilization techniques is [22] a(w, u) + c(v; w, u) + X where P e is the mesh Péclet number, defined as P e = vh/2a.h represents the mesh size parameter and v the modulus of the convective velocity.The problem with this type of stabilization is that, on one side, this value has been defined for one-dimensional problems and linear finite elements in order to represent exactly the problem solution.
Within the PGD framework, in the enrichment stage we look for an improved representation of the essential field in the form or, equivalently, u n+1 (x, y) = u n (x, y) + R(x) • S(y).
The test function will then be given by where we have omitted, for clarity, the obvious dependence of R in x and S in y.
After applying such an approximation, the weak form of the problem given by Eq. ( 148) is solved by using the fixed-point algorithm (extensively used throughout the present paper) in which the R and S functions are sought iteratively.
For instance, assuming R(x) known, the resulting expression will be After some lengthy, yet simple, algebra and after taking into account that second derivatives of functions R and S vanish if we approximate them by standard, firstorder, finite elements, we obtain for the term A, For the term B we obtain The impact of the choice of τ is being evaluated and PGD solutions compared with fully 2D SUPG stabilized solutions.

-Separated representation of Molecular Dynamics
The molecular dynamic problems are based on the resolution of a set of Newton's equations for a given number of particles denoted here by N .The system is described by a set of discrete positions varying during time x 1 (t), x 2 (t), ..., x N (t).The dynamics of each particle depends on the position of the surrounding particles through an interaction potential as well as on its own position in the case of an existence of an external field generating a force on each particle.The general equations are then written in the next form: m i is the mass and f i is the force applied on the i th particle which is non linearly dependent on the system state.These N equations can be gathered into a compact notation involving a vector x = (x 1 , ..., x N ) T .
where the vectorial discrete operator L contains all forces and masses contributions.
The main difficulty associated to the solution of this equations is associated to the high number of iteration steps for which we have to solve high number of discrete equations related to each particle.In the context of the PGD we can look for the solution directly as a product of time functions T j (t) by some configurations states of the particles given by a position coordinates vector X j .X j could be interpreted as a discrete significant mode of the position vector during the time evolution that has been to be calculated 'a priori'.If we look for a general solution that writes and if we assume known the first n couples of configuration states and time functions (n < R) then the new couple Y, S(t) deriving from the enrichment process must satisfy The residual minimization of section 3.4 could be applied onto this equation.The main difficulty is related to the updating of the discrete operator L that accounts for high non linearity.This requires a complete reconstitution of the time-space solution for each enrichment stage (if we consider the simplest updating of the non linearity described in section 3.2).Even if this technique allows to find the solution with an 'a priori' estimation it costs, unfortunately, R times more expensive than the direct simulation.

-New dreams
The possibility of solving efficiently multidimensional models including hundreds of coordinates [4] opens new possibilities in that concerns the solution of models that are in nature multidimensional (quantum chemistry, statistical mechanics, biology, computational finance, ...) but also many others (standard in nature) that could be transformed into highly multidimensional models.
We could expect that some non-linear problems could be linear in a higher dimensional spaces, making possible its solution in such spaces and then come back to the initial one solving a non-linear algebraic problem.Other possibility that we recently explored lies in the solution of some simple models that exhibit bifurcations in theirs solution.The treatment of these models needs some appropriate technique able to identify the bifurcation points and to follow each solution branch.Our idea was to immerse these physics in higher dimensional spaces, such that the solution can represent simultaneously all the solution richness.Thus, in the buckling of a column in 2D, we obtained a solution that after reaching the critical load exhibit all the possible solutions simultaneously.The price to be paid is obviously the consideration of many extra-coordinates.Thus, if a solver able to address highly multidimensional models is available, new and some times unimaginable possibilities appears suddenly!

Conclusions
In this paper we presented the ability of the Proper Generalized Decomposition -PGD -for solving highly multidimensional models.This technique operates by constructing a separated representation of the solution, such that, the solution complexity scales linearly with the dimension of the space in which the model is defined, instead the exponentially-growing complexity characteristic of mesh based discretization strategies.
The PGD makes possible the efficient solution of models defined in multidimensional spaces, as the ones encountered in quantum chemistry, kinetic theory description of complex fluids, genetics (chemical master equation), financial mathematics, ... but also those, classically defined in the standard space and time, to which we can add new extra-coordinates (parametric models, ...) opening numerous possibilities (optimization, inverse identification, real time simulations, ...).
Despite the first success in the treatment of the models described in this paper, this numerical proposal is too recent to conclude about its limits, ... Thus, the analysis should continue for identifying new opportunities but also its limitations.
)τ R(u)dΩ = (w, s) + (w, t) ΓN = l(w).(143) where P(w) is some operator applied to the test functions and R(u) = L(u) − s is the residual of the equation.In the SUPG method, P(w) = v • ∇w.The exact nodal values for 1D convection-diffusion equation are obtained for a value
1 , • • • C M and select the minimum value, i.e. kop = k j such that C j = min{C 1 , • • • C M }.Thus, we could compute the global minimum, instead of a simple minimum (not necessary the global one) obtained by using standard minimization techniques.-Evolvingdomains Some times we are interested in solving models defined in a domain that is evolving in time, i.e. the partial differential equations are defined in Ωx(t).There are probably many ways to treat this kind of models, but a simple procedure consists in writing the problem in the reference domain Ω ref x