The calculation of theoretical (synthetic) seismograms in seismology is not done for its own sake. In most branches the ultimate goal is to compare the theoretical calculations with observations. In this data fitting process (the inverse problem) seismograms are calculated for earth models through iterative procedures or extensive searches of the model space. So far this model search was carried out using predominantly linear or linearized approaches. Arguably the most elegant method to solve inverse problems is the probabilistic approach whereby a probability density function is calculated that is descriptive of the information available accounting for observations and uncertainties in theory and observations. Even though we are aware that with complete 3D solutions we are far away from being able to employ Monte Carlo type approaches to solve seismic inverse problems on a large scale our goal is to investigate whether – for example through combination of techniques mentioned above – we can start to go beyond linearized inversion using full waveforms. This will involve the use of clever global search algorithms such as simulated annealing, genetic algorithms, evolutionary programming or neighbourhood algorithm.