Network analysis has played a key role in knowledge discovery and data mining. In many real-world applications in recent years, we are interested in mining multilayer networks, where we have a number of edge sets called layers, which encode different types of connections and/or time-dependent connections over the same set of vertices. Among many network analysis techniques, dense subgraph discovery, aiming to find a dense component in a network, is an essential primitive with a variety of applications in diverse domains. In this paper, we introduce a novel optimization model for dense subgraph discovery in multilayer networks. Our model aims to find a stochastic solution, i.e., a probability distribution over the family of vertex subsets, rather than a single vertex subset, whereas it can also be used for obtaining a single vertex subset. For our model, we design an LP-based polynomial-time exact algorithm. Moreover, to handle large-scale networks, we also devise a simple, scalable preprocessing algorithm, which often reduces the size of the input networks significantly and results in a substantial speed-up. Computational experiments demonstrate the validity of our model and the effectiveness of our algorithms.