  We introduce the first bold diagrammatic Monte Carlo approach to deal with
polaron problems at finite density non-perturbatively, i.e., by including
vertex corrections to high orders. Using Holstein model on a square lattice as
a prototypical example, we demonstrate that our method is capable of providing
accurate results in the thermodynamic limit in all regimes from renormalized
Fermi-liquid to single polarons, across the non-adiabatic region where Fermi
and Debye energies are of the same order of magnitude. By accounting for vertex
corrections the accuracy of theoretical description is increased by orders of
magnitude relative to the lowest-order self-consistent Born approximation
employed in most studies. We also find that for electron-phonon coupling
typical for real materials, the quasiparticle effective mass increases and the
quasiparticle residue decreases with increasing the system density.
