We have derived and implemented analytical gradients for broken-symmetry unrestricted density functional calculations (BS-UDFT) with removal of spin contamination by Yamaguchi's approximate spin projection method. Geometry optimizations with these analytical gradients (AGAP-opt) yield results consistent with those obtained with the previously available numerical gradients (NAP-opt). The AGAP-opt approach is found to be more precise, efficient, and robust than NAP-opt. It allows full geometry optimizations for large open-shell systems. We report results for three types of organic diradicals and for a binuclear vanadium(II) complex to demonstrate the merits of removing the spin contamination effects during geometry optimization (AGAP-opt vs BS-UDFT) and to illustrate the superior performance of the analytical gradients (AGAP-opt vs NAP-opt). The results for the vanadium(II) complex indicate that the AGAP-opt method is capable of handling pronounced spin contamination effects in large binuclear transition metal complexes with two magnetic centers.