Identifying key differences between linear stochastic estimation and neural networks for fluid flow regressions

Neural networks (NNs) and linear stochastic estimation (LSE) have widely been utilized as powerful tools for fluid-flow regressions. We investigate fundamental differences between them considering two canonical fluid-flow problems: (1) the estimation of high-order proper orthogonal decomposition coefficients from low-order their counterparts for a flow around a two-dimensional cylinder, and (2) the state estimation from wall characteristics in a turbulent channel flow. In the first problem, we compare the performance of LSE to that of a multi-layer perceptron (MLP). With the channel flow example, we capitalize on a convolutional neural network (CNN) as a nonlinear model which can handle high-dimensional fluid flows. For both cases, the nonlinear NNs outperform the linear methods thanks to nonlinear activation functions. We also perform error-curve analyses regarding the estimation error and the response of weights inside models. Our analysis visualizes the robustness against noisy perturbation on the error-curve domain while revealing the fundamental difference of the covered tools for fluid-flow regressions.

Recent boom of neural networks (NNs) has expanded into fluid mechanics 1,2 . Although these have still been staying at the fundamental level, NNs have indeed shown their great potentials for data estimation [3][4][5] , control 6 , reduced-order modeling 7 , and turbulence modeling 8 . Despite the recent enthusiastic trends on the use of NNs, we still have to rely largely on linear theories because their generalizability and transparency are, for the moment, superior to NN-based methods. This indicates that we may be able to obtain some clues from the relationships between NNs and linear theories to develop more advanced NN-based techniques for fluid flow analyses.
Motivated above, of particular interest here is an analogy between NNs and linear methods. Milano & Koumoutsakos 9 used an autoencoder (AE) to perform low-dimensionalizaton for the Burgers equation and a turbulent channel flow. They also reported that the operation inside the multi-layer perceptron (MLP) based AE with linear activation functions is equivalent to that of proper orthogonal decomposition (POD) 10 . More recently, Murata et al. 11 investigated the analogy between a convolutional neural network (CNN)-based AE and POD. The strength of the nonlinear activation function used inside NNs was also demonstrated. From the perspective on state estimation, Nair & Goza 12 have recently compared the MLP, Gappy POD, and linear stochastic estimation (LSE) for the POD coefficient estimation from local sensor measurements of laminar wake of a flat plate.
To clarify the fundamental difference between linear methods and NNs, we here compare abilities of LSE, MLP, and CNN by considering two canonical fluid flow regression problems whose complexities are different with each other; 1. estimation of high-order POD coefficients from low-order counterparts of a flow around a cylinder 13 , and 2. state estimation from information on the wall in a turbulent channel flow 14 , as illustrated in Fig. 1. In particular, we seek key attributes for their estimations by focusing on influences of biases inside NNs, optimization methods, and robustness against noisy inputs.

