We introduce a new numerical technique -- bosonic auxiliary-field Monte Carlo (bAFMC) -- which allows to calculate the thermal properties of large lattice-boson systems within a systematically improvable semiclassical approach, and which is virtually applicable to any bosonic model. Our method amounts to a decomposition of the lattice into clusters, and to an Ansatz for the density matrix of the system in the form of a cluster-separable state -- with non-entangled, yet classically correlated clusters. This approximation eliminates any sign problem, and can be systematically improved upon by using clusters of growing size. Extrapolation in the cluster size allows to reproduce numerically exact results for the superfluid transition of hardcore bosons on the square lattice, and to provide a solid quantitative prediction for the superfluid and chiral transition of hardcore bosons on the frustrated triangular lattice.