The dynamics of magnetization relaxation in ferrofluids are studied with statistical-mechanical theory and Brownian dynamics simulations. The particle dipole moments are initially perfectly aligned, and the magnetization is equal to its saturation value. The magnetization is then allowed to decay under zero-field conditions toward its equilibrium value of zero. The time dependence is predicted by solving the Fokker-Planck equation for the one-particle orientational distribution function. Interactions between particles are included by introducing an effective magnetic field acting on a given particle and arising from all of the other particles. Two different approximations are proposed and tested against simulations: a first-order modified mean-field theory and a modified Weiss model. The theory predicts that the short-time decay is characterized by the Brownian rotation time τB, independent of the interaction strength. At times much longer than τB, the asymptotic decay time is predicted to grow with increasing interaction strength. These predictions are borne out by the simulations. The modified Weiss model gives the best agreement with simulation, and its range of validity is limited to moderate, but realistic, values of the dipolar coupling constant.