Transfer learning for metamaterial design and simulation

1.1 Contributions of this work

In this work we demonstrate the tremendous potential benefits of transfer learning for AEM tasks, motivating its broader adoption in the AEM community. To study the benefits of transfer learning, we use it to develop surrogate DNN-based models for the multi-scale problem of scattering from electrically large random metasurface arrays – a problem commonly encountered in applications such as metamaterial lens design [28], [29], disordered metasurfaces [30], [31], and Huygens metasurfaces [32]. We begin by generating a large number of simulations for a source task, and then explore several different target tasks with varying problem features – the base material, geometric shape, incident angle, and polarization of the excitation. We evaluate the accuracy of DNN-based surrogate models on each target task as a function of the quantity of data used to train the model, and perform a direct quantitative comparison to the accuracy of source tasks with and without transfer learning. Our results indicate that, using transfer learning, we can often dramatically reduce the quantity of training data that is needed to achieve a desired prediction accuracy for new problems, e.g., often over an 80 % reduction, and as high as 90 %. Crucially, we also find that transfer learning is never detrimental over our diverse target problem settings, suggesting that there is little risk in utilizing it.

To enable the benefits of transfer learning, we propose a geometric parameterization that makes transfer learning straightforward across all of the proposed tasks, and allows us to leverage existing successful DNN-based models from the computer vision literature. We therefore represent the geometric structure of the AEMs as a 2D array – similar to an image – where each entry in the array corresponds to a geometric parameter of the unit-cell, as illustrated in Figure 1. This structure allows us to adopt convolutional DNN model architectures from the computer vision community, which can take advantage of the spatial structure of the underlying metasurface when making predictions. This parameterization also means that little or no modification is required to utilize pre-trained model parameters from the source task, making it easy to use pre-trained weights for each target task.

This paper is organized as follows. First, we present the analysis of using a residual neural network structure to solve a base metamaterial array scattering problem, which serves as the source problem from which we transfer knowledge. Second, we briefly introduce the discrete dipole approximation (DDA) method for efficient data generation. Next, we explore five target tasks involving various condition changes, including element material, shape, and excitation. Finally, we perform a quantitative analysis on the effectiveness of transfer learning in terms of error reduction as a function of dataset size, as well as data reduction for given levels of accuracy. For each target task we also make a comparison between the transfer learning approach and for the case where we use random initialization and do not use any transfer knowledge.

2.1 Array scattering problem settings

We aim to solve the scattering problem for different random arrays under similar, yet distinct, physical conditions. Specifically, one class of problems, or equivalently one task, examines scattering from a random metasurface array consisting of elements with distinct morphology but with fixed periodicity and excitation fields. Another problem class involves array scattering from the base problem, but with different excitation fields. In each task instance, the input variables are the geometric parameters of different random arrays, denoted as g, and the output variables are the far-field patterns, i.e., the electromagnetic responses of the array, denoted as s.

2.2 Deep convolutional neural network

The scattering patterns from the metamaterial array are determined by the specific spatial arrangement of single elements in the array subject to different excitations and other environmental conditions. Therefore to make our method as general as possible, i.e. transferable to as many different types of metamaterials and physical scenarios as possible, we represent the geometry space g as a numeric matrix. This numeric matrix g is of the size C × N × M, where C represents the total number of tunable geometric parameters for each single element, e.g. some critical length or radius, and N, M are the 2D array sizes. Therefore, large arrays of metamaterial elements can be characterized in a compact and scalable fashion even if the elements are three-dimensional with fine features. The use of a CNN also permits preservation of the spatial relationship between elements in the input numeric matrix, thereby allowing the convolutional filters to learn element interactions.

While the geometric layout can be readily represented using matrices, the random array far-field radiation pattern is generally a function of 3D spherical angles, as depicted in Figure 1. For simplicity, and without loss of generality, we consider only one representative slice in the 3D radiation pattern along a constant forward scattered azimuthal angle. The output far-field patterns |E| are sampled along the forward scattering direction with varying elevation (polar) angles from 0 (North Pole) to 180 (South Pole) degrees, represented by a 181 × 1 positive number vector, so that s ∈ R + , 181 .

Due to the large magnitude variations of the far-field patterns, we formulated a mixed linear-log loss, which we found to be more effective, in comparison to the traditional MSE loss (applied in either the linear or log scale.) We define our training loss as:

