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 Euclidean Monte Carlo algorithms due to the sign-problem. I will 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, infinite-dimensional U(1) gauge field Hilbert space, i.e. without truncating the U(1) gauge group to some finite subgroup. 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, I will present several benchmarking cases, in particular against an existing Euclidean Monte Carlo simulation at zero chemical potential where the sign-problem is absent. Moreover, in limiting cases where the exact ground state is known I will show that the ansatz is able to capture this behavior. Finally, a sign-problem affected regime will be studied by probing density-induced phase transitions.