Recently, Kapustin, Willett and Yaakov have found, by using localization techniques, that vacuum expectation values of Wilson loops in ABJM theory can be calculated with a matrix model. We show that this matrix model is closely related to Chern-Simons theory on a lens space with a gauge supergroup. This theory has a topological string large N dual, and this makes possible to solve the matrix model exactly in the large N expansion. In particular, we find the exact expression for the vacuum expectation value of a 1/6 BPS Wilson loop in the ABJM theory, as a function of the 't Hooft parameters, and in the planar limit. This expression gives an exact interpolating function between the weak and the strong coupling regimes. The behavior at strong coupling is in precise agreement with the prediction of the AdS string dual. We also give explicit results for the 1/2 BPS Wilson loop recently constructed by Drukker and Trancanelli