We generalize the method of computing functional determinants with a single excluded zero eigenvalue developed by McKane and Tarlie to differential operators with multiple zero eigenvalues. We derive general formulas for such functional determinants of $r\times r$ matrix second order differential operators O with $0 < n \leqslant 2r$ linearly independent zero modes. We separately discuss the cases of the homogeneous Dirichlet boundary conditions, when the number of zero modes cannot exceed r, and the case of twisted boundary conditions, including the periodic and anti-periodic ones, when the number of zero modes is bounded above by 2r. In all cases the determinants with excluded zero eigenvalues can be expressed only in terms of the n zero modes and other $r-n$ or $2r-n$ (depending on the boundary conditions) solutions of the homogeneous equation $O h=0$ , in the spirit of Gel'fand–Yaglom approach. In instanton calculations, the contribution of the zero modes is taken into account by introducing the so-called collective coordinates. We show that there is a remarkable cancellation of a factor (involving scalar products of zero modes) between the Jacobian of the transformation to the collective coordinates and the functional fluctuation determinant with excluded zero eigenvalues. This cancellation drastically simplifies instanton calculations when one uses our formulas.