Unraveling the complexity of non-collinear magnetism in materials with strongly correlated electrons is a considerable task that requires developing and applying state of the art theories and computational methods. Using the Coury model Hamiltonian, which includes spin and orbital degrees of freedom and generalizes the collinear Stoner Hamiltonian, we derive an extension of the collinear LSDA+U approximation to non-collinear magnetic configurations and explore the magnetic ground state of the archetypal spin-orbit correlated oxide UO2. We show that parameterizing a non-collinear LSDA+U model requires only one parameter U, as opposed to the difference between the Hubbard and Stoner parameters U −J found in an earlier derivation based on a collinear model Hamiltonian. To find the magnetic ground state of UO2 in the non-collinear configuration space, we combine LSDA+U with a spin adiabatic occupation matrix approach, involving the construction of a magnetic energy surface that follows the adiabatic evolution of the occupation matrix as a function of the spin canting angle. Our results show that the strong spin-orbit coupling (SOC) is the key factor stabilizing the so-called 3k spin ordered magnetic ground state of UO2. Using a relativistic atomic Hamiltonian we find that the SOC strength is colossal, 1.49 eV per uranium atom, the largest value ever found in relativistic materials. This unusually strong SOC implies that the spin and orbital degrees of freedom are virtually inseparable. As a result, to derive and quantify spin-spin interactions it is necessary to adopt the pseudospin picture. By constructing an extended effective multipolar pseudospin Hamiltonian, we prove the significance of octupolar and dipole-dipole exchange couplings in establishing the 3k magnetic phase, consistent with the non-collinear spin arrangement, and associated with the non-canted orbital ordering of uranium f-orbitals. Importantly, our study reveals that despite the strong spin-lattice interactions in action in UO2, the cooperative Jahn-Teller instability does not contribute to the onset of the non-collinear 3k state, which remains the most favorable ground state even in the undistorted cubic lattice. Finally, we discuss the role of prefactor U in the LSDA+U scheme and provide evidence that the choice of this parameter has a substantial quantitative influence on the predicted properties of the oxide, in particular the magnetic exchange interactions and, perhaps trivially, on the band gap: the value of U computed fully ab initio by the constrained random phase approximation, U=3.46 eV, delivers a band gap of 2.11 eV in good agreement with experiment, and a balanced account of the other relevant energy scales.