Model reduction by separation of variables: a comparison between Hierarchical Model reduction and Proper Generalized Decomposition

Hierarchical Model reduction and Proper Generalized Decomposition both exploit separation of variables to perform a model reduction. After setting the basics, we exemplify these techniques on some standard elliptic problems to highlight pros and cons of the two procedures, both from a methodological and a numerical viewpoint.


Introduction
This paper is meant as a first attempt to compare two procedures which share the idea of exploiting separation of variables to perform model reduction, albeit with different purposes. Proper Generalized Decomposition (PGD) is essentially employed as a powerful tool to deal with parametric problems in several fields of application [3,14,23]. Parametrized models characterize multi-query contexts, such as parameter optimization, statistical analysis or inverse problems. Here, the computation of the solution for many different parameters demands, in general, a huge computational effort, and this justifies the development of model reduction techniques.
For this purpose, projection-based techniques, such as Proper Orthogonal Decomposition (POD) or Reduced Basis methods, are widely used in the literature [11]. The idea is to project the discrete operators onto a reduced space so that the problem can be solved rapidly in the lower dimensional space. PGD adopts a completely different way to deal with parameters. Here, parameters are considered as new independent variables of the problem, together with the standard space-time ones [5]. Although the dimensionality of the problem is inevitably increased, PGD transforms the computation of the solution for new values of the parameters into a plain evaluation of the reduced solution, with striking computational advantages. Hierarchical-Model (HiMod) reduction has been proposed to improve one-dimensional (1D) partial differential equation (PDE) solvers for problems defined in domains with a geometrically dominant direction, like slabs or pipes [6,17]. The main applicative field of interest is hemodynamics, in particular the modeling of blood flow in patient-specific geometries. Purely 1D hemodynamic models completely drop the transverse dynamics, which, however may be locally important (e.g., in the presence of a stenosis or an aneurism). HiMod aims at providing a numerical tool to incorporate the transverse components of the 3D solution into a conceptually 1D solver. To do this, the driving idea is to discretize main and transverse dynamics in a different way. The latter are generally of secondary importance and can be described by few degrees of freedom using a spectral approximation, in combination, for instance, with a finite element (FE) discretization of the mainstream.
The parametric version of HiMod (namely, HiPOD) is a more recent proposal [4,13]. On the other hand, PGD is not so widely employed in a non-parametric setting, despite its original formulation [12]. Nevertheless, for the sake of comparison, in this paper we consider the non-parametric as well as the parametric versions of both the HiMod and PGD approaches. The goal is to begin a preliminary comparative analysis between the two methodologies, to highlight the respective weaknesses and strengths. The main limit of PGD remains its inability to deal with non-Cartesian geometries without losing the computational benefits arising from the separability of the spatial coordinates. HiMod turns out to be more flexible from a geometric viewpoint. On the other hand, PGD turns out to be extremely effective for parametric problems thanks to the explicit expression of the PGD solution in terms of the parameters, while HiPOD can be classified as a projection-based method with all the associated drawbacks. In perspective, the ultimate goal is to merge HiMod with PGD to emphasize the good features and mitigate the intrinsic limits of the two methods taken alone.

