The theory of Wilson loops for gauge theories with unitary gauge groups is formulated in the language of symmetric functions. The main objects in this theory are two generating functions, which are related to each other by the involution that exchanges an irreducible representation with its conjugate. Both of them contain all information about the Wilson loops in arbitrary representations as well as the correlators of multiply-wound Wilson loops. This general framework is combined with the results of the Gaussian matrix model, which calculates the expectation values of 1/2-BPS circular Wilson loops in N = 4 Super-Yang-Mills theory. General, explicit, formulas for the connected correlators of multiply-wound Wilson loops in terms of the traces of symmetrized matrix products are obtained, as well as their inverses. It is shown that the generating functions for Wilson loops in mutually conjugate representations are related by a duality relation whenever they can be calculated by a Hermitian matrix model. * localization techniques [8,9] have provided a wealth of new solutions, in particular for highly symmetric loop configurations in supersymmetric gauge theories. A paradigmatic case is the case of 1 2 -BPS circular Wilson loops in N = 4 Super-Yang-Mills theory with gauge group U(N ) or SU(N ). On the one hand, the holographic dual fully captures the planar approximation in the limit of large 't Hooft coupling λ [10-15], and a lot of effort has been dedicated to obtain corrections in 1/λ [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. On the other hand, localization reduces the calculation of the Wilson loops to the solution of a Gaussian matrix model [30][31][32][33][34], which, in principle, is exact in both, λ and N and provides an easier path to an asymptotic 1/N -expansion [35][36][37][38][39][40][41]. Wilson loops in N = 4 Super-Yang-Mills theory with fewer symmetries have been studied, e.g., in [42][43][44][45][46].Localization can be applied more generally in N = 2 Super-Yang-Mills theories. Interested readers are referred to the recent paper [47] and references therein.The purpose of this paper is to further develop a recent result [48], which relates the connected correlators of multiply-wound Wilson loops in N = 4 Super-Yang-Mills theory [40] to the exact solution of the corresponding Hermitian matrix model. This relation was worked out explicitly in [40] for the connected n-loop correlators with n ≤ 4 and was generalized in [48] to any n by recognizing the combinatorical pattern. In [48], it was conjectured that the relation has a deeper, group-theoretical, origin. It will be shown here that this is indeed true. To achieve this goal, it turns out to be most natural and effective to employ the framework of symmetric functions. Symmetric functions have already been used in the matrix model solution of the 1/2-BPS Wilson loops in general representations [34]. However, the generality of this framework does not seem to be widely appreciated. In fact, it allows to define the generating function(s) for ...