The possibility to simulate the properties of many-body open quantum systems with a large number of degrees of freedom is the premise to the solution of several outstanding problems in quantum science and quantum information. The challenge posed by this task lies in the complexity of the density matrix increasing exponentially with the system size. Here, we develop a variational method to efficiently simulate the non-equilibrium steady state of Markovian open quantum systems based on variational Monte Carlo and on a neural network representation of the density matrix. Thanks to the stochastic reconfiguration scheme, the application of the variational principle is translated into the actual integration of the quantum master equation. We test the effectiveness of the method by modeling the two-dimensional dissipative XYZ spin model on a lattice.Open quantum systems have evolved into a major field of studies in recent years. Focus of these studies are the characterization of emergent phenomena and dissipative phase transitions [1-23], as well as the ongoing debate about whether quantum computing schemes are still hard to simulate classically -and thus achieve quantum supremacy -when in presence of some degree of noise-induced decoherence [24][25][26][27].Assuming a Markovian interaction with the environment, the dynamics of open quantum systems is governed by the quantum master equation in Lindblad form [28]. Only few models within this description admit an analytical solution [29,30]. The quest for efficient numerical methods to simulate the dynamics and the asymptotic steady state resulting from the Lindblad master equation, is a research field that is still in its infancy. Many recent tools have been developed following in the footsteps of well established numerical methods for the simulations of closed, Hamiltonian quantum systems. In particular, matrix-product state and tensor network schemes [31][32][33][34], a real-space renormalization approach [35], cluster mean-field [19], and other ad-hoc approximation schemes [23,36] have recently emerged.A groundbreaking progress in the numerical simulation of both the ground state and the dynamics of closed quantum systems has recently been made with the introduction of the neural-network variational ansatz [37][38][39][40][41][42][43], which efficiently represents highly correlated quantum states and whose parameters are easily optimized by means of the variational Monte Carlo (VMC) method. Recently, a self-adjoint and positive semidefinite parametrization of the density matrix, in terms of a neural network has been introduced [44].The steady state of an open quantum system can be characterized by a variational principle [31,33,45,46], whereby the dissipative part of the real-time dynamics, under quite general conditions, drives the system towards a unique steady state, in analogy to the imaginary-time Schrödinger equation that leads to the ground state of Hamiltonian systems.In this Letter, we present a VMC approach to simulate the non-equilibrium steady state (NESS)...