The HiMod approach
Hierarchical Model reduction proved to be an efficient and reliable method to deal with phenomena charaterized by dominant dynamics [10]. In general, the computational domain itself exhibits an intrinsic directionality. We assume where Ω 1D ⊂ R denotes the supporting fiber aligned with the main stream, while γ x ⊂ R d−1 is the transverse fiber at x ∈ Ω 1D , parallel to the transverse dynamics.
For the sake of simplicity, we identify Ω 1D with a straight segment, (x 0 , x 1 ). We refer to [15,18] for the case where Ω 1D is curvilinear. From a computational viewpoint, the idea is to exploit a map, Ψ : Ω →Ω, transforming the physical domain, Ω, into a reference domain,Ω, and to make explicit computations inΩ only. Typically,Ω coincides with a rectangle in 2D, with a cylinder with circular section in 3D. To define Ψ, for each x ∈ Ω 1D , we introduce the map, ψ x : γ x →γ d−1 , from fiber γ x to the reference transverse fiber,γ d−1 , so that the reference domain coincides withΩ = x ∈Ω 1D {x} ×γ d−1 . The supporting fiber is preserved by map Ψ, which modifies the lateral boundaries only.
We consider now the (full) problem to be reduced. Due to the comparative purposes of the paper, we focus on a scalar elliptic equation, and, in particular, on the associated weak formulation, where V ⊆ H 1 (Ω), a(·, ·) : V ×V → R is a continuous and coercive bilinear form and F(·) : V → R is a continuous linear functional. To provide the HiMod formulation for problem (1), we introduce the hierarchical reduced space k=1 denotes a modal basis of functions orthogonal with respect to the L 2 (γ d−1 )-scalar product. Index m sets the hierarchical level of the HiMod space, being V m ⊂ V m+1 , for any m. Concerning V h 1D , we adopt here a standard FE space, although any discrete space can be employed (see, e.g., [18], where an isogeometric discretization is used). Functions in V h 1D have to include the boundary conditions on {x 0 } × γ x 0 and {x 1 } × γ x 1 ; analogously, the modal functions have to take into account the boundary data along the horizontal sides. In Sect. 4 further comments are provided about the selection of the modal basis and of the modal index m. The HiMod formulation for problem (1) thus reads To ensure the well-posedness of formulation (3) and the convergence of the HiMod approximation, u HiMod m , to the full solution, u, we endow the HiMod space with a conformity and a spectral approximability hypothesis, and we introduce a standard density assumption on the discrete space V h 1D (see [17] for all the details). The HiMod solution can be fully characterized by introducing a basis, {θ l } N h l=1 , for the space V h 1D . Actually, each modal coefficient,ũ k , of u HiMod m can be expanded in terms of such a basis, so that, we obtain the modal representation The actual unknowns of problem (3) become the mN h coefficients {ũ k,l } m, N h k=1,l=1 . With reference to the Poisson problem, −∆u = f , completed with full homogeneous Dirichlet boundary data, the corresponding HiMod formulation, after exploiting (4) in (3) and picking v m (x, y) = θ i (x)ϕ j (ψ x (y)) with i = 1, . . . , N h and j = 1, . . . , m, reduces to the system of mN h 1D equations in the mN h unknowns Information associated with the transverse dynamics are lumped in the coefficients {r a,b jk }, so that the HiMod system is solved on the supporting fiber, Ω 1D . Collecting the HiMod unknowns, by mode, in the vector u HiMod we can rewrite the HiMod system in the compact form According to (5), for each modal index j, between 1 and m, the nodal index, i, takes the values 1, . . . , N h . Thus, HiMod reduction leads to solve a system of order mN h , independently of the dimension of the full problem (1).

The PGD approach
To perform PGD, we have to introduce on problem (1) a separability hypothesis with respect to both the spatial variables and the data [5,22]. Thus, domain Ω ⊂ R d coincides with the rectangle Ω x × Ω y if d = 2, with the parallelepiped Ω x × Ω y × Ω z (total separability) or with the cylinder Ω x × Ω y (partial separability) if d = 3, for Ω x , Ω y , Ω z ⊂ R and Ω y ⊂ R 2 , being y = (y, z). In the following, we focus on partial separability, since it is more suited to match HiMod reduction with PGD. Analogously, we assume that the generic problem data, d = d(x, y, z), can be written as d = d x (x)d y (y). The separability is inherited by the PGD space of Ω x and Ω y , respectively. In general, W x h and W y h are FE spaces, although, a priori, any discretization can be adopted. It turns out that W m is a tensor function space, being . Index m plays the same role as in the HiMod reduction, setting the level of detail for the reduced solution (see Sect. 4 for possible criteria to choose m). PGD exploits the hierarchical structure in W m to build the generic function w m ∈ W m . In particular, w m is computed as where w x k and w y k are assumed known for k = 1, . . . , m − 1, so that the enrichment functions, w x m and w y m , become the actual unknowns. To provide the PGD formulation for the Poisson problem considered in Sect. 2, we exploit representation (8) for the PGD approximation, u PGD m , and we pick the test function as X(x)Y (y), with X ∈ W x h and Y ∈ W y h . The coupling between the unknowns, u x m and u y m , leads to a nonlinear problem, which is tackled by means of the Alternating Direction Strategy (ADS) [5]. The idea is to look for u x m and u y m , separately via a fixed point procedure. We introduce an auxiliary index to keep trace of the ADS iterations, so that, at the p-th ADS iteration we compute u where the separability of f is exploited (the dependence on the independent variables, x and y, is omitted to simplify notation). Successively, we compute u y, p m , after setting u x m to u x, p m and choosing function X as to u x, p m in the test function, so that we obtain, for any Y ∈ W y h , The algebraic counterpart of (9) and (10) is obtained by introducing a basis, . Thanks to these expansions and to the arbitrariness of X and Y , we can rewrite (9) and (10) as and h are the stiffness and mass matrices associated with xand y-variables, with (11) and (12) are solved at each ADS iteration, so that the computational effort characterizing PGD is the one associated with the solution of two systems of order N x h and N y h , respectively, for each ADS iteration. When a certain stopping criterion is met (see the next section for more details), ADS procedure yields vectors u x m and u y m which identify the enrichment functions u x m and u y m .

