In this contribution we investigate the inclusive double quarkonium production in $pp$ collisions at the Run 2 LHC energies. The cross sections are estimated using the Non-Relativistic QCD (NRQCD) factorization formalism, taking into account the color singlet and color octet contributions. We present predictions for the total cross sections for the $J/\Psi J/\Psi$ and $\Upsilon \Upsilon$ production, as well as for the rapidity and transverse momentum distributions. Our results indicate that these distributions are dominated by the color singlet contribution, with the color octet one being important only at large $p_T$. Finally, we compare our results with the predictions for the double quarkonium production by the double parton scattering mechanism.