We develop a theoretical and computational framework to study polarons in semiconductors and insulators from first principles. Our approach provides the formation energy, excitation energy, and wavefunction of both electron and hole polarons, and takes into account the coupling of the electron or hole to all phonons. An important feature of the present method is that it does not require supercell calculations, and relies exclusively on electron band structures, phonon dispersions, and electron-phonon matrix elements obtained from calculations in the crystal unit cell. Starting from the Kohn-Sham (KS) equations of density-functional theory, we formulate the polaron problem as a variational minimization, and we obtain a nonlinear eigenvalue problem in the basis of KS states and phonon eigenmodes. In our formalism the electronic component of the polaron is expressed as a coherent superposition of KS states, in close analogy with the solution of the Bethe-Salpeter equation for the calculation of excitons. We demonstrate the power of the methodology by studying polarons in LiF and Li2O2. We show that our method describes both small and large polarons, and seamlessly captures Frohlich-type polar electron-phonon coupling and non-Frohlich coupling to acoustic and optical phonons. To analyze in quantitative terms the electron-phonon coupling mechanisms leading to the formation of polarons, we introduce spectral decompositions similar to the Eliashberg spectral function. We validate our theory using both analytical results and direct calculations on large supercells. This study constitutes a first step toward complete ab initio many-body calculations of polarons in real materials.