Regression methods
Multi-layer perceptron. Let us first introduce a multi-layer perceptron (MLP) 15 . The MLP, which mimics the neurons in human's brain, has widely been utilized in physical sciences 16 . In the present paper, the MLP is used for low-dimensional fluid flow regressions, i.e., POD coefficient estimation, as an example of NN. The MLP is an aggregate of the minimum units called perceptrons. A linear superposition of input data q with weights W www.nature.com/scientificreports/ are passed inside a perceptron while biases b being added, and then a nonlinear activation function φ is applied such that where l denotes a layer index. Weights among all edges W ij are optimized with back-propagation 18 so as to minimize a loss function E. The present MLP has hidden units of 4-8-16-8, while the number output nodes is 4 (3rd to 6th POD modes). The number of input nodes varies depending on considered cases, whose details will be provided in Sect. 3.1. In the simplest case, the present MLP M attempts to output a out = {a 3 , a 4 , a 5 , a 6 } from 2 inputs a in = f (a 1 , a 2 ) , and the problem setting regarding weights w m (representing both W and b ) inside the MLP can be represented as We use the L 2 error norm as a loss function. Note that penalization terms, e.g., Lasso and Ridge penalties, are not considered in the present loss function because of the difficulty and cost in tuning the hyperparameter for regularization 17 . Convolutional neural network. One of the issues associated with the MLP is that the number of edges inside it may explode due to its fully-connected structure when handling high-dimensional data such as fluid flows. To overcome this issue in dealing with fluid flow problems, a convolutional neural network (CNN) 19 has widely been accepted as a good candidate 3,20 . We capitalize on the combination of two-and three-dimensional CNNs for the state estimation task in the present study. The convolutional layer, which is a fundamental operation inside CNNs, extracts spatial features from input data using a filter operation, where C = floor(H/2) , K is the number of filters in a convolution layer, and b (l) m is the bias. This filter operation achieves efficient data handling for two-or three-dimensional flow fields 2 . For the three-dimensional CNN, the three-dimensional convolution is performed similarly to Eq. (3). In the present paper, the size H is set to 5. In the CNN, the output from the filtering operation is passed through an activation function φ , analogous to the operation inside an MLP (Eq. (1)). Filters h are optimized with back-propagation to minimize a loss function.
We use the combination of two-and three-dimensional CNNs to output a three-dimensional turbulent state u from streamwise wall-shear stress τ x,wall (details will be provided in Sect. 3.2). The optimization problem of weights w c (representing both h and b ) for the CNN C can be expressed as Again, the aim of this study is the comparison between LSE and NN. To achieve a fair comparison, it is ideal to consider the same amount of weights for both LSE and NNs. To this end, we perform singular value decomposition (SVD) for the LSE weights w l to align the number of weights contained inside the CNNs and the LSE.  www.nature.com/scientificreports/ The details for the weight design for the LSE will be provided in the next section. The operation can be expressed as w l = U ŴV T , where U ∈ R n input ×n rank and V T ∈ R n rank ×n output are singular vectors, and Ŵ(∈ R n rank ×n rank ) is a diagonal matrix whose diagonal elements are singular values. Since the rank of LSE weights in the present study is n rank = rank(w l ) = 193 according to our preliminary test, the present number of LSE weights based on SVD reduction n w LSE,SVD is 197,632 (= n input × n rank = 1024 × 193) . Capitalizing on this SVD-based weight reduction, we are now able to obtain the target number of weights of approximately 197,000 to determine the CNN parameters. In this study, our CNN contains 196 608 weights with the parameters in Table 1.

Linear stochastic estimation.
For comparison to the NNs, we use linear stochastic estimation (LSE) 14,21 .
In this study, we express target data Q ∈ R n data ×n output (output) as a linear map w l ∈ R n input ×n output with respect to input data P ∈ R n data ×n input such that Q = Pw l , where n data represents the number of training snapshots, n input represents the number of input attribute, and n output represents the number of output attribute. Analogous to the optimization for NNs, the linear map w l can be obtained through a minimization such that Note that penalization terms are not considered in the present loss function for the fair comparison to the NNs in terms of weight updating. We also emphasize that the LSE is analytically optimized solving Eq. (5) while the NNs are numerically optimized through back-propagation. Hence, we can also compare the NNs and the LSE regarding the difference of the optimization approaches. The optimized weights w l can then be applied to test data.

