An orientation distribution is a necessary input in any crystal plasticity simulation. The computational time involved in crystal plasticity simulations scales linearly with the number of crystal orientations in the input distributions. Reducing the number of crystal orientations in representing the input orientation distributions quantitatively is a critical and necessary requirement for performing computationally efficient crystal plasticity simulations of deformation processes. A procedure for the compaction of orientation distribution functions (ODFs) relying on a spectral representation using series of generalized spherical harmonics (GSH) basis functions was recently developed. Linear fitting of the spectral representation of an ODF containing a compact set of weighted orientations towards a full-size ODF containing many crystal orientations was in the core of the procedure. This paper advances the compaction procedure by replacing the linear with a nonlinear optimization in Matlab for which a suitably defined error, gradient, and Hessian matrix are derived to allow for more efficient, accurate, and greater compactions. The utility of the new procedure is to allow for not only fitting the weights of a compacted set of orientations but to allow for optimizing equally weighted orientations or a combination of optimizing weights and orientations. The new compaction procedure has been successfully applied to compactions of large ODFs of cubic, hexagonal, and orthorhombic polycrystalline metals. In doing so, the evolution of texture, twinning, and stress-strain are predicted at large plastic strains with compact ODFs to agree with the corresponding full size ODFs using crystal plasticity models. In closing, guidance for effective texture compaction trading off the accuracy and computational gains are provided.