Most of the neural networks (NNs) learned via state-of-the-art machine learning techniques are black-box models. For a widespread success of machine learning in science and engineering, it is important to develop new NN architectures to effectively extract high-level mathematical knowledge from complex datasets. Motivated by this understanding, this paper develops a new NN architecture called the Gumbel-Max Equation Learner (GMEQL) network. Different from previously proposed Equation Learner (EQL) networks, GMEQL applies continuous relaxation to the network structure via the Gumbel-Max trick and introduces two types of trainable parameters: structure parameters and regression parameters. This paper also proposes a two-stage training process with new techniques to train structure parameters in both online and offline settings based on an elite repository. On 8 benchmark symbolic regression problems, GMEQL is experimentally shown to outperform several cutting-edge machine learning approaches.