We introduce and explain the matrix model description of Chern-Simons theory and present new results on the unitary matrix models that describe U(N) Chern-Simons theory on S³. We show that they are integral representations of symmetric polynomials, which establishes a relationship with (1/2)-BPS Wilson loops in N=4 supersymmetric Yang-Mills theory and with free fermions.