Non-perturbative approaches to strongly interacting quantum field theories commonly provide access to Euclidean correlation functions in terms of numerical observations. From these imaginary-time data, real-time properties can in principle be computed by numerically solving a Fredholm integral equation, the Källén-Lehmann spectral representation. By extending standard Gaussian process regression to inference from indirect observations, one obtains a powerful non-parametric algorithm for the probabilistic treatment of such ill-conditioned linear inverse problems. Asymptotic constraints can then be incorporated analytically in the spectral reconstruction by applying Mercer's theorem. In this talk, I will introduce the general approach and present two applications: 1. Estimation of bound state masses from resonant interaction channels of gauge-fixed correlation functions. 2. Determination of non-perturbative S-matrix building blocks for scattering processes.