where h(x) = ln(1 + x). The log-scale term is used to decrease the loss for the scattering pattern expressed in a dB scale. We found for the hyperparameter α that α = 5 is a suitable choice for the far-field scattering problem in this paper, which balances the loss between the linear and logarithm scales. For the test set we use the mean absolute relative error (MARE), defined as:

This metric ensures that the reported loss is independent of the absolute magnitude level of the far-field signal, making it suitable for comparisons across different problems.

CNNs can be constructed with different architectures, which refer to the order in which data is processed (e.g., its connectivity) and the specific processing that takes place. In this work, we use a residual convolutional neural network (ResNet), which has been found highly effective for computer vision problems [33] by exploiting the local information between close neighbors – similar to the input data we consider here. Using a 2-dimensional matrix encoding of our material design, we can leverage state-of-the-art CNN architectures from the computer vision literature. Furthermore, we hypothesize that the 2-D ResNet convolution operations resemble the far-field inference procedure employed for Green’s functions, and therefore, may be well-suited for modeling far-field radiation patterns. Our ResNet consists of one base convolution layer, followed by six residual blocks with (64, 64, 128, 128, 256, 256) convolution channels. Each residual block in the ResNet architecture consists of two identically sized convolutional layers connected by a skip connection. Each convolutional (Conv) layer is followed by a batch normalization (BN) layer and a ReLU activation. The base convolution kernel has a size of 7 × 7, and the remaining convolutional layers have a kernel size of 3 × 3. These residual blocks are followed by one max pooling layer of 2 × 2, and one fully-connected layer to yield the 181 × 1 vector output. If sufficient training data are collected and used to train a CNN model such that s = f ̂ ( g ; Θ ) predicts accurate responses, the f ̂ surrogate network model can be generalized to solve the scattering problem of any new random array within the training domain, which is not in the training dataset, without expensive full-wave computation. In this context, the CNN model serves as a fast solver for the random array scattering problem. Note that this trained CNN surrogate model f ̂ is tied to the problem with all the physics conditions fixed, (the element material, shape, and incident fields), and all implicitly embedded in the network parameters.

2.3 Transfer learning using fine-tuning

As mentioned, the trained CNN model is accurate only for the dataset on which it was trained, and therefore cannot be used directly to solve a new task. An entirely new dataset D _tgt for the target problem needs to be generated, and a new CNN must be trained using D _tgt . However, if the new task shares some common underlying physics, transfer learning can leverage the information incorporated in the trained CNN for the new task. We begin with an initial modeling problem, termed the source problem, for which a large training dataset exists, denoted D _src . We then train a ResNet model on D _src , resulting in a trained CNN with weights θ _src that achieve highly accurate predictions for (g, s) ∈ D _src . In the target problem, denoted D _tgt , we vary some aspects of the experimental scenario of the base problem, such as the geometry or excitation. In the traditional approach, we would need to collect a large quantity of data again to obtain accurate model parameters for D _tgt ; however, we hypothesize that initializing the gradient descent training procedure with the model parameters θ _src will lead to improved performance compared to randomly initializing the model parameters – the procedure typically used when training CNNs for a new task.

2.4 Discrete dipole approximation method for scattering

To generate the dataset in a practically realizable time, a semi-analytic method called the discrete dipole approximation (DDA) is used as an alternative to commercial full-wave simulation software. A finite-sized metasurface array may consist of dielectric or metallic elements and may be of large electrical size. It is computationally expensive or even impossible to simulate the entire array depending on the electrical length and desired wavelength range, which specifies the grid size and therefore the total problem size. The local slowly-varying approximation (SVA) is usually applied to these types of array problems by solving for the electromagnetic response of a single metamaterial element with periodic boundary conditions, thereby modeling an infinite identical array [34]. The local periodic response is then used to approximate the actual response for the element when embedded in a random array. The SVA is valid if the response fields from element to element are smooth, and the change in the coupling between elements, due to the random arrangement, is negligible. However, this method may result in poor accuracy since the local response of one element can be drastically different from the periodic response of the same element with different neighboring elements. The problem we explore here cannot be studied with SVA and therefore we use the polarizability retrieval method together with the discrete dipole approximation method to extract the response of the elements embedded in a random array. The DDA method, along with the retrieval method, enables accurate and fast large-area metasurface simulation.

Many sub-wavelength metasurfaces can be accurately approximated by an electric dipole p ̄ and a magnetic dipole m ̄ . The polarizability tensor α ̄ ̄ relates p ̄ , m ̄ with the local field E, H by

where for example,

