A fuzzy geometry is a certain type of spectral triple whose Dirac operator crucially turns out to be a finite matrix. This notion was introduced in [J. Barrett, J. Math. Phys. 56, 082301 (2015)] and accommodates familiar fuzzy spaces like spheres and tori. In the framework of random noncommutative 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 noncommutative 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.