We report on recent progress in testing the factorization formalism of non-relativistic quantum chromodynamics (NRQCD) at next-to-leading order (NLO) for J/psi yield and polarization. We demonstrate that it is possible to unambiguously determine the leading color-octet (CO) long-distance matrix elements (LDMEs) in compliance with the velocity scaling rules through a global fit to experimental data of unpolarized J/psi production in pp, p (p) over bar, ep, gamma gamma, and e(+)e(-) collisions. Three data sets not included in the fit, from hadro-production and from photo-production in the fixed-target and colliding-beam modes,are nicely reproduced. The polarization observables measured in different frames at DESY HERA and CERN LHC reasonably agree with NLO NRQCD predictions obtained using the LDMEs extracted from the global fit, while measurements at the FNAL Tevatron exhibit severe disagreement. We demonstrate that the alternative LDME sets recently obtained, with different philosophies, in two other NLO NRQCD analyses of J/psi yield and polarization also fail to reconcile the Tevatron polarization data with the other available world data.