The polarizability tensor α ̄ ̄ can accurately describe both the material parameters and geometric features of a metamaterial element in a periodic environment. Therefore we may build up a library of α ̄ ̄ ’s for a series of desired metamaterial unit-cells which we intend to use for a given problem, as depicted in Figure 2. The polarizability retrieval method entails performing several scattering parameters (S-parameters) simulations to retrieve α ̄ ̄ [35]. A 4 × 4 S-parameter matrix, corresponding to 2 ports each with 2 orthogonal polarization states (h,v), of a single metamaterial element with periodic boundary conditions is found through full wave simulations. Note that the full wave simulation is only used in this step for each given single element. The inverse retrieval operator F is then applied on the S-parameters to get the 4 × 4 in-plane sub-tensor of the full 6 × 6 polarizability tensor, given by α ̄ ̄ x y = F [ S ] 4 × 4 . To retrieve the full 6 × 6 polarizability tensor, one or two additional simulations may be required where the element is rotated about a principle axis. This technique is applied for every metamaterial unit-cell geometry and frequency of interest in advance to generate the {geometry; polarizability} library which can then be configured in any arrangement to approximate an array of actual metamaterial elements.

Figure 2:

Flowchart of the DDA solver and the neural network proxy. (a) Single element simulations are used to obtain a library of polarizabilities for different elements. (b) Using the polarizability library, an array of elements can be mathematically replaced by a matrix of the corresponding polarizabilities. (c) The scattering from a complex array is reduced to a system of linear equations, where the coefficients are determined by the exact elements and the environment. The right-hand side of the boxed equation is the excitation E i inc . Any problem change can be directly modified in this system of linear equations.

With a polarizability tensor determined for every desired metamaterial element, the random metasurface array scattering problem is then reduced to solving a linearized integral equation system which, for an arbitrary jth metamaterial element in the array, is described by,

where E ₀, H ₀ are the incident electric and magnetic fields at the jth element location, and G ̄ ̄ is the Green’s function for the environment. The superscripts on the Green’s function stand for electric field (e), magnetic field (h), electric dipole (p), and magnetic dipole (m). The above scattering equation holds for every element in the metamaterial array, leading to a linear equation system. Then Eq. (5) can be solved by multiple standard methods, for instance, generalized minimal residual method (GMRES) to obtain p ̄ j , m ̄ j for every element in the array.

The far-field pattern is the coherent summation of the radiated fields from all the electric and magnetic dipoles with the far-field approximation, i.e. E ( θ , ϕ ) = ∑ i E i p + E i m . For an individual electric dipole the free space far-field radiation is given by,

and for a magnetic dipole by,

where Z 0 = μ 0 / ϵ 0 , ϵ ₀ and μ ₀ are the permittivity and permeability of free-space, respectively, n ̂ is the surface normal of the array, and k is the magnitude of the free-space wavevector. We note that in the case where dipoles are embedded within a substrate, the radiation is within a dielectric and therefore possesses an analytically solution for the case of planar stratified substrates [36].

3.1 Tasks and datasets

All six problems are shown in Table 1, (including the original base problem), and consist of the computed scattered far-fields of a 2D metamaterial array of 50 × 50 randomly chosen elements in free space. As mentioned, each problem differs in either the type of radiation illuminating the array, or in the material and geometry of the array elements. Different instances of these problems, namely the samples in the dataset, are of different individual element sizes with other conditions fixed. The base problem involves calculating the scattered far-field of an array of silicon spheres under a TE polarized plane wave incident at 13° from the surface normal at a wavelength of 11 μm. The dataset for each problem consists of 50,000 training data in total and another 5000 testing data, all generated by the DDA method described here. The geometry parameters for each metamaterial element in the array are randomly drawn from a uniform distribution spanning from 0.45 μm to 1.1 μm, which represents the radius for a sphere array or the half length of an array consisting of metallic cross elements. We have chosen transfer tasks which we expect will have varying similarity to the base problem, and therefore which will yield varying degrees of success. The first three target problems shown in Table 1 each keep constant the metasurface geometry and material, but allow the state of external radiation to vary. In tasks 4 and 5 the base material is changed to PEC and aluminum (Al), respectively, and in task 5 both the material and shape of the metasurface element changes, thus representing a two-factor change from the base problem.

