In this study, an efficient stochastic gradient-free method, the ensemble neural networks (ENN), is developed. In the ENN, the optimization process relies on covariance matrices rather than derivatives. The covariance matrices are calculated by the ensemble randomized maximum likelihood algorithm (EnRML), which is an inverse modeling method. The ENN is able to simultaneously provide estimations and perform uncertainty quantification since it is built under the Bayesian framework. The ENN is also robust to small training data size because the ensemble of stochastic realizations essentially enlarges the training dataset. This constitutes a desirable characteristic, especially for real-world engineering applications. In addition, the ENN does not require the calculation of gradients, which enables the use of complicated neuron models and loss functions in neural networks. We experimentally demonstrate benefits of the proposed model, in particular showing that the ENN performs much better than the traditional Bayesian neural networks (BNN). The EnRML in ENN is a substitution of gradient-based optimization algorithms, which means that it can be directly combined with the feed-forward process in other existing (deep) neural networks, such as convolutional neural networks (CNN) and recurrent neural networks (RNN), broadening future applications of the ENN.3 the quality of text that has been machine-translated (Papineni et al., 2002;Reiter, 2018). However, it is difficult to build a loss function based on this evaluation criterion since it is not differentiable. However, this will no longer pose a problem if we can find a gradient-free optimization method.Considering the aforementioned problems, a salient question is: are there any alternatives for the optimization method in a neural network that are able to perform uncertainty analysis and perform well with a small dataset, but do not rely on derivative calculations?These obstacles are encountered in numerous engineering fields, such as petroleum engineering. Uncertainty quantification is critical because underground geological parameters are highly heterogeneous. High-dimension models are always solved based on a small dataset due to the expensive and time-consuming data collection. It is also difficult to identify gradients of a target variable with respect to model parameters because the corresponding physical models are highly nonlinear and too complicated to solve analytically. In response to these problems, the ensemble randomized maximum likelihood algorithm (EnRML) is proposed by Gu and Oliver (2007) in the field of history matching in petroleum engineering. History matching is an inverse modeling method, which adjusts a model of a reservoir until it closely reproduces its past behavior (Oliver et al., 2008;Stordal & Naevdal, 2018). It should be mentioned that the word "ensemble" here indicates a different meaning from that in ensemble averaging (Naftaly et al., 1997). In the former, it means the ensemble of realizations generated from the same model, rather tha...