We have reached a high level of sophistication in the computational solution of partial differential equations (PDEs). They allow to model complex physical processes accurately in order to guide engineering design or to predict weather and climate. However, in practical applications the coefficients, the initial conditions or the geometry are typically uncertain or known only partially. Thus, it is crucial to calibrate PDE models on measured data or to assimilate such data into the models. But even high performance PDE tools can quickly reach their limits, when the solution varies on many scales and a large number of parameters in the model need to be calibrated. Yet, the accuracy achievable with fundamentally grounded PDE models is unsurpassable by purely data-driven machine learning approaches, especially in situations were accurate data is scarce or difficult to obtain. The most common way to overcome these challenges is to use surrogates or reduced order models both in the numerical modelling of the PDEs as well as in the modelling of the data likelihoods in Bayesian inference and learning. In this talk, I will argue that surrogates can help to tame or break the curse of dimensionality in Bayesian inference and then describe in more detail a highly effective, localised model-order reduction technique for multiscale PDEs based on generalised finite element methods (GFEM). This highly parallelisable approach localises the surrogate construction via a partition of unity and then identifies the optimal local approximation spaces as solutions of generalised eigenproblems. A new comprehensive convergence theory gives rigorous guarantees for the efficacy of the model reduction for a large range of important PDE models, including Helmholtz, Maxwell, convection-diffusion or Stokes.