In Figure 3(a) we show a schematic of the ResNet used in our study. The geometry of the metasurface array is represented here as a colormap matrix with the colorbar showing the radius of each sphere in the array. The matrix then feeds into the ResNet, which we depict here in AlexNet style [37]. As mentioned, the ResNet is trained on a dataset obtained from our DDA solver and thus the output predicts the far field scattering. To validate the accuracy of our DDA modeling approach we perform a direct comparison to a commercial computational electromagnetic solver. However, CEMS is limited in the maximum array size which can be studied; therefore here we reduce the metasurface array size to 30 × 30. On the right side of Figure 3(a) we show a slice of the far-field radiation taken along ϕ = 180° and for polar angles from 0° ≤ θ ≤ 180°, where the ground truth is shown as a dashed blue curve and the ResNet prediction is shown as the solid red curve.

Figure 3:

Dataset and random training results: (a) the base problem is the scattering resulting from a 50 × 50 array, consisting of silicon (ϵ _r = 11.3) spheres. The spheres are grouped in a 5 × 5 subarrays with identical radii. The input matrices are shown with the radii ranging from 0.45 µm to 1.05 µm. The forward scattering slice in the farfield pattern is plotted in polar coordinates in dB scale in the original unit V/m of E field. This slice is at the azimuthal angle of 180° and the elevation angle is ranging from 0 to 180°. Two random sample points from the test dataset are shown for all 6 problems. Ground truths are computed using the DDA method and compared to the output from a trained model using a total of 50,000 data. (b) All the trainable parameters are transferred from the trained base network for the initialization of the neural network for any new problem. Random networks represent the baseline where all the neural network parameters are randomly initialized. (c) The farfield slices of samples from all the five target problems are shown.

Figure 3(b) presents a schematic depicting the transfer learning approach used for our study. The top AlexNet diagram in Figure 3(b) represents the base problem from which we transfer the parameters of the trained ResNet to the new problem under study – middle AlexNet diagram shown. We also compare the transfer learning approach to training a Resnet from scratch, which is shown in the bottom of Figure 3(b). In Figure 3(c) we compare the output of our DDA trained ResNet to computational simulations obtained from CST Microwave Studio for all target problems listed in Table 1. Here the network predictions are shown as solid red curves, while the ground truth simulations are shown as the dashed blue curves. All the neural network predictions are from the best model trained from random intializations using a total of 50,000 data.

3.2 Comparison of training from random initialization and transfer learning

To evaluate the effectiveness of transfer learning approach for the aforementioned five target tasks, two training schemes are implemented and compared for each problem. The first approach involves training a neural network with randomly initialized parameters, while the second is the transfer learning approach which involves fine-tuning a neural network initialized with the best model from the base problem – see Table 1. In the experiments conducted, the transferred neural network parameters are obtained from the best training results of the base problem, which used a total of 50,000 data points. Identical datasets of varying sizes are used for both training schemes, ranging from {500} to {1, 000} with a step size of 100, {1, 000} to {10, 000} with a step size of 1000, and {10, 000} to {50, 000} with a step size of 5000. The optimal model for each case is determined by searching over a grid of learning rates ([0.1, 0.01, 0.001, 0.0001]) and regularizers ([0.01, 0.001, 0.0001, 0.00001]). Each experiment is trained five times using the best hyperparameter combination to account for randomness. The resulting best models for both training schemes at different dataset sizes are then evaluated using a test dataset of size 5000. Although the mixed loss shown in Eq. (1) is employed during the training phase, we use a linear MARE loss in the testing phase, and for reporting and visualization purposes. Moreover, the MARE test loss is a normalized measure which enables a comparison among different problem types.

Figure 4 shows the test MARE of the transfer learning (orange) and random (blue) models as a function of data set size for each target problem. Here the solid curve is the average of the five aforementioned experiments, and the associated shaded area represents the 95 % confidence interval. We find that the transfer learning scheme consistently results in lower test loss for all target problems compared to random initialization. Further, the test MARE difference between transfer learning and random is a maximum for low dataset sizes which gradually diminishes for increasing dataset size. Notably, we do find that the advantage of transfer learning varies among the five target tasks, which can be intuitively observed through the gap between the two loss curves. The benefits of transfer learning are good for the two angle variation cases shown in Figure 4(b) and (c), while the TM polarization case shows the most significant improvement.

Figure 4:

Results of transfer learning and training from scratch for all target problems. (a) The training curves of transferred models and random models using 25,000 and 50,000 new data for the new problem. (b)–(f) The averaged test error of the transferred models and the random models as a function of dataset size. Five repetitive training are used to compute the mean error (solid line) and the 95 % confidence interval (shaded regions).

3.3 Quantitative analysis of the advantage of transfer learning

