The structure and constitution of a material dictate the properties of the waves it can support, which are thoroughly described by the material's dispersion relation. A material's dispersion relation is the relationship between the frequencies and wavevectors of all linear elastic waves that can travel in the material. Dispersion relations contain information about bandgaps (frequencies at which waves experience spatial attenuation), group velocity (the velocity with which energy can be transported by waves in the material), and phase velocity (the velocity at which a wave peak will move in the material). Bandgap information can be leveraged to engineer filtering and reflection, and group velocity and phase velocity information can be leveraged to achieve refraction. These fundamental behaviors can in turn be utilized to design complex wave propagation devices such as wave guides [1][2][3], wave lenses [4][5][6], wave multiplexers [1,2,7], and more.
Current methods to calculate dispersion relations are based on numerically computing the eigenvalues of the elastic wave equation [8]. Computational limitations are especially restrictive when these expensive forward models are used to solve the inverse problem of design [9][10][11][12][13]. As a result, scientists and engineers are often forced to simplify design workflows by limiting the design space to be relatively narrow, considering only partial dispersion information and/or relatively basic design objectives [1,2,4,9,[12][13][14].
Our method aims to increase the accessibility of full dispersion relations to alleviate the need for simplifications such as partial dispersion computations. Full dispersion relations are defined over a region of reciprocal space (wavevector space) called the irreducible Brillouin zone (IBZ). Due to computational costs, dispersion relations are often computed only along one-dimensional slices of these regions called irreducible Brillouin zone contours (IBZ contours). However it has been recently shown that using this simplification can lead to errors [15]. Further, plane symmetry groups without mirror symmetries do not have a uniquely defined IBZ contour. For unit cells belonging to these plane symmetry groups, dispersion analysis along a contour will yield different results depending arbitrarily on the choice of contour. As a second example, a common shortcut when designing for group velocity is to compute the dispersion relation at only a handful of wavevectors or frequencies of interest. This is useful for achieving desired wave behavior in specific directions of propagation [12], or at specific frequencies of interest [14]. Although these methods are sufficient for designing simple wave devices, the design of next-generation, multi-functional wave devices will require the efficient computation of full dispersion relations.
Techniques have been developed specifically to efficiently solve dispersion problems such as Bloch mode synthesis [16,17], spectral X-FEM [18], specialized mixed-variational formulations [19], reduced Bloch mode expansion [20], reduced Bloch operator finite element methods [21], and wave finite element methods (WFEM) [22,23]. However, even these specially tailored methods executed by a highly parallelized and optimized implementation on a multi-CPU multi-GPU cluster admit the need for more efficient ways to explore the design space [13]. To perform design optimization on a 48 3 voxel domain, a distributed computing setup using 16 NVIDIA GTX-780 GPUs takes 6.5 min per optimization iteration, even when computing the dispersion relation only along the IBZ contour rather than over the full IBZ [13]. Fundamentally, each of [16][17][18][19][20][21][22][23] reduces the cost of dispersion computations by reducing the number of effective degrees of freedom involved in the model, thereby reducing the size of the eigenvalue problem that needs to be solved at each wavevector point. In contrast, our surrogate model reduces the total number of eigenvalue problems that need to be solved. Due to this construction, our surrogate model can seamlessly integrate with any of these aforementioned techniques to produce an efficiency improvement that compounds with that of the physics-based technique being paired with our model. Surrogate models, often involving machine learning, have been proposed in a variety of settings to reduce the computational cost of purely physics-based methods. The applications range from multi-scale modeling [24][25][26][27], to dynamics [28], statics [29], turbulence [30], optics [31], and many others. For dispersion relation computations specifically, a fixed-resolution convolutional neural network (CNN) surrogate model has been proposed [32]. This CNN model is shown to vastly outperform the traditional densely connected neural network, suggesting that the convolutional architecture is more appropriate for these types of problems. Because CNNs can make predictions so quickly once the model is trained, it seems like a very attractive option. However, neural networks often require large amounts of training data, which are expensive to compute (hundreds of thousands of precomputed dispersion relations [32]) and require a lengthy training period (13 min on a GPU reported just for the 1D problem [32]). Unless the trained model is subsequently used to make predictions for more dispersion relations than were included in the training set, the time used to generate the training data outweighs the time saved by using the model. Further, the CNNs architecture locks the user into a fixed unit cell resolution and fixed dispersion relation resolution. If a user desires to change the resolution of the discretization of the unit cell or IBZ, a new model architecture needs to be trained on a new training set with the corresponding resolution. Recently, there has been a push towards surrogate models that exhibit resolution independence, meaning that the model can train on data of one spatio-temporal resolution, and make predictions at another resolution [33][34][35]. In the spirit of this effort, the surrogate model we propose exhibits resolution independence both in terms of unit cell discretization and IBZ discretization.
We propose using Gaussian process regression [36] to obtain principled surrogate models of full dispersion relations. In particular, we propose a data-driven covariance function that only requires a small training set of dispersion relations (tens to thousands) and allows for an efficient acquisition function. In comparison to deep learning-based surrogate models, our model achieves higher accuracy with less training data and is invariant to unit cell representation.
The key features of our model are:
1. The model is independent of unit cell representation. For example, if the model learned covariance information from a dataset of parametrically defined unit cells, it can still be used to make predictions for pixelated designs, and vice versa. 2. The model requires a relatively low amount (tens to thousands) of training data samples to make accurate predictions. The model learns from data by measuring covariance information, which takes less than 20 s even for reasonably large training set sizes. 3. The model uses a physics-based method to sample the dispersion relation at various points in the IBZ. This physics-based method is a modular component of the model that can be swapped for other, faster physics-based methods for dispersion computations. Because of this, the efficiency gains provided by our model directly augment the efficiency gains of the physics-based model. 4. The model not only supplies predictions, but also offers uncertainty estimates. 5. The model is differentiable, allowing for its use in gradient-based topology optimization.
In this paper, we first discuss the traditional methods currently used for the computation of dispersion relations and gradientbased optimization of dispersion relations. We then go on to provide a brief overview of Gaussian process regression (GPR). Next, we describe the details of our model and demonstrate how it would be applied to one-dimensional problems. Finally, we apply our model to two-dimensional unit cells and evaluate its performance against traditional methods as well as another existing surrogate model.
To compute the dispersion relation of a periodic structured material, the physical domain of analysis is reduced to a single unit cell. The frequency domain dynamics of linear-elastic periodic solids are governed by the harmonic elastic wave equation with Bloch-Floquet periodic boundary conditions. The harmonic elastic wave equation looks for wave mode solutions that obey linear elasticity, and assumes that these solutions will be temporally sinusoidal. The Bloch-Floquet boundary conditions enforce that the solution is periodic (since the domain is periodic), but with a complex phase offset. The complex phase offset is used to make sure we are finding modal wave solutions with a variety of different wavelengths (not just wavelengths that are integer divisors of the unit cell's dimensions).
The harmonic elastic wave equation in 𝜔 dimensions is
and the Bloch-Floquet periodic boundary conditions are
The density, 𝜀(𝜗), and the material stiffness tensor, 𝜕 𝜍𝛻ℵℶ (𝜗), define the linear-elastic material properties of the unit cell, 0. The wavevector, ⊲, and the eigenfrequency, 𝜛, inversely describe the length and time scales of the wave. The displacement field, 𝜚(𝜗), describes the vibrational mode by which the wave travels. Finally, the lattice vectors, ℷ ℸ , define the lattice of the periodic material, and 𝛚 is the imaginary unit. The equation is a generalized eigenvalue problem with 𝜚(𝜗) and 𝜛 2 as the eigenvector and eigenvalue respectively. To compute the dispersion relation, the eigenvalues of this equation must be determined as a function of ⊲, ϑ⊲ ϖ 1 (letting 1 denote the irreducible Brillouin zone).
With the exception of a few simple cases, analytical solutions for dispersion relations are not obtainable. For this reason, numerical solutions are often used to discretize the generalized eigenvalue problem (Eq. ( 1)). Although there are many numerical methods to compute dispersion relations, one of the most well-known (though not necessarily the most efficient) is the Finite Element Method. Within the framework of FEM, the dispersion relation can be computed by assembling the stiffness and mass matrices, 2(⊲), and 3(⊲), and solving the generalized eigenvalue problem for each ⊲ in a discretization of the IBZ:
There are several other specially tailored methods that have been shown to be more efficient than the finite element method when solving this particular PDE [16][17][18][19][20][21][22][23]. We note that the surrogate model we have developed can be used in conjunction with any of these more highly-specialized methods.
Gradient-based optimization is a powerful design tool that can be used in many contexts (including wave propagation [9,12,13]) to obtain material designs that exhibit desired behaviors. Designs are represented by a design vector, 4, containing parameter values that uniquely determine the geometry and constitutive properties of the unit cell design. The goal of design optimization is to determine a 4 that results in the desired wave behavior, often achieved by maximizing or minimizing an objective function (which encodes design objectives) via gradient-based methods. To use these gradient-based methods, gradients of the objective function are needed, which in turn depend on gradients of quantities embedded in the dispersion relation. In this section, we quickly step through a few known equations commonly used in optimization problems surrounding wave behavior. In this context, 2 and 3 are not only functions of the wavevector, but also of the design vector, 4. Derivations of these equations can be found in Appendix C.
To modify the frequency, 𝜛, at a specific wavevector by incrementally changing the design variable, the following sensitivity equation is useful. Assume all eigenvectors are normalized by the mass matrix; i.e., 𝜚 ω 3𝜚 = 1 where (⋛) ω denotes conjugate transpose. For ease of reading, assume here that 4 is a scalar. These equations can easily be extended to vector 4 by repeating the following computations for each entry of 4. The gradient of one frequency in the dispersion relation with respect the design variable is given by
To modify the group velocity, 𝜑𝜛 𝜑⊲ , at a specific wavevector, the sensitivity computation is more complicated and more computationally expensive. The group velocity can be computed analytically by
The design sensitivity of the group velocity can be computed analytically by
where 5 = 𝜑2 𝜑⊲ ε 𝜛 2 𝜑3 𝜑⊲ , and 𝜑𝜚 𝜑4 can be computed by solving the linear system:
)
)
Not only is it computationally expensive to compute dispersion relations, but it is also computationally expensive to compute the design sensitivities of dispersion relations. Lu et al. report sensitivity analysis accounting for nearly half the time consumed by one iteration of optimization [13]. Later, we will explain how our surrogate model gives multiple options for computing sensitivities, opening the possibility for further computational savings. We also provide explicit equations for how to compute these gradients when using the surrogate model.
Gaussian process regression (GPR) is a probabilistic regression method that constructs a predictive distribution for an unknown function, 7 (𝜗), by combining observations of the function with prior information [36]. Because the function is unknown, we think of 7 (𝜗) as being a random function distributed according to a Gaussian process ∱∲(8(𝜗), ℵ(𝜗, 𝜗 ϱ )) defined by mean function 8(𝜗) and covariance function ℵ(𝜗, 𝜗 ϱ ).
Evaluation of this function at any finite number of points is described by a joint Gaussian distribution. Let 8 and ℵ operate row-wise when given matrix input, with 8(9) yielding the mean vector corresponding to the points indicated by the rows of 9, and ℵ(9, 9 ϱ ) yielding the covariance matrix for the points in 9 and 9 ϱ . With matrix 9 ω and vector 7 ω denoting the set of points for which we would like to make evaluations, the prior distribution is
where ∳ denotes a normal distribution. Given a dataset . = {(𝜗 𝜍 , , 𝜍 )}, 𝜍 = 7 (𝜗 𝜍 ) + <, < ς ∳ (0, ℏ 2 7 )} > 𝜍=1 of (possibly noisy) observations of 7 (𝜗), the predictive distribution is formed by conditioning the prior distribution on the observations. With 9 and , denoting the observations, and ⋆ denoting the identity matrix, the predictive distribution is
Numerically, Eq. ( 10) is often implemented in two steps. The first step is to find the model weights, ≨, which are independent of the prediction locations, by solving the linear system in Eq. (11). Subsequently, predictions can be made at arbitrary locations 9 ω using the pre-computed model weights ≨ in Eq. ( 10) rather than solving the system anew.
Because GPR yields a predictive distribution, we not only infer a prediction as the mean of the distribution, but also a measure of uncertainty as the variance of the distribution.
In later sections, we perform eigenfunction analysis of our proposed covariance function. Eigenfunction analysis is often used to assess the convergence and expressive power of a Gaussian process regressor. A function 𝐴(𝜗) is an eigenfunction of the covariance function ℵ(𝜗, 𝜗 ϱ ) with eigenvalue 𝐵, with respect to measure 6(𝜗) if it satisfies
These eigenfunctions form a basis which spans the space of all possible functions that can be attained by the mean function of the Gaussian process regressor, which is called the hypothesis space. A Gaussian process regressor can exactly represent any function in the span of the eigenfunctions with nonzero eigenvalues. A covariance function with an infinite number of nonzero eigenvalues is said to be nondegenerate, while one with only a finite number of nonzero eigenvalues is said to be degenerate.
A Gaussian process regressor with a nondegenerate covariance function has almost unlimited expressive power since its hypothesis space is the span of an infinite eigenbasis. In the limit of infinite observations, the model can exactly fit any function that can be represented by the generalized fourier series ⌉ ∇ 𝜍=1 ℷ 𝜍 𝐴 𝜍 (𝜗). In other words, for a GP regressor with a nondegenerate covariance function, the posterior mean function will converge to the expected value of the true underlying distribution as the number of observations goes to infinity (assuming the true underlying distribution is reasonably well-behaved).
In contrast, a Gaussian process regressor with a degenerate covariance function has limited expressive power, since its eigenfunctions span only a finite dimensional hypothesis space. For a GP regressor with a degenerate covariance function, in the limit of infinite observations, the posterior mean function will instead converge to the best 6-weighted 𝐷 2 approximation of the expected value of the true underlying distribution within the span of the eigenfunctions. Practically, this means that although a Gaussian process regressor with a degenerate covariance function may not converge to the exact solution in the limit of infinite regression points, it can still make a highly accurate prediction if its eigenbasis is comprehensive enough for the desired application.
In this section, we give the specifics of our surrogate model, and demonstrate its application to relatively simple one-dimensional dispersion computations. This toy problem is meant only to serve as an example -the substance of our results can be found in the following section. Traditional dispersion computations require the solution of an eigenvalue problem at each wavevector in a discretization of the irreducible Brillouin zone (IBZ). The goal of our surrogate model is to predict the eigenvalue solutions for each wavevector after solving the eigenvalue problems at only a small subset of these wavevectors. For most unit cells, which have many physical degrees of freedom, it is common in dispersion analysis to only consider the smallest > eigenvalues. In our surrogate model framework, we use > separate GPR models to fit the > eigenvalue bands, {𝜛 𝜍 (⊲)} > 𝜍=1 . The structures used in this section for demonstration are periodic one-dimensional bars. They can only transmit waves along their axis and have relatively simple dispersion relations. We randomly generate unit cell designs which have smoothly varying constitutive properties, except at randomly selected interfaces where the property values are allowed to be discontinuous (details in Appendix A.1). These interfaces lead to wave scattering behavior in the material -a mechanism well-known to cause band gaps in the dispersion relation. Examples of these unit cells and their dispersion relations are shown below in Fig. 1.
Our model uses a zero mean function (8(𝜗) = 0). This is common in GPR, and does not mean that the predictive distribution will have a zero mean. Nonzero mean functions are often used when the function to be fit has an underlying trend that can be written explicitly (such as a signal that, on average, increases linearly with time but has unpredictable variations around that average). Because dispersion relations have no such trends that can be explicitly written, we leave our mean function as zero. This leaves the entirety of the model's capacity to be carried by the covariance function.
The covariance function is usually the most crucial component in any Gaussian process model because it encodes the statistical coupling between function values at each point in the domain, defining underlying patterns and relationships in the data. Often, Gaussian process models use simple isotropic covariance functions such as the squared exponential or the Matérn covariance functions. Isotropic covariance functions are functions only of the distance between points in input space, meaning that they can only encode assumptions about the underlying function globally. These functions offer a convenient starting point for basic applications of Gaussian process modeling; however, they lack the expressiveness that can be offered by more thoughtful constructions of the covariance function. We propose a non-isotropic covariance function that encodes richer covariance information in a principled way that varies throughout the domain.
Our surrogate model uses a dataset of precomputed dispersion relations to construct its covariance function. Given a dataset 𝐸 = {(⊲, 𝜛) 𝜍 } of dispersion relations evaluated at a fixed set of wavevectors 1 𝐸 , we measure the covariance between the frequency values at all pairs of wavevectors in 1 𝐸 . This covariance computation is extremely fast, especially when compared to the training process of neural networks. For example, to compute the covariance matrix for 10,000 dispersion relations with 10,000 points evaluated in each dispersion relation, only 6 s on a 9th generation Intel Core i7-9750H are required. Once the covariance matrix is computed, we use cubic interpolation of the measured covariance values to define the covariance function, ℵ ∂ 1 𝐹 -1 𝐺 R 𝐹-𝐺 mapping any two lists of wavevectors to the appropriate covariance matrix. Each eigenvalue band is described by its own separate covariance function. Although this data-driven covariance function is degenerate (i.e., in the limit of infinite regression points, the estimated dispersion relation converges to the best solution spanned by the finite eigenbasis rather than the exact numerical solution), it achieves higher accuracy than a data-agnostic covariance function given a limited number of regression points. Convergence to the exact numerical solution is not necessary in practical design problems, since even the exact solution does not perfectly represent real-world physics.
To demonstrate how this works for the one-dimensional case, we compute the covariance between frequency values at each pair of wavenumbers for the dataset shown in the right-most panel of Fig. 1. The results of this computation are shown in Fig. 2. For this demonstration, the dataset is small, with only 5 dispersion relations, and coarse, with the domain discretized into only 10 points. Even so, we can see that we have already captured critical information about the type of function we are trying to fit. For example, we see that the variance is 0 at ⊲ = 0, and that the frequency at this wavenumber has no covariance with the frequencies at any other wavenumber. This makes sense because, for the first eigenvalue band, this frequency is always zero, representing a rigid body mode, regardless of the rest of the dispersion relation. Conversely, the variance is maximum at ⊲ = 𝐻, and the frequency at this wavenumber is covariant with the frequency at many other wavenumbers, especially those nearby. In later sections, when working with 2D structures, there are two wavevector components to consider. In these cases, the covariance function should be thought of as a continuous scalar-valued function of four scalar variables, and covariance interpolation must happen in this four-dimensional space. Regardless, for all computational operations other than interpolation, standard matrix operations can still be applied after reshaping the four-dimensional covariance array into a matrix.
With the covariance and mean functions defined, our Gaussian process prior is defined. We can visualize our prior by sampling it. With the exception that the samples are allowed to assume negative values, they have characteristics similar to the dispersion bands (see the first panel of Fig. 3). This shows that our covariance function yields a prior that encodes an appropriate hypothesis space, even before sampling any regression points.
A.C. Ogren et al.
To fit the model to a specific dispersion relation, we supply it with a set of points (regression points) sampled from that dispersion relation. Because each regression point has a computational cost, it is beneficial to try to minimize the number of regression points needed, and maximize the information that each regression point contributes. An acquisition function is a function that tells the model how to select the locations of regression points. Here, we propose a greedy acquisition function that sequentially determines the most informative next point to sample. The acquisition function we have constructed uses the covariance function to evaluate the posterior variance over a discretization of the IBZ. The variance can be computed by evaluating the diagonal entries of the posterior covariance, 𝜕 𝐹𝐼𝐽𝐾⊳𝐿𝜍𝐼𝐿 , given by
where 1 ω are the wavevectors constituting the discretization of the IBZ, 1 . are the wavevectors representing the locations of the optimal regression points, and ⋆ denoting the identity matrix.
Starting with no regression points (1 . is empty, and the entire second term is zero), we consider the wavevector in 1 ω with the highest variance to be the most informative next wavevector to sample, so it becomes the next regression point location. We append this wavevector to 1 . , re-evaluate 𝜕 𝐹𝐼𝐽𝐾⊳𝐿𝜍𝐼𝐿 , and repeat. This process can be iterated until a desired number of regression point locations is obtained.
When this process is finished, 1 . represents an ordered list of optimal regression points. During the fitting process (described in Section 5.4), subsets of 1 . can be selected by taking the first > elements if the computational budget does not permit the solution of the phononic eigenvalue problem for every wavevector in 1 . . Because our covariance function is degenerate, there is a finite limit on the number of regression points that can be effectively used by the model, and 1 . will not be able to grow beyond this size (the inversion in Eq. ( 13) will run into conditioning issues). This limit can be quite large, and depends on the size and resolution of the dataset that was used to construct the covariance function.
Conveniently, our acquisition function is entirely independent of the observed values of the regression points. This means the locations (wavevectors) of the optimal regression points can be determined ahead of time without actually making observations of the regression points themselves, allowing for the parallel computation of regression point values (frequencies) later. Additionally, this means that the final 1 . is dependent only on the covariance function, and does not change as a function of the specific dispersion relation that the model is being fit to.
We visualize how our acquisition function works for a simple one-dimensional problem in the bottom row of panels of Fig. 3. Without performing eigenvalue computations, we iteratively evaluate the posterior variance, and append the point of maximum variance to the ordered list of regression points. Fig. 3. Top row: samples generated from the posterior (predictive) distribution with 0, 1, 2, and 10 regression points sampled (note that when 0 points are sampled, the posterior is the prior). As more points are sampled, the posterior distribution becomes much tighter, as can be seen by the decrease in variation of samples drawn from the posterior distribution. Bottom row: the variance of the posterior distribution with 0, 1 2, and 10 points sampled. The first three panels of the bottom row demonstrate the first three sequential steps of the algorithm used to determine the order in which points are sampled.
With the mean, covariance, and acquisition functions determined, we can fit the GP model to a specific dispersion relation. To fit the model, we supply it with a set of regression points sampled from that dispersion relation. Our acquisition function determines the locations (wavevectors) of these points ahead of time. To perform the fit, we need to observe the eigenfrequency solutions associated with these wavevectors by solving the corresponding eigenvalue problem. Because the locations of the regression points have been determined ahead of time, these computations can be done in parallel. After observing regression points (1 . , 𝑀 . ), we compute the model weights ≨ by solving the linear system
where 1 . and 𝑀 . are the wavevectors and frequencies of the regression points, ℏ 2 𝜛 is the assumed variance of the noise on the regression point measurements, and ⋆ is the identity matrix. Because our regression points are generated by a well-chosen finite element model, it is reasonable to assume values of ℏ 2 𝜛 as low as 0, if assuming the eigenvalue solutions are exact, or 10 ε16 to represent the uncertainty associated with using double precision floating point data format. We find that using ℏ 2 𝜛 = 0 works well for our purposes, and does not produce significantly different behavior from using ℏ 2 𝜛 = 10 ε16 . Once the model is fit to the regression points, the predictive distribution can be computed (Eq. ( 15)). Because the predictive distribution is Gaussian, the maximum a posteriori estimate of the full dispersion relation is the mean of the predictive distribution. Further, the diagonal of the covariance matrix of the predictive distribution is the variance of each point being predicted. The variance gives us a metric for the model's predictive uncertainty at each wavevector in 1 ω . The predictive distribution is
where 1 ω are the wavevectors at which we would like to predict the frequency values of the dispersion band, and 𝑀 ω are the predicted frequencies. In the top row of Fig. 3, we draw samples from this predictive distribution to show its evolution as more regression points are used. When fewer regression points are used, the predictive distribution is quite broad, indicating uncertainty in the prediction. However, as more regression points are used, the predictive distribution becomes extremely tight, indicating much more certainty.
So far, we have detailed how to make predictions with our model. However, to use our model in gradient-based design optimization, we must be able to supply efficient design sensitivity computations. In this subsection, we describe how to compute gradients for frequency optimization and group velocity optimization.
Typically, computing the design sensitivity of the frequency at one wavevector requires not only the eigenvalue solution, but also the eigenvector solution (Eq. ( 4)). However, our model predicts only the eigenvalue, not the eigenvector. There are at least two ways to amend this.
First, we could simply use the surrogate model to infer the dispersion relation, and use this inference to determine the wavevector at which we need to optimize the frequency. Then we could solve the eigenvalue problem corresponding to that single wavevector, and compute the design sensitivity using Eq. ( 4). This is a great option if we only need to compute the design sensitivity of the frequency at a single wavevector, but do not know a priori at which wavevector (e.g. for bandgap optimization).
Alternatively, if the eigenvalue problems are particularly expensive, or if there are many wavevectors at which we desire to optimize the frequency, we can use a trick to avoid solving more eigenvalue problems. By employing the chain rule we can compute the gradient of any set of frequencies without performing additional eigenvalue computations. To do this, we first compute the gradient of the frequency at any wavevector with respect to the frequencies measured at the regression points. This gradient is defined by the GPR model. Then, since we do know the eigenvector at the regression points, we compute the gradient of the frequency at the regression point with respect to the design variables. This second gradient is defined by Eq. ( 4). We then multiply these two quantities to obtain the design sensitivity of any frequencies we are interested in optimizing. That is,
where
and
where 𝑁 . is a matrix containing the eigenvectors of the regression points, (⋛) ω denotes conjugate transpose, and the colon subscript indicates every element along the corresponding axis of the array (e.g. (𝑁 . ) ∂,𝜍 represents the 𝜍th column of 𝑁 . ). Note that although this method does not require further eigenvalue solutions, it does require traditional sensitivity computations for a potentially large number of regression points. In practice, these two costs should be weighed against each other.
Frequency is not the only subject of dispersion optimization. Group velocity is also a common optimization objective [12]. Computing the design sensitivity of group velocity again requires eigenvector information at the wavevector of interest. It should be noted that traditional sensitivity computations for group velocities are more expensive than traditional sensitivity computations for frequencies because they require solving an extra system of equations. As with frequency optimization, there are two methods that we can use which will be advantageous in different situations.
First, we can use the surrogate model to infer the dispersion relation, and identify a wavevector of interest. Then we can solve the eigenvalue problem at that wavevector, and use Eqs. ( 6), (7) to compute the group velocity design sensitivities.
Alternatively, we can avoid group velocity sensitivity computations altogether by employing the chain rule. Using the surrogate model, the analytical expression for the predicted group velocity, 𝑂 ω , is
where 𝑃 is the gradient of the covariance function with respect to its first input. Since the covariance function is fixed, this can be estimated numerically once, and stored for further use. Subsequently, we apply the chain rule to write the sensitivity of the predicted group velocity with respect to the design variables:
where
and
𝜑4 is again given by Eq. ( 18). Using the second method requires traditional frequency sensitivity computations for a potentially large number of regression points. However, it gracefully avoids the need for any traditional group velocity sensitivity computations, which require solving a rather large linear system. This second method is advantageous to use when we desire to optimize group velocity at a large number of wavevector points.
In this section, we apply our surrogate model to dispersion relations of two-dimensional materials and evaluate its performance against traditional methods and an existing deep-learning-based surrogate model.
When modeling materials in two dimensions, the evaluation of the dispersion relation becomes more complicated and more computationally costly than in one dimension. This is for two reasons. First, a two-dimensional material requires more degrees of freedom to be accurately modeled. Second, wave propagation in a 2-D material must be described by a wavevector, rather than a wavenumber, since wave propagation may occur in any combination of the two directions. This means that one must sweep over a 2-D parameter space of wavevectors, solving an eigenvalue problem at each point in this 2-D space. These issues are further exacerbated when working with 3-D materials.
For this reason, when considering 2-D and 3-D materials, engineers and scientists often evaluate the dispersion relation only on the IBZ contour or a ''virtual IBZ contour'' rather than the full IBZ. The definition of the IBZ contour is dependent on the symmetries of the unit cell, and is defined as the boundary of the IBZ for symmetry groups with reflection or glide-reflection symmetries. For symmetry groups without reflection or glide-reflection symmetries, the IBZ contour is not defined since the IBZ can be arbitrarily rotated, giving an infinite number of possible choices for the IBZ contour, and an infinite number of possible resulting dispersion band diagrams.
For example, Fig. 4 shows a unit cell lacking symmetry. Unit cells lacking symmetry belong to the 𝐹1 plane symmetry group (or wallpaper group). The IBZs of unit cells in the 𝐹1 symmetry group can be rotated arbitrarily, so a unique IBZ contour is not defined. As a substitute, the authors of [15] use a ''virtual IBZ contour''. Although the construction of this contour is entirely artificial, we include it because many people in the acoustics and phononics field are accustomed to seeing dispersion information in this format. Even for unit cells that have well-defined IBZ contours, the commonly used contour shortcut yields only partial information about wave propagation through the medium. The goal of our surrogate model is to increase the accessibility of dispersion information over the full IBZ.
In Fig. 4, we show the unit cell along with the first three bands of its dispersion relation evaluated along the ''virtual IBZ contour'', and over the full IBZ. We also show the inference of the surrogate model over these domains. ). The side length of the unit cell is ℷ = .01 m. This unit cell has no symmetry, making it part of symmetry group 𝐹1 [15]. The IBZ of this symmetry group can be arbitrarily rotated, so a well-defined IBZ contour does not exist [15]. However, since so many are accustomed to visualizing dispersion information in this format, we have included the dispersion relation evaluated on the ''virtual IBZ contour'' as defined in [15]. Here, we again show the true dispersion relation along with the inferred dispersion relation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In this section, we consider 2D unit cell designs of the 𝐹1 symmetry group (i.e. no symmetry) using a plane stress model. To generate random unit cell designs, we sample a Gaussian process with a 2-D periodic covariance function. To construct this covariance function, we start with a 1-D periodic covariance function,
where ℏ 7 controls the overall variance of the Gaussian process, ⋝ controls the lengthscale of the randomness, and 𝐹 is the period.
In our case, we set 𝐹 to be the side length, ℷ, of the unit cell. We form the 2-D periodic covariance function as the product of the 1-D periodic covariance function applied to each direction:
where 𝜗 and , are the horizontal and vertical coordinates of locations in the unit cell. We emphasize that the Gaussian process used to generate unit cell designs has no direct connection with the Gaussian process used by the surrogate model for inference. Beyond the periodic covariance parameters (which control the variance and the lengthscale of features in the unit cell), there are several other parameters that define the unit cell. The lattice parameter ℷ dictates the side length of the square unit cell. The resolution of the unit cell is determined by the integer > 𝐹𝜍𝜗 , which defines the number of pixels both vertically and horizontally. The constitutive properties of each pixel are bounded by 𝑆 8𝜍ℸ , 𝑆 8ℷ𝜗 , 𝜀 8𝜍ℸ , 𝜀 8ℷ𝜗 , 𝑇 8𝜍ℸ , and 𝑇 8ℷ𝜗 , corresponding to elastic modulus, mass density, and Poisson's ratio respectively. The number of unique materials in each design is controlled by > 8ℷ𝐾 (for example, a bimaterial unit cell would be defined by > 8ℷ𝐾 = 2). Further details of unit cell design generation are given in Appendix A.2. Examples of how these unit cells might look are provided in Section 6.6.
In this section, we visualize the covariance and acquisition functions for the two-dimensional problem. In Fig. 5, we show the evolution of the acquisition function and the resulting variance fields for a covariance function constructed from a dataset of 2000 dispersion relations. 1 As we compute the order of regression point sampling, we can observe that points that are not near to each other can be highly covariant. For example, the covariance function is taking advantage of the periodicity of the dispersion relation. In the example shown in Fig. 5, the unit cells in the dataset used to construct the covariance function belong to the p1 plane symmetry group (they have no symmetry), which is a superset of all other plane symmetry groups. Due to the construction of this symmetry group's IBZ and the periodicity of the dispersion relation, the solution along the boundary of the IBZ must be symmetric across the line ⊲ 𝜗 = 0. Indeed, we see from the (1, 1) and (2, 1) panels that a single regression point on the far left side of the IBZ causes a ''cold spot'' in the variance field on the right side of the IBZ in the region that is ''close'' to the regression point when considering the periodicity of the dispersion surface. This indicates that the model gains insight about the function value on the right side of the domain from information observed on the left. Our covariance function is able to implicitly learn crucial patterns such as periodicity and symmetry from the data, without the need for a user to manually encode these assumptions. Note that when the model uses 0 regression points, the posterior distribution is identical to the prior distribution. Also, we observe that as more regression points are used, the samples become more similar to each other, and to the true dispersion relation in the right column. Further, note that the posterior mean function is the final hypothesis of the surrogate model, not a randomly drawn posterior sample.
We now show examples of our surrogate model infering the dispersion relation for two-dimensional structures using a sparse set of regression points in Fig. 6. In each case, the model uses a dataset of 2000 dispersion relations (the same dataset used to generate Fig. 5) and 50 regression points (the same regression points shown in the bottom left panel of 5). The dataset was constructed with the following parameters: ℏ 7 = 1, ⋝ = 0.5, ℷ = 1 m, > 𝐹𝜍𝜗 = 16, > 8ℷ𝐾 = 2, 𝑆 8𝜍ℸ = 200 MPa, 𝑆 8ℷ𝜗 = 200 GPa, 𝜀 8𝜍ℸ = 800 kg+m 3 , 𝜀 8ℷ𝜗 = 8000 kg+m 3 , 𝑇 8𝜍ℸ = 0, 𝑇 8ℷ𝜗 = 0.5. The wavevector resolution of each dispersion relation in the dataset is 101 -51. In Fig. 6 we use 50 regression points to achieve a relatively accurate prediction; however, in practice, the number of regression points can be adjusted fluidly to strike the right balance of accuracy versus computational cost depending on the application.
In this section, to examine the relationship between the size of the covariance function training set and the expressivity of the model, we examine the spectra of covariance functions based on datasets of different sizes. Recall from Section 4.2 that the number of nonzero eigenvalues of the covariance function corresponds to the size of the covariance function's eigenbasis, which is the basis the model uses to fit to the regression points. A covariance function with a more comprehensive eigenbasis results in a model with more expressive power, and a model that can make use of a larger number of regression points.
For a range of logarithmically spaced dataset sizes, we compute the spectrum of the corresponding covariance function, and plot the spectra in Fig. 7. As the size of the dataset increases, the number of nonzero eigenvalues increases, meaning that the covariance function's eigenbasis is expanding to be more expressive. Eventually, however, the addition of more dispersion relations to the dataset has a decreasing effect on the spectrum of the covariance function. For the largest dataset size (20,000 dispersion relations) we plot the eigenfunctions corresponding to the 10 largest eigenvalues of the covariance function in Fig. 8. These eigenfunctions Because we plot in log scale, eigenvalues that are zero are not plotted. As expected from basic statistics, we observe that datasets with larger numbers of dispersion relations produce covariance matrices with more nonzero eigenvalues. This indicates that covariance functions generated from larger datasets will have larger eigenbases, and therefore have more expressive power. However, we observe that the change in spectrum appears to slow down at some point, and we get diminishing returns after about 1700 dispersion relations. We can expect this ''saturation'' effect to happen at different points in different scenarios, and will depend on many factors including the resolution of the dispersion relations in the datasets and the method of generating unit cell designs. are the first ten elements of the eigenbasis, and can be thought of as the ten most prominent and fundamental function shapes that the model uses to fit a dispersion relation.
The expressivity of the model (i.e. the broadness of the hypothesis space) depends on the dataset used to construct the covariance function. In general, larger datasets will allow for more expressivity. This is not surprising, since other data-driven methods such as neural networks require larger amounts of data to characterize more complex functions. However, there is a difference in how the model complexity is determined. For a neural network, the hypothesis space is determined by the network architecture, which is typically determined before training the model. If too little training data is supplied, the hypothesis space of the network remains large, but the model's performance will be low since the training data was not sufficient to narrow down the hypothesis space. In contrast, the expressivity of our model's data-driven covariance function is determined by the dataset itself. With our method of covariance function construction, a small dataset will automatically result in a smaller hypothesis space, and a larger dataset will automatically result in a larger hypothesis space, without the need to adjust complicated model hyperparameters (such as network architecture).
In the following subsections, we compare the performance of the GPR model against two other models in the context of the twodimensional problem. The first model is linear interpolation, representing the traditional model for computing dispersion relations. 2The second model is an existing convolutional neural network surrogate model from the literature [32].
To evaluate the performance of the GPR model, we consider the accuracy of the model as a function of the size of the covariance training set and the number of regression points for both in-distribution and out-of-distribution scenarios. We generate three datasets, one for training the covariance function, one for testing in-distribution performance, and one for testing out-of-distribution performance. The training dataset contains dispersion relations of 5000 unit cells generated as described in Section 6.
To measure the error of the model, we use the 𝐷 2 norm and the 𝑈 1 norm, defined in Eqs. (24), (25). The 𝐷 2 norm allows us to assess how well the model infers frequency values, and the 𝑈 1 norm to assess how well the model infers group velocity values. Errors are evaluated at the resolution of the test sets, which are about five times more resolved along each of the two wavevector directions, as described in Section 6.6.1. After evaluating model error on each eigenvalue band of each dispersion relation in the test set, we average over eigenvalue bands, and take the 0.95 quantile over dispersion relations to get a sense for worst-case error. The norms are given as
where 𝜑 ≨ denotes the multi-index partial differential operator, }≨} = ⌉ 𝜍 ≨ 𝜍 , and 𝜔 is some domain over which we are interested in measuring the function. In our case, 𝜔 will be the IBZ, and the function 7 will be the error of the model.
To compare against the linear interpolation model, we consider the number of regression points used by each model and the error that each model incurs in its prediction. In this subsection, because of the 1000shift in elastic modulus of the out-of-distribution test set, we divide the out-of-distribution errors by ⟩ 1000 so that they can be readily compared to the in-distribution errors. 3 We find that the GPR model always uses fewer regression points than the corresponding linear interpolation model achieving the same accuracy. At most, in this particular study, the GPR model uses 7.75 times fewer regression points than the equivalent linear interpolation model. Because almost all of the computational cost is due to regression point eigenvalue computations, this can be seen as a speed up of approximately the same factor. The results are shown below in Fig. 10.
A comparison of the GPR model against the traditional linear interpolation model. Performance is evaluated on two datasets: in-distribution and outof-distribution. In both cases, the GPR model performs substantially better than the linear interpolation model. In the 𝑈 1 plot for in-distribution performance, the GPR model is shown to achieve the same error as the linear interpolation model while using up to 7.75fewer regression points.
We compare our model performance to the performance of the CNN surrogate model proposed by Finol et al. [32]. Please refer to Appendix B for details about the CNN architecture and training. Note that the CNN uses the material properties of the unit cell in the input array. This is in contrast to the GPR model, which does not directly observe the material properties of the unit cell, but instead observes regression points of the dispersion relation.
In this subsection, we allow the GPR models to observe only 50 regression points, and limit both the CNN and GPR models training set sizes to range from only 50 to 5000 dispersion relations. For each training set size, we train the CNN for a maximum of 30 epochs with early stopping based on a validation set. Training the CNN on the training set of 5000 dispersion relations takes approximately 24 minutes/epoch on an NVIDIA GeForce RTX 3080 Ti GPU. For each training set size, we also compute a new covariance function for the GPR, taking only 1 s for the largest training set size (5000). Because training a neural network depends on the randomly initialized weights of the network, we repeat the training 10 times with different random initializations to get a sense for the variation in results. After measuring the 𝐷 2 and 𝑈 1 norms of the error, we average over all 10 eigenvalue bands, and 3 This exactly transforms the errors to what their values would be if the modulus range were 1000smaller, since the dispersion relation of a unit cell whose modulus is increased by 1000will experience an upshift in frequency at a factor of ⟩ 1000, and because both the linear interpolation and GPR models' predictions (and therefore errors) are linearly proportional to the frequency values of the regression points.
take the 95% quantile error over unit cells. Additionally, for the CNN, we measure the standard deviation of this error over the ten trials and plot the ± one standard deviation spread. The results (Fig. 11) show that the GPR model achieves errors significantly lower than the CNN for all training set sizes. This is because we are using relatively small training sets, and neural networks typically require more data to train effectively. For reference, Finol et al. used a minimum of 20,000 and up to 100,000 dispersion relations to train their CNN. In general, it is a nontrivial task to adjust the architecture (and therefore expressivity) of a neural network to be appropriate for the amount of training data that you have. This process generally involves a substantial amount of trial-and-error and hyperparameter tuning, which can be costly in terms of time and effort. Note that in this section, we do not adjust the out-of-distribution errors by the factor of ⟩ 1000 because this transformation would not be valid for the CNN, so in-distribution and out-of-distribution results cannot be directly compared against each other. For both in-distribution and out-of-distribution, we see the CNN gradually improve with training set size in the 𝐷 2 norm, but not in the 𝑈 1 norm. This makes sense, because the loss function that the CNN uses (mean squared error) is not sensitive to mismatches in gradients (𝑈 1 error).
Interest in the design of multi-functional wave devices and the high computational cost of dispersion relation computations motivate the need for more efficient methods of dispersion relation evaluation, including surrogate models. We have introduced a new GPR-based surrogate model for the accelerated computation of dispersion relations over the full irreducible Brillouin zone.
The GPR model builds its covariance function from data and is able to learn from datasets that are 2-3 orders of magnitude smaller than datasets required by the demonstrated deep learning approaches to date. Instead of involving a typical deep learning training process, the model uses only a simple covariance calculation to capture the data trends. The low requirement for data, along with the avoidance of a typical training process, massively reduces the overhead cost to set up the model, allowing the model to be more easily repurposed for new design tasks.
Because the model never directly observes the parameters governing the unit cell design, it can seamlessly train on unit cells represented in one format (e.g. pixelated/voxelated) and subsequently generalize its predictions to unit cells of other representations (e.g. parametric). This inherent quality makes the model versatile, emphasizing its potential to be effectively applied across a wide spectrum of design scenarios.
Finally, the model can supply analytically computed gradients for use with gradient-based design optimization, along with uncertainty estimates -both highly useful tools for design tasks. For optimization problems where deviation from the solution of the variational method cannot be tolerated, the surrogate model can simply be used for the first several iterations of optimization and then handed-off to a slower, but more exact dispersion evaluation method. Further, the model can work in tandem with existing, highly specialized variational methods for dispersion computations. We hope that by increasing the computational accessibility of where ω denotes conjugate transpose.
By taking a gradient with respect to a generic variable 4, we obtain Eq. (C.3). For now, consider 4 to be a scalar. To generalize to vector 4 (possibly as a design vector or a wavevector), simply repeat these computations for each entry of 4. The product rule is used here.
There are two things that can be done with this equation. First, if we only want the gradient of the eigenvalue with respect to a design vector, we can left-multiply Eq. (C.3) by 𝜚 ω . Then, applying Eq. (C.2), we can see that
Alternatively, if the gradient of the eigenvector 𝜚 is also desired, Eq. (C.3) can be rearranged to
If 4 ϖ R and 2, 3 ϖ R ℸ , then this system has ℸ equations but 2ℸ unknowns (or in general if 4 ϖ R 𝐶 , and 2, 3 ϖ R ℸ , then this system has ℸ𝐶 equations but ℸ(𝐶 + 1) unknowns). This is because eigenvectors are only determined up to their scale. To provide the extra ℸ equations required to solve, we can use an eigenvector normalization equation. It is standard to normalize eigenvectors by their mass matrix.
By differentiating, we arrive at
If our eigenvector 𝜚 were purely real, we could say:
This would make Eq. (C.7) easy to rewrite so that the unknowns on the LHS only receive operations from the left. Then we could solve
However, in this case, our Bloch-Floquet periodic boundary conditions cause our eigenvector 𝜚 to be complex, so in general
𝜑4 ) ω 3𝜚) ω = 𝜚 ω 3 𝜑𝜚 𝜑4 (C.9) and therefore ( 𝜑𝜚 𝜑4 ) ω 3𝜚 = (𝜚 ω 3 𝜑𝜚 𝜑4 ) ω (C.10) By substituting Eq. (C.10) into Eq. (C.7), we can arrive at (𝜚 ω 3 𝜑𝜚 𝜑4 ) ω + 𝜚 ω 3 𝜑𝜚 𝜑4 = ε𝜚 ω 𝜑3 𝜑4 𝜚 (C.11) In other words, 2 Re ⟨ 𝜚 ω 3 𝜑𝜚 𝜑4 ⟩ = ε𝜚 ω 𝜑3 𝜑4 𝜚 (C.12) And there is no constraint on the imaginary component. In other words, Im ⟨ 𝜚 ω 3 𝜑𝜚 𝜑4 ⟩ = anything (C.13)
To solve this, we can break up the imaginary and real equations, and write them individually. This turns 3 complex equations into 5 real equations. The reason we do not obtain 6 real equations from the 3 complex equations (as one may expect) is that Eq. (C.11) actually only provides one real equation (note that of Eqs. (C.12), (C.13), only Eq. (C.12) tells us something useful). Luckily, because our eigenvectors are complex, they not only need a constraint in magnitude, but also in complex angle. Note that if 𝜚 is an eigenvector, then 𝜚.⊳ 𝜍𝑊 is also an eigenvector for any ., 𝑊 ϖ R. In other words, the eigenvector can be multiplied by any complex scalar, and still be an eigenvector. For this reason, we not only need to constrain the magnitude of 𝜚 (as we have done with Eq. (C.6)), we also need to constrain the complex angle of the entries of 𝜚. To do this, we can simply enforce that the first entry of 𝜚 is purely real. This can be written as )
)
)
Now, let us say we want to optimize the group velocity of a wave. The group velocity is the gradient of the frequency, 𝜛, with respect to the wavevector, ⊲. As with 4, it is helpful to think of ⊲ as a scalar initially, then generalize to vector ⊲ by repeating the computation for each component of the wavevector. First, we know from Eq. (C.4) that the gradient of the eigenvalue, 𝜛 2 , with respect to the wavevector, ⊲, for a given branch and wavevector can be determined by
Note that Eq. (C.17) gives 𝜑𝜛 2 𝜑⊲ , and does not directly give the group velocity, which is 𝜑𝜛 𝜑⊲ . To obtain the group velocity, we can use the chain rule. Consider two functions 7 (⊲) and 𝑃(⊲), such that 7 (⊲) = 𝑃(⊲) 2 . 𝜑𝜛 𝜑⊲ = 1 2𝜛 𝜑6 𝜑⊲ (C.19) Using Eqs. (C.17), (C.19), we can write an analytical equation to compute the group velocity. 𝜑𝜛 𝜑⊲ = 1 2𝜛 𝜚 ω ( 𝜑2 𝜑⊲ ε 6 𝜑3 𝜑⊲ )𝜚 (C.20) Note that this same equation can be used to compute the design sensitivity of the eigenfrequency, 𝜛, (in contrast to Eq. (C.4), which computes the gradient of the eigenvalue, 6). 𝜑𝜛 𝜑4 = 1 2𝜛 𝜚 ω ( 𝜑2 𝜑4 ε 6 𝜑3 𝜑4 )𝜚 (C.21)
By taking the gradient of Eq. (C.17) with respect to a design vector, we can aim to obtain the design sensitivity that would be used for a gradient-based optimization of group velocity.
First, to simplify notation, define the matrix 5 as follows.
5 = 𝜑2 𝜑⊲ ε 6 𝜑3 𝜑⊲ (C.22) Substituting into Eq. (C.20) shows 𝜑𝜛 𝜑⊲ = 1 2𝜛 𝜚 ω 5𝜚 (C.23) Next, take the gradient with respect to a design vector, 4. 𝜑 𝜑4 ⌋ 𝜑𝜛 𝜑⊲ { = 𝜑 2 𝜛 𝜑4𝜑⊲ = 𝜑 𝜑4 ⟨ 1 2𝜛 ⟩ 𝜚 ω 5𝜚 + 1 2𝜛 𝜑𝜚 𝜑4 ω 5𝜚 + 1 2𝜛 𝜚 ω 𝜑5 𝜑4 𝜚 + 1 2𝜛 𝜚 ω 5 𝜑𝜚 𝜑4 (C.24)
Next, let us note how to find some of the quantities needed to use Eq. (C.24).
• 𝜑𝜚 𝜑4 can be found by solving the linear system given by Eq. (C.16).
As for computing the derivatives of the system matrices 2, 3, we can recall that they are constructed by applying wavevectordependent transformation matrices to design-dependent full system matrices 2 7 , 3 7 . 2 = 𝐸 ω 2 7 𝐸 3 = 𝐸 ω 3 7 𝐸 (C.26) 𝜑2 𝜑⊲ = 𝜑𝐸 𝜑⊲ ω 2 7 𝐸 + 𝐸 ω 2 7 𝜑𝐸 𝜑⊲ 𝜑3 𝜑⊲ = 𝜑𝐸 𝜑⊲ ω 3 7 𝐸 + 𝐸 ω 3 7 𝜑𝐸 𝜑⊲ (C.27) 𝜑2 𝜑4 = 𝐸 ω 𝜑2 7 𝜑4 𝐸 𝜑3 𝜑4 = 𝐸 ω 𝜑3 7 𝜑4 𝐸 (C.28) 𝜑 2 2 𝜑4𝜑⊲ = 𝜑𝐸 𝜑⊲ ω 𝜑2 7 𝜑4 𝐸 + 𝐸 ω 𝜑2 7 𝜑4 𝜑𝐸 𝜑⊲ 𝜑 2 3 𝜑4𝜑⊲ = 𝜑𝐸 𝜑⊲ ω 𝜑3 7 𝜑4 𝐸 + 𝐸 ω 𝜑3 7 𝜑4 𝜑𝐸 𝜑⊲ (C.29)
The dataset was constructed with the following parameters: ℏ 7 = 1, ⋝ = 0.5, ℷ = 1 m, > 𝐹𝜍𝜗 = 16, > 8ℷ𝐾 =
2, 𝑆 8𝜍ℸ = 200 MPa, 𝑆 8ℷ𝜗 = 200 GPa, 𝜀 8𝜍ℸ = 800 kg+m3 , 𝜀 8ℷ𝜗 = 8000 kg+m3 , 𝑇 8𝜍ℸ = 0, 𝑇 8ℷ𝜗 = 0.5. The wavevector resolution of each dispersion relation in the dataset is 101 -51.
For the traditional linear interpolation model, one simply evaluates the dispersion relation on a grid of evenly spaced points and linearly interpolates if evaluation is needed at points that do not belong to the grid.
For example, if a dataset of dispersion relations has a resolution in reciprocal space of 201 -201, the covariance matrix is of size 40 401 -40 401, occupying 13.06 GB of memory when stored in double precision. This means that if the covariance between even two bands at this resolution is considered, the resulting covariance matrix will be 80 802 -80 802, occupying
52.24 GB of memory.