Lead halide perovskites have emerged as highly efficient solar cell materials. However, to date, the most promising members of this class are polymorphs in which a wide-band-gap delta phase competes with the photoactive perovskite alpha form and the intrinsic physical interactions that stabilize one phase over the other are currently not well understood. Classical molecular dynamics simulations based on suitably parametrized force fields (FF) enable computational studies over broad temperature (and pressure) ranges and can help to identify the underlying factors that govern relative phase stability at the atomic level. In this article, we present a force-matching-based approach for the automatized generation of polarizable (pol) as well as nonpolarizable (npol) FFs from high-level reference data and apply it to the all-inorganic lead halide material CsPbI3 as a prototype system exhibiting a delta/alpha polymorphism. These force-matched npol and pol FFs have been determined based on extensive reference data from first-principles molecular dynamics simulations over a wide range of temperatures. While both FFs are able to describe the perovskite as well as the nonperovskite delta phase, finer structural details, as well as the relative phase stability, are better reproduced with the polarizable version. A comparison of these ab initio-derived interatomic potentials allows direct insight into the physical origin of the interactions that govern the interplay between the two competing phases. It turns out that explicit polarization is the essential factor that stabilizes the strongly anisotropic delta phase over the high-symmetry (cubic) perovskite alpha phase at lower temperatures. This fundamental difference between alpha and delta phases appears universal and might thus also hold for other perovskite compounds with delta/alpha polymorphism providing rational guidance for synthetic efforts to stabilize the photoactive perovskite phase at room temperature.