  Lattice gauge theories coupled to fermionic matter account for many
interesting phenomena in both high energy physics and condensed matter physics.
Certain regimes, e.g. at finite fermion density, are difficult to simulate with
traditional Monte Carlo algorithms due to the so-called sign-problem. We
present a variational, sign-problem-free Monte Carlo method for lattice gauge
theories with continuous gauge groups and apply it to (2+1)-dimensional compact
QED with dynamical fermions at finite density. The variational ansatz is
formulated in the full gauge field basis, i.e. without having to resort to
truncation schemes for the $U(1)$ gauge field Hilbert space. The ansatz
consists of two parts: first, a pure gauge part based on Jastrow-type ansatz
states (which can be connected to certain neural-network ansatz states) and
secondly, on a fermionic part based on gauge-field dependent fermionic Gaussian
states. These are designed in such a way that the gauge field integral over all
fermionic Gaussian states is gauge-invariant and at the same time still
efficiently tractable. To ensure the validity of the method we benchmark the
pure gauge part of the ansatz against another variational method and the full
ansatz against an existing Monte Carlo simulation where the sign-problem is
absent. Moreover, in limiting cases where the exact ground state is known we
show that our ansatz is able to capture this behavior. Finally, we study a
sign-problem affected regime by probing density-induced phase transitions.
