We introduce a new numerical technique, the bosonic auxiliary-field Monte Carlo method, which allows us 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 nonentangled, 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 us to reproduce numerically exact results for the superfluid transition of hard-core 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.