Abstract: Shear wave splitting is perhaps the most unambiguous signature of the effect of anisotropic materials on the propagation of seismic waves. It has been used extensively to study anisotropy in the Earth, at global scales from the inner core to the tectonics of the uppermost mantle and crust, and at smaller scales for imaging deformation in hydrocarbon reservoirs. Well-established techniques exist for measuring shear wave splitting in a single (three-component) seismogram and more recently these have been extended to treat shear wave splitting in a tomographic fashion: determining non-uniform anisotropic models using large data sets of splitting measurements. Here, I propose an extension to a recent shear wave splitting tomography methodology which incorporates the data analysis into the inversion itself. This methodology uses a non-linear neighbourhood algorithm inversion to explore the parameter space defined by an anisotropic model consisting of a number of uniform domains. Each candidate model is assessed by applying the splitting it predicts to the entire data set. This approach is computationally expensive, but is highly amenable to parallelization. I apply the methodology to three simple synthetic cases to demonstrate the utility of the method. Finally, I apply the approach to the problem of inferring two-layer anisotropy from SKS splitting, which is a commonly attempted problem in global seismology. This uses data from the seismic station EKTN, where two-layer splitting has been previously inferred. This highlights some of the inherent trade-offs with such studies, and emphasizes the need to incorporate extra information to resolve these. This method is applicable to shear wave anisotropy analysis in a broad range of settings from global to reservoir scale.