Results
Example 1: POD coefficient of two-dimensional cylinder wake at Re D = 100. As presented in Fig. 1a, we first aim to estimate high-order POD coefficients a out = {a 3 , a 4 , a 5 , a 6 } of a two-dimensional cylinder wake at Re D = 100 from information of low-order counterparts a in = f (a 1 , a 2 ) such that a out = F 1 (a in ) , where F 1 denotes a model for this purpose. The LSE and the MLP are used as the model F 1 . Flow snapshots are generated using a two-dimensional direct numerical simulation (DNS). The governing equations are the incompressible Navier-Stokes equations, where u and p denote the velocity vector and the pressure, respectively. All quantities are non-dimensionalized using the fluid density, the free-stream velocity, and the cylinder diameter. The size of the computational domain is (L x , L y ) = (25.6, 20.0) , and the cylinder center is located at (x, y) = (9, 0) . The grid spacing and the time step are respectively x = y = 0.025 and t = 2.5 × 10 −3 , while imposing the no-slip boundary condition on the cylinder surface using an immersed boundary method 22 . The number of grid points used for the DNS is (N x , N y ) = (1024, 800) . For the POD, the vorticity field ω around the cylinder is extracted as a domain of 192). We then take the POD for the collected snapshots to decompose the flow field q as q = q 0 + M i=1 a i ϕ i , where ϕ denotes a POD basis, a is the POD coefficient, q 0 is the temporal average of the flow field, and M represents the number of POD modes. For training the present MLP and LSE, we use 5000 snapshots. For comparison with LSE, we do not divide the training data for MLP into training and validation. We also consider additional 5000 snapshots for the assessment. We here compare the LSE with the linear MLP and the nonlinear MLP with ReLU activation function 23 . The ReLU is known as a good candidate to prevent a vanishing gradient issue. We consider three patterns of input a in = f (a 1 , a 2 ) for the LSE and the linear MLP; the input f (a 1 , www.nature.com/scientificreports/ nonlinear MLP. Since Loiseau et al. 13 reported that the high-order coefficients a out can be represented using the quadratic expression of a 1 and a 2 due to its triadic interaction, this analysis enables us to check two viewpoints; 1. whether the nonlinear function inside MLP works to capture such nonlinear interactions, and 2. whether the linear models can also be utilized if a proper input including the essential combinations of nonlinear terms is given. Let us demonstrate the estimation of a out = {a 3 , a 4 , a 5 , a 6 } from only the information of first-order coefficients a 1st = {a 1 , a 2 } such that a out = F 1 (a 1st ) , as shown in Fig. 2a. The nonlinear MLP shows its clear advantage against the LSE and the linear MLP for both coefficient maps. The L 2 error norm ε = ||a out,ref − a out,est || 2 /||a out,ref || 2 for each case are 1.00 (LSE), 1.00 (linear MLP), and 0.0119 (nonlinear MLP), respectively. This suggests that the nonlinear activation function plays an important role in estimation. Noteworthy here, however, is that this nonlinearity can be recovered by giving a proper input, i.e., a in = {a 2nd , a 3rd } , even if we only use the linear methods, as presented in Figs. 2b,c. The reasonable estimation for {a 3 , a 4 } can be achieved utilizing the input up to the 2nd order term a 2nd , while that for {a 5 , a 6 } requires the 3rd order term a 3rd with both the LSE and the linear MLP. This trend is analogous to the observation by Loiseau et al. 13 , as introduced above. The LSE outperforms the linear MLP with the high-order coefficient inputs in this example, as shown in Fig. 2c; however, they will show a significant difference in terms of noise robustness as discussed later.
We then compare the LSE and the MLP in terms of the availability of training data. The dependence of the L 2 error on the number of training snapshots is examined in Fig. 3. Based on the results in Fig. 2, we choose the third-order coefficients a in = a 3rd for the linear models and the first-order coefficients for the nonlinear MLP, as the input for the models. The LSE shows its advantage over the MLP models when the number of training snapshots is limited. This is because the degree of freedom in the MLP is larger than that in the LSE. Note that the MLP will show a clear advantage against the LSE in terms of noise robustness, which reveals the fundamental difference between the linear MLP and the LSE.
We here consider the Gaussian white noise defined by the signal-to-noise ratio (SNR), SNR = σ 2 data /σ 2 noise , where σ 2 data and σ 2 noise are the variances of input data and noise, respectively. The behaviors for noisy inputs are summarized in Fig. 4. As the linear models, we use the LSE and the linear MLP with a in = a 3rd . For comparison,  www.nature.com/scientificreports/ we also monitor the nonlinear MLP with a in = a 1st . The response of LSE is much more sensitive than that of the covered MLPs. This is caused by two considerable reasons: one is the influence of biases contained in the MLPs, as expressed in Eq. (1), while the other being the difference of optimization methods. Hereafter, we will seek which has the main contribution for the noise robustness. For the investigation of the main contribution as mentioned above, we consider the LSE and three types of MLPs as follows: 1. LSE model: the same LSE model as that used above. 2. Linear MLP model with bias M 1 : the same linear MLP as that used above 3. Linear MLP model without bias M 2 : the biases are removed from the model M 1 to investigate the influence on the bias. 4. Shallow linear MLP model without bias M 3 : the MLP with a single hidden layer is prepared to align the number of weights with the LSE so that the difference of optimization methods can be assessed.
The dependence of the L 2 error norm on the noise magnitude of each model is summarized in Fig. 5a. There is almost no difference among the covered MLP models. This suggests that the bias and the number of layers do not contribute to the noise robustness so much. The other observation here is that the shallow linear MLP is still more robust than the LSE, even though the model structures are identical with each other. To examine this point, we visualize the weights inside the LSE and the shallow linear MLP in Fig. 5b. The weights for the secondorder term input ( a 1 a 2 , a 2 1 , a 2 2 ) are optimized to the same values with each other, while that for the first-order input term input exhibits a significant difference. This is caused by the difference of the optimization methods. Moreover, this point is examined by visualizing an error surface for the input of a 1 and a 3 1 with the output of a 5 , as shown in Fig. 5d. The reason for the choice of this input-output combination is that the output a 5 is one of the most sensitive components to the noise for the LSE, as presented in Fig. 5c. The optimized solutions are different to each other, which is likely caused by the difference in optimization methods. What is notable here is that the noise addition drastically changes the error-surface shape of the LSE, while that of the MLP changes only slightly. This difference can be quantified using the mean structural similarity index measure (MSSIM) 24 . We apply this to the elevation of each error surface, i.e., without E r ∈ R M×N and with noise E n ∈ R M×N , where M, N are the number of samples on the error surface for each weight axis. MSSIM is mathematically expressed as This measurement can assess a similarity between two images E r ∈ R M×N and E n ∈ R M×N by considering their mean µ and standard deviation σ . To obtain MSSIM, the SSIM in a small window of two images e r ∈ R m×n and e n ∈ R m×n , where M ′ = M − m + 1 and N ′ = N − n + 1 , is computed and its average is taken over the image. As the constant values C 1 and C 2 in Eq. (7), we set {C 1 , C 2 } = {0.16, 1.44} following Wang et al. 24 . As presented in Fig. 5d, the MSSIM of the shallow linear MLP reports 0.879 while that of LSE is 0.428, which indicates that the deformation of the error surface is substantially larger in the LSE. Due to this large deformation of the error surface, the optimum point of the LSE is pushed up vertically in the error space of Fig. 5d. This indicates that the weights obtained by the LSE in an analytical manner guarantee the global optimal solution over the training data; however, this solution may not be optimal from the viewpoint of noise robustness. On the other hand, the MLP provides a noise-robust solution, although it is not the exact global optimum over the training data since the MLP weights are a numerical solution obtained through back-propagation.
Considering the characteristics of LSE which fits training data, we then add noisy data of SNR = 0.05 to training data for both LSE and the shallow linear MLP (without bias), as summarized in Fig. 6. As shown in Fig. 6a, .  Example 2: velocity field in a minimal turbulent channel flow at Re τ = 110. We then perform the comparison of NNs and LSE for a more complex problem. Let us consider the estimation of a velocity field u in a minimal turbulent channel flow at Re τ = 110 from the streamwise wall shear stress input τ x , as illustrated in Fig. 1b, such that u = F 2 (τ x ) , where F 2 denotes a model for example 2. As already mentioned, we use LSE and CNN as the model F 2 . The training data is generated using a three-dimensional DNS which numerically solves the incompressible Navier-Stokes equations, where u and p represents the velocity vector and pressure, respectively 25,26 . The quantities are nondimensionalized with the channel half-width δ and the friction velocity u τ . The computational domain is (L x , L y , L z ) = (πδ, 2δ, 0.5πδ) with the number of grid points of (N x , N y , N z ) = (32, 64, 32) . The grid is arranged uniformly in the x and the z directions, while nonuniformly in the y direction. The time step is t + = 0.0385 , where the subscript + denotes the wall units. We use 10,000 snapshots to train the models. For equivalent comparison with LSE, we do not divide the training data into training and validation. We also prepare additional 2700 snapshots for the assessment. The channel flow fields estimated by the models are assessed in Fig. 7a. Note again that the number of weights inside the CNN and the LSE is almost the same with each other as explained in Sect. 2.2. We here also consider the linear CNN (i.e., the same CNN structure, but with the linear activation function) to examine whether the similar observation to the linear MLP of the cylinder example can be found or not even for a CNN whose filter operation is different from a fully-connected MLP. The estimated fields are visualized using the second invariant of the velocity gradient tensor Q + of −0.005 . The field reconstructed by the LSE shows the similar behavior to the reference DNS data qualitatively, e.g., the amount of vortical structure; however, it should be emphasized that the LSE field merely provides turbulent-like structure, which is different from that of the DNS.
To investigate this point, we also visualize the x − z sectional streamwise velocity distributions in Fig. 7b. Notably, the nonlinear CNN outperforms the LSE and the linear CNN, which is not intuitive from the observation with the Q isosurfaces. Especially, the advantage of the nonlinear method can be clearly found at y + = 30.1 . We also present the time-averaged Reynolds shear stress −u ′ v ′ in Fig. 7c. The LSE curve looks to be in reasonable agreement in its shape (although overestimated), despite its high-error level as stated above. With the observation for the reasonable reconstruction of LSE in Fig. 7a, it implies that the flow field estimated by the LSE model is similar to the DNS data in a time-ensemble sense, although it does not match for each instantaneous field.
We also investigate the dependence of estimation on the y position in Fig. 7d. Analogous to previous studies 14,27 , it is hard to estimate the velocity field in the region away from the wall because of the lack of (8) www.nature.com/scientificreports/ correlation between the wall shear stress τ x and the velocity field u away from the wall. Within the range that can be estimated, the nonlinear CNN presents the better estimation than the linear methods in terms of both the L 2 error and the statistics. We further assess the difference among the linear and nonlinear methods focusing on noise robustness, as summarized in Fig. 8. Analogous to the noise investigation for the POD coefficient estimation, the noisy input for the the wall shear stress is designed with the SNR. The L 2 error of the LSE explodes rapidly with the noise addition, which is the same behavior as that in Fig. 4a. In contrast, the CNNs are still able to reconstruct the large scale structure despite that the LSE cannot even with 1/SNR = 0.05 . These findings also make us to suspect the role of optimization methods.
Let us then discuss the contribution of optimization methods for noise robustness. We here skip the influence on bias uses since there is no significant effects according to our preliminary test, which is akin to the trend with the cylinder example. For the particular demonstration here, we arrange a shallow linear CNN, in which there are n output filters and each filter has a shape of (n input , 1) . Note that this shallow CNN is distinct against the CNNs used above, which were based on the original SVD-based weight reduction expressed in Sect. 2.2. The use of the shallow linear CNN enables us to observe the weight sensitivity against noisy inputs while comparing to the LSE directly, thanks to its filter shape.
The error surfaces of the LSE and the shallow linear CNN are visualized in Figs. 9a-c. Note that we here choose six points for the visualization. The error used in the error surface is arranged by the streamwise velocity u at an arbitrary point in the x − z cross section at y + = 15.4 . In the case of the shallow linear CNN, the noise has little influence on the error surface. On the other hand, the error surfaces of the LSE drastically change their shape with the existence of noise. These trends can also be found with the MSSIM. It implies that the CNN can also obtain noise robustness thanks to the gradient method, while it is hard to obtain with the LSE, which is analogous to the observation with the POD coefficient estimation.

Conclusions
Fundamental differences between neural networks (NNs) and linear methods were investigated considering canonical fluid flow regression problems: 1. the estimation of POD coefficients of a flow around a cylinder, and 2. the state estimation from wall measurements in a turbulent channel flow. We compared linear stochastic estimation (LSE) with multi-layer perceptron (MLP) and convolutional neural network (CNN). For both regression problems, efficacy of nonlinear function can be observed. We also found that a linear model could surrogate a nonlinear model by giving an appropriate combination of inputs under the consideration of nonlinear relationship. This enables us to expect that the combination of nonlinear activation functions and proper inputs can further enhance the prediction capability of the model, which is similar to the observation in several previous studies [28][29][30] . In addition, the linear NNs were more robust against noise than the LSE, and the reason for this was revealed by visualizing the error surface. The error surface told us that the difference in optimization methods has a significant contribution to the noise robustness.
Although we observed the strength of nonlinear NNs from several perspectives, we should note that the learning process of NNs can also be unstable depending on a problem setting since it is founded on a gradient method. This implies that we may not reach a reasonable valley of a weight surface, especially when the problem has a multimodality and the availability of training data is limited 31 . In this sense, the LSE can provide us a stable solution in a theoretical manner. Hence, we may be able to unify these characteristics for further improvement of NN learning pipelines, e.g., transfer learning 32,33 , so that a learning process of NN can be started from a reasonable solution while achieving a proper noise robustness. Furthermore, it should also be emphasized that we can introduce a loss function associated with a priori knowledge from physics since both NNs and LSE are based on a minimization manner in terms of weights. Inserting a physical loss function may be one of the considerable paths towards practical applications of both methods [34][35][36] . www.nature.com/scientificreports/