<p>We introduce a highly efficient fully Bayesian approach for anisotropic multidimensional smoothing. The main challenge in this context is the Markov chain Monte Carlo (MCMC) update of the smoothing parameters as their full conditional posterior comprises a pseudo-determinant that appears to be intractable at first sight. As a consequence, existing implementations are computationally feasible only for the estimation of two-dimensional tensor product smooths, which is, however, too restrictive for many applications. In this paper, we break this barrier and derive closed-form expressions for the log-pseudo-determinant and its first and second order partial derivatives. These expressions are valid for arbitrary dimension and very fast to evaluate, which allows us to set up an efficient MCMC sampler with derivative-based Metropolis–Hastings (MH) updates for the smoothing parameters. We derive simple formulas for low-dimensional slices and averages to facilitate visualization and investigate hyperprior sensitivity. We show that our new approach outperforms previous suggestions in the literature in terms of accuracy, scalability and computational cost and demonstrate its applicability through an illustrating temperature data example from spatio-temporal statistics.</p>