{'fraction_non_alphanumeric': 0.04326725905673274, 'fraction_numerical': 0.014490772385509228, 'mean_word_length': 4.488433162435913, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Grassmannian diffusion maps based surrogate modeling via geometric harmonics

In this paper, a novel surrogate model based on the Grassmannian diffusion maps (GDMaps) and utilizing geometric harmonics is developed for predicting the response of engineering systems and complex physical phenomena. The method utilizes the GDMaps to obtain a low-dimensional representation of the underlying behavior of physical/mathematical systems with respect to uncertainties in the input parameters. Using this representation, geometric harmonics, an out-of-sample function extension technique, is employed to create a global map from the space of input parameters to a Grassmannian diffusion manifold. Geometric harmonics is also employed to locally map points on the diffusion manifold onto the tangent space of a Grassmann manifold. The exponential map is then used to project the points in the tangent space onto the Grassmann manifold, where reconstruction of the full solution is performed. The performance of the proposed surrogate modeling is verified with three examples. The first problem is a toy example used to illustrate the development of the technique. In the second example, errors associated with the various mappings employed in the technique are assessed by studying response predictions of the electric potential of a dielectric cylinder in a homogeneous electric field. The last example applies the method for uncertainty prediction in the strain field evolution in a model amorphous material using the shear transformation zone (STZ) theory of plasticity. In all examples, accurate predictions are obtained, showing that the present technique is a strong candidate for the application of uncertainty quantification in large-scale models.

1. Introduction. Surrogate models (aka emulators or metamodels) have become an important tool for uncertainty quantification because they afford a computationally efficient means of approximating (often complex) input-output relations generated from high-fidelity computational models. Surrogate models are constructed by training (or learning) a mathematical model from a finite set of input-output observations referred to as the training set. Given a new set of input parameters, where the solution of the model is unknown, the surrogate can then be used to predict the solution at minimal cost. Surrogate models are typically classified as either intrusive or non-intrusive [45,60]. Intrusive methods, such as those based on Galerkin schemes, typically provide good convergence [28,3], but their complexity limits their flexibility [10] and requires the development of new numerical schemes and, hence, entirely new solvers. In contrast, non-intrusive methods are trained from realizations of a deterministic solver (i.e. they can leverage pre-existing numerical methods and solvers) at selected sample points generated from the uncertain parameters. Moreover, surrogate models often provide an improvement over statistical approaches based on Monte Carlo simulation (MCS). Although MCS is a versatile and non-intrusive method, it typically offers slow convergence with the number of samples [26]. To circumvent this limitation, quasi Monte Carlo methods [9], adaptive sampling techniques [59,58], and other intelligent sampling techniques leveraging variance reduction techniques [49,21,18,5] are employed. Increasingly, these enhanced sampling methods are being leveraged to improve the training efficiency of surrogate models; enabling surrogates to be developed from far fewer training data.
Among the most widely-used surrogate models in UQ are polynomial chaos expansions (PCE) [28] and Gaussian process (GP) regression (Kriging) [43,55]. PCE was originally proposed by Wiener [64] based on the projection of random solutions onto a basis of Hermite polynomials, which are orthogonal with respect to the Gaussian measure. Ghanem et al. [28] introduced the Stochastic Galerkin projection for PCE, an intrusive method that requires the formulation of a system of algebraic equations for distinct classes of problems. The generalized PCE (gPCE) [67] provides improved convergence and enhanced flexibility by utilizing polynomial bases from the Wiener-Askey scheme. Moreover, gPCE can be used in the construction of non-intrusive surrogate based on collocation schemes aiming at enhanced versatility and convergence properties. On the other hand, a GP approximates the input-output relation as a Gaussian stochastic process completely described by its mean and correlation function. Learning a GP thus consists of determining the mean and correlation functions. However, the computational complexity of the GP training process with n observations is on the order O(n 3 ), whereas the memory requirements scale with O(n 2 ).
It is well known that the relationship between data dimensionality and model fidelity has a strong influence on the computational performance of UQ for large-scale models. In this regard, dimension reduction techniques have become attractive due to their ability to represent high-dimensional data in a low-dimensional and more informative space (manifold) [52,23,12,20]. Dimension reduction techniques can be classified into linear and nonlinear methods. Linear methods include principal component analysis (PCA) [27], locality preserving projections (LPP) [34], linear regression [54], and singular value decomposition [62]. On the other hand, nonlinear methods are useful for constructing nonlinear maps between the high-and low-dimensional spaces. Nonlinear methods include Isomaps [63], locally linear embedding (LLE) [56,16], Kernel PCA [44], and diffusion maps (DMaps) [15,13]. It is worth noting that several methods based on DMaps have been proposed in the literature such as the method proposed by Soize and Ghanem [61] for sampling a random vector whose probability distribution is constrained to an Euclidean manifold. They used DMaps for discovering the underlying structure of the dataset. Moreover, DMaps has been used in the development of surrogate models either based on local polynomial interpolations and the Nyström out-of-sample extension [37], or based on neural networks and the Laplacian pyramids [38]. Further, a subspace extension of Diffusion Maps, the Grassmannian Diffusion Maps (GDMaps), was introduced by dos Santos et al. [19]. In this technique, a model dimension hyper-reduction is achieved by combining a pointwise linear dimensionality reduction technique, projecting the high-dimensional data onto a low-dimensional Grassmann manifold, and a multipoint nonlinear dimension reduction using Diffusion maps which reveals the intrinsic structure of the data on the manifold. Recently, Kontolati et al. [42] employed GDMaps and PCE in the development of a manifold learning based method for UQ in systems describing complex processes.
In this work, we leverage the GDMaps to construct surrogate models for very high-dimensional systems. More specifically, a set of low-dimensional coordinates embedding the high-dimensional model data on the Grassmann manifold is obtained via Grassmannian diffusion maps (GDMaps). This low-dimensional diffusion space serves as a connecting space between the parameter space and the Grassmann manifold, where full solution reconstruction can be performed. To connect these spaces, maps are constructed using the idea of out-sampling extension [8,7,48,46]. Inspired by the the Nyström extension [65,35], Coifman and Lafon [14] introduced a scheme, referred to as geometric harmonics (GH) for extending empirical functions only available at few locations. They demonstrated that this process relates the function complexity and its extension, with important implications to the construction of the lifting and restriction operators. In this paper, GH is employed in the construction of an out-of-sample extension to create maps between spaces of interest. A global GH surrogate is constructed between the parameter and diffusion spaces, while local GH surrogates are constructed between the diffusion space and the Grassmann manifold via the tangent space, a flat inner-product space allowing local exponential mapping onto the Grassmannian, where the predicted model response can be constructed. Moreover, GH is also used to create a map between the parameter space and the space of singular values used in the contruction of the model prediction.
This paper is organized as follows. Section 2 discusses important background on the Grassmann manifold. Section 3 provides an overview of the Grassmannian diffusion maps technique. Section 4 introduces GH as an out-of-sample function extension technique. Section 5 contains a detailed description of the surrogate modeling approach, referred to as Grassmannian-Geometric Harmonics Maps (Grassmannian-GHMaps), developed herein. In section 6, three examples are provided. The first is a simple example used to explain the proposed method in a manner that is conceptually understandable and easy to visualize. The second example uses the GDMaps surrogates to predict the electric potential for a dielectric cylinder in homogeneous electric field, and is used to perform error analysis on the proposed method. The third example develops a surrogate model for the evolution of the strain field of amorphous solids under simple shear using the shear transformation zone (STZ) theory of plasticity. Finally, concluding remarks are provided in Section 7. Furthermore, the algorithms presented in this paper have been implemented in UQpy (Uncertainty Quantification with python) [53], a general purpose Python toolbox for modeling uncertainty in physical and mathematical systems.
2. Grassmann manifold. The concepts presented in this section are essential for the development of a surrogate model based on the Grassmannian diffusion maps. Let us begin by defining a p-plane as a p-dimensional subspace, and a p-frame as a coordinate system that spans that subspace. Based on these two concepts one can define two important manifolds, the Stiefel and the Grassmann manifold. The Stiefel manifold V(p, n) is the set of all p-frames in R n such that V(p, n) = {X ∈ R n×p : X X = I p }, where I p ∈ R p×p is the identity matrix and X ∈ R n×p is an orthonormal matrix [2]. The Grassmann manifold (or Grassmannian) G(p, n) is the set of p-planes in R n , where a point is given by X = span (Ψ), with Ψ ∈ V(p, n) [68]. It is worth noting that X is identified as an equivalence class of n × p matrices under orthogonal transformation of the Stiefel manifold [68,69,47]. Therefore, a point on the Grassmann manifold is represented by an orthonormal matrix Ψ ∈ R n×p (the Stiefel representation).
2.1. Tangent Space: Exponential and logarithmic maps. The Grassmann manifold is a smooth and continuously differentiable manifold, which enables numerous mathematical operations such as differentiation and optimization [22,1]. Given that the Grassmann manifold is smooth and continuously differentiable, one can define a trajectory γ(z), z ∈ [0, 1], known as geodesic, defining the shortest path between two points, γ(0) = X 0 = span(Ψ 0 ) and γ(1) = X 1 = span(Ψ 1 ), on the manifold G(p, n) [22]. The derivative of this trajectory at any point X (represented by Ψ) defines the tangent space (T X G(p, n)), which is given by the set of all tangent vectors Γ such that T X G(p, n) = {Γ ∈ R n×p : Γ Ψ = 0}. Therefore, given a tangent space T X0 G(p, n) at X 0 , one can map Γ 1 onto the Grassmannian point γ(1) = X 1 represented by Ψ 1 via the exponential map where Q ∈ R n×n is an orthogonal matrix satisfying the following expressions.
After appropriate manipulation, one can obtain the following expression.
Consequently, one can write the logarithmic map from the Grassmannian to the tangent space, log X : G(p, n) → T X G(p, n) as 2.2. Grassmannian distance. The properties of the Grassmann manifold further afford a notion of distance between points on it. Many definitions of distance exist and can be expressed in terms of the principal angles between subspaces [31]. One can easily see that the cosine of the principal angles θ i ∈ [0, π/2] between two subspaces X = span(Ψ x ) and Y = span(Ψ y ) can be computed from the singular values of Ψ T x Ψ y =ŪSV T , whereŪ ∈ O(k),V ∈ O(l), and S = diag(σ 1 , σ 2 , . . . , σ p ), with p = min(k, l). Thus, the principal angles are computed as θ i = cos −1 (σ i ) [50]. A well-known and commonly used distance is the geodesic distance, d G(p,n) (X , Y), corresponding to the distance along the geodesic curve γ(z) parameterized by z ∈ [0, 1], and given by d G(p,n) (X , Y) = s 2 [66,68,29], where s = (θ 1 , θ 2 , . . . , θ p ) is the vector of principal angles. See [19] for more definitions of distances on the Grassmann manifold.

Karcher mean.
Consider a set of points (subspaces) on G(p, n). The Riemannian center of mass of these points is known as Karcher mean, µ G(p,n) , and corresponds to the point Y ∈ G(p, n) that minimizes the cost function σ 2 G(p,n) : G(p, n) → R + [39,30] given by where dP (X ) = ρ(X )dG(p, n) is a probability measure on the Grassmann manifold with probability density function ρ(X ). Thus, the Karcher mean µ G(p,n) = span(M) can be computed by solving the optimization problem One can easily notice a similarity between Eq. (2.6) and the variance of a continuous random variable, which gives a notion of dispersion around the mean. The Karcher variance σ 2 G(p,n) and the Karcher mean µ G(p,n) can therefore be interpreted as the mean and variance of a continuous random variable on G(p, n). With this interpretation, the Karcher mean can be estimated from a discrete set of points on the Grassmann manifold S = {X 1 , . . . , X N } ⊂ G(p, n), by solving the following minimization

Grassmannian kernels.
Kernel-based dimensionality reduction techniques such as the conventional diffusion maps depend on the appropriate definition of a real-valued positive semi- In this regard, the Gaussian kernel (Eq. 2.9) is perhaps the most widely used.
where X i and X j are the high-dimensional data and is a length-scale parameter. However, the Gaussian kernel is not suitable to represent the underlying subspace structure of datasets. On the other hand, Grassmannian kernels are endowed with this feature, which is advantageous in the analysis of high-dimensional data. A Grassmannian kernel is defined as a real symmetric map k : G(p, n) × G(p, n) → R embedding the Grassmann manifold into a reproducing kernel Hilbert space. Moreover, a Grassmannian kernel is invariant to the choice of basis and is positive semidefinite. Several families of Grassmannian kernels, with different characteristics, are proposed in the literature (see [31,32,33]). However, the most popular kernels are the Binet-Cauchy and Projection kernels. The Binet-Cauchy kernel is constructed by embedding the Grassmann manifold G(p, n) into a projective space P ( p R n ). Considering two subspaces X = span(Ψ x ) and Y = span(Ψ y ), the Binet-Cauchy kernel is given by or equivalently in terms of principal angles [31,33] (2.11) Another kernel frequently used in kernel-based methods on the Grassmann manifold is the projection kernel. It is constructed based on the projection embedding Π : G(p, n) → R n×n such that Π (Ψ) = ΨΨ T . This kernel is defined as Further discussion of Grassmannian kernels and their specific use for Grassmannian diffusion maps can be found in dos Santos et al. [17].
3. Grassmannian diffusion maps. The Grassmannian diffusion maps (GDMaps) [19] is a nonlinear dimension reduction technique that uses diffusion maps to learn the low-dimensional structure of a dataset on the Grassmann manifold. GDMaps is a two-stage dimension reduction. The first dimensional reduction is a pointwise projection of the elements of a dataset onto the Grassmann manifold. Next, a connected graph is created across the data on the Grassmann manifold where a random walk is performed to embed the data into a low-dimensional Euclidean space. This dimension reduction is illustrated in Fig. 1 and the procedure is detailed herein.
Given a set of high-dimensional data S X = {X i , . . . , X N }, with X i ∈ R n×m possessing a lowrank structure, one can project X i onto a Grassmann manifold via Singular Value Decomposition Next, considering a positive semi-definite Grassmannian kernel k : G(p, .) × G(p, .) → R (i.e. the Projection or Binet-Cauchy kernel), we construct the kernel matrices Next, we compile to composite kernel matrix as either We use this composite kernel to define a random walk over the data on the manifold, denoted is the joint probability distribution of the input parameters, and P is the transition probability matrix. The matrix P = [P ij ] is constructed by first building the following diagonal matrix D such that the stationary distribution of the random walk is given by Next, by normalizing the kernel k ij as follows the transition matrix P can be constructed as The eigendecomposition of P yields a set of eigenvectors Φ = [φ 0 , . . . , φ N ] and their respective eigenvalues Λ = {λ 0 , . . . , λ N }. Thus, every element X i of S X has a representation on a lowdimensional Euclidean space defined by the points 4. Geometric Harmonics. The construction of a mapping between the reduced (i.e., diffusion space) and ambient (i.e., parameter space, response space) space has been discussed in the literature by several authors [52,12,23,11,24]. In this regard, the construction of the lifting (from the reduced space to the ambient space) and restriction (from the ambient space to the reduced space) operators relies on the extension of empirical functions only known at specific locations of the domain, also known as out-of-sample extension.
For example, the Nyström extension, which is very closely related to GP regression [4], is commonly used to construct the restriction operator within the conventional diffusion maps framework [52,11]. In particular, given a new sample X k in the dataset S X , the diffusion coordinates can be obtained as follows where P (·, ·) is the transition matrix. The Nyström extension has some disadvantages such as the computational complexity of the required diagonalization. Moreover, it can become ill-conditioned [4], so instead we leverage a variation on the Nyström extension, known as geometric harmonics [14], for constructing both the lifting and restriction operators. The geometric harmonics aim to extend a function f defined on a set Ω to a set Ω such that Ω ⊂ Ω. This extension depends on the selection of an appropriate positive semi-definite kernel k(·, ·) (such as the Gaussian kernel in Eq. (2.9)), which defines a unique reproducing kernel Hilbert space H of functions in Ω. Therefore, restricting k(·, ·) to Ω one can define an operator K : where dµ is a measure with dµ < +∞, ω ∈ Ω, and ω ∈ Ω. A lemma presented by Coifman and Lafon [14] shows that the adjoint operator K * : H → L 2 (Ω, dµ) is in fact the restriction operator. Moreover, as this operator is self-adjoint and compact, its eigendecomposition exists and one can write the Geometric Harmonics as which can be summarized by the following expressions Therefore, the mechanization of the extension algorithm is given by two main steps [14]. First, f is projected onto Second, the extension Ef is given by Details on implementation of GH can be found in Algorithm 4.1. In this algorithm, one can observe that the sets S x and S y as well as S * x and S * y are composed of column vectors. Therefore, if the Euclidean space of interest is the space of matrices, one can transform an element of this space into vectors by stacking its columns. Further, one can easily observe that the restriction operator between a point on the Grassmann manifold G(p, n) and the Grassmannian diffusion space exists since appropriate kernels can be defined on the Grassmann manifold, see Section 3. On the other hand, there is no guarantee that the inverse map (lifting operator) exists due to the orthogonality constraints of the Grassmann manifold. However, assuming that the Grassmann manifold is locally approximated by a flat inner-product space (i.e. the tangent space T X0 G(p, n)) constructed at X 0 , one can define a local lifting operator from the Grassmannian diffusion space defined on a set Ω ∈ R k to T X0 G(p, n). It is then straightforward to apply the exponential mapping to project the extended sampling points in T X0 G(p, n) onto G(p, n) (see Section 2.1).

Grassmannian-Geometric Harmonics Maps.
Consider the random vector Θ ∈ R Q having joint probability distribution F Θ (Θ 1 , . . . , Θ Q ) as the input parameters to a model M(·). One can obtain samples Θ i as elements of a set where Π is the parameter space. For each element of S Θ , the model M(·) (e.g., finite element model) produces a high-dimensional response X i ∈ R n×m such that With the set S X and assuming that X i has a low-rank structure, we begin by projecting X i onto a Grassmann manifold. This operation is performed via singular value decomposition (SVD), as presented in Section 2. Thus, one can decompose X i as Selecting an appropriate Grassmannian kernel [19], we next construct a connected graph on the sets S U = {U 1 , . . . , U N } ⊂ G(p, n) and S V = {V 1 , . . . , V N } ⊂ G(p, m) and apply the procedure of Section 3 to determine the new coordinates (Grassmannian diffusion coordinates) embedding the data on the Grassmann manifolds G(p, n) and G(p, m) into a low-dimensional Euclidean space (Grassmannian diffusion space). Once the Grassmannian diffusion coordinates ξ = {ξ 1 , . . . , ξ N } ∈ ∆ are obtained, we construct a global map (surrogate) using geometric harmonics (see Section 4) between Π and ∆ (GH 0 : Π → ∆), although local maps can also be constructed if necessary. For this mapping, the Gaussian kernel k(Θ i , Θ j ) (Eq. 2.9) is used because we want to construct a map between Euclidean spaces. A second map between Π and Σ (GH 1 : Π → Σ) is constructed using GH as well. With both maps constructed, we can predict the dimension reduced response of M(·) for any new set of input parameters. In other words, considering that a new set of input parameters Θ * is sampled from F Θ , we estimate the corresponding Grassmannian diffusion coordinates in ∆ by ξ * = GH 0 (Θ * ). Simultaneously, we estimate the singular values by S * = GH 1 (Θ * ).
Once the estimated coordinate ξ * in ∆ is obtained, it is necessary to expand this reduced dimension solution from the Grassmannian diffusion manifold back to the full-dimensional solution in the ambient space. The first step of this decoding is to define a mapping between ∆ and the Grassmann manifolds G(p, n) and G(p, m). To achieve this, we define a series of local GH lifting operators Λ 0 : ∆ → Tμ u G(p, n) and Λ 1 : ∆ → Tμ v G(p, m) whereμ u andμ v are the reference points on the Grassmann manifold where the tangent spaces are constructed. The local operators are determined by identifying the k nearest neighbors to the point ξ * . These points define the vicinity of ξ * on the Grassmann manifold, and the tangent space can be constructed either around their Karcher mean (see Section 2.3); or around the nearest neighbor of ξ * on the Grassmann manifold, which is a computationally efficient method since no optimization is performed. These local data are then used to construct the GH lifting operator. Using these lifting operators, we obtain the points Γ * u = Λ 0 (ξ * ) and Γ * v = Λ 1 (ξ * ), where Γ * u ∈ R m×p and Γ * u ∈ R m×p represent points on their respective tangent spaces. Next, we apply the exponential map to obtain the corresponding points on G(p, n) and G(p, m) (see Section 2.1) as U * = expμ u (Γ * u ) and V * = expμ v (Γ * v ), where U * = span (U * ) ∈ G(p, n) and V * = span (V * ) ∈ G(p, m). Finally, the solution X * for the new set of input parameters Θ * can be predicted by the following matrix product Next, two algorithms are presented summarizing this method. Algorithm 5.1 describes the construction of the maps between the spaces of interest (training), and Algorithm 5.2 shows how to predict the response using the constructed maps. Moreover, the proposed surrogate modeling approach is illustrated conceptually in Fig. 2. Compute the thin Singular Value Decomposition:
6. Examples. In this section, three examples are considered to demonstrate the versatility of the proposed surrogate modeling approach. We begin with a toy example in which structured points on the Grassmann manifold can be easily visualized as points on the unit sphere. In the second example, the electrical potential field of an infinitely long dielectric cylinder suspended in a homogeneous electric field is predicted considering that the cylinder's radius r 0 and the strength of the electric field E ∞ are random variables. The third example considers the evolution of the strain field in an amorphous solid under simple shear using the shear transformation zone (STZ) theory of plasticity.
The projection kernel in Eq. (2.12) is adopted in all examples presented in this section, and the kernel composition by the Hadamard product is considered. Moreover, the accuracy of the predicted solutions is evaluated by using the entry-wise relative error for matrices (Eq. (6.1)) and the relative error in a L 2 -norm (Frobenius for matrices) sense (Eq. (6.2)).
6.1. Structured data on the unit sphere in R 3 . Consider the following set of equations, x = |r|sin(t)cos(s), y = rsin(t)sin(s), z = |r|cos(t); such that r is uniformly distributed in the interval [−2, 2], t is uniformly distributed in the interval [−π/2, π/2], and s = sin −1 cos(t) 2 . We draw N = 3, 000 sample pairs (r, t) ∈ Π to obtain a collection of N points constrained on two cone-like structures in R 3 as presented in Fig. 3a, with the colors representing the magnitude x 2 + y 2 + z 2 . In effect, we have a model (Eq. (6.3)) that maps two random variables onto a surface in R 3 . Each point is represented by a column vector X i = [x i , y i , z i ] T , which together compose the set S X = {X 1 , . . . , X N }. One can readily see that these points can be projected onto the Grassmann manifold G (1, 3), which is the unit sphere S 3 . A point X i on G(1, 3) is given by the unit vector obtained from the normalization of X i such that X i = X i / X i 2 , which reveals two inverted teardrop shaped structures on the sphere as illustrated in Fig. 3b. Applying diffusion maps to these points on the Grassmann manifold, we see that a well-defined parametrization is obtained, as revealed by the Grassmannian diffusion coordinates in Fig. 3c.
Geometric harmonics is used to create a map GH 0 : Π → ∆ from the parameter space Π to the Grassmannian diffusion manifold ∆. Therefore, it can be considered as a manifold learning technique, where the position on the Grassmannian diffusion manifold (Fig. 3c) can be predicted for any point Θ in the parameter space Π. To verify the accuracy of this learning process, we draw 3,000 additional samples Θ (Fig. 4). One can easily see in Fig. 4 that the trained GH can reliably predict the shape of the Grassmannian diffusion manifold.
We obtain a new parameter vector Θ * = (r, t) = (1, 1.3) by sampling Π, and the GDMaps-based surrogate model is used to predict the vector X * . Using the map GH 1 : Π → Σ one can predict the first two nontrivial diffusion coordinates as ξ * = [−1.3657×10 −2 , −6.5931×10 −3 ], represented by the red star in Fig. 5a. Moreover, we observe that S * = r = 1 determines the magnitude of the predicted point in R 3 . Selecting the k = 3 nearest neighbors of ξ * in Algorithm 5.2, we predict the point in the ambient space (i.e. on the cone structure) corresponding to Θ * by mapping ξ * onto the tangent space T , constructed in the closest neighbor of ξ * (one can also use the Karcher mean alternatively), for posterior projection onto the Grassmann manifold G (1, 3), as illustrated by the red star in Figure 5b. This point coincides very closely with with true point denoted by the black dot. From the the point Next, we draw N = 10, 000 pairs (r, t) ∈ Π to assess the overall performance of the proposed surrogate modeling technique. Using this new set of input parameters we predict N points on the cone-like structure and compare them with the exact points corresponding to the set of input parameters. Figure 6 shows the predicted cone-like structures from these 10,000 surrogate model evaluations. Comparing with Figure 3, we can see that the points closely match the true structure.
To assess the overall quality of the predictions, the marginal probability density functions (PDF) for each dimension (x, y, z) are estimated using the kernel density estimation (KDE) and shown in Fig. 7 for both the true samples and the surrogate predictions. The PDFs for the surrogate predictions match those of the true samples very closely.
6.2. Dielectric cylinder in homogeneous electric field. In this example, we study variations in the electrical potential of an infinitely long dielectric cylinder suspended in a homogeneous electric field, resulting from uncertainty in the input parameters. The problem is defined over a rectangular domain Ω = [−1, 1] × [−1, 1] with the embedded cylinder domain D c = x = (x, y) | x 2 + y 2 ≤ r 0 , where r 0 is the cylinder's radius. We assume Dirichlet boundary conditions, Γ D , on the left and right boundaries and Neumann boundary conditions Γ N on the top and bottom boundaries.
The electric potential u(x) in Ω can be computed by solving the Laplace equation where n denotes the outer normal unit vector. The permittivity ε (x) is given by and u * , which is also the analytical solution to the problem, is given by where E ∞ is the strength of the homogeneous electric field. We consider as random variables the cylinder's radius r 0 , and the strength of the electric field E ∞ . These variables are uniformly distributed as defined in Table 1. The GDMaps based surrogate model presented in this paper is composed of three distinct maps based on geometric harmonics: (1) A global map between the parameter space and the Grassmannian diffusion manifold (GH 0 : Π → ∆); (2) a global map between the parameter space and the singular values space (GH 1 : Π → Σ); and (3) a local map between the Grassmannian diffusion manifold and the tangent space of a region of the Grassmann manifold Λ : ∆ → T . Moreover, an additional map is built for the projection of the points in the tangent space onto the Grassmann manifold. Considering this sequence of mappings, one can expect that errors can propagate within the framework pipeline, negatively affecting the response prediction. Therefore, the identification of the source of errors is relevant for increasing the accuracy of the predicted outcomes. Therefore, the model presented in this example is used to demonstrate how the errors can be reduced by an appropriate selection of parameters for the construction of maps based on GH.
For the error analysis presented in this section, let's assume a training set of input parameters Θ = (r 0 , E ∞ ) of N = 300 samples and the corresponding electrical potential fields X = M(Θ) ∈ R m×n over the domain Ω discretized into (n × m) = (300 × 300) = 90, 000 mesh points. In this example, the Grassmann manifolds G(p, n) = G(74, 300) and G(p, m) = G(74, 300) are sufficient to encode the geometric structure of the electrical potential fields X. Using Algorithm 5.1, we obtain points on the Grassmannian diffusion manifold; and according to the parsimonious representation of the diffusion maps [20], only the first non-trivial diffusion coordinate is relevant to encode the geometrical information of the underlying physical phenomenon. Therefore, a dimension reduction from n × m = 90, 000 to r = 1 is achieved using Grassmannian diffusion maps.
To assess the ability of the global map GH 0 : Π → ∆ to learn the relationship between points in the parameter space and the structure of the data on the Grassmannian diffusion manifold (∆), we start by sampling 1,000 additional points from Π to be mapped in ∆ using GH 0 . In Fig. 8, which shows both training data and predicted points in the Grassmannian diffusion coordinates, we see that GH 0 has adequately learned the shape of the data on the Grassmannian diffusion manifold using Gaussian kernel with length-scale equal to G = 1, and retaining q = 50 eigenvalues in the GH framework. However, these parameter, together with the number of samples in the training set, can influence the accuracy of the geometric harmonics maps. Next, the influence of ( G ) and q on the accuracy of both GH 0 : Π → ∆ and GH 1 : Π → Σ is analyzed.
We start by plotting the decay of the eigenvalues for GH 0 for different values of G . It is clear from Fig. 9 that as G increases, the eigenvalues tends to decay quicker. This behavior will have a strong influence on the prediction error of GH 0 and GH 1 because the construction of the matrix B in Algorithm 4.1 depends on the reciprocal of the corresponding eigenvalues. Therefore, if a larger number of very small eigenvalues are retained, their reciprocal could lead to large errors and numerical instabilities. Thus, the selection of G is inherently connected with the number of eigenvalues one should retain for the construction of B. This analysis is presented in Fig. 10a for GH 0 , and in Fig. 10b for GH 1 where we show the average error in the GH predictions for different combinations of G and q. In both cases, we see that for large values of G , a smaller number of eigenvectors and their respective eigenvalues should be retained in the construction of the matrix B, which has a direct influence in its rank. Selecting a larger G along with a high q will introduce large errors. Next, the influence of the size of the training set N is investigated. In this regard, the initial 300 samples and their respective diffusion coordinates are split into a training and a testing set to which the predicted Grassmannian diffusion coordinates can be compared. Considering N train varying from 20 to 200, and keeping G = 1 and q = 50 constants, the average error in the L 2 -norm sense as a function of N train is presented in Fig. 11a, for GH 0 , and in Fig. 11b, for GH 1 . In both cases the error reduces, although at a limited rate, after a certain value of N due to the fact that a residual error remains due to the selected values for G and q.
Next, we investigate the accuracy of the local maps, Λ 0 (·) and Λ 1 (·), between ∆ (Grassmannian diffusion space) and the tangent spaces Tμ u and Tμ v , respectively. The exact diffusion coordinates for each of the 300 training points are used, and the prediction of the corresponding points on the tangent spaces are obtained using the local maps with different numbers of neighbors (k) used to construct the local geometric harmonics maps. The probability density functions (PDFs) for the error L2 are estimated using kernel density estimation for k =3, 5, and 10 (closest neighbors points) in Figure 12a. We clearly see that the errors induced by these local maps are minimal and that they are not strongly influenced by k in this specific problem, because a point is predicted in a region close to the reference point where the tangent space is constructed on. Note also that the length-scale parameter for the local maps is taken to be 0.25 times the square value of the median of the pairwise distances of the k neighbors.
The cumulative error associated with the full process is analyzed by drawing 1,000 sample  points from the parameter space and computing posterior error estimation (error L2 ) of the predicted response. The estimated PDF of the error is shown in Fig. 12b, where the mean is equal to 1.8489 × 10 −3 and the standard deviation is equal to 2.7598 × 10 −3 . This reveals that the overall errors in the prediction solutions are very small compared to their true solutions, even considering An important theoretical hypothesis about the behavior of amorphous solids (e.g., metallic glasses) says that irreversible plastic deformation is mediated through atomic rearrangements in small clusters of atoms known as shear transformation zones (STZs) [25]. Consequently, amorphous materials subject to large shear stresses are often prone to the formation of shear bands due to the rearrangements of STZs in localized regions. It has been proposed that one can connect these large-scale plastic deformations to an effective temperature [6] (6.7) where U c and S c are the potential energy and entropy of the configurational degrees of freedom under the assumption that both the total energy U and total entropy S are separated into kinetic and configurational components, i.e. U = U c + U k and S = S c + S k , respectively. This effective temperature provides a measure of the degree of structural disorder (characterizing the density of STZs) and can be dedimensionalized as where k B = 1.38 × 10 −23 is the Boltzman factor, E z is the STZ formation energy. Given a spatially varying initial effective temperature field χ on a material domain, the STZ theory defines two coupled equations to model the evolution of plastic strain in the material. The first describes a plastic flow rule that relates the plastic rate of deformation tensor D pl to the effective temperature as: Note that this flow rule is monotonic with respect to σ/σ y , with σ = |σ 0 | given as the magnitude of the deviatoric shear stress σ 0 = σ = 1 3 1Tr(σ). Therefore, plastic deformation does not occur when σ/σ y < 1.
The second equation describes the evolution of χ as Other parameters are defined in Table 2.  [57,51] is used to solve the system of Eqs. (6.9-6.10) . This method utilizes an Eulerian finite-difference method under quasistatic conditions. As mentioned previously, the STZ theory assumes that the effective temperature has a spatial distribution that influences the material response. Therefore, the evolution of Eqs. (6.9-6.10) depends on the initial χ field. This field is assumed to be Gaussian [36,30,41], therefore it can be characterized by the mean µ χ and coefficient of variation c χ . Herein, it is assumed that these parameters are both uncertain with uniform distributions as described in Table 3. The associated correlation structure and length-scale are inferred from molecular dynamics simulations [36,40]. In the simulations, a simple shear up to 50% strain is imposed to a simulation box of size 400Å× 400Å. A grid of size 32 × 32 is considered in the discretization, where each element has a size of 12.5Å× 12.5Å. Therefore, each snapshot of this simulation is given by a matrix X i ∈ R 32×32 . We obtain N = 196 samples of the pair Θ = (µ χ , c χ ) via stratified sampling to train a surrogate model for full evolution of the plastic strain field. The evolution of the plastic strain field for a given pair Θ = (µ χ , c χ ) is presented in Fig. 14 as a sequence of 101 snapshots of size L x and L y at discrete values of the imposed shear strain¯ . A matrix Y ∈ R 1024×101 is then constructed for a given pair Θ = (µ χ , c χ ), where each column of Y correspond to the vectorized snapshot of the plastic strain field. Figure 14: Evolution of plastic strain with a snapshot of the strain field at a given strain level extracted.
In this problem, two Grassmann manifolds given by G(p, n) (left manifold) and G(p, m) (right manifold) are associated with the left and right singular vectors of the matrices Y ; where n = 1024, m = 101, and p = 10 suffices to encode the geometric structure of each data point. Using GDMaps, we obtain a set of 196 Grassmannian diffusion coordinates embedding the high-dimensional data into a low-dimensional Euclidean space as shown in Fig. 15 with r = 3.
Once the surrogate model is trained using the 196 samples obtained previously, we can predict the full evolution of the plastic strain field for any pair Θ = (µ χ , c χ ). Considering a representative case with Θ = (530.1748, 0.0792), the simulated and predicted evolution of the plastic strain field, as well as their relative error, are presented in Fig. 16 for five different levels of imposed strain (0%, 12.5%, 25%, 37.5%, 50%). The error in the L 2 -norm sense for this plastic strain field is equal to 2.0473 × 10 −3 . Next, considering 100 additional samples we compute the mean and standard deviation of the plastic strain fields at the different levels of imposed strain (25%, 30%, 35%, 40%, 45%, 50%) as presented in Figs. 17, and 18. These figures include the statisti- cal characterization obtained by using the numerical model and the surrogate model, along with the relative errors. From these results, we see that the surrogate model developed herein can predict the uncertain response of a complex model with high-accuracy, by taking advantage of low-dimensional subspace structure of the problem to reduce the computational burden associated with running high-fidelity models for UQ.

0%
12.5% 25% 37.5% 50% 7. Conclusions. This paper introduces a fully data-driven surrogate model for uncertainty quantification of high-dimensional models of complex physical/engineering systems. It takes advantage of the low-dimensional representation of high-dimensional input/output data obtained via Grassmannian diffusion maps to create a set of geometric harmonics based maps. A global map is constructed to predict the Grassmannian diffusion coordinates corresponding to any new element in the set of input parameters with good accuracy. Once the Grassmannian diffusion coordinates corresponding to a new set of input parameters are predicted, the k-nearest neighbors points in the Grassmannian diffusion space and their associated points on the Grassmann manifold are utilized to estimate, via geometric harmonics, a local map from the Grassmannian diffusion space to a tangent space. Next, the exponential map project the point onto the Grassmann manifold, a result used to predict the high-fidelity solution of the problem.
The method developed herein used the descriptive power of the Grassmannian diffusion maps and the computational performance of geometric harmonics to provide an efficient and accurate prediction of the solution of complex systems described by algebraic equations and partial/ordinary differential equations. Three examples were considered to evaluate the performance of this technique. The first one consisted of a toy example to demonstrate the ability of the technique to predict data with complex geometry using spectral methods in a way that is easy to understand and visualize.
In the second example, the performance of the surrogate modeling developed herein was verified in a physical model (i.e., electric potential of a cylinder in homogeneous electric field) with highdimensional response, also considering discontinuities in the system response. It was demonstrated that some parameters such as the length-scale parameter and the number of retained eigenvalues, both for the Gaussian kernel used in the geometric harmonics framework; as well as the amount of data in the training set, are important quantities affecting the accuracy of the presented technique. The third problem analyzed in this paper evaluated the plastic deformation of amorphous solids using the shear transformation zone (STZ) theory of plasticity. The uncertainty was imposed in the mean and coefficient of variation of the initial nondimensionalized effective temperature field (χ). In this case, the evolution of the strain field with the strain level is also taken into consideration, and the uncertainties of the plastic strain field are predicted accurately. In all the cases considered herein, a good accuracy was identified in the predicted solutions in comparison with the exact ones. The method proves advantageous due to its computational performance and ability to make reliable predictions for high-dimensional responses considering a highly sparse set of points.