We study lattice QED in $2+1$ dimensions coupled to Wilson fermions in the Hamiltonian formulation and develop a sign-problem-free variational Monte Carlo algorithm. Using a continuous basis to describe the gauge degrees of freedom and fermionic Gaussian states for the matter fields, we construct a variational ansatz for the ground state. The parameters of the ansatz are optimized using gradient descent and/or imaginary time evolution based on the time-dependent variational principle. For a lattice with $N$ matter sites, the computational complexity for sampling a gauge configuration is $\mathcal{O}(N^2\log N)$, and the complexity for computing the energy per gauge configuration is $\mathcal{O}(N^4)$. Comparing our method to results from exact diagonalization on small lattices, we demonstrate good agreement. In addition, we present preliminary results for the ground state for two fermion flavors, $N_f=2$, in the presence of a non-zero chemical potential, a regime in which the conventional action-based Monte Carlo methods suffer from the sign problem.