  We report on recent progress in testing the factorization formalism of
nonrelativistic 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 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\bar{p}$, $ep$,
$\gamma\gamma$, and $e^+e^-$ collisions.Three data sets not included in the
fit, from hadroproduction and from photoproduction 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 from the FNAL Tevatron exhibit severe disagreement. We
demonstrate that alternative LDME sets recently obtained in two other NLO NRQCD
analyses of $J/\psi$ yield and polarization, with different philosophies, also
fail to reconcile the Tevatron polarization data with the other available world
data.
