Summary: A computational method is developed for evaluating the plastic strain gradient hardening term within a crystal plasticity formulation. While such gradient terms reproduce the size effects exhibited in experiments, incorporating derivatives of the plastic strain yields a nonlocal constitutive model. Rather than applying mixed methods, we propose an alternative method whereby the plastic deformation gradient is variationally projected from the elemental integration points onto a smoothed nodal field. Crucially, the projection utilizes the mapping between Lie groups and algebras in order to preserve essential physical properties, such as orthogonality of the plastic rotation tensor. Following the projection, the plastic strain field is directly differentiated to yield the Nye tensor. Additionally, an augmentation scheme is introduced within the global Newton iteration loop such that the computed Nye tensor field is fed back into the stress update procedure. Effectively, this method results in a fully implicit evolution of the constitutive model within a traditional displacement-based formulation. An elemental projection method with explicit time integration of the plastic rotation tensor is compared as a reference. A series of numerical tests are performed for several element types in order to assess the robustness of the method, with emphasis placed upon polycrystalline domains and multi-axis loading.