Accurate modelling of the primary beam is an important but difficult task in radio astronomy. For high dynamic range problems such as 21cm intensity mapping, small modelling errors in the sidelobes and spectral structure of the beams can translate into significant systematic errors. Realistic beams exhibit complex spatial and spectral structure, presenting a major challenge for beam measurement and calibration methods. In this paper series, we present a Bayesian framework to infer per-element beam patterns from the interferometric visibilities for large arrays with complex beam structure, assuming a particular (but potentially uncertain) sky model and calibration solution. In this first paper, we develop a compact basis for the beam so that the Bayesian computation is tractable with high-dimensional sampling methods. We use the Hydrogen Epoch of Reionization Array (HERA) as an example, verifying that the basis is capable of describing its single-element E-field beam (i.e. without considering array effects like mutual coupling) with a relatively small number of coefficients. We find that 32 coefficients per feed, incident polarization, and frequency, are sufficient to give percent-level and ∼10 per cent errors in the mainlobe and sidelobes respectively for the current HERA Vivaldi feeds, improving to $\sim 0.1{{\%}}$ and $\sim 1{{\%}}$ for 128 coefficients.