title: Kernel-based active subspaces with application to computational fluid dynamics parametric problems using the discontinuous Galerkin method
authors: Francesco Romor Marco Tezzele Andrea Lario Gianluigi Rozza
text: 
  
   
  
  ]Francesco Romor[francesco.romor@sissa.it]
  ]Marco Tezzele[marco.tezzele@sissa.it]
  ]Andrea Lario[andrea.lario@sissa.it]
  ]Gianluigi Rozza[gianluigi.rozza@sissa.it]
  
  
  Mathematics Area, mathLab, SISSA, via Bonomea 265, I-34136
    Trieste, Italy
  
  Kernel-based Active Subspaces with application to CFD problems
    using Discontinuous Galerkin method
      [
      
  ====================================================================================================
  
  
  
  
  Nonlinear extensions to the active subspaces method have brought remarkable
  results for dimension reduction in the parameter space and response
  surface design. We further develop a kernel-based nonlinear method. In particular
  we introduce it in a broader mathematical framework that contemplates also the
  reduction in parameter space of multivariate objective functions. The
  implementation is thoroughly discussed and tested on more challenging
  benchmarks than the ones already present in the literature, for which dimension
  reduction with active subspaces produces already good results. Finally, we show
  a whole pipeline for the design of response surfaces with the new methodology in
  the context of a parametric CFD application solved with the Discontinuous
  Galerkin method.
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Nowadays, in many industrial settings the simulation of complex
  systems requires a huge amount of computational power. Problems
  involving high-fidelity simulations are usually large-scale, moreover
  the number of solutions required increases with the number of
  parameters. In this context, we mention optimization tasks, inverse problems,
  optimal control problems, and uncertainty
  quantification; they all suffer from the curse of dimensionality, that
  is, in this case, the computational time grows exponentially with the
  dimension of the input parameter space. Data-driven reduced order
  methods (ROM) <cit.> have been developed to deal with such costly outer loop
  applications for parametric PDEs, but the limit for high dimensional parameter spaces remains.
  
  One approach to alleviate the curse of dimensionality is to identify
  and exploit some notion of low-dimensional structure of the
  model or objective function that maps the inputs to the outputs of interest. A possible linear input coordinate transformation technique
  is the Sliced Inverse Regression (SIR) <cit.> approach
  and its extensions <cit.>. Sharing some characteristics with
  SIR, there is the active subspaces (AS)
  property[Some authors refer to the
  active subspaces method, we prefer to employ the term active subspaces property,
  as suggested by Constantine.] <cit.> which, in the last years, has
  emerged as a powerful linear data-driven technique to construct ridge
  approximations using gradients of the model function. AS has been
  successfully applied to quantify uncertainty in the numerical
  simulation of the HyShot II scramjet <cit.>,
  and for sensitivity analysis of an integrated hydrologic
  model <cit.>.
  Reduction in parameter space has been coupled with model order
  reduction techniques <cit.> to enable more complex
  numerical studies without increasing the computational load. We
  mention the use of AS in cardiovascular applications with
  POD-Galerkin <cit.>, in nonlinear structural
  analysis <cit.>, in nautical and naval
  engineering <cit.>, coupled with POD with interpolation for structural
  and computational fluid dynamic (CFD) analysis <cit.>, and with Dynamic Mode
  Decomposition in <cit.>. Applications in
  automotive engineering within a
  multi-fidelity setting can be found in <cit.>, for
  turbomachinery see <cit.>, while
  for results in chemistry see <cit.>.
  Advances in efficient global design optimization with surrogate
  modeling are presented in <cit.> and applied to the shape design of the N+2
  Supersonic Passenger Jet. Applications to enhance optimization methods
  have been developed in <cit.>. AS has also been successfully
  used to reduced the memory consumption of highly parametrized systems
  such as artificial neural networks <cit.>.
  
  Possible extensions and variants of the active subspaces property are
  the Local Active Subspace method <cit.>, the Active
  Manifold method <cit.> which reduces the problem to the
  analysis of a 1D manifold by traversing the level sets of the model
  function at the expense of high online costs, the shared Active Subspace
  method <cit.>, the active subspaces property for
  multivariate functions <cit.>, and more recently an
  extension of AS to dynamical systems <cit.>. Another
  method is Nonlinear Level set Learning
  (NLL) <cit.> which exploits RevNets to reduce the
  input parameter space with a nonlinear transformation.
  
  The search for low dimensional structures is also investigated in machine
  learning with manifold learning algorithms. In this context the Active
  Subspaces methodology can be seen as a supervised dimension reduction
  technique along with Kernel Principal Component Analysis
  (KPCA) <cit.> and Supervised Kernel
  Principal Component Analysis
  (SKPCA) <cit.>. Other methods in the context of
  kernel-based ROMs are <cit.>.  In <cit.> a
  non-linear extension of the active subspaces property based on
  Random Fourier Features <cit.> is
  introduced and compared with machine learning manifold learning
  algorithms for the construction of Gaussian process regressions
  (GPR) <cit.>.
  
  From the preliminary work <cit.> in the context of supervised dimension
  reduction algorithms in machine learning, we develop the kernel-based active subspaces
  (KAS) method. The novelties of our contribution are the following:
  
      
    * regarding the AS theoretical background, we provide an upper bound of the ridge
      approximation error (<ref>) for
      vector-valued objective functions and for a wide collection of probability
      distributions (see Assumption <ref>).
      
    * we extend kernel-based AS to
      vector-valued model functions and develop a detailed algorithmic procedure
      for the optimization of the feature map. We also test different spectral
      measures (see Equation (<ref>) for the definition), differently from <cit.> where only the Gaussian
      measure is employed.
      
    * the application to several test problems of increasing
      complexity. In particular, we mainly test KAS on problems where the active
      subspace is not present or the behaviour is not linear, differently
      from <cit.>, where the comparison is made with KPCA and
      its variants on datasets with linear trends in the reduced parameter space,
      apart from the hyperparaboloid test case that we have also included among
      our toy problems.
      
    * the KAS method is finally applied to a computational fluid
      dynamics problem and compared with the standard AS technique. We study
      the evolution of fluid flow past a NACA 0012 airfoil in a duct
      composed by an initialization channel and a chamber. The motion is
      modelled with the unsteady incompressible Navier-Stokes
      equations, and discretized with the Discontinuous Galerkin
      method (DG) <cit.>. Physical and geometrical parameters
      are introduced and sensitivity analysis of the lift and drag coefficients with respect to
      these parameters is provided.
  
  
  The work is divided as follows: in <ref> we briefly
  present the active subspaces property of a model function with a focus
  on the construction of Gaussian process response surfaces.
  Then, <ref> illustrates the novel method called kernel-based
  active subspaces for both scalar and vector-valued model
  functions.
  Several tests to compare AS and KAS are provided in <ref>
  where we start from scalar functions with radial symmetry, we analyze an epidemiology model and a
  vector-valued output generated from a stochastic elliptic PDE. A
  parametric CFD test case for the study of the flow past a NACA airfoil
  using the Discontinuous Galerkin method is presented in
  <ref>. Finally, we outline some perspectives and future
  studies in <ref>.
  
  
  
  § ACTIVE SUBSPACES FOR PARAMETER SPACE REDUCTION
  
  
  
  Active Subspaces (AS) approach proposed in <cit.>
  and developed in <cit.> is a technique for
  dimension reduction in parameter space. In brief AS are defined as the leading
  eigenspaces of the second moment matrix of the model function's
  gradient (for scalar model functions) and constitutes a global
  sensitivity index <cit.>. In the context of ridge
  approximation, the choice of the active subspace corresponds to the
  minimizer of an upper bound of the mean square error obtained through
  Poincaré-type inequalities <cit.>. After
  performing dimension reduction in the parameter space through AS, the
  method can be applied to reduce the computational costs of different
  parameter studies such as inverse problems, optimization tasks and
  numerical integration. In this work we
  are going to focus on the construction of response surfaces with
  Gaussian process regression.
  
  
  The quantities related to the input space are:
  
  
  
    * m∈ℕ the dimension of the input space,
  
    * (Ω, ℱ, P) the probability space,
  
    * :(Ω, ℱ,P)→ℝ^m, the absolutely continuous random vector representing the parameters,
  
    * ρ:ℝ^m→ℝ, the probability density of 𝐗 with support 𝒳⊂ℝ^m.
  
  The quantities related to the output are:
  
  
    * d∈ℕ the dimension of the output space,
  
    * V=(ℝ^d, R_V) the Euclidean space with metric
  R_V∈ℳ(d× d)[In this work with ℳ(m × n) we denote the set of real matrices with
  m rows and n columns.] and norm
  
      ‖𝐱‖^2_R_V=𝐱^TR_V𝐱,
  
  
    * f:𝒳⊂ℝ^m→ V, the quantity/function
  of interest, also called objective function in optimization tasks.
  
  
  
  Let ℬ ( ℝ^m ) be the Borel σ-algebra of
  ℝ^m. We will consider the Hilbert space L^2 (ℝ^m, ℬ (ℝ^m), ρ ;  V), of the measurable functions f:(ℝ^m,
  ℬ(ℝ^m),ρ)→(ℝ^d,R_V)
  such that
  
      ‖ f‖^2_L^2:=∫_𝒳‖ f()‖^2_R_V dρ()≤∞;
  
  and the Sobolev space H^1(ℝ^m, ℬ (ℝ^m), ρ ;  V) of measurable functions f:(ℝ^m,
  ℬ(ℝ^m),ρ)→(ℝ^d,R_V) such that
  
      ‖ f‖^2_H^1:=‖ f‖^2_L^2+‖∇ f‖^2_L^2=‖ f‖^2_L^2+|f|^2_H^1≤∞
  
  where ∇ f is the weak derivative of f, and ‖∇
  f‖_L^2=:|f|_H^1.
  
  We briefly recall how dimension reduction in parameter space is
  achieved in the construction of response surfaces. The first step
  involves the approximation of the model function with ridge
  approximation. We will follow <cit.>
  for a review of the method.
  
  The ridge approximation problem can be stated in the following way:
  
  
  Let ℬ(ℝ^m) be the
  Borel σ-algebra of ℝ^m. Given
  r∈ℕ, r≪ d and a tolerance ϵ≥ 0, find the
  profile h:(ℝ^m, ℬ(ℝ^m),
  ρ)→ V and the r-rank projection
  P_r:ℝ^m→ℝ^m such that
  
  
      𝔼_P[‖
      f()-h(P_r)‖^2_R_V]≤ϵ^2.
  
  
  
  In particular we are interested in the minimization problem
  
  
      _P_r ∈ℳ(m× m) 𝔼_P [ ‖
      f()-h(P_r) ‖^2_R_V ],
  
  where h=𝔼_ρ[ f |σ(P_r) ] is the conditional
  expectation of f under the distribution ρ given the
  σ-algebra σ(P_r). The range of the projector P_r,
  ℝ^r∼Im(P_r)⊂ℝ^m, is the reduced
  parameter space. The kernel of the projector P_r,
  ℝ^m-r∼Im(P_r)⊂ℝ^m, is the inactive subspace. The
  existence of h is guaranteed by the Doob-Dynkin
  lemma <cit.>. The function h is
  proven to be the optimal profile for each fixed P_r, as a
  consequence of the definition of the conditional expectation of a
  random variable with respect to a σ-algebra.
  
  Dimension reduction is effective if the
  inequality (<ref>) is satisfied for a specific
  tolerance. The choice of r is certainly of central
  importance. The dimension of the reduced parameter space can be chosen a
  priori for a specific parameter study (for example r-dimensional
  regression), it can be chosen in order to satisfy the
  inequality (<ref>) or it is determined to guarantee a
  good accuracy of the numerical method used to evaluate it [Corollary
  3.10, <cit.>].
  
  Dividing the left term of the inequality (<ref>) with 𝔼_ρ[‖
  f()-𝔼_ρ[f(])‖^2_R_V] we obtain the Relative Root Mean
  Square Error (RRMSE) and since it is a normalized quantity, we will use it to
  make comparisons between different models
  
  
      RRMSE = √(𝔼_P[‖
          f()-h(P_r)‖^2_R_V]/𝔼_P[‖
          f()-𝔼_P[f()]‖^2_R_V]).
  
  We remark that P_r is not unique.
  It can be shown that if h is the optimal profile, then P_r is not
  uniquely defined and can be chosen arbitrarily from the set
  { Q_r : ℝ^m →ℝ^m |   Q_r
   =   P_r}, see [Proposition 2.2, <cit.>].
  
  The following lemma is the key ingredient in the proof of the existence of an
  active subspace. It is inherently linked to probability Poincaré inequalities of
  the kind
  
  
      ∫_𝒳‖ h(𝐱)‖_L^2^2 dρ()≤ C_P(𝒳, ρ) ∫_𝒳‖∇ h(𝐱)‖_L^2^2 dρ(),
  
  for zero-mean functions in the Sobolev space h∈
  H^1(𝒳), where C_P(𝒳, ρ) is the Poincaré
  constant dependent on the domain 𝒳 and on the probability
  density functions (p.d.f.), ρ. We need to make the following
  assumption to prove the next lemma and the next theorem.
  
  
  The probability density function ρ:𝒳→ℝ belongs to one of the following classes:
  
  
    * 𝒳 is convex and bounded, ∃δ, D>0: 0<δ≤‖ρ(𝐱)‖_L^∞≤ D<∞ ∀∈𝒳,
  
    * ρ()∼exp(-V()) where V:ℝ^m→ (-∞,∞] ,V∈𝒞^2 is α-uniformly convex,
  
      𝐮^THess(V())𝐮≥ α‖𝐮‖^2_2,  ∀, 𝐮∈ℝ^m
  
  where Hess(V()) is the Hessian of V().
  
    * ρ()∼exp(-V()) where V is a convex function. In this case we require also f Lipschitz continuous.
  
  
  In particular the uniform distribution belongs to the first class, the
  multivariate Gaussian distribution 𝒩(m, Σ) to the
  second with α=1/(σ_max(Σ)) and the exponential and
  Laplace distributions to the third. A complete analysis of the various
  cases is done in <cit.>.
  
  
  Let (Ω, ℱ, P) be a probability space, :(Ω,
  ℱ,P)→ℝ^m an absolutely continuous
  random vector with probability density function ρ belonging to
  one of the classes from the Assumption <ref>. Then the
  following inequality is satisfied
  
      𝔼_ρ[(h-𝔼_ρ[h|σ (P_r)])^2|σ(P_r)]≤ C_P(P_r, ρ) 𝔼_ρ[‖(I-P_r^T)∇ h‖^2_2|σ (P_r)]
  
  for all scalar functions h∈ H^1(𝒳) and for all r-rank
  orthogonal projectors, P_r, where C_P(P_r, ρ) is the Poincaré constant depending on P_r and on the p.d.f. ρ.
  
  A summary of the values of the Poincaré constant in relationship with the choice of the probability density function ρ is reported in <cit.>.
  
  In the next theorem the projection P_r will depend on the output function f, so also the Poincaré constant C_P(P_r, ρ) will depend in fact on f.
  
  We introduce the following notation for the matrix that substitutes
  the uncentered covariance matrix of the gradient ∇ f in the case of the
  application of AS to scalar model
  functions <cit.>
  
      H=∫_𝒳 (D_f())^T R_V(ρ) (D_ f())   dρ ().
  
  where D_ f() ∈ℳ(d× m) is the
  Jacobian matrix of f. The matrix R_V(ρ) depends on the class which ρ belongs to, see Appendix <ref>.
  
  
  
  
  Under the hypothesis <ref>, let f∈ H^1
  (ℝ^m, ℬ (ℝ^m), ρ ;  V) and let the
  p.d.f. ρ satisfy Lemma <ref> and
  Assumption <ref>. Then the solution P_r of the ridge
  approximation problem <ref> is the orthogonal
  projector to the eigenspace of the first r-eigenvalues of H ordered by magnitude
  
      Hv_i=λ_i v_i    ∀ i∈{1,…,m},    P_r=∑_j=1^rv_j⊗ v_j,
  
  with r∈ℕ chosen such that
  
  
      𝔼_ρ[‖
      f-h‖^2_R_V]≤ C(C_P, τ) (∑_i=r+1^mλ_i)^1/1+τ≤ϵ^2.
  
  with C(C_P, τ) a constant depending on τ>0 related to the choice of
  ρ and on the Poincaré constant from
  lemma <ref>, and h=𝔼_ρ[ f
  |σ(P_r) ] is the conditional expectation of f given the σ-algebra
  generated by the random variable P_r∘.
  .
  
  
  
  This theorem summarizes the results from Proposition 2.5 and Proposition 2.6
  of <cit.>, and from Lemma 3.1, Lemma 4.2, Lemma 4.3, Lemma 4.4
  and Theorem 4.5 of <cit.>. The proof is expanded in
  Appendix <ref>.
  
  The eigenspace span{v_1,…, v_r}⊂^m
  is the active subspace and the remaining eigenvectors generate the
  inactive subspace
  span{v_r+1,…,v_m}⊂^m. The
  condition f∈ L^2 (^m, ℬ (^m), ρ ;  V) is
  necessary for f to satisfy the error bound (<ref>).
  
  For the explicit procedure to compute the active subspace given its dimension
  r see Algorithm <ref>: from W_1 and W_2 we define the approximations of
  the projector P_r with P_r = W_1W_1^T.
  
  
  
  
  
  
  
  
   §.§ Response surfaces
  
  
  The term response surface refers to the general procedure of finding
  the values of a model function f for new inputs without directly
  computing it but exploiting regression or interpolation from a
  training set {_i, f(_i)}. The procedure for
  constructing a Gaussian process response is reported in
  Algorithm <ref>, while in Algorithm <ref> we
  show how to exploit it o predict the model function
  at new input parameters.
  
  Directly applying the simple Monte Carlo method with N samples we get a
  reduced approximation of f as
  
  
      (h_ϵ∘ P_r) () = 𝔼_ρ[ f |
            σ(P_r) ] ≈1/N∑^N_i=1 f(P̂_r
           + (I_d - P̂_r ) 𝐘_i ) =: ĥ_ϵ, N
          (P̂_r ),
  
  where we have made explicit the dependence of the optimal profile
  h_ϵ on ϵ,
  𝐘_1,…,𝐘_N are independent and
  identically distributed samples of 𝐘∼ρ, and
  P̂_r is an approximation of P_r obtained with the simple Monte Carlo
  method from H, see Algorithm <ref>. An intermediate approximation error is obtained
  employing the Poincaré inequality and the central limit theorem for
  the Monte Carlo approximation
  
  
      𝔼_P[ (f()-ĥ_ϵ, N (P̂_r))^2
          ] ≤ C_1 (1+N^-1/2)^2
          (λ_n+1+ … + λ_m ),
  
  where C_1 is a constant, and λ_n+1, …, λ_m are
  the eigenvalues of the inactive subspace of H [Theorem 4.4, <cit.>].
  
  In practice ĥ_ϵ, N(P̂_r) is approximated with a
  regression or an interpolation such that a response surface
  ℛ satisfying
  𝔼_ρ[(ĥ_ϵ, N (P̂_r)-ℛ(P̂_r)^2)
  ] ≤ C_2 δ is built, where C_2 is a constant, and
  δ depends on the chosen method. An estimate for the successive
  approximations
  
      f() ≈h_ϵ(P_r)≈ĥ_ϵ, N(P̂_r)≈ℛ_ϵ,N,δ(P̂_r),
  
  is given by
  
  
      𝔼_P   [(f()-ℛ(P̂_r))^2] 
         ≤ C_1(1+N^-1/2)^2(
            τ(λ_1+…+λ_n)^1/2+(λ_n+1+…+λ_m)^1/2)^2+C_2λ
  
  where dist(Im(P_r), Im(P̂_r))
  ≤τ, and λ_i are the eigenvalues of H [Theorem 4.8, <cit.>].
  
  In our numerical simulations we will build the response surface
  ℛ with Gaussian process
  regression (GPR) <cit.>.
  
  
  
  
  
  
  
  
  
  § KERNEL-BASED ACTIVE SUBSPACES EXTENSION
  
  
  
  
  Keeping the notations of section <ref>, :(Ω,
  ℱ, P) →^m is the absolutely continuous
  random vector representing the m-dimensional inputs with density
  ρ:𝒳⊂^m →, and
  f:𝒳⊂^m → (V, R_V) is the
  model function that we assume to be continuously differentiable and
  Lipschitz continuous.
  
  One drawback of sufficient dimension reduction with AS applied to
  ridge approximation is that if a clear linear trend is missing,
  projecting the inputs as P_r represents a loss of accuracy on the
  approximation of the model f that may not be compensated even by
  the choice of the optimal profile
  h∘ P_r=𝔼_ρ [f | σ(P_r)].
  In order to overcome this,
  non-linear dimension reduction to one-dimensional parameter space
  could be achieved discovering a curve in the space of parameters that
  cuts transversely the level sets of f, this variation is presented
  in <cit.> as Active Manifold. Another approach could
  consist in finding a diffeomorphism ϕ that reshapes the level
  sets such that subsequently applying AS dimension reduction to the new
  model function f∘ϕ=f
  could be more profitable:
  
      𝒳⊂ℝ^m rϕ[swap]drf   ϕ(𝒳)⊂ℝ^m
          df
          V
  
  
  Unfortunately constructing the Active Manifold or finding the right
  diffeomorphism ϕ could be a complicated matter. If we renounce to
  have a backward map and we weaken the bond of the method with the
  model, we can consider an immersion ϕ from the space of
  parameters 𝒳 to an infinite-dimensional Hilbert space
  ℍ obtaining
  
      𝒳⊂ℝ^m rϕ[swap]drf   ϕ(𝒳)⊂ℍdf
          V
  
  This is a common procedure in machine learning in order to increase
  the number of features <cit.>. Then AS is applied
  to the new model function
  f:ϕ(𝒳)⊂ℍ→ V with
  parameter space ϕ(𝒳)⊂ℍ. A response surface can be built with <ref> remembering to replace every occurrence of the inputs  with their images ϕ(). A synthetic scheme of the procedure is represented in <ref>.
  
  
  
  
  
  
  
  
  
  
  
  In practice we consider a discretization of the infinite-dimensional
  Hilbert space ^D ≃ℍ with D>m. Dimension reduction with AS results in the choice of a r-rank
  projection in the much broader set of r-rank projections in
  ℍ.
  
  Since for AS only the samples of the
  Jacobian matrix of the model function are employed, we can ignore the
  definition of the new map
  f: ϕ(𝒳)⊂ℍ→ (V, R_V)
  and focus only on the computation of the Jacobian matrix of
  f with respect to the new input variable
  :=ϕ(). The uncentered covariance
  matrix becomes
  
  
      H    = ∫_ϕ(𝒳) [
                            (D_f)^T() 
                            ] R_V  [
                            (D_f)()  ]  
                            dμ()
         = ∫_𝒳 [
                                                  (D_f)^T(ϕ())
                                                   ] R_V  [
                                                  (D_f)(ϕ())
                                                   ]  dℒ_(),
  
  where μ:=ϕ_#(ℒ_) is the pushforward
  probability measure of ℒ_ (the law of
  probability of ) with respect to the map ϕ. Simple
  Monte Carlo can be applied sampling from the distribution ρ in
  the input space 𝒳
  
  
      H    = ∫_𝒳 [
                            (D_f)^T(ϕ())
                             ] R_V  [
                            (D_f)(ϕ())
                             ]  dℒ_()
         ≈1/M∑_i=1^M  [(D_f)^T(ϕ(_i))
         ]  R_V  [ (D_f)(ϕ(_i))  ] .
  
  The gradients of f with respect to the new input variable
   are computed from the known values D_f with
  the chain rule.
  
  
  The application of the chain rule to the composition of
  functions f∘ϕ:ℝ^m→ℍ→
  V is applicable if f is defined in an open set
  U⊃ϕ(𝒳).
  If ϕ is non singular and also injective the new input space is a
  m-dimensional submanifold of ℍ. If ϕ is also smooth
  there exists a smooth extension of
  f:ϕ(𝒳)⊂ℍ→ V onto the whole
  domain ℍ,
  see Proposition 1.36 from <cit.>.
  
  
  If the Hilbert space ℍ has finite dimension ℍ∼ℝ^D this procedure leaves us with an underdetermined linear system to solve
  for D_f
  
      D_f(ϕ())Dϕ()=D_f(),
        
         D_f(ϕ())=D_f()(Dϕ())^†,
  
  where ^† stands for the right Moore-Penrose inverse of the matrix Dϕ() with rank r, that is
  
      (Dϕ())^†=V Σ^†U^T,
  
  with the usual notation for the singular value decomposition (SVD) of Dϕ()
  
      Dϕ() = U Σ V^T,
  
  and Σ^†∈ℳ(r× r) equal to the diagonal matrix
  with the inverse of the singular values as diagonal elements. As anticipated if
  f is smooth enough and ϕ is an embedding, so that Dϕ has full rank, the previous system has an
  unique solution. The most crucial part is the evaluation of the gradients
  D_f() from the input output couples, when they are not available
  analytically or from the adjoint method applied to PDEs models:
  different approaches are present in the literature, like local polynomial regressions and
  Gaussian process regression on the whole domain to approximate the gradients;
  both are available in the ATHENA package <cit.>. For an estimate
  of the ridge approximation error due to inaccurate gradients see <cit.>.
  
  Finally, we remark that in the AS method we approximate the random variable  as
  
      P_r  = 𝐯_1(𝐯_1·)+…+𝐯_r
      (𝐯_r·),
  
  with {𝐯_i}⊂^m the active
  eigenvectors, whereas with KAS the reduced input space is contained in
  ℋ
  
      P_r  = 𝐯_1 (𝐯_1 ·ϕ())+…+𝐯_r
      (𝐯_r ·ϕ()),
  
  with {𝐯_i}⊂ℋ the active eigenvectors of
  KAS. In this case the model is enriched by the non-linear feature map
  ϕ.
  
  
  
  
  
   §.§ Choice of the Feature Map
  
  The choice for the map ϕ is linked to the theory of Reproducing
  Kernel Hilbert Spaces (RKHS) <cit.>, and it is
  defined as
  
      =ϕ()=√(2/D) σ_f cos(W+𝐛),
      cos(W+𝐛):=1/√(D)(cos(W[1, :] · +b_1),…, cos(W[D, :] · +b_D))^T
  
  where σ_f is an hyperparameter corresponding to the empirical
  variance of the model, W∈ℳ(D× m) is the
  projection matrix whose rows are sampled from a probability
  distribution μ on ^m and 𝐛∈^D is a bias term whose components are
  sampled independently and uniformly in the interval [0, 2π]. We
  remark that its Jacobian can be computed analytically as
  
      ∂ z^j/∂
        x^i=-√(2/D) σ_f sin(∑_k=1^D W_ik_k+𝐛_k) W_ij,
  
  for all i ∈{1, …, m}, and for all j ∈{1, …, D}.
  
  We remark that in order to guarantee the correctness of the procedure for evaluating
  the gradients we have to prove that the feature map is injective and
  non singular. In general however the feature map (<ref>) cannot
  not be injective due to the periodicity of the cosine but at least it is almost
  surely non singular if the dimension of the feature space is high enough.
  
  The feature map (<ref>) is not the only effective
  immersion that provides a kernel-based extension of the active
  subspaces. For example an alternative is the following composition of
  a linear map with a sigmoid
  
      ϕ()=C/1+α  e^ -W,
  
  where C is a constant, α is an hyperparameter to be tuned, and
  W ∈ℳ (D, m) is, as before, a matrix whose rows are
  sampled from a probability distribution on ^m.
  
  Other choices involve the use of Deep Neural Networks to learn the
  profile h and the projection function P_r of the ridge approximation
  problem <cit.>.
  
  The tuning of the hyperparameters of the spectral measure consists
  in a global optimization problem where the dimension of the domain can
  vary between 1 and the dimension of the input space m. The object
  function to optimize is the relative root mean square error (RRMSE)
  
      RRMSE(Y_test,
      T_test) = √(∑^N_i=1(t_i-y_i)^2/∑^N_i=1(t_i-y)^2),
  
  where T_test=(t_i)_i∈{1, …, N} are the
  predictions obtained from the response surface built
  with KAS and associated
  to the test set, Y_test=(y_i)_i∈{1, …, N} are the targets
  associated to the test set, and y is the mean value of the
  targets. We implemented a logarithmic grid-search, see Algorithm <ref>, making use of
  the SciPy library <cit.>. Another choice could be
  Bayesian stochastic optimization implemented in the open-source library GPyOpt <cit.>.
  
  The tuning of the hyperparameters of the spectral measure chosen is
  the most computationally expensive part of the procedure. We report the
  computational complexity of the  algorithms introduced to have a better understanding of
  the additional cost implied by the implementation of response
  surface design with KAS. Let us assume that the number of random Fourier features D, the
  number of input, output, and gradient samples M, and the dimension of the
  parameter space m, are ordered in this manner D>M>m, as is usually the case,
  and that the quantity of interest f
  is a scalar function. The cost of computing an active subspace is O(Mm^2), that is the cost of the
  SVD of the gradients matrix
  dY used to get the active and inactive eigenvectors in
  Algorithm <ref>. The cost of the training of a response surface with
  Gaussian process regression in Algorithm <ref> depends on
  the cost of minimization of the log-likelihood: each evaluation of the
  log-likelihood involves the
  computation of the determinant and the inverse of the regularized Gram matrix K(θ)+σ
  I_M, that is O(M^3). Finally, the cost for the evaluation of the
  kernel-based active subspace is associated to the SVD of dY that is O(DM^2) in
  Algorithm <ref>, and to the resolution of the overdetermined linear system to
  obtain the gradients dY, that is M times O(Dm^2) since it is related to the evaluation of the pseudo-inverse of Dϕ. So, the
  computational complexity for the response surface design with AS and GPR is
  O(n_GPRM^3), while for the response surface design with KAS and GPR is
  O(n_grid-search
  n(DM^2/n^2+M/nDm^2+n_GPRM^3/n^3)), where
  n_GPR is the maximum number of steps of the
  optimization algorithm used to minimize the log-likelihood,
  n_grid-search is the number of hyperparameter instances γ∈
  G to try in Algorithm <ref>, and n is the number of batches in the
  n-fold cross validation procedure. In particular, for each grid search
  hyperparameter the main cost is associated to the GPR training since
  n_GPR usually satisfy Dn<n_GPRM, when the optimizer chosen
  is L-BFGS-B from SciPy <cit.>, accounting also for the number of restarts of the
  optimizer: in the numerical tests we performed the number of restarts of the training of
  the GPR is problem-dependent but always less than 10. In
  general, the number n_grid-search depends on the chosen application,
  and the multiplicative factor between the computational complexity of the
  response surface design procedure with KAS or AS is lower than 3n_grid-searchn.
  
  
  
  
  
  
   §.§ Random Fourier Features
  
  The motivation behind the choice for this map from Equation (<ref>)
  comes from the theory on Reproducing Kernel Hilbert Spaces. The infinite-dimensional Hilbert
  space (ℍ, ⟨·,·⟩) is assumed to be a
  RKHS with real shift-invariant kernel
  k:𝒳×𝒳→ℝ with k(0)=1
  and feature map ϕ.
  
  In order to get a discrete approximation of
  ϕ:𝒳⊂ℝ^m →ℍ, random
  Fourier features are employed <cit.>. Bochner's theorem <cit.> guarantees the existence of a
  spectral probability measure μ such that
  
      k(𝐱, 𝐲) = ∫_ℝ^m e^iω·(𝐱-𝐲) dμ(ω).
  
  From this identity we can get a discrete approximation of the scalar
  product ⟨·,·⟩ with Monte Carlo method,
  exploiting the fact that the kernel is real
  
      ⟨ϕ(),ϕ() ⟩=k(,
      )≈1/D∑_i=1^Dcos(ω_i· + b_i)
      cos(ω_i · +b_i)=^T, 
      =1/√(D)(cos(ω_1 · +b_1),…, cos(ω_D · +b_D)),
  
  and from this relation we obtain the approximation
  ϕ≈𝐳. The sampled vectors
  {ω_i}_i=1,…,D are called random Fourier
  features. The scalars {b_i}_i=1,…,D are bias terms
  introduced since in the approximation we have excluded some trigonometric terms from the following initial expression
  
      1/D∑_i=1^D( cos(ω_i ·𝐱)
      cos(ω_i ·𝐲)-sin(ω_i ·𝐱)
      sin(ω_i ·𝐲)).
  
  
  Random Fourier features are frequently used to approximate
  kernels. We consider only spectral probability measures which have a
  probability density, usually named spectral density. In the approximation of the
  kernel with random Fourier features, under some
  regularity conditions on the kernel, an explicit
  probabilistic bound depending on the dimension of the feature space
  D can be proved <cit.>.  This technique is used to scale up Kernel Principal
  Component Analysis <cit.> and Supervised Kernel Principal Component
  Analysis <cit.>, but in the case of kernel-based AS
  the resulting overdetermined linear system employed to compute the
  Jacobian matrix of the new model function increases in dimension
  instead.
  
  The most famous kernel is the squared exponential kernel also called Radial Basis Function kernel (RBF)
      
      k_RBF(, 𝐲)=exp( -‖-𝐲‖^2/2l^2),
  
      where l is the characteristic length-scale. The spectral density
      is Gaussian 𝒩(0, 1/4π^2l^2):
      
      S(ω)=(2π l^2)^D/2exp(-2π^2l^2ω^2).
  
  
  Thanks to Bochner's theorem to every probability distribution that
  admits a probability density function corresponds a stationary
  positive definite kernel. So having in mind the
  definition of the feature map ϕ from
  Equation (<ref>), we can choose any probability
  distribution for sampling the
  random projection matrix W∈ℳ(D× m)
  without focusing on the corresponding kernel since it is not needed by
  the numerical procedure.
  
  After the choice of the spectral measure the corresponding
  hyperparameters have to be tuned. This is linked to the choice of the
  hypothesis models in machine learning and it is usually carried out for
  the hyperparameters of the employed kernel. From the choice of the
  kernel and the corresponding hyperparameters some regularity
  properties of the model are implicitly
  assumed <cit.>.
  
  
  
  
  § BENCHMARK TEST PROBLEMS
  
  
  In this section we are going to present some benchmarks to prove
  the potential gain of KAS over standard linear AS, for both scalar and
  vectorial model functions. In particular we test KAS on radial
  symmetric functions, with 2-dimensional and 8-dimensional
  parameter spaces, on the approximation of the reproduction number
  R_0 of the SEIR model, and finally on a vectorial output function
  that is the solution of a Poisson problem.
  
  One dimensional response surfaces are built following the algorithm described in
  <ref>. The tuning of the hyperparameters
  of the feature map is carried out with a logarithmic grid-search and
  5-fold cross validation for the Ebola test case, while for the other
  cases we employed Bayesian stochastic optimization implemented
  in <cit.> with 3-fold cross validation. The score function
  chosen is the Relative Root Mean Square Error (RRMSE). The spectral measure
  for each test case is chosen by brute force among the Laplace, Gaussian, Beta
  and multivariate Gaussian distributions. The number of Fourier features is not
  established based on a criterion but we have seen experimentally that above a certain
  threshold the number of features is high enough to at least reproduce the accuracy of
  the AS method. Since the most sensitive part to the final accuracy of the
  response surface is the tuning
  of the hyperparameters of the spectral measures, we suggest to choose an
  affordable number of features between 1000 and 2000, and focus on the tuning of
  said hyperparameters instead.We remark that the number of samples employed
  is problem dependent: some heuristics to determine it can be found
  in <cit.>, but the crucial point is that additional
  training samples with respect to the ones used for the AS method are not needed for the
  novel KAS method. Moreover, the CPU time for the hyperparameters tuning
  procedure is usually negligible with respect to the time required to obtain
  input-output pairs from the numerical simulation of PDEs models: in our
  applications the tuning procedure's computational time is in the order of
  minutes (usually around 10-15 minutes for most testcases), while for
  the CFD application of Section <ref> it is in the order of days and
  for the stochastic elliptic partial differential equation of
  Subsection <ref> it is in the order of hours. We also remark that the
  tuning Algorithm <ref>, the GPR training restarts, and the choice
  of the spectral measure can be easily
  parallelized.
  
  For the radial
  symmetric and Ebola test cases the inputs are sampled from a uniform
  distribution with problem dependent ranges.For the stochastic
  elliptic partial differential case the inputs are the coefficients of
  a Karhunen-Loève expansion and are sampled from a normal
  distribution. All the computations regarding AS and KAS are done using
  the open source Python package called ATHENA <cit.>.
  
  
  
  
   §.§ Radial symmetric functions
  
  Radial symmetric functions represent a class of model functions for
  which AS is not able to unveil any low dimensional behaviour. In fact
  for these functions any rotation of the parameter space produce the
  same model representation. Instead Kernel-based AS is able to overcome
  this problem thanks to the mapping onto the feature space.
  
  We present two benchmarks: an 8-dimensional hyperparaboloid defined as
  
      f:[-1, 1]^8 ⊂^8 → ,         
      f() = 1/2‖‖^2 ,
  
  and the surface of revolution in ^3 with generatrix g(x)=
  sin(x^2)
  
      f:[-3, 3]^2 ⊂^2 → ,         
      f() = g(‖‖)=sin(‖‖^2) .
  
  The gradients are computed analytically.
  
  For the hyperparaboloid we use N_s = 500 independent, uniformly
  distributed training samples in [-1, 1]^8, while for the sine case the
  training samples are N_s = 800 in [-3, 3]^2. In both cases the test samples are 500. The feature space has dimension
  1000 for both the first and the second case. The spectral
  distribution chosen is the multivariate normal with
  hyperparameter a uniform variance λ I_d, and a product of
  Laplace distributions with γ and b as
  hyperparameters, respectively. The tuning is carried out with 3-fold
  cross validation. The results are summarized in <ref>.
  
  
  
  Looking at the eigenvalues of the uncentered covariance matrix of the
  gradients H for the hyperparaboloid case
  in <ref>,
  we can clearly see how the decay for AS is almost absent, while using
  KAS the decay after the first eigenvalue is pronounced, suggesting the
  presence of a kernel-based active subspace of dimension 1.
  
  
  
  The one-dimensional sufficient summary plots, which are f()
  against W_1^T — in the AS case — or against W_1^T ϕ
  () — in the KAS case —, are shown in <ref> and
  <ref>, respectively. On the left panels we present the
  Gaussian process
  response surfaces obtained from the active subspaces reduction, while
  on the right panels the ones obtained with the kernel-based AS
  extension. As we can see AS fails to properly reduce the parameter
  spaces, since there are no preferred directions over which the model
  functions vary the most. The KAS approach, on the contrary, is able to
  unveil the corresponding generatrices. This results in a reduction of the RMS by
  a factor of at least 3 (see <ref>).
  
  
  
  
  
  
  
   §.§ SEIR model for the spread of Ebola
  
  In most engineering applications the output of interest presents a
  monotonic behaviour with respect to the parameters. This means that,
  for example, the increment in the inputs produces a proportional response
  in the outputs. Rarely the model function has a radial symmetry, and
  in such cases the parameter space can be divided in subdomains, which
  are analyzed separately. In this section we are going to present a
  test case where there is no radial symmetry, showing that, even in
  this case the kernel-based AS presents better performance with respect
  to AS.
  
  For the Ebola test case[The dataset was taken from
    <https://github.com/paulcon/as-data-sets>.], the output of
  interest is the basic reproduction
  number R_0 of the SEIR model, described in <cit.>,
  which reads
  
      R_0 =β_1 +β_2ρ_1 γ_1/ω +
        β_3/γ_2ψ/γ_1+ ψ,
  
  with parameters distributed uniformly in Ω⊂^8. The
  parameter space Ω is an hypercube defined by the lower and upper bounds
  summarized in <ref>.
  
  
  
  We can compare the two one-dimensional response surfaces obtained with
  Gaussian process regression. The training samples are N_s = 800, and
  we use 1000 features. As spectral measure we use again the multivariate
  gaussian distribution 𝒩(0, Σ) with
  hyperparameters the elements of the diagonal of the covariance
  matrix. The tuning is carried out with 5-fold cross validation. Even
  in this case, the KAS approach results in smaller RMS with respect to
  the use of AS (around 60% less), as reported in <ref>.
  In <ref> we report the comparison of the two approaches
  over an active subspace of dimension 1.
  
  
  
  
  
  
   §.§ Elliptic Partial Differential Equation with random coefficients
  
  
  In our last benchmark we apply the kernel-based AS to a
  vectorial model function, that is the solution of a Poisson
  problem with heterogeneous diffusion coefficient. We refer
  to <cit.> for an application, on the
  same problem, of the AS approach.
  
  We consider the following stochastic Poisson problem on the square 𝐱=(x, y)∈Ω := [0,
  1]^2:
  
      -∇· (κ ∇ u)=1,    ∈Ω, 
      
      u = 0,    ∈∂Ω_top∪∂Ω_bottom,
      
      u = 10 y(1-y),    ∈∂Ω_left,
      𝐧·∇ u = 0,    ∈∂Ω_right,
  
  with homogeneous Neumann boundary condition on the right side of the domain, that is ∂Ω_right, Neumann boundary conditions on the left side of the domain, that is ∂Ω_left, and Dirichlet boundary conditions on the remaining part of ∂Ω.
  The diffusion coefficient κ:(Ω, 𝒜, P)×Ω→, where
  𝒜 is a σ-algebra, is such that log(κ) is a Gaussian random field, with covariance function C(,) defined by
  
      C(, ) = exp(-‖ - ‖^2/β^2),  ∀,∈Ω,
  
  where β=0.03 is the correlation length. This random field is approximated with the truncated Karhunen-Loève decomposition
  
      κ(s, ) ≈exp(∑_i=0^m X_i(s)
        γ_iψ_i () )    ∀ (s,
      ) ∈Ω×Ω,
  
   where (X_i)_i∈ 1,…, m are independent standard normal
   distributed random variables, and (γ_i,
   ψ_i)_i∈ 1,…, d are the eigenpairs of the
   Karhunen-Loève decomposition of the zero-mean random field κ.
  
  In our simulation the domain Ω is discretized with a triangular
  unstructured mesh 𝒯 with 3194 triangles. The parameter
  space has dimension m=10. The simulations are carried out with the finite
  element method (FEM) with polynomial order one, and for each
  simulation the parameters (X_i)_i=1,… m are sampled from a
  standard normal distribution.  The solution u is evaluated at
  d=1668 degrees of freedom, thus (V, R_V)≈(ℝ^d,
  S+M) where the metric R_V is approximated with the sum of the
  stiffness matrix S∈ℝ^d×ℝ^d and the
  mass matrix M∈ℝ^d×ℝ^d. This sum is a
  discretization of the norm of the Sobolev space H^1(Ω). The
  number of features used in the KAS procedure is D=1500, the number
  of different independent simulations is M=1000.
  
  Three outputs of interest are considered. The first target function
  f:^m → is the mean value of
  the solution at the right boundary ∂Ω_right, which reads
  
      f(𝐗)=
      1/|∂Ω_right|∫_∂Ω_rightu(s) ds,
  
  and it is used to tune the feature map minimizing the RRMSE of the
  Gaussian process regression, as described in Algorithm <ref>. A
  summary of the results for the first output is reported in
  <ref>. The plots of the regression are reported in
  <ref>. Even in this case both from a qualitative
  and a quantitative point of view, the kernel-based approach achieves
  the best results.
  
  
  
  The second output we consider is the solution function
  
      f:^m→(V, R_V)≈(^d, S),     
      f()=u∈^d.
  
  This output can be employed as a
  surrogate model to predict the solution u given the parameters
  𝐗 that define the diffusion coefficient instead of carrying
  out the numerical simulation. The surrogate model should be
    constructed over the span of the modes identified by the chosen
    reduction strategy, after projecting the data. AS and KAS modes are
  distinguished but can detect some common regions of interest as shown in
  <ref>.
  
  The third output is the evaluation of the solution at a specific degree of freedom
  with index î, that is
  
      f:ℝ^m→ℝ,    
      f()=u_î∈,
  
  in this case the dimension of the input space is
  m=100. Since we use a Lagrangian basis in the finite element formulation and the polinomial order is 1, the node of the mesh associated to the chosen degree of freedom has coordinates [0.27, 0.427]∈Ω. Qualitatively we can see from <ref> that the AS
  modes locate features in the domain which are relatively more regular
  with respect to the KAS modes. To obtain this result we increased the dimension of the input space, otherwise not even the AS modes could locate properly the position in the domain Ω of the degree of freedom.
  
  In the second and third case the diffusion coefficient is given by
  
      κ() = exp(∑_i=1^D 𝐯_j[i]
        ψ_j() )    ∀ (s,
      ) ∈Ω×Ω,
  
  where 𝐯_j∈ℝ^D, j∈{1, …,D}, is the
  j-th active eigenvector from the KAS procedure and the functions
  Ψ:=(ψ_1,…,
  ψ_D) are defined by
  
      Ψ = ϕ(Ψ),
  
  where ϕ is the feature map defined in
  Equation (<ref>) with the projection matrix W and bias
  b, and Ψ:=(γ_1ψ_1,…,
  γ_mψ_m).
  
  The gradients of the three outputs of interest considered are evaluated with the adjoint method.
  
  
  
  
  
  
  § A CFD PARAMETRIC APPLICATION OF KAS SOLVED WITH THE DG METHOD
  
  
  
  We want to test the kernel-based extension of the active subspaces
  in a computational fluid dynamics context. The lift and drag
  coefficients of a NACA 0012 airfoil are considered as model
  functions.
  Numerical simulations are carried out with different input parameters for
  quantities that describe the geometry and the physical conditions of
  the problem. The evolution of the model is protracted until a periodic
  regime is reached.  Once the simulation data have been
  collected, sensitivity analysis is performed searching for an active
  subspace and response surfaces with GPR are then built from the
  application of AS and KAS techniques.
  
  The fluid motion is modelled through the unsteady
  incompressible Navier-Stokes equations approximated through the
  Chorin-Temam operator-splitting method implemented in HopeFOAM <cit.>.
  HopeFOAM is an extension of OpenFOAM <cit.>, an open source
  software for the solution of complex fluid flows problems, to variable
  higher order  element method and it adopts a Discontinuous Galerkin Method,
  based on  the formulation proposed by Hesthaven and
  Warburton <cit.>.
  
  The Discontinuous Galerkin method (DG) is a high-order method, which
  has appealing features such as the low artificial viscosity and a
  convergence rate which is optimal also on unstructured grids, commonly
  used in industrial frameworks.
  In addition to this, DG is naturally suited for the solution of problems
  described by conservative governing equations (Navier Stokes equations,
  Maxwell's equations and so on) and for parallel computing.
  All these properties are due to the fact that, differently from
  formulations based on standard finite elements, no continuity is imposed
  on the  cell boundaries and neighboring elements only exchange a common
  flux. The major drawback of DG is its high computational cost with respect
  to continuous Galerkin methods, due to the need of evaluating fluxes
  during each time step and the presence of extra degree of freedoms in
  correspondence of the elemental edges.
  
  Nowadays efforts are aimed at applying the DG in problems which involve
  deformable domains <cit.> and at improving the
  computational efficiency of the DG adopting techniques based on
  hybridization methods, matrix-free implementations, and massive
  parallelization <cit.>.
  
  
  
  
   §.§ Domain and mesh description
  
  
  The domain Ω of the fluid dynamic simulation is a
  two-dimensional duct with a sudden area expansion and a NACA 0012
  airfoil is placed in the largest section.
  The inflow ∂Ω_I is placed at the beginning of the
  narrowest part of  the duct, and here the fluid velocity is set constant
  along all the inlet boundary. The outlet is placed on the
  right hand side and it is denoted with  ∂Ω_O.
  We refer with
  ∂Ω_W := ∂Ω\{∂Ω_O∪∂Ω_I}
  to the boundaries of the airfoil and to the walls of the duct, where
  no slip boundary conditions are applied.
  The horizontal lengths of the sections of the channels are 0.6 m and
  1.35 m, respectively. The vertical length of the duct after the area
  expansion is 0.4 m, while the width of the first one depends on two
  distinct parameters.
  The airfoil has a chord-length equal to 0.1 m but its position with
  respect to the duct and its angle of attack are described by
  geometric parameters. Further details about the geometric
  parameterization of the geometry are provided in the following
  section. A proper triangulation is designed with the aid of the
  gmsh <cit.> tool and the domain is discretized with 4445
  unstructured elements.
  
  The evaluation of adimensional
  magnitudes, commonly used for characterizing the fluid flow field,
  requires the definition of some reference magnitudes. For
  the problem at hand we consider the equivalent diameter of the
  channel in correspondence of the inlet as the reference lengthscale,
  while the reference velocity is the one imposed at the inlet.
  
  
  
  
   §.§ Parameter space description
  
  We chose 7 heterogeneous parameters for the model: 2 physical, and
  5 geometrical which describe the width of the channel and the
  position of the airfoil. In <ref> are reported the
  ranges for the geometrical and physical parameters of the
  simulation. U is the first component of the initial velocity, ν
  is the kinematic viscosity, x_0 and y_0 are the horizontal and
  vertical components of the translation of the airfoil with respect to its reference position (see <ref>), α is the
  angle of the counterclockwise rotation and the center of rotation is
  located right in the middle of the airfoil, y^+ and y^- are the
  module of the vertical displacements of the upper and lower side of
  the initial conduct from a prescribed position.
  
  
  In <ref> are reported different configurations
  of the domain for the minimum and maximum values of the parameters
  α, x_0, y_0, and the minimum opening of the
  channel.
  
  
  We have considered only the counterclockwise rotation of the airfoil
  for symmetrical reasons. The range of the Reynolds number varies from
  400 to 2000, still under the regime of laminar flow.
  
  
  
  
  
   §.§ Governing equations
  
  The CFD problem is modeled through the incompressible Navier-Stokes
  and the open source solver HopeFOAM <cit.> has been employed
  for solving this set of equations <cit.>.
  
  Let Ω⊂^2 be the two-dimensional domain introduced in
  <ref>, and let us consider the incompressible
  Navier-Stokes equations. Omitting the dependence on
  (, t) ∈Ω×^+ in the first two
  equations for the sake of compactness, the governing equations are
  
      ∂_t   + (·∇)  =-∇ p+νΔ               ∈Ω, 
      ∇·  = 0     ∈Ω, 
       ( ,0) =  _0,    p( ,0)=0     ∈Ω, 
       ( ,t) =  _0,   𝐧·∇ p( ,t)=0
                 ∈∂Ω_I, 
       ( ,t) = 0,   𝐧·∇ p( ,t)=0    ∈∂Ω_W, 
      𝐧·∇ ( ,t)=0,    p( ,t)=1    ∈∂Ω_O,
  
  where p is the scalar pressure field, =(u,v) is the velocity
  field, ν is the viscosity constant and _0 is the initial
  velocity. In conservative form, the previous equations can be rewritten
  as
  
      ∂_t 𝐮+∇·ℱ=-∇ p + νΔ𝐮,
      ∇·𝐮=0,
  
  with the flux ℱ given by
  
      ℱ=[𝐅_1,𝐅_2]=[[ u^2  uv;  uv v^2 ]].
  
  From now on, in order to have a more compact notation, the
   advection term is  written as 𝒩(𝐮)=∇·ℱ(𝐮).
  
  For each timestep the procedure is broken into three stages
  accordingly to the algorithm proposed by Chorin and adapted for a DG
  framework by Hesthaven et al. <cit.>: the solution
  of the advection dominated conservation law component, the pressure
  correction weak divergence-free velocity projection, and the viscosity
  update. The non-linear advection term is treated explicitly in time
  through a second order Adams-Bashforth method <cit.>, while the diffusion
  term implicitly. The Chorin algorithm is reported in
  Algorithm <ref>.
  
  In order to recover the Discontinuos Galerkin formulation, the
  equations introduced by the Chorin method are projected onto the
  solution space by introducing a proper set of test functions and then
  the variables are approximated over each element as a linear
  combination of local shape functions.
  The DG does not impose the continuity of the solution between
  neighboring elements and therefore it requires the adoption of methods
  for the evaluation of the flux exchange between neighboring
  elements. In the present work the convective fluxes are treated
  accordingly to the Lax-Friedrichs scheme, while the viscous
  ones are solved through the Interior Penalty
  method <cit.>.
  
  
  
  
  The aerodynamic quantities we are interested in are the lift and drag
  coefficients in the incompressible case computed from the quantities
  𝐮, p, ν, A_ref, and 𝐮_0
  with a contour integral along the airfoil Γ as
  
      f=
          ∮_Γp𝐧 -
          ν(∇𝐮+∇𝐮^T)𝐧 d𝐬.
  
  The vector 𝐧 is
  the outward normal along the airfoil surface. The circulation in
  Γ is affected by both the pressure and stress distributions
  around the airfoil. The
  projection of the force along the horizontal and vertical directions
  gives the drag and lift coefficients respectively
  
      C_D=f·𝐞_1/1/2|𝐮_0|^2A_ref,
  
  
      C_L=f·𝐞_2/1/2|𝐮_0|^2A_ref,
  
  where the reference area A_ref is the chord of the
  airfoil times a length of 1 m. For the aerodynamic analysis of the fluid flow past an airfoil see <cit.>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
            
  
  
  
  
  
   §.§ Numerical results
  
  In this section a brief review of the procedure and some details about
  the numerical method and the computational domain will be presented
  along the results obtained. For what concerns the DG the polynomial
  order chosen is 3. The total number of degrees of freedom is
  133350. Small variations on the mesh are present in each of the 285
  simulations due to the different configurations of the domain.
  Each simulation is carried out until a periodic behaviour is reached
  and for this reason the final times range between 3.5 and 5 s,
  depending on the specific configuration.
  The integration time intervals are variable and they are updated at
  the end of each step in order to satisfy the CFL condition.
  The 7 physical and geometrical parameters of the simulation are
  sampled uniformly from the intervals in <ref>. In
  total we consider a dataset of 285 samples.
  
  With the purpose of qualitatively visualizing the results, 4 different
  simulations are reported in <ref> for the
  module of the velocity field and the scalar
  pressure field, respectively, both evaluated at the last time instant.
  These simulations were chosen from the 285 collected in order to show
  significant differences in the evolution of the fluid flow. In
  <ref> are reported the corresponding
  parameters. Depending on the position of the airfoil and the other
  physical parameters, different fluid flow patterns can be qualitatively
  observed.
  
  
  
  
  
  
  The lift (C_L) and drag (C_D) coefficients are evaluated when
  stationary or periodic regimes are reached, starting from the values
  of pressure and viscous stresses evaluated on the nodes close to the
  airfoil. After this sensitivity analysis is carried out. First the AS
  method is applied. The gradients necessary for the application of the
  AS method are obtained from the Gaussian process regression of the
  model functions C_L and C_D on the whole parameters'
  domain. The eigenvalues of the uncentered covariance matrix for the
  lift and drag coefficients suggest the presence of a one-dimensional
  active subspace in both cases.
  
  The plots of the first active eigenvector components are useful as
  sensitivity measures, see <ref>. The greater
  the absolute value of a component is, the greater is its influence on
  the model function. We observe that the lift coefficient is influenced
  mainly by the vertical position of the airfoil and the angle of
  attack, while the drag coefficient depends mainly on the initial
  velocity, and secondarily on the viscosity and on the angle of
  attack.
  
  As one could expect from physical considerations, the angle of attack
  affects both drag and lift coefficients, while the viscosity, which
  governs the wall stresses, is relevant for the evaluation of the C_D.
  The vertical position of the airfoil with respect to the symmetric
  axis of the section of the duct after the area expansion also greatly
  affects both coefficients, and this is mainly due to the fact that the
  fluid flow conditions change drastically between the core, where the
  speed is higher, and the one close to the wall of the duct, where the
  speed tends to zero. On the other hand, the horizontal translation has
  almost no impact on the results, given the regularity of the fluid
  flow along the x-axis for the considered range of x_0.
  Moreover, the non-symmetric behaviour of the upper and lower
  parameters which determine the opening of the channel is due to the
  non-symmetric choice of the range considered for the angle of attack.
  
  
  
  The KAS method was applied with 1500 features. In order to compare
  the AS and KAS methods 5-fold cross validation was implemented. The
  score of cross validation is the relative root mean square error
  (RRMSE) defined in <ref>.
  
  
  
  
  
  The Gaussian process regressions for the two methods are shown in
  <ref> for the lift coefficient, and in
  <ref> for the drag coefficient. They were obtained as
  a single step of 5-fold cross validation with one fifth of the 285
  samples used as test set. The spectral distribution of the feature map
  is the Gaussian distribution for the lift, and the Beta for the drag, respectively.
  The RRMSE mean and standard deviation from 5-fold cross validation,
  are reported for different active dimensions in
  <ref>. The feature map from Equation (<ref>) was
  adopted. The hyperparameters of the spectral distributions were tuned
  with logarithmic grid-search with 5-fold cross validation as
  described in Algorithm <ref>.
  
  Regarding the drag coefficient, the relative gain using the KAS method reaches
  the 19.2% on average when employing the Beta spectral measure for the
  definition of the feature map. The relative gain of the one
  dimensional response surface built with GPR from the KAS method is
  7% on average for the lift coefficient. This result could be
  due to the higher noise in the evaluation of the C_L. In this case
  the relative gain increases when the dimension of the
  response surface increases to 2 with a gain of 14.6%. A
  slight reduction of the AS RRMSE relative to the drag
  coefficient is ascertained when increasing the dimension of the response surface.
  
  
  
  
  
  § CONCLUSIONS AND PERSPECTIVES
  
  
  In this work we presented a new nonlinear extension of the active
  subspaces property that introduces Kernel-based Active Subspaces (KAS). The
  method exploits random Fourier features to find active
  subspaces on high-dimensional feature spaces. We
  tested the new method over 5 different benchmarks of increasing
  complexity, and we provided pseudo-codes for every aspects of the
  proposed kernel-extension. The
  tested model functions range from scalar to vector-valued. We also
  provide a CFD application discretized by the Discontinuous Galerkin
  method. We compared the kernel-based active subspaces to the standard
  linear active subspaces and we observed in all the cases an increment
  of the accuracy of the Gaussian response surfaces built over the reduced
  parameter spaces.
  The most interesting results regard the possibility to apply the KAS method
  when an active subspace does not exist. This was shown
  for radial symmetric model functions.
  
  Future developments will involve the study of more efficient
  procedures for tuning the hyperparameters of the spectral
  distribution. Other possible advances could be done finding an
  effective back-mapping from the targets to the actual parameters in
  the full original space. This could promote the implementation of
  optimization algorithms or other parameter studies enhanced by the
  kernel-based active subspaces extension.
  
  
  
  § APPENDIX - PROOF DETAILS
  
  
  
  In this section we provide an expanded version of the proof of <ref>.
  
  
  The proof is remodeled from <cit.>, and it is developed in five steps:
  
  
    * Since R_V∈ℳ(d,d) is symmetric positive definite
    there exists a basis of eigenvectors
    (𝐰_i)_i∈{1,…, d} and a corresponding set of positive
    eigenvalues (β_i)_i∈{1,…, d} such that
  
  
      R_V=∑^d_i=1β_i 𝐰_i⊗𝐰_i.
  
  
    * Let us define the ridge approximation error as
  
      𝐞=‖𝐟-𝐡∘
      P_r‖_L^2(ℝ^m, ℬ(ℝ^m),
        ρ;V) = 𝔼_P[ ‖(
        𝐟(𝐗)-𝐡(P_r(𝐗))‖^2_R_V].
  
  Then we can decompose the error analysis for each component employing
  the spectral decomposition (<ref>)
  
      𝔼_P[
      ‖𝐞(𝐗)‖^2_R_V]    = 𝔼_P[tr((R_V𝐞(𝐗)) ⊗𝐞 (𝐗)) ]=
       
         =∑^d_i=1β_i 𝔼_P[tr(((𝐰_i⊗𝐰_i)𝐞(𝐗))⊗𝐞(𝐗))]=
        
         =∑^d_i=1β_i 𝔼_P[(𝐰_i·𝐞
      (𝐗)) tr(𝐰_i⊗𝐞(𝐗))]=
        
         =∑^d_i=1β_i 𝔼_P[(𝐰_i·𝐞(𝐗))^2],
  
  so we can define
  e_i(𝐗)=𝐰_i·𝐞(𝐗)=f_i(𝐗)-h_i(P_r(𝐗)), ∀
  i∈{1, …, d} and treat each component separately.
  
    * The next step involves the application of
    Lemma <ref> to the scalar functions
    f_i(𝐗)-h_i(P_r(𝐗)), ∀ i∈{1,
    …, d}
  
      𝔼_P[(f_i(𝐗)-h_i(P_r()))^2]
            =
          𝔼_P[𝔼_P[(f_i
          ()-h_i(P_r ()))^2|σ(P_r)]]
        
         ≤𝔼_P[C_p(ρ, P_r(𝐗))
                       𝔼_P[‖(I-P_r^T)∇
                       f_i(𝐗)‖^2_2|σ(P_r)]]
        
         ≤𝔼_P[C_p(ρ,
                       P_r(𝐗))]^1/p 𝔼_P[‖(I - P_r^T)∇
                       f_i(𝐗)‖^2_2]^1/q ,
  
  where we used the Hölder inequality with indexes (p, q)=(∞,
  1) when ρ belongs to the first and second classes of
  Assumption <ref>, and (p, q)=(τ+1/τ, 1+τ)
  when ρ belongs to the third class.
  
  Then we can bound
  𝔼_P[(C_p,
    P_r(𝐗))]^1/p with a constant
  C(C_p(ρ, P_r(𝐗))) which depends on the class of
  ρ (see Lemma 3.1, Lemma 4.2, Lemma 4.3, Lemma 4.4 and Theorem 4.5
  of <cit.>) as follows
  
      𝔼_P[C_p(ρ,
                      P_r(𝐗))]^1/p 𝔼_P[‖(I - P_r^T) ∇
                      f_i(𝐗)‖^2_2]^1/q≤
         ≤ C(C_p(ρ, P_r(𝐗))) tr(𝔼_P[(I - P_r^T)∇ f_i ()(∇ f_i())^T (I -P_r) ]^1/q)= 
         = C(C_p(ρ, P_r(𝐗))) tr((I - P_r^T)
      𝔼_P[∇ f_i()(∇ f_i())^T](I - P_r))^1/q.
  
  
  
    * The spectral decomposition (<ref>) is
    employed again and the covariance matrix H is introduced in the
    last equation
  
      𝔼_P[
      ‖𝐞(𝐗)‖^2_R_V]≤
         ≤∑^d_i=1β_i C(C_p(ρ, P_r(𝐗))) tr((I - P_r^T)𝔼_P[((∇𝐟(𝐗))^T𝐰_i)⊗((∇𝐟(𝐗))^T𝐰_i)](I - P_r))^1/q= 
         =C(C_p(ρ, P_r(𝐗))) tr((I - P_r^T)𝔼_P[(∇𝐟(𝐗))^T(∑^d_i=1β_i^q 𝐰_i⊗𝐰_i)∇𝐟(𝐗)](I - P_r))^1/q= 
         =C(C_p(ρ, P_r(𝐗))) tr((I - P_r^T)𝔼_P[(∇𝐟(𝐗))^TR_V(ρ)∇𝐟(𝐗)](I - P_r))^1/q= 
         =C(C_p(ρ, P_r(𝐗))) tr((I - P_r^T)H(I - P_r))^1/q,
  
  where R_V(ρ) is the original metric matrix if ρ belongs to the first or second class of Assumption <ref> and is equal to
  
      ∑^d_i=1β_i^1+τ 𝐰_i⊗𝐰_i ,
  
  if ρ belongs to the third class.
  
    * Finally the bound in the statement of the theorem is recovered
    solving the following minimization problem with classical model
    reduction arguments employing singular value decomposition (SVD)
  
      P_r = _P_r∈𝒪(m, m)tr((I - P_r^T)H(I - P_r)).
  
  
  
  
  
  
  § ACKNOWLEDGEMENTS
  
  This work was partially supported by an industrial Ph.D. grant sponsored
  by Fincantieri S.p.A. (IRONTH Project), by MIUR (Italian ministry for
  university and research)
  through FARE-X-AROMA-CFD project, and partially funded by European
  Union Funding for Research and Innovation — Horizon 2020 Program — in
  the framework of European Research Council Executive Agency: H2020 ERC
  CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methods with
  Applications in Computational Fluid Dynamics” P.I. Professor Gianluigi Rozza.
  
  abbrvurl
  
  
  10
  
  OpenFOAM
  <https://openfoam.org/>.
  
  gmsh
  Gmesh, a three-dimensional finite element mesh generator with built-in pre- and
    post-processing facilities.
  <http://gmsh.info/>.
  
  gpyopt2016
  GPyOpt: A bayesian optimization framework in python.
  <http://github.com/SheffieldML/GPyOpt>, 2016.
  
  HopeFOAM
  Hopefoam extension of OpenFOAM.
  <https://github.com/HopeFOAM/HopeFOAM>, Last release 15th
    September 2017.
  
  aguiar2018dynamic
  I. P. Aguiar.
  Dynamic Active Subspaces: a Data-Driven Approach to Computing
    Time-Dependent Active Subspaces in Dynamical Systems.
  Master's thesis, University of Colorado Boulder, 2018.
  
  arnold1982interior
  D. N. Arnold.
  An interior penalty finite element method with discontinuous
    elements.
  SIAM journal on numerical analysis, 19(4):742–760, 1982.
  https://doi.org/10.1137/0719052 doi:10.1137/0719052.
  
  barshan2011supervised
  E. Barshan, A. Ghodsi, Z. Azimifar, and M. Z. Jahromi.
  Supervised principal component analysis: Visualization,
    classification and regression on subspaces and submanifolds.
  Pattern Recognition, 44(7):1357–1371, 2011.
  
  berlinet2011reproducing
  A. Berlinet and C. Thomas-Agnan.
  Reproducing kernel Hilbert spaces in probability and
    statistics.
  Springer Science & Business Media, 2011.
  
  bobrowski2005functional
  A. Bobrowski.
  Functional analysis for probability and stochastic processes: an
    introduction.
  Cambridge University Press, 2005.
  
  Bridges2019ActiveMA
  R. A. Bridges, A. D. Gruber, C. Felder, M. E. Verma, and C. Hoff.
  Active Manifolds: A non-linear analogue to Active Subspaces.
  In Proceedings of International Conference on Machine Learning,
    2019.
  
  brunton2019data
  S. L. Brunton and J. N. Kutz.
  Data-driven science and engineering: Machine learning,
    dynamical systems, and control.
  Cambridge University Press, 2019.
  
  constantine2015active
  P. G. Constantine.
  Active subspaces: Emerging ideas for dimension reduction in
    parameter studies, volume 2 of SIAM Spotlights.
  SIAM, 2015.
  
  constantine2017global
  P. G. Constantine and P. Diaz.
  Global sensitivity metrics from active subspaces.
  Reliability Engineering & System Safety, 162:1–13, 2017.
  
  constantine2014active
  P. G. Constantine, E. Dow, and Q. Wang.
  Active subspace methods in theory and practice: applications to
    kriging surfaces.
  SIAM Journal on Scientific Computing, 36(4):A1500–A1524, 2014.
  
  constantine2015exploiting
  P. G. Constantine, M. Emory, J. Larsson, and G. Iaccarino.
  Exploiting active subspaces to quantify uncertainty in the numerical
    simulation of the HyShot II scramjet.
  Journal of Computational Physics, 302:1–20, 2015.
  
  cook2005sufficient
  R. D. Cook and L. Ni.
  Sufficient dimension reduction via inverse regression: A minimum
    discrepancy approach.
  Journal of the American Statistical Association,
    100(470):410–428, 2005.
  
  cui2020active
  C. Cui, K. Zhang, T. Daulbaev, J. Gusak, I. Oseledets, and Z. Zhang.
  Active subspace of neural networks: Structural analysis and universal
    attacks.
  SIAM Journal on Mathematics of Data Science, 2(4):1096–1122,
    2020.
  https://doi.org/10.1137/19M1296070
    doi:10.1137/19M1296070.
  
  demo2021hull
  N. Demo, M. Tezzele, A. Mola, and G. Rozza.
  Hull Shape Design Optimization with Parameter Space and Model
    Reductions, and Self-Learning Mesh Morphing.
  Journal of Marine Science and Engineering, 9(2):185, 2021.
  https://doi.org/10.3390/jmse9020185
    doi:10.3390/jmse9020185.
  
  demo2019cras
  N. Demo, M. Tezzele, and G. Rozza.
  A non-intrusive approach for reconstruction of POD modal
    coefficients through active subspaces.
  Comptes Rendus Mécanique de l'Académie des Sciences,
    DataBEST 2019 Special Issue, 347(11):873–881, November 2019.
  https://doi.org/10.1016/j.crme.2019.11.012
    doi:10.1016/j.crme.2019.11.012.
  
  demo2020asga
  N. Demo, M. Tezzele, and G. Rozza.
  A Supervised Learning Approach Involving Active Subspaces for an
    Efficient Genetic Algorithm in High-Dimensional Optimization Problems.
  SIAM Journal on Scientific Computing, 43(3):B831–B853, 2021.
  https://doi.org/10.1137/20M1345219
    doi:10.1137/20M1345219.
  
  diaz2018modified
  P. Diaz, P. Constantine, K. Kalmbach, E. Jones, and S. Pankavich.
  A modified SEIR model for the spread of Ebola in Western Africa and
    metrics for resource allocation.
  Applied Mathematics and Computation, 324:141–155, 2018.
  
  gazdag1976time
  J. Gazdag.
  Time-differencing schemes and transform methods.
  Journal of Computational Physics, 20(2):196–207, 1976.
  
  ghoreishi2019adaptive
  S. F. Ghoreishi, S. Friedman, and D. L. Allaire.
  Adaptive Dimensionality Reduction for Fast Sequential Optimization
    With Gaussian Processes.
  Journal of Mechanical Design, 141(7):071404, 2019.
  
  gpy2014
  GPy.
  GPy: A Gaussian process framework in Python.
  <http://github.com/SheffieldML/GPy>, since 2012.
  
  guo2018reduced
  M. Guo and J. S. Hesthaven.
  Reduced order modeling for nonlinear structural analysis using
    Gaussian process regression.
  Computer Methods in Applied Mechanics and Engineering,
    341:807–826, 2018.
  
  heas2020generalized
  P. Héas, C. Herzet, and B. Combes.
  Generalized kernel-based dynamic mode decomposition.
  In ICASSP 2020-2020 IEEE International Conference on Acoustics,
    Speech and Signal Processing (ICASSP), pages 3877–3881. IEEE, 2020.
  
  hesthaven2016certified
  J. S. Hesthaven, G. Rozza, and B. Stamm.
  Certified Reduced Basis Methods for Parametrized Partial
    Differential Equations.
  Springer, 2016.
  
  hesthaven2007nodal
  J. S. Hesthaven and T. Warburton.
  Nodal discontinuous Galerkin methods: algorithms, analysis, and
    applications.
  Springer Science & Business Media, 2007.
  
  jefferson2015active
  J. L. Jefferson, J. M. Gilbert, P. G. Constantine, and R. M. Maxwell.
  Active subspaces for sensitivity analysis and dimension reduction of
    an integrated hydrologic model.
  Computers & geosciences, 83:127–138, 2015.
  
  ji2019quantifying
  W. Ji, Z. Ren, Y. Marzouk, and C. K. Law.
  Quantifying kinetic uncertainty in turbulent combustion simulations
    using active subspaces.
  Proceedings of the Combustion Institute, 37(2):2175–2182,
    2019.
  https://doi.org/10.1016/j.proci.2018.06.206
    doi:10.1016/j.proci.2018.06.206.
  
  ji2018shared
  W. Ji, J. Wang, O. Zahm, Y. M. Marzouk, B. Yang, Z. Ren, and C. K. Law.
  Shared low-dimensional subspaces for propagating kinetic uncertainty
    to multiple outputs.
  Combustion and Flame, 190:146–157, 2018.
  
  kevrekidis2015kernel
  I. Kevrekidis, C. Rowley, and M. Williams.
  A kernel-based method for data-driven Koopman spectral analysis.
  Journal of Computational Dynamics, 2(2):247–265, 2015.
  
  kundu2012fluid
  P. K. Kundu, I. M. Cohen, and D. R. Dowling, editors.
  Fluid Mechanics (Fifth Edition).
  Academic Press, Boston, 2012.
  https://doi.org/https://doi.org/10.1016/C2009-0-63410-3
    doi:https://doi.org/10.1016/C2009-0-63410-3.
  
  li1991sliced
  K.-C. Li.
  Sliced inverse regression for dimension reduction.
  Journal of the American Statistical Association,
    86(414):316–327, 1991.
  
  li2007sparse
  L. Li.
  Sparse sufficient dimension reduction.
  Biometrika, 94(3):603–613, 2007.
  
  li2019a
  Z. Li, J.-F. Ton, D. Oglic, and D. Sejdinovic.
  Towards a unified analysis of random Fourier features.
  In Proceedings of the 36th International Conference on Machine
    Learning, volume 97 of Proceedings of Machine Learning Research, pages
    3905–3914. ICML, 2019.
  
  lukaczyk2015surrogate
  T. W. Lukaczyk.
  Surrogate Modeling and Active Subspaces for Efficient
    Optimization of Supersonic Aircraft.
  PhD thesis, Stanford University, 2015.
  
  lukaczyk2014active
  T. W. Lukaczyk, P. Constantine, F. Palacios, and J. J. Alonso.
  Active subspaces for shape optimization.
  In 10th AIAA multidisciplinary design optimization conference,
    page 1171, 2014.
  
  meneghetti2021dimensionality
  L. Meneghetti, N. Demo, and G. Rozza.
  A dimensionality reduction approach for convolutional neural
    networks.
  arXiv preprint arXiv:2110.09163, Submitted, 2021.
  
  mika1999kernel
  S. Mika, B. Schölkopf, A. J. Smola, K.-R. Müller, M. Scholz, and
    G. Rätsch.
  Kernel PCA and de-noising in feature spaces.
  In Advances in neural information processing systems, pages
    536–542, 1999.
  
  mohri2018foundations
  M. Mohri, A. Rostamizadeh, and A. Talwalkar.
  Foundations of machine learning.
  MIT press, 2018.
  
  nguyen2009implicit
  N. C. Nguyen, J. Peraire, and B. Cockburn.
  An implicit high-order hybridizable discontinuous galerkin method for
    linear convection–diffusion equations.
  Journal of Computational Physics, 228(9):3232–3254, 2009.
  
  palaci2018gaussian
  M. Palaci-Olgun.
  Gaussian process modeling and supervised dimensionality reduction
    algorithms via stiefel manifold learning.
  Master's thesis, University of Toronto Institute of Aerospace
    Studies, 2018.
  
  parente2020generalized
  M. T. Parente, J. Wallin, B. Wohlmuth, et al.
  Generalized bounds for active subspaces.
  Electronic Journal of Statistics, 14(1):917–943, 2020.
  
  pazner2017stage
  W. Pazner and P.-O. Persson.
  Stage-parallel fully implicit Runge–Kutta solvers for discontinuous
    Galerkin fluid simulations.
  Journal of Computational Physics, 335:700–717, 2017.
  https://doi.org/10.1016/j.jcp.2017.01.050
    doi:10.1016/j.jcp.2017.01.050.
  
  quarteroni2014reduced
  A. Quarteroni and G. Rozza.
  Reduced order methods for modeling and computational reduction,
    volume 9.
  Springer, 2014.
  https://doi.org/10.1007/978-3-319-02090-7
    doi:10.1007/978-3-319-02090-7.
  
  rahimi2008random
  A. Rahimi and B. Recht.
  Random features for large-scale kernel machines.
  In Advances in neural information processing systems, pages
    1177–1184, 2008.
  
  romor2021multi
  F. Romor, M. Tezzele, M. Mrosek, C. Othmer, and G. Rozza.
  Multi-fidelity data fusion through parameter space reduction with
    applications to automotive engineering.
  arXiv preprint arXiv:2110.14396, Submitted, 2021.
  
  romor2020athena
  F. Romor, M. Tezzele, and G. Rozza.
  ATHENA: Advanced Techniques for High dimensional parameter spaces to
    Enhance Numerical Analysis.
  Software Impacts, 10:100133, 2021.
  https://doi.org/10.1016/j.simpa.2021.100133
    doi:10.1016/j.simpa.2021.100133.
  
  romor2021las
  F. Romor, M. Tezzele, and G. Rozza.
  A local approach to parameter space reduction for regression and
    classification tasks.
  arXiv preprint arXiv:2107.10867, Submitted, 2021.
  
  morhandbook2020
  G. Rozza, M. Hess, G. Stabile, M. Tezzele, and F. Ballarin.
  Basic Ideas and Tools for Projection-Based Model Reduction of
    Parametric Partial Differential Equations.
  In P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, W. H. A.
    Schilders, and L. M. Silveira, editors, Model Order Reduction,
    volume 2, chapter 1, pages 1–47. De Gruyter, Berlin, Boston, 2020.
  https://doi.org/10.1515/9783110671490-001
    doi:10.1515/9783110671490-001.
  
  rozza2018advances
  G. Rozza, M. H. Malik, N. Demo, M. Tezzele, M. Girfoglio, G. Stabile, and
    A. Mola.
  Advances in Reduced Order Methods for Parametric Industrial Problems
    in Computational Fluid Dynamics.
  In R. Owen, R. de Borst, J. Reese, and P. Chris, editors, ECCOMAS ECFD 7 - Proceedings of 6th European Conference on Computational
    Mechanics (ECCM 6) and 7th European Conference on Computational Fluid
    Dynamics (ECFD 7), pages 59–76, Glasgow, UK, 2018.
  
  russi2010uncertainty
  T. M. Russi.
  Uncertainty quantification with experimental data and complex
    system models.
  PhD thesis, UC Berkeley, 2010.
  
  salmoiraghi2016advances
  F. Salmoiraghi, F. Ballarin, G. Corsi, A. Mola, M. Tezzele, and G. Rozza.
  Advances in geometrical parametrization and reduced order models and
    methods for computational fluid dynamics problems in applied sciences and
    engineering: Overview and perspectives.
  ECCOMAS Congress 2016 - Proceedings of the 7th European Congress
    on Computational Methods in Applied Sciences and Engineering, 1:1013–1031,
    2016.
  https://doi.org/10.7712/100016.1867.8680
    doi:10.7712/100016.1867.8680.
  
  scholkopf1998nonlinear
  B. Schölkopf, A. Smola, and K.-R. Müller.
  Nonlinear component analysis as a kernel eigenvalue problem.
  Neural computation, 10(5):1299–1319, 1998.
  
  scholkopf2002learning
  B. Schölkopf, A. J. Smola, F. Bach, et al.
  Learning with kernels: support vector machines, regularization,
    optimization, and beyond.
  MIT press, 2002.
  
  seshadri2018turbomachinery
  P. Seshadri, S. Shahpar, P. Constantine, G. Parks, and M. Adams.
  Turbomachinery active subspace performance maps.
  Journal of Turbomachinery, 140(4):041003, 2018.
  https://doi.org/10.1115/1.4038839
    doi:10.1115/1.4038839.
  
  shahbazi2005explicit
  K. Shahbazi.
  An explicit expression for the penalty parameter of the interior
    penalty method.
  Journal of Computational Physics, 205(2):401–407, 2005.
  https://doi.org/10.1016/j.jcp.2004.11.017
    doi:10.1016/j.jcp.2004.11.017.
  
  sriperumbudur2017approximate
  B. Sriperumbudur and N. Sterge.
  Approximate kernel PCA using random features: computational vs.
    statistical trade-off.
  arXiv preprint arXiv:1706.06296, 2017.
  
  tezzele2018combined
  M. Tezzele, F. Ballarin, and G. Rozza.
  Combined parameter and model reduction of cardiovascular problems by
    means of active subspaces and POD-Galerkin methods.
  In D. Boffi, L. F. Pavarino, G. Rozza, S. Scacchi, and C. Vergara,
    editors, Mathematical and Numerical Modeling of the Cardiovascular
    System and Applications, volume 16 of SEMA-SIMAI Series, pages
    185–207. Springer International Publishing, 2018.
  https://doi.org/10.1007/978-3-319-96649-6_8
    doi:10.1007/978-3-319-96649-6_8.
  
  tezzele2018model
  M. Tezzele, N. Demo, M. Gadalla, A. Mola, and G. Rozza.
  Model order reduction by means of active subspaces and dynamic mode
    decomposition for parametric hull shape design hydrodynamics.
  In Technology and Science for the Ships of the Future:
    Proceedings of NAV 2018: 19th International Conference on Ship & Maritime
    Research, pages 569–576. IOS Press, 2018.
  https://doi.org/10.3233/978-1-61499-870-9-569
    doi:10.3233/978-1-61499-870-9-569.
  
  tezzele2018ecmi
  M. Tezzele, N. Demo, A. Mola, and G. Rozza.
  An integrated data-driven computational pipeline with model order
    reduction for industrial and applied mathematics.
  In M. Günther and W. Schilders, editors, Novel Mathematics
    Inspired by Industrial Challenges, number 38 in Mathematics in Industry.
    Springer International Publishing, 2022.
  https://doi.org/10.1007/978-3-030-96173-2_7
    doi:10.1007/978-3-030-96173-2_7.
  
  tezzele2019marine
  M. Tezzele, N. Demo, and G. Rozza.
  Shape optimization through proper orthogonal decomposition with
    interpolation and dynamic mode decomposition enhanced by active subspaces.
  In R. Bensow and J. Ringsberg, editors, Proceedings of MARINE
    2019: VIII International Conference on Computational Methods in Marine
    Engineering, pages 122–133, 2019.
  
  tezzele2020enhancing
  M. Tezzele, N. Demo, G. Stabile, A. Mola, and G. Rozza.
  Enhancing CFD predictions in shape design problems by model and
    parameter space reduction.
  Advanced Modeling and Simulation in Engineering Sciences,
    7(40), 2020.
  https://doi.org/10.1186/s40323-020-00177-y
    doi:10.1186/s40323-020-00177-y.
  
  tezzele2021multi
  M. Tezzele, L. Fabris, M. Sidari, M. Sicchiero, and G. Rozza.
  A multi-fidelity approach coupling parameter space reduction and
    non-intrusive POD with application to structural optimization of passenger
    ship hulls.
  arXiv preprint arXiv:2206.01243, Submitted, 2022.
  
  tezzele2018dimension
  M. Tezzele, F. Salmoiraghi, A. Mola, and G. Rozza.
  Dimension reduction in heterogeneous parametric spaces with
    application to naval engineering shape design problems.
  Advanced Modeling and Simulation in Engineering Sciences,
    5(1):25, Sep 2018.
  https://doi.org/10.1186/s40323-018-0118-3
    doi:10.1186/s40323-018-0118-3.
  
  tripathy2019deep
  R. Tripathy and I. Bilionis.
  Deep active subspaces: A scalable method for high-dimensional
    uncertainty propagation.
  In ASME 2019 International Design Engineering Technical
    Conferences and Computers and Information in Engineering Conference.
    American Society of Mechanical Engineers Digital Collection, 2019.
  
  tripathy2016gaussian
  R. Tripathy, I. Bilionis, and M. Gonzalez.
  Gaussian processes with built-in dimensionality reduction:
    Applications to high-dimensional uncertainty propagation.
  Journal of Computational Physics, 321:191–223, 2016.
  https://doi.org/10.1016/j.jcp.2016.05.039
    doi:10.1016/j.jcp.2016.05.039.
  
  2020SciPy-NMeth
  P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy,
    D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright,
    S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman,
    N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson,
    C. Carey, İ. Polat, Y. Feng, E. W. Moore, J. Vand erPlas,
    D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero,
    C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van
    Mulbregt, and S. . . Contributors.
  SciPy 1.0: Fundamental Algorithms for Scientific Computing in
    Python.
  Nature Methods, 2020.
  https://doi.org/https://doi.org/10.1038/s41592-019-0686-2
    doi:https://doi.org/10.1038/s41592-019-0686-2.
  
  vohra2019active
  M. Vohra, A. Alexanderian, H. Guy, and S. Mahadevan.
  Active subspace-based dimension reduction for chemical kinetics
    applications with epistemic uncertainty.
  Combustion and Flame, 204:152–161, 2019.
  https://doi.org/10.1016/j.combustflame.2019.03.006
    doi:10.1016/j.combustflame.2019.03.006.
  
  warner1983foundations
  F. W. Warner.
  Foundations of differentiable manifolds and Lie groups,
    volume 94.
  Springer Science & Business Media, 1983.
  
  weller1998
  H. G. Weller, G. Tabor, H. Jasak, and C. Fureby.
  A tensorial approach to computational continuum mechanics using
    object-oriented techniques.
  Computers in physics, 12(6):620–631, 1998.
  https://doi.org/10.1063/1.168744
    doi:10.1063/1.168744.
  
  williams2006gaussian
  C. K. Williams and C. E. Rasmussen.
  Gaussian Processes for Machine Learning.
  Adaptive Computation and Machine Learning series. MIT press
    Cambridge, MA, 2006.
  
  wu2009localized
  Q. Wu, S. Mukherjee, and F. Liang.
  Localized sliced inverse regression.
  In Advances in neural information processing systems, pages
    1785–1792, 2009.
  
  zahm2020gradient
  O. Zahm, P. G. Constantine, C. Prieur, and Y. M. Marzouk.
  Gradient-based dimension reduction of multivariate vector-valued
    functions.
  SIAM Journal on Scientific Computing, 42(1):A534–A558, 2020.
  https://doi.org/10.1137/18M1221837
    doi:10.1137/18M1221837.
  
  zahr2016adjoint
  M. J. Zahr and P.-O. Persson.
  An adjoint method for a high-order discretization of deforming domain
    conservation laws for optimization of flow problems.
  Journal of Computational Physics, 326:516–543, 2016.
  
  zhang2019learning
  G. Zhang, J. Zhang, and J. Hinkle.
  Learning nonlinear level sets for dimensionality reduction in
    function approximation.
  In Advances in Neural Information Processing Systems, pages
    13199–13208, 2019.
  
  
  
  
  