Spectropolarimetric inversion techniques, known as Zeeman Doppler imaging (ZDI), have become the standard tools for reconstructing surface magnetic field maps of stars. Accurate and efficient uncertainty quantification of such magnetic field maps is an open problem in current research, and the high dimensionality of the spherical-harmonic magnetic field parameterization makes inference inherently difficult. We propose a probabilistic machine learning framework for stellar surface magnetic field reconstruction using a gradient-based Metropolis-adjusted Langevin algorithm. By efficient implementation in JAX, our framework allows for reliable uncertainty quantification of the global stellar magnetic field topology. We test the proposed scheme on the bright, massive star Tau Scorpii, and show that our approach enables accurate computation of the posterior magnetic field distribution with fast convergence.