We construct a hierarchy of loop equations for invariant circular ensembles. These are valid for general classes of potentials and for arbitrary inverse temperatures $ {\rm Re}\,\beta>0 $ and number of eigenvalues $ N $. Using matching arguments for the resolvent functions of linear statistics $ f(\zeta)=(\zeta+z)/(\zeta-z) $ in a particular asymptotic regime, the global regime, we systematically develop the corresponding large $ N $ expansion and apply this solution scheme to the Dyson circular ensemble. Currently we can compute the second resolvent function to ten orders in this expansion and also its general Fourier coefficient or moment $ m_{k} $ to an equivalent length. The leading large $ N $, large $ k $, $ k/N $ fixed form of the moments can be related to the small wave-number expansion of the structure function in the bulk, scaled Dyson circular ensemble, known from earlier work. From the moment expansion we conjecture some exact partial fraction forms for the low $ k $ moments. For all of the forgoing results we have made a comparison with the exactly soluble cases of $ \beta = 1,2,4 $, general $ N $ and even, positive $ \beta $, $ N=2,3 $.