A fuzzy geometry is a certain type of finite-dimensional spectral triple whose Dirac operator crucially turns out to be a (finite) matrix. In the framework of random noncommutatative geometry, we use Barrett's characterization of Dirac operators of fuzzy geometries in order to systematically compute the spectral action $S(D)= \mathrm{Tr} f(D)$ for $2n$-dimensional fuzzy geometries. In contrast to the original Chamseddine-Connes spectral action, we take a polynomial $f$ with $f(x)\to \infty$ as $ |x|\to\infty$ in order to obtain a well-defined path integral that can be stated as a random matrix model with action of the type $S(D)=N \cdot \mathrm{tr}\, F+\textstyle\sum_i \mathrm{tr}\,A_i \cdot \mathrm{tr} \,B_i $, being $F,A_i $ and $B_i $ noncommutative polynomials in $2^{2n-1}$ complex $N\times N$ matrices that parametrize the Dirac operator $D$. For arbitrary signature---thus for any admissible KO-dimension---formulas for 2-dimensional fuzzy geometries are given up to a sextic polynomial, and up to a quartic polynomial for 4-dimensional ones, with focus on the octo-matrix models for Lorentzian and Riemannian signatures. The noncommutatative polynomials $F,A_i $ and $B_i$ are obtained via chord diagrams and satisfy: independence of $N$; self-adjointness of the main polynomial $F$ (modulo cyclic reordering of each monomial); also up to cyclicity, either self-adjointness or anti-self-adjointness of $A_i $ and $B_i $ simultaneously, for fixed $i$. Collectively, this favors a free probabilistic perspective for the large-$N$ limit we elaborate on.