Random noncommutative geometry can be seen as a Euclidean path-integral approach to the quantization of the theory defined by the Spectral Action in noncommutative geometry (NCG). With the aim of investigating phase transitions in random NCG of arbitrary dimension, we study the non-perturbative Functional Renormalization Group for multimatrix models whose action consists of noncommutative polynomials in Hermitian and anti-Hermitian matrices. Such structure is dictated by the Spectral Action for the Dirac operator in Barretts spectral triple formulation of fuzzy spaces.The present mathematically rigorous treatment puts forward coordinate-free language that might be useful also elsewhere, all the more so because our approach holds for general multimatrix models. The toolkit is a noncommutative calculus on the free algebra that allows to describe the generator of the renormalization group flow -- a noncommutative Laplacian introduced here -- in terms of Voiculescus cyclic gradient and Rota-Sagan-Stein noncommutative derivative. We explore the algebraic structure of the Functional Renormalization Group Equation and, as an application of this formalism, we find the $beta$-functions, identify the fixed points in the large-$N$ limit and obtain the critical exponents of $2$-dimensional geometries in two different signatures.