HiMod reduction versus PGD
Both HiMod reduction and PGD exploit the separation of variables and, according to [5], belong to the a priori approaches, since they do not rely on any solution to the problem at hand. Nevertheless, we can easily itemize features which distinguish the two techniques. The most relevant ones concern the geometry of Ω, the selection of the transverse basis and of the modal index, and the numerical implementation of the two procedures. Pros and cons of the two methods are then here highlighted.

Domain geometry
HiMod reduction and PGD advance precise hypotheses on the geometry of the computational domain.
According to the HiMod approach, Ω is expected to coincide with a fiber bundle and to be mapped into the reference domain,Ω, by a sufficiently regular transformation. Actually, map Ψ is assumed differentiable, while map ψ x is required to be a C 1 -diffeomorphism, for all x ∈ Ω 1D [17]. These hypotheses introduce some constraints, in particular, on the lateral boundary of Ω which, e.g., cannot exhibit kinks. Additionally, geometries of interest in many applications, such as bifurcations or, more in general, networks are ruled out from the demands on ψ x and Ψ. An approach based on the domain decomposition technique is currently under investigation as a viable way to deal with such geometries. The isogeometric version of HiMod (i.e., the HIgaMod approach) will play a crucial role in view of HiMod simulations for the blood flow modeling in patient-specific geometries [18].
The constraints introduced by PGD on the geometry of Ω are more restrictive. The separability hypothesis leads to consider essentially only Cartesian domains. This considerably reduces the applicability of PGD to practical contexts. Some techniques are available in the literature to overcome this issue. For instance, in [9] a generic domain is embedded into a Cartesian geometry, while in [7] the authors introduce a parametrization map for quadrilateral domains.
Overall, HiMod reduction exhibits a higher geometric flexibility with respect to PGD, in its straightforward formulation. As discussed in Sect. 5, this limitation can be removed when considering a parametric setting.

