We present an textit{ab initio} theory for superconductors, based on a unique mapping between the statistical density operator at equilibrium, on the one hand, and the corresponding one-body reduced density matrix $gamma$ and the anomalous density $chi$, on the other. This new formalism for superconductivity yields the existence of a universal functional $mathfrak{F}_beta[gamma,chi]$ for the superconductor ground state, whose unique properties we derive. We then prove the existence of a Kohn-Sham system at finite temperature and derive the corresponding Bogoliubov-de Gennes-like single particle equations. By adapting the decoupling approximation from density functional theory for superconductors we bring these equations into a computationally feasible form. Finally, we use the existence of the Kohn-Sham system to extend the Sham-Schluter connection and derive a first exchange-correlation functional for our theory. This reduced density matrix functional theory for superconductors has the potential of overcoming some of the shortcomings and fundamental limitations of density functional theory of superconductivity.