A recent reformulation [1] of the problem of Coulomb gases in the presence of a dynamical dielectric medium showed that finite temperature simulations of such systems can be accomplished on the basis of completely local Hamiltonians on a spatial lattice by including additional bosonic fields. For large systems, the Monte Carlo algorithm proposed in Ref. [1] becomes inefficient due to a low acceptance rate for particle moves in a fixed background multiboson field. We show here how this problem can be circumvented by use of a coupled particle-multiboson update procedure that improves acceptance rates on large lattices by orders of magnitude. The method is tested on a one-component plasma with neutral dielectric particles for a variety of system sizes.