Non-Hermitian (NH) quantum many-body systems exhibit a rich array of physical phenomena, including NH skin effects and exceptional points, that remain largely inaccessible to existing numerical techniques. In this Letter, we investigate the application of variational Monte Carlo and neural network wave function representations to examine their ground-state properties. Due to the breakdown of the Rayleigh-Ritz variational principle in NH settings, we develop a self-consistent symmetric optimization framework based on variance minimization with a dynamically updated energy estimate. Our approach respects the biorthogonal structure of left and right eigenstates, and is further strengthened by exploiting system symmetries and pseudo-Hermiticity. Leveraging this tool, we probe and report accurate NH physical observables for a two-dimensional transverse-field Ising model with a complex longitudinal field, spanning both parity-time symmetric and broken phases. Lastly, we show, through extensive numerical evidence, that our method offers a scalable and flexible computational tool to investigate NH quantum many-body systems, beyond the reach of conventional numerical techniques such as the density-matrix renormalization group algorithm.