In this talk I'll discuss problems governed by the frequency domain wave (Helmholtz) equation, where the coefficients of the PDE (e.g. wave speed and material density) are parametrized random fields, and one wants to compute expected values or other moments of quantities derived from the PDE solutions. For the (simpler) random diffusion problem, methods based on (multilevel) Quasi-Monte Carlo methods can effectively handle problems with very high stochastic dimension. However, in the Helmholtz case, the (increasing) frequency acts as a `bad parameter', causing difficulty in the theory of the PDE itself, in the error estimates for QMC methods, and in the implementation of sampling methods in general. I'll start with a brief overview of some relevant theory for the heterogeneous Helmholtz problem and for the application of QMC methods to this problem. Then I'll discuss two pieces of recent work which aim to avoid some of the difficulties described. First we show that, if the Helmholtz system matrix corresponding to a given parameter choice is LU factorized, then the factors can be used as a preconditioner for other `nearby' problems (in parameter space) and the resulting iterative method can be frequency robust. We give a numerical illustration where QMC estimates for a 2D Helmholtz problem are computed (with frequency-independent iteration numbers and accuracy) using exact factorizations at about 1% of the QMC points needed for the estimator. In the second work we show how the application of `hybrid numerical-asymptotic methods' together with a new Filon-Smolyak method for multidimensional oscillatory integration leads to a new method for Helmholtz UQ in which the accuracy actually improves as frequency grows. At present this latter algorithm is only worked out for 1D Helmholtz problems with random refractive index. The talk contains joint work with Euan Spence, Owen Pembery, Frances Kuo, Dirk Nuyens, Ian Sloan, Zhizhang Wu and Zhiwen Zhang.