Discovering interaction effects on a response of interest is a fundamental problem faced in biology, medicine, economics, and many other scientific disciplines. In theory, Bayesian methods for discovering pairwise interactions enjoy many benefits -- including coherent uncertainty quantification, the ability to incorporate background knowledge, and desirable shrinkage properties. In practice, however, Bayesian methods are often computationally intractable for even moderate-dimensional problems. Our key insight is that many hierarchical models of practical interest admit a Gaussian process representation such that a posterior over all O(p^2) interactions need never be maintained explicitly, only a vector of O(p) kernel hyper-parameters. This implicit representation allows us to run MCMC over model hyper-parameters in time and memory linear in p per iteration. On datasets with a variety of covariate and parameter behaviors such as sparsity, we show that: (1) our method improves running time by orders of magnitude over naive applications of MCMC, (2) that our method offers improved Type I and Type II error relative to state-of-the-art LASSO-based approaches, and (3) that our method offers improved computational scaling in high dimensions relative to existing Bayesian and LASSO-based approaches.