We consider Gaussian distributions on certain Riemannian symmetric spaces. In contrast to the Euclidean case, it is challenging to compute the normalization factors of such distributions, which we refer to as partition functions. In some cases, such as the space of Hermitian positive definite matrices or hyperbolic space, it is possible to compute them exactly using techniques from random matrix theory. However, in most cases which are important to applications, such as the space of symmetric positive definite (SPD) matrices or the Siegel domain, this is only possible numerically. Moreover, when we consider, for instance, high-dimensional SPD matrices, the known algorithms for computing partition functions can become exceedingly slow. Motivated by notions from theoretical physics, we will discuss how to approximate the partition functions in the large $N$ limit: an approximation that gets increasingly better as the dimension of the underlying symmetric space (more precisely, its rank) gets larger. We will give formulas for leading order terms in the case of SPD matrices and related spaces. Furthermore, we will characterize the large $N$ limit of the Siegel domain through a singular integral equation arising as a saddle-point equation.