Heavy quarkonium production at high transverse momentum ($p_T$) in hadronic collisions is explored in the QCD factorization approach. We find that the leading power in the $1/p_T$ expansion is responsible for high $p_T$ regime, while the next-to-leading power contribution is necessary for the low $p_T$ region. We present the first numerical analysis of the scale evolution of coupled twist-2 and twist-4 fragmentation functions (FFs) for heavy quarkonium production and demonstrate that the QCD factorization approach is capable of describing the $p_T$ spectrum of hadronic $J/psi$ production at the LHC.