  In this talk, we present the QCD factorization formula for heavy quarkonium
production at large $p_T$ with factorized leading-power and next-to-leading
power contributions in the $1/p_T$ expansion. We show that the leading order
analytical calculations in this QCD factorization approach can reproduce
effectively the full next-to-leading order numerical results derived using
non-relativistic QCD (NRQCD) factorization formalism. We demonstrate that the
next-to-leading power contributions are crucial to the description of the
channels that are the most relevant for the rate as well as polarization of
$J/\psi$ production at current collider energies.
