Here we develop a theoretical and practical framework for the tomographic inversion of shear wave splitting intensity measurements for anisotropic structure in the upper mantle using a model space search approach. Treating the anisotropic scatterers as a first order perturbation to the background isotropic state, we implement the Born approximation to compute the integral sensitivity kernels in a finite frequency framework. We implement a parametrization of the anisotropy based on insights from olivine elasticity and fabric development that involves three parameters (corresponding to the azimuth and dip of the anisotropic symmetry axis, plus a strength parameter). Previous work on finite-frequency shear wave splitting tomography has implemented a linearization technique to invert splitting intensity data for the spatial distribution of anisotropic scatterers. The inverse problem, however, is strongly non-linear in terms of several of the involved parameters (those that describe the orientation of the symmetry axis), and their variation is not of first order. Therefore, in the case of a realistic upper mantle where anisotropic structure varies in a complicated manner, a linearization technique may not be adequate. To ameliorate these problems, we implement a model space search approach (specifically, a Markov chain Monte Carlo with Gibbs sampling algorithm) to the tomographic inversion of splitting intensity data. This approach allows for the visualization of posterior probability distributions for anisotropic parameters in the inversion. We perform a suite of synthetic resolution tests to demonstrate the reliability of our method, using a station distribution from an actual deployment of a dense seismic network. These resolution tests show that anisotropic structure may be resolved up to a length scale of roughly 50 km with teleseismic SKS waves for station spacing of 10-15 km.