  We introduce a novel method of efficiently simulating the non-equilibrium
steady state of large many-body open quantum systems with highly non-local
interactions, based on a variational Monte Carlo optimization of a matrix
product operator ansatz. Our approach outperforms and offers several advantages
over comparable algorithms, such as an improved scaling of the computational
cost with respect to the bond dimension for periodic systems. We showcase the
versatility of our approach by studying the phase diagrams and correlation
functions of the dissipative quantum Ising model with collective dephasing and
long-ranged power law interactions for spin chains of up to $N=100$ spins.