We comprehensively examine the advantages of the transfer learning approach by comparing its performance using two metrics. In the first case we measure the relative reduction in error between the transfer learning case and the random case, assuming a fixed quantity of target-task training data is available – denoted r _e . We also examine the relative reduction in training data required in the transfer learning case to achieve a desired level of prediction error, as compared to the random case – denoted r _D . To compare the relationship between the dataset size and error, we fit the data from Figure 4 using a third order polynomial to obtain a smooth e = g(|D|), where |D| denotes the size of a dataset D. Using the fitted functions both r _e and r _D can be calculated and results are shown in Figure 5. The error reduction ratio is a vertical cut of the data presented in Figure 4(b)–(f) and defined as,

where e ^TL is the test MARE using transfer learning at a dataset size of D, and e ^RI is the test MARE using random initializations also for D. In Figure 5(a) we see that for small data sizes r _e is optimal and, for example at D = 300, we find that r _e ≤ 10 % for Tasks 3–5, and r _e ≈ 65 % for Tasks 1 and 2 – see Table 1. For all dataset sizes explored we find that transfer learning always leads to improved accuracy on the test set in comparison to training from random initialization, i.e. r _e is less than unity. However, we find that the advantage of transfer learning gradually diminishes for increasing dataset size, and as we approach D ≈ 50, 000 we find r _e → 1.

Figure 5:

Gain of transfer learning: (a) the ratio of the test loss error of transferred training and that of training from scratch using the same amount of data. (b) The reduction in required dataset size using transfer learning to reach a fixed level of accuracy.

The data reduction ratio r _D explores the data required for the transfer learning approach to achieve a specific test MARE compared to the data for random initialization to achieve the same test MARE. Therefore r _D are horizontal cuts of the data presented in Figure 4(b)–(f), and defined as,

where D ^TL is the dataset size for transfer learning and D ^RI is the dataset size for random initialization. In Figure 5(b) we plot r _D as a function of new sample simulations. If we take a test MARE of 3 % as a point of comparison, we find that the data can be reduced by over a factor of 2000 compared to random initialization for the TM Task 3 case, whereas we find r _D values ranging from 94 % to approximately 22 % for the PEC Sphere (Task 4) to the 0° (Task 1), respectively. For all Tasks explored we find that a significant reduction in data is achieved for transfer learning compared to random initialization. The exact numbers of dataset sizes or MARE are compared in Table 2 at two critical metrics, i.e. MARE = 0.03 and data of D = 1, 000 points. For this comparison, the data of Figure 4(b)–(f) have been interpolated and, as observed from the data presented in Table 2, we find that to achieve a test MARE of 3 %, various tasks require drastically different amounts of data. For example, Task 3 needs only D = 3, whereas Tasks (4, 5, 2, 1) required D = (776, 924, 13691, 17,353). We note that for TM Task 3 case a third order polynomial fit – rather than the interpolated data – is used for Table 2, as denoted by the asterisk.

The large variance in data required to achieve a certain error level presented in Table 2 hints that some task problems share more or less underlying physics with the base problem. To elucidate this connection, we perform a similarity analysis on the output far field scattering for each of the transfer Task problems. Figure 6(a) shows the results of a principal component analysis (PCA) performed on the farfield scattering patterns of 181 by 1 data points. The analysis utilized a total 60,000 data points, with each Task – including the base Task – contributing 10,000. As evident in the feature-wise basis, there is some overlap between the TM Task 3 (orange) and Base Task 0 (brown) data points in the principal components’ space. Other explored target Tasks (blue, purple, green and red) have centroids of varying distance to the Base Task. However, it is crucial to remember that the PCA operates on the output (i.e., the farfield scattering patterns). Therefore, close proximity in the output space (codomain) does not guarantee the same in the input space (domain). In other words, the two overlapping points from the TM and base problems may not correspond to the same input geometry features. None-the-less we find that the centroid distance of the Target task clusters in Figure 6(a) correlates relatively well with the critical transfer learning metrics presented in Table 2. As shown in Figure 6(b) we plot the PCA centroid distance and r _D to reach a test MARE = 0.03 (open circle symbols) and r _e for D = 1000 (open triangle symbols).

Figure 6:

Simularity analysis: (a) PCA representation of the far field scattering data from the six problems. (b) The correlation of the PCA centroid distance from (a) and data reduction ratio (circle symbols) at 3 % MARE and the error ratio (triangle symbols) at D = 1000 for each of the target tasks.