{'fraction_non_alphanumeric': 0.02805594287404403, 'fraction_numerical': 0.010309713947690877, 'mean_word_length': 4.78114313629702, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hyperparameter Search using Genetic Algorithm for Surrogate Modeling of Geophysical Flows

The computational models for geophysical flows are computationally very expensive to employ in multi-query tasks such as data assimilation, uncertainty quantification, and hence surrogate models sought to alleviate the computational burden associated with these full order models. Researchers have started applying machine learning algorithms, particularly neural networks, to build data-driven surrogate models for geophysical flows. The performance of the neural network highly relies upon its architecture design and selection of other hyperparameters. These neural networks are usually manually designed through trial and error to maximize their performance. This often requires domain knowledge of the neural network as well as the problems of interest. This limitation can be addressed by using an evolutionary algorithm to automatically design architecture and select optimal hyperparameters of the neural network. In this paper, we apply the genetic algorithm to effectively design the long short-term memory (LSTM) neural network to build the non-intrusive reduced order model of the sea-surface temperature field.


I. Introduction
There are numerous tasks in engineering design, control, and climate modeling that require solving a partial differential equation (PDE) based forward model for many ensembles. For example, in sequential data assimilation, an ensemble of models is run to approximate the covariance matrix [1]. As the desired level of accuracy increases, the resolution of the spatiotemporal numerical discretization also increases. This can cause a significant increase in the computational resources and can become the bottleneck in the outer loop of the design and forecast cycle. For example in computational fluid dynamics (CFD), the three-dimensional simulation of turbulent flows [2] is often used to obtain an exact solution to the Navier-Stokes equation but is less frequently used in related tasks such as shape optimization [3]. Therefore, there is a sustained interest in developing a reduced-order model (ROM) that is computationally much cheaper than the full order model (FOM) and also provides the solution with a sufficient level of accuracy [4][5][6].
In recent years, there is a growing interest in developing non-intrusive ROMs for physical systems [7][8][9]. The non-intrusive refers to the use of only data to construct the ROM. The non-intrusive ROMs are particularly attractive for systems where the perfect model is unknown, but there is an abundance of data available for that system. This situation is commonplace for geophysical flows. The geophysical system model can be imperfect due to a variety of reasons such as approximation of subgrid-scale processes due to insufficient grid resolution, uncertainty in model parameters or incorrect structure of the model itself [10]. However, there is an abundance of data from remote sensing, satellite observations, and experimental measurements for geophysical flows. As a result, there have been several works to harness this data for the effective forecast of geophysical processes [11,12].
In this work, we focus on the use of linear dimensionality reduction technique followed by the machine learning algorithm to evolve the latent space [13][14][15]. Specifically, we use proper orthogonal decomposition (POD) to identify the latent space of the FOM, and then the long short-term memory (LSTM) neural network is utilized to model the evolution of latent space. One of the main advantages of this surrogate modeling technique is that it is purely data-driven and hence is particularly appealing to geophysical flows. The archival data collected from remote sensing and in-situ measurements can be exploited to build the surrogate model and then the surrogate model can be employed for the forecast task or in the data assimilation cycle. As the new data becomes available, the surrogate model can be retrained using the transfer learning approach [16] to improve the forecast performance.
Apart from ROM, the use of the deep neural network has significantly increased for many scientific applications [17][18][19]. One of the main challenges with neural networks is that their performance is highly dependent on their architecture [20]. Additionally, the neural network has a large number of hyperparameters that are problem-dependent. Usually, the neural networks are manually designed through trial and error and this procedure can be time-consuming. Moreover, rich domain knowledge about the data and neural networks is required to achieve a promising performance. There are methods like grid search [21] or a random search [22] to find a good combination of hyperparameters. But these methods are not scalable as the dimension of the search space goes up. Hence, there is a growing interest in automating the design of deep neural network architecture and hyperparameter search that will allow the user with no specialized knowledge of the neural network to apply to their problem of interest [23]. In the present study, we use the genetic algorithm (GA) to optimize the architecture design and hyperparameters of the LSTM neural network. Readers are referred to a recent survey by Liu et al. [24] on most recent evolutionary computation-based neural architecture search methods in light of the core components, to systematically discuss their design principles and justifications on the design. The LSTM network is used within the surrogate model design of the NOAA sear-surface temperature (SST) data set.
The rest of the paper is structured as follows. Firstly, we describe the non-intrusive ROM methodology and data pre-processing in Section II. Then, the details of the proposed algorithm are provided in Section III. Section IV discusses the experimental results and their analysis. Finally, in Section V, conclusions and future works are outlined.

A. Proper orthogonal decomposition
We use the proper orthogonal decomposition (POD) to extract the dominant modes from the data. We collect the data snapshots u 1 , u 2 , . . . , u ∈ R at different time instances. Here, is the spatial degree of freedom which is equal to the total number of grid points, and is the number of snapshots. In POD, we construct a set of orthonormal basis functions that optimally describes the field variable of the system in 2 -norm. The snapshot data matrix A is formed as where the mean-subtracted (or anomaly) fieldsũ ∈ R are computed as follows whereū is the time-average of the solution field. Once the snapshot data matrix is formed, we use singular value decomposition (SVD) to compute the left and right singular vectors of the matrix A. In matrix form, the SVD can be written as where W ∈ R × , ∈ R × , and V ∈ R × . The W and V contains the left and right singular vectors which are identical to the eigenvectors of AA and A A, respectively. Also, the square of singular vales are equal to the eigenvalues, i.e., = 2 . The vectors w k (also the eigenvectors of AA ) are the POD basis functions and we denote them as in this text. The POD basis functions are orthonormal (i.e., , = , where is the Kronecker delta) and are computed in an optimal manner in the 2 sense [25,26]. The state of the dynamical system can be approximated using these POD basis functions as follows where is the number of retained basis functions such that << and are the time-dependent modal coefficients. The POD basis functions minimize the mean-square error between the field variable and its truncated representation. The number of retained modes is usually decided based on their energy content. Using these retained modes, we can form the POD basis set = { } =1 to build the ROM.

B. Model-free evolution of the ROM
The model-free prediction of the evolution of ROMs has a promising application, especially for geophysical dynamical systems. For a variety of tasks such as data assimilation, uncertainty quantification of geophysical flows, the ensemble of the forward model needs to be run and this can be computationally expensive. Instead of evolving the forward model, the surrogate model based on ROMs can be utilized. Also for many geophysical dynamical systems, the mechanistic description of the dynamical system is unavailable or incomplete due to coarse grid resolution. However, there is an abundance of data obtained from local and satellite observations for many past decades. In recent years, data-driven approaches exploiting the recurrent neural network (RNN) are proven successful for the model-free prediction of spatiotemporal chaotic dynamical systems [27,28].
The time-dependent modal-coefficients given in Equation 4 are obtained by projecting the mean subtracted field on the POD basis as follows where angle parenthesis denotes the inner product of two functions. Therefore, if we can model the prediction of the evolution of modal coefficients we can reconstruct the solution field in the future. One of the popular approaches of time-series prediction is the RNN and has been recently applied to several studies on non-intrusive ROMs for physical systems [13,[29][30][31]. This is motivated by the fact that RNN takes into account the past history of the system for future state prediction and the availability of historical data to train the neural network. The long short-term memory (LSTM) neural network is one of the most successful variants of the RNNs that solves the vanishing gradient problem [32]. In this study, we consider an input sequential data matrix X and the output sequential data matrix Y. Each sample of the input training matrix X, i.e., X is constructed as ( ) 1 , . . . , ( ) ; . . . ; ( − +1) 1 , . . . , ( − +1) and the corresponding output sample in output sequential data matrix Y, i.e., Y is ( +1) 1 , . . . , ( +1) , where corresponds to the time index, and is for the th sample. The parameter is called the lookback time-windows. Since we are using the information of only past temporally consecutive states as the input, the LSTM can capture dependencies up to previous time steps.

C. Data preparation
Model-free prediction is particularly important for real-world data which is a manifestation of many multiscale processes that cannot be approximated accurately with any of the model-based approaches. The dynamics of the flow in real-life is governed by the coupling with external systems, and this coupling might be unknown or cannot be modeled. Therefore, the real-world data is not well-behaved and non-intrusive methods are well suited for these complex flows. To this end, we investigate the application of non-intrusive ROMs for the NOAA Optimum Interpolation Sea Surface Temperature (SST) V2 data set * . This data set consists of the weekly average sea surface temperature (SST) data snapshots on a one-degree resolution grid and is produced by combining satellite and local measurements. The seasonal fluctuations give rise to a strong periodic structure to the temperature field in this data set.
Each snapshot of the data set has the dimension of 180 × 360 along the latitudes and longitudes with one-degree grid resolutions. The data points corresponding to land are removed by using the mask operation and the surrogate model is built for the flattened data corresponding to only the ocean surface. This data is available from October 22, 1981, to June 30, 2018 (i.e., total 1,914 snapshots).
Once the data prepossessing is done, then the snapshot data matrix A is formed using the mean-subtracted anomalies of the temperature field. The number of modes for the ROM, R, is selected based on a relative information content (RIC) formula as given below The RIC represents the fraction of information (variance) of the total data that can be recovered using basis functions. Figure 1 shows the RIC percentage for the NOAA SST data set. We fix the number of retained modes to be = 8 that captures around 92% (i.e., 9.182 × 10 1 ) of the variance of the data and these modes are sufficient to capture the seasonal trends in the SST data set. It can also be seen that the increase in the number of modes after 8 modes gives a very small increase in the total information variance, as these modes are mainly responsible for capturing small-scale fluctuations. * https://psl.noaa.gov/ (Accessed 2020)

III. Genetic Algorithm for LSTM Network
The performance of the neural network is highly dependent on the design of the neural network architecture and the selection of other hyperparameters such as activation function, optimizer, weight initialization, etc. Similar to other hyperparameter optimization studies [33,34], we encode the entire LSTM architecture using different types of small architectures. Figure 2 shows different types of small architectures that act as the building blocks of the main LSTM network. In each of these small architectures, skip-connection is employed as it allows training a deep neural network without the vanishing gradient problem [35,36]. The final LSTM network is constructed by ordering these building blocks. The number of LSTM cells and the number of building blocks are the two additional parameters related to the design of the neural network architectures that need to be optimized using the genetic algorithm. The other hyperparameters to optimize are the type of the optimizer (for example, SGD and Adam), the learning rate of the optimizer, type of initialization for weights and biases such as RandomUniform and Xavier initialization, and the type of activation function (for example, ReLU and tanh). Algorithm 1 lists the different steps of the proposed framework. We start by initializing the population in a given size. The LSTM network is randomly constructed for each individual from the predefined building blocks and other hyperparameters are assigned based on the assignment. The neural network is evaluated using the three-fold cross-validation on the training dataset. Once the fitness (i.e., the validation mean squared error in this study) is evaluated for each individual in the population, they are sorted based on their fitness values. From this population, a certain percentage of elite individuals (i.e., the individuals with better fitness values) are retained and some individuals are selected using the tournament selection for offspring generation. The offspring are generated from the selected parents using crossover and mutation operators. The new population is formed by combining the retained elite individuals from the old population and the newly generated offspring. The counter for the number of generations is increased by one and this procedure is continued for the specified number of generations.

Algorithm 1 The proposed genetic algorithm
Input: A set of predefined building blocks, population size, maximum number of generations, dataset for training, elite percentage to retain. Output: Discovered LSTM architecture along with other hyperparameters 1: 0 ← Initialize the population with the given population size.

2:
← 0 3: for = 0 to ℎ do 4: Evaluate the fitness of each individual within the population ;

5:
← Retain the elite population with the maximum value of the fitness; 6: ← Select the population for offspring generation using tournament selection; 7: ← Generate offspring from the selected population using the proposed crossover and mutation operations 8:

9:
= + 1 10: end for 11: return the individual with the best fitness form Each individual in the population is defined using 6 integer numbers and 1 real number. These numbers represent the hyperparameter of the LSTM network. The upper and lower limit of these hyperparameters and the type of distribution used for generating these numbers are reported in Table 1. The first parameter is the type of the building block of the network and four different building blocks as shown in Figure 2 that can be used to design the LSTM network. The second and third parameters are related to how deep the LSTM network will be and how wide the neural network within the LSTM cell will be. The other parameters such as activation function, optimizer, learning rate, and weight initialization are related to the training of the LSTM network. During the evaluation, the LSTM network corresponding to each individual is trained and the validation mean squared error (MSE) is assigned as the fitness for that individual. The population is sorted based on the value of the fitness in ascending order (lower MSE means better fitness). The next step is the selection of parents for the offspring generation. We use tournament selection to select parents. Algorithm 2 shows the procedure of tournament selection. In tournament selection, a few individuals are selected randomly from the population and the best individual is selected for offspring generation. The details of generating the offspring are provided in Algorithm 3. The offspring population is generated in two stages. The first stage is the crossover operation (lines 2 -13), and we utilize the uniform crossover operator. In the uniform crossover operation, two individuals are selected from the population. Then for every gene of the individual, a Learning rate of the optimizer 10 −4 10 −2 Lognormal Real number between 10 −4 to 10 −2 random number is generated, and if this number is greater than the crossover probability then the corresponding gene at the position of two individuals is swapped. During the mutation, the gene at the position is changed if the random number is greater than the mutation probability.

IV. Numerical Experiments
In this Section, we provide the result of the GA in optimizing the architecture design and hyperparameters of the neural network. Then, we show the performance of the optimized LSTM network in emulating and forecasting the SST. The NOAA SST data set is available from October 22, 1981, to June 30, 2018, which corresponds to 1,914 snapshots. We utilize 70% of the data selected randomly from the first 1,500 snapshots as the training data set. During the hyperparameter optimization of the LSTM network, the individual within the population is trained for 100 epochs and three-fold cross-validation is used to avoid overfitting. The lookback time window of the LSTM network is set at = 8, which is found to capture the temporal correlation in modal coefficients. The fitness of each individual is the mean of the mean squared error (MSE) for the held-out data set from three folds. The batch size is fixed at 64 for the training. The population size for this experiment over the number of generations is set at 20. Figure 3 displays the evolution trajectory of the fitness of individuals within the population with generations of the GA. The boxplot is utilized to show the statistics of the population at each generation. Figure 3 also reports the median and minimum MSE of the population at each generation. As the evolution progresses, the MSE on the validation data set decreases. The height of the box indicates the variance of the fitness of the population at each generation, and we can observe that the variance is substantially low after the first generation. The MSE decreases sharply from the first generation to the second generation, which can be attributed to the random initialization of the population at the beginning of the GA. For the problem investigated in this study, the five to ten generations seem enough to find the best set of hyperparameters.
The optimal neural network architecture discovered by the genetic algorithm is shown in Figure 4. Once the best architecture and other hyperparameters are discovered using the GA, the best LSTM network is trained for 1200 epochs with a batch size of 64. During the deployment of the trained network, the auto-regressive method is used to forecast the modal coefficients. The initial condition for modal coefficients for the first 8 time steps (equal to the lookback of the LSTM) is provided. This information is used to predict the modal coefficients at the 9th time step. Then the modal coefficients from the 2nd time step to the 9th time step are used to predict the modal coefficients at the 10th time step. This procedure is repeated up to the final time step, i.e., the 1,914th-time step. Since we use the auto-regressive deployment of the trained LSTM network, only the initial condition corresponding to the LSTM lookback is required. After a few time steps, only the prediction of the LSTM is used to forecast the modal coefficients in the future. In Figure 5, the true modal coefficients and the predicted modal coefficients are shown. We can see that there is a very good agreement between the true and predicted modal coefficients, especially for the first few modes which are responsible for capturing large-scale fluctuations and seasonal patterns. Figure 6 depicts the true and reconstructed temperature field at two different times. The temperature field is ← + 1 13: end while 14: for each individual in do 15: ← Randomly generate a number between [0, 1]; 16: if < then 17: ← Randomly choose a point in ; 18: Change the gene at the point of ; 19: end if 20: end for 21: return Offspring population  reconstructed using Equation 4, where the average temperature field is computed using the first 1,500 snapshots. The large patterns in the temperature field are accurately captured with the LSTM based ROM which demonstrates the potential of data-driven ROMs for geophysical flows. In Figure 7, the 2 -norm of the difference between the true and predicted temperature field is shown. We can observe that the error grows toward the final time, which could be due to the inaccurate prediction of modal coefficients and the inability of the POD basis functions to capture the spatial pattern with high accuracy in the extrapolated time region (i.e., beyond 1500). One of the remedies to this behavior with the LSTM error piling up is to train and deploy the network in a non-auto-regressive manner [37].

V. Concluding Remarks
We presented a genetic algorithm (GA) for automatic architecture search and hyperparameter optimization of the long short-term memory (LSTM) neural network for the task of surrogate modeling of geophysical flows. We have successfully demonstrated the surrogate model built using the optimized LSTM is able to forecast the sea-surface temperature (SST) field. To facilitate the training of deeper neural networks, we adopt an encoding strategy in which the LSTM network is designed using small building blocks that incorporate skip connections. The optimized LSTM network can predict the seasonal variation in the SST field accurately without any instability issue for a sufficiently long period.
We observe that the discrepancy between the true and predicted modal coefficients is higher for the forecast period compared to the training zone. This issue can be addressed either using the non-auto-regressive deployment or transfer learning approach where the LSTM is retrained as the new data becomes available. In the present study, we assume that the number of LSTM cells in each hidden layer is constant. This constraint can be removed and the LSTM network can be optimized with the variable-length encoding strategy for each individual in the population. The variable-length encoding strategy will let us design the neural network architecture that employs a different number of hidden units in each of the LSTM cells and we plan to investigate this in our future studies. Another challenge with the neural architecture search is the computational resources required to train the neural network, and this challenge can be addressed using asynchronous computation.