[en] We advocate a method to improve systematically the self-consistent harmonic approximation (or the gaussian approximation), which has been employed extensively in condensed matter physics and statistical mechanics. Such a method was previously applied to the IIB matrix model, a conjectured nonperturbative definition of type-IIB superstring theory in ten dimensions. Remarkably the dominance of four-dimensional space-time in the partition function was suggested from calculations up to the 3rd order. Recently this calculation has been extended to the 5th order, and the same conclusion has been obtained. Here we apply this gaussian expansion method to the bosonic version of the IIB matrix model, where Monte Carlo results are available, and demonstrate the convergence of the method by explicit calculations up to the 7th order. More generally we study matrix models obtained from dimensional reduction of SU(N) Yang-Mills theory in D dimensions, where the D=10 case corresponds to the bosonic IIB matrix model. Convergence becomes faster as D increases, and for D≥10 it is already achieved at the 3rd order. (author)