AbstractFor M a simple surface, the nonlinear statistical inverse problem of recovering a matrix field from discrete, noisy measurements of the SO(n)‐valued scattering data CΦ of a solution of a matrix ODE is considered (n ≥ 2). Injectivity of the map Φ ↦ CΦ was established by Paternain, Salo, and Uhlmann in 2012. A statistical algorithm for the solution of this inverse problem based on Gaussian process priors is proposed, and it is shown how it can be implemented by infinite‐dimensional MCMC methods. It is further shown that as the number N of measurements of point evaluations of CΦ increases, the statistical error in the recovery of Φ converges to 0 in L2(M)‐distance at a rate that is algebraic in 1/N and approaches for smooth matrix fields Φ. The proof relies, among other things, on a new stability estimate for the inverse map CΦ → Φ.Key applications of our results are discussed in the case n = 3 to polarimetric neutron tomography. © 2020 The Authors. Communications on Pure and Applied Mathematics published by Wiley Periodicals LLC