A current objective of low-energy nuclear theory is to build nonempirical nuclear energy density functionals (EDFs) from underlying internucleon interactions and many-body perturbation theory (MBPT). The density matrix expansion (DME) of Negele and Vautherin is a convenient method to map highly nonlocal Hartree-Fock expressions into the form of a quasi-local Skyrme functional with density-dependent couplings. In this work, we assess the accuracy of the DME at reproducing the nonlocal exchange (Fock) contribution to the energy. In contrast to the scalar part of the density matrix for which the original formulation of Negele and Vautherin is reasonably accurate, we demonstrate the necessity to reformulate the DME for the vector part of the density matrix, which is needed for an accurate description of spin-unsaturated nuclei. Phase-space-averaging techniques are shown to yield a significant improvement for the vector part of the density matrix compared to the original formulation of Negele and Vautherin. The key to the improved accuracy is to take into account the anisotropy that characterizes the local momentum distribution in the surface region of finite Fermi systems. Optimizing separately the DME for the central, tensor, and spin-orbit contributions to the Fock energy, one reaches a few-percent accuracy over a representative set of semi-magic nuclei. With such an accuracy at hand, one can envision using the corresponding Skyrme-like energy functional as a microscopically constrained starting point around which future phenomenological parametrizations can be built and refined.