To be able to write this into a matrix form, we need a priori to escape from the many-body picture. To do so we are going to suppose that for a reason due to the physics of the problem, the spin-down fermions have no dynamics (the spin is represented by $\sigma$, we are facing a Fermi-Hubbard model). This means hat we can replace $\hat n_{\downarrow}$ by its mean value $n_{\downarrow}$.
From this we see that the total Hamiltonian $\hat H$ can be decomposed in the following way:
$$ T^{(\uparrow)} = -\sum_{<i,j>} \hat c_{i\uparrow} \hat c^\dagger_{j\uparrow} + H.c. $$
$$ V^{(\uparrow)} = \sum_i \hat n_{i\uparrow} V_i.$$
Their spin-down counterpart are simply (sorry for the redundancy):
$$ T^{(\downarrow)} = -\sum_{<i,j>} \hat c_{i\downarrow} \hat c^\dagger_{j\downarrow} + H.c. $$
$$ V^{(\downarrow)} = \sum_i \hat n_{i\downarrow} V_i $$
$$ U^{(\uparrow)} = U \sum_i \hat n_{i\uparrow} n_{i\downarrow}.$$
This allows us to separate the spin-up and spin-down Hamiltonians:
$$ \hat H^{(\uparrow)} = T^{(\uparrow)} + V^{(\uparrow)} + U^{(\uparrow)}$$
$$ \hat H^{(\downarrow)} = T^{(\downarrow)} + V^{(\downarrow)}. $$
Note that $\hat n_{i\downarrow}$ has lost its hat. By writing this we assume that there is no back-action from the spin-up fermions on the spin-down fermions. Here we just used the assumption of no dynamics for the spin-down fermions.
From this we can first diagonalize $\hat H^{(\downarrow)}$ and pick our favorite eigenstate or linear combination of eigensates, take the absolute squared value (which is the density $n_{\downarrow}$) and plug it into $\hat H^{(\uparrow)}$. We then diagonalize $\hat H^{(\uparrow)}$ if we want to know its spectrum. Or we use $\hat H^{(\uparrow)}$ to compute the time-evolution of the spin-up fermions in presence of spin-down fermions as a background.