Modeling of the transverse dynamics
In the HiMod expansion, y-components, ϕ k (ψ x (y)), are selected before starting the model reduction. This choice, although coherent with an a priori approach, introduces a constraint on the dynamics that can be described, so that hints about the solution trend along the transverse direction can be helpful to select a representative modal basis. In the original proposal of the HiMod procedure, sinusoidal functions are employed according to a Fourier expansion [6,17]. This turns out to be a reasonable choice when Dirichlet boundary conditions are assigned on the lateral surface, Γ lat = {x} × ∂γ x , of Ω. Legendre polynomials, properly modified to include the homogeneous Dirichlet data and orthonormalized, are employed in [17] as an alternative to a trigonometric expansion. Nevertheless, Legendre polynomials require high-order quadrature rules to accurately compute coefficients {r a,b jk }.
In [1], the concept of educated modal basis is introduced to impose generic boundary conditions on Γ lat . The idea is to solve an auxiliary Sturm-Liouville eigenvalue problem on the transverse reference fiberγ d−1 , to build a basis which automatically includes the boundary values on Γ lat . The eigenfunctions of the Sturm-Liouville problem provide the modal basis. A first attempt to generalize the educated-HiMod reduction to three-dimensional (3D) cylindrical geometries is performed in [10], where the Navier-Stokes equations are hierarchically reduced to model the blood flow in pipes. This generalization is far from being straightforward due to the employment of polar coordinates. To overcome this issue, we are currently investigating the HIgaMod approach [18], which allows us to define the transverse basis as the Cartesian product of 1D modal functions, independently of the considered geometry. Additionally, we remark that any modal basis can be precomputed on the transverse reference fiber before performing the HiMod reduction, thanks to the employment of map Ψ. This considerably simplifies computations. When applying PGD, y-components are unknown as the ones associated with x. This leads to the nonlinear problems (9)-(10), thus loosing any advantage related to a precomputation of the HiMod modal basis. On the other hand, PGD does not constrain the transverse dynamic to follow a prescribed (e.g., sinusoidal) analytical shape as HiMod procedure does. The educated-Himod reduction clearly is out of this comparison, since the modal basis strictly depends on the problem at hand.
Finally, we observe that HiMod modes are orthonormal with respect to the L 2 (γ d−1 )-norm. This property is not ensured by PGD.
Concerning the selection of the modal index m in (2) and (7), as a first attempt, both HiMod reduction and PGD resort to a trial-and-error approach, so that the modal index is gradually increased until a check on the accuracy of the reduced solution is satisfied. For instance, in [6,17] a qualitative investigation of the contour plot of the HiMod approximation drives the choice of m. Concerning PGD, the check on the relative enrichment is usually employed, with TOL E a user-defined tolerance [5]. An automatic selection of index m can yield a significant improvement. In [19,21], an adaptive procedure is proposed for HiMod, based on an a posteriori modeling error analysis. In particular, the estimator in [19] is derived in a goal-oriented setting to control a quantity of interest, and exploits the hierarchical structure (i.e., the inclusion V m ⊂ V m+d , ∀m, d ∈ N + ) typical of a HiMod reduction. A similar modeling error analysis is performed in [2] for PGD, although no adaptive algorithm is here set to automatically pick the reduced model. Paper [21] generalizes the a posteriori analysis in [19] to an unsteady setting, providing the tool to automatically select m together with the partition T h along Ω 1D and the time step. Finally, HiMod allows to tune the modal index along the domain Ω, according to the local complexity of the transverse dynamics. In particular, m can be varied in different areas of Ω or, in the presence of very localized dynamics, in correspondence with specific nodes of the partition T h . We refer to these two variants as to piecewise and pointwise HiMod reduction, in contrast to a uniform approach, where the same number of modes is adopted everywhere [16,20]. This flexibility in the choice of m is currently not available for PGD. Adaptive strategies to select the modal index are available for the three variants of the HiMod procedure [19,21].

Computational aspects
From a computational viewpoint, HiMod reduction and PGD lead to completely different procedures. Indeed, for a fixed value of m, we have to solve the only system (6) of order mN h when applying HiMod, in contrast to PGD which demands a multiple solution of systems (11) subdivide Ω 1D into 285 subintervals. We set the PGD discretization along y as well as the PGD and the HiMod index m in order to ensure the same accuracy, TOL, on the reduced approximations with respect to a reference FE solution, computed on a 2500 × 500 structured mesh. In particular, for TOL = 8 · 10 −3 , we have to subdivide interval (0, 1) into 20 uniform subintervals, and to set m to 6 and to 9 in the PGD and the HiMod discretization, respectively. Sinusoidal functions are chosen for the HiMod modal basis. The ADS iterations are controlled in terms of the relative increment, as u with TOL FP = 10 −2 . Fig. 1 shows the reduced approximations (which are fully comparable with the FE one, here omitted). The contourplots are very similar. The coarse PGD y-discretization justifies the slight roughness of the PGD contourlines.
Another distinguishing feature between HiMod and PGD is the domain discretization. Indeed, HiMod requires only the partition T h along Ω 1D , independently of the dimension of Ω. No discretization is needed in the y-direction, although we have to carefully select the quadrature nodes to compute coefficients {r a,b jk }. This task becomes particularly challenging when dealing with polar coordinates [10]. With PGD to benefit of the computational advantages associated with a 1D discretization, we are obliged to assume the full separability of Ω; actually, a partial separability demands a 1D partition for Ω x , and a two-dimensional partition of Ω y . As explained in Sect. 5, non-Cartesian domains require a 3D discretization of Ω.
Finally we analyze the interplay between the enrichment and the ADS iterations in the PGD reduction. We investigate the possible relationship between TOL FP in (14) and TOL E in (13), to verify if a small tolerance for the fixed point iteration improves the accuracy of the PGD approximation, thus reducing the number of enrichment steps. To do this, we adopt the same test case used above. Table 1 gathers the number of ADS iterations, #IT FP , the number, m, of enrichment steps, and the CPU time1 (in seconds) demanded by the PGD procedure, for two different values of TOL E and three different choices of TOL FP . In particular, in column #IT FP we specify the number of ADS iterations required by each enrichment step. As expected, there exists a link between the two tolerances, namely, when a higher accuracy constrains the fixed point iteration, a smaller number of enrichment steps is performed to ensure the accuracy TOL E . Table 1 Quantitative analysis for PGD in terms of fixed point iterations and enrichment steps

HiMod reduction and PGD for parametrized problems
The actual potential of PGD becomes more evident when considering a parametric setting, i.e., when problem (1) is replaced by the formulation with µ a parameter, which may represent any data of the problem, e.g., the coefficients of the considered PDE, the source term, a boundary value or the domain geometry. The technique adopted by PGD to deal with the parametric dependence in (15) is very effective. Parameter µ is considered as an additional independent variable which varies in a domain Ω µ [5]. Thus, the PGD space (7) changes into the new one with W µ h a discretization of the space L 2 (Ω µ ; R Q ), being Q the length of vector µ. Generalizing the enrichment paradigm in (8), at the m-th step of the PGD approach applied to problem (15) we have to compute three unknown functions, u x m , u y m and u µ m , by picking the test function as ∫ Ω µ∇u · ∇v + b · ∇u dΩ with b = [2.5, 0] T and µ the parameter to be varied in Ω µ = [1,5], The problem is completed with mixed boundary conditions, namely a homogeneous Dirichlet data on Γ up ∪ Γ down , the non-homogeneous Dirichlet condition, u = u in with u in = y(1 − y), on Γ in and a homogeneous Neumann value on Γ out = {3} × (0, 1). We apply the PGD reduction for m = 2, and we uniformly subdivide Ω x , Ω y , Ω µ , being N x h = 150, N y h = 50, N µ h = 500. The tolerance in (14) is set to 10 −2 . Fig. 2 compares the PGD approximation for µ = 1 and µ = 2.5 with a reference full solution coinciding with a linear FE approximation computed on a 300 × 100 structured mesh. The qualitative matching between the corresponding solutions is significant. From a quantitative viewpoint, the L 2 (Ω)-norm of the relative error associated with the PGD approximation does not vary significantly by increasing m, whereas a slight error reduction is detected by increasing µ.
The parametric counterpart of the HiMod reduction, known as HiPOD, merges HiMod with POD [4,13]. HiPOD pursues a different goal with respect to PGD. Indeed, for a new value, µ * , of the parameter, PGD provides an approximation for the full solution u(µ * ), while HiPOD approximates the HiMod solution associated with µ * . The offline/online paradigm of POD is followed also by HiPOD. The peculiarity is that the offline step is now performed in the HiMod setting to contain the computational burden typical of this stage and by relying on the good properties of HiMod in terms of reliability-versus-accuracy balance. Thus, we choose P different values, µ = µ i with i = 1, . . . , P, for parameter µ, and we collect the HiMod approximation for the corresponding problem (15) into the response matrix, S = u HiMod m (µ 1 ), u HiMod m (µ 2 ), . . . , u HiMod m (µ P ) ∈ R mN h ×P , according to representation (5). Successively, we define the null-average matrix and we apply the Singular Value Decomposition (SVD) to V, so that V = ΦΣΨ T , where Φ ∈ R (mN h )×(mN h ) and Ψ ∈ R P×P are the unitary matrices of the left-and of the right-singular vectors of V, respectively while Σ = diag (σ 1 , . . . , σ ρ ) ∈ R (mN h )×P denotes the pseudo-diagonal matrix of the singular values of V, being σ 1 ≥ σ 2 ≥ · · · ≥ σ ρ ≥ 0 and ρ = min(mN h , P) [8]. The POD basis is identified by the first l left singular vectors, φ i , of V, so that the reduced POD space is V l POD = span{φ 1 , . . . , φ l }, with dim(V l POD ) = l and l mN h . In the numerical assessment below, value l coincides with the smallest integer such that σ 2 l < ε, with ε a prescribed tolerance. The online phase of HiPOD approximates the HiMod solution to problem (15) for a new value, µ * , of the parameter by exploiting the POD basis instead of solving system (6). This is performed via a projection step. After assembling the HiMod stiffness matrix and right-hand side, A HiMod m (µ * ) and f HiMod m (µ * ), associated with the new value of the parameter, we solve the POD system of order l where A POD (µ * ) = (Φ l POD ) T A HiMod m (µ * ) Φ l POD and f POD (µ * ) = (Φ l POD ) T f HiMod m (µ * ) denote the POD stiffness matrix and right-hand side, respectively with Φ l POD = [φ 1 , . . . , φ l ] ∈ R (mN h )×l the matrix collecting the POD basis vectors. The HiMod solution is thus approximated by vector Φ l POD u POD (µ * ) ∈ R mN h , i.e., after solving a system of order l instead of mN h . Overall, HiPOD requires to solve P linear systems of order mN h during the offline phase, additionally to a system of order l in the online phase.
To check the performances of HiPOD, we adopt the test case used above for PGD, for the same values of the parameters, µ * = 1 and µ * = 2.5. The reference solution is the corresponding HiMod approximation computed by using m = 15 sinusoidal functions in the y-direction, and a linear FE discretization along the mainstream based on a uniform subdivision of Ω 1D into 50 subintervals. The same HiMod discretization is adopted to build the response matrix. Concerning the HiPOD approximation, we pick P = 100 by uniformly sampling the interval [1,5], and we select ε = 2.5 · 10 −15 . This choice sets the dimension of the POD space to l = 8, so that we have to solve a system of order 8 instead of 750. The contour plots in Fig. 3 qualitatively compare the HiMod solution with the HiPOD approximation for l = 1. The correspondence between the two approximations is good despite a single POD mode is employed (in l = 1 l = 4 µ * = 1 4.06e-02 2.53e-06 µ * = 2.5 2.74e-03 1.11e-07 random 7.19e-03 2.65e-07 l = 6 l = 8 µ * = 1 1.79e-09 5.58e-12 µ * = 2.5 4.05e-10 1.21e-13 random 2.97e-10 3.66e-13 Fig. 3 Contour plots: comparison between the reference HiMod solution (left) and the HiPOD approximation with l = 1 (right), for µ * = 1 (top) and µ * = 2.5 (bottom). Table: relative error between HiMod and HiPOD solutions with respect to the L 2 (Ω)-norm.
such a case, system (18) reduces to a scalar equation). We do not provide the HiPOD approximations for l = 8 since they qualitatively coincide with the corresponding HiMod solution. The left panels can be additionally compared with the FE solutions in Fig. 2 to verify the reliability of the HiMod procedure. Finally, the table in Fig. 3 gathers the L 2 (Ω)-norm of the relative error between HiMod and HiPOD solutions, for four different POD bases and for three choices of the viscosity (1, 2.5 and the average over a sampling of 30 random values of µ). The error monotonically decreases for larger and larger values of l, independently of the choice for µ. If we compare the values for µ = 1 and for µ = 2.5 (one of the endpoints and the midpoint of the sampling interval, respectively), we notice a higher accuracy (of about one order of magnitude) for the latter choice. This is rather standard in projection-based reduced order modeling [11]. Concerning the computational saving in terms of CPU time, HiPOD method requires on average O(10 −3 )[s] to be compared with O(10)[s] demanded by HiMod, resulting in a speedup of 10 4 .
Although PGD and HiPOD are not directly comparable due to the different purpose they pursue, we highlight the main pros and cons of the two methods. The explicit dependence of the approximation on the parameters makes PGD an ideal tool to efficiently deal with parametric problems. For any new parameter, a direct evaluation yields the corresponding PGD approximation. On the other hand, HiPOD suffers of the drawbacks typical of the projection-based methods. The main bottleneck is the assembling of the HiMod arrays involved in A POD (µ * ) and f POD (µ * ). When PGD is applied to parametric problems, we recover the possibility to deal with any geometric domain. In such a case, a partial separability is applied to the problem, so that the space independent variables are kept together whereas parameters are separated. This approach clearly looses the computational advantages due to space separability. On the contrary, HiPOD inherits the geometric flexibility of the HiMod reduction, without giving up the spatial dimensional reduction of the problem.