In this paper we discuss the extension of Gaussian beam analysis to multi-beam quasi-optical systems. The normal approach taken in analyzing the propagation of multiple beams is to ray trace the optical axes of the individual beams and then superimpose the beam evolution of the on-axis beam. This is a rather inelegant and cumbersome process. We show how an efficient paraxial description of off-axis beams using Gaussian beam modes can be developed. The computational feasibility of the beam mode description clearly depends on finding an on-axis mode set in terms of which any off-axis beam can be approximated by a modal sum of modest size. We show how the optimum choice is made, and illustrate the power of the approach with an example.