We use localization techniques to compute the expectation values of supersymmetric Wilson loops in Chern-Simons theories with matter. We find the path-integral reduces to a non-Gaussian matrix model. The Wilson loops we consider preserve a single complex supersymmetry, and exist in any theory, though the localization requires superconformal symmetry. We present explicit results for the cases of pure Chern-Simons theory with gauge group U (N), showing agreement with the known results, and ABJM, showing agreement with perturbative calculations. Our method applies to other theories, such as Gaiotto-Witten theories, BLG, and their variants.