Cost of solving a matrix equation using the FFT
If $m \times m$ symmetric matrix $\rm H^{\top} H$ is circulant, then its spectral decomposition is
$$\rm H^{\top} H = Q \Lambda Q^*$$
where the eigenvalues of $\rm H^{\top} H$ are given by the Discrete Fourier Transform (DFT) of the first row of $\rm H^{\top} H$, and where the eigenvalue matrix is
$$\rm Q =  \frac{1}{\sqrt m} \, F_m$$
where $\rm F_m$ is the $m \times m$ Fourier matrix. Hence,
$$\rm H^{\top} H + I_m = Q \Lambda Q^* + I_m = Q \Lambda Q^* + Q Q^* = Q \, (\Lambda + I_m) \, Q^*$$
and
$$\rm  \left( H^{\top} H + I_m \right)^{-1} U = Q \, (\Lambda + I_m)^{-1} \, Q^* U = \frac 1m \, F_m \, (\Lambda + I_m)^{-1} \, F_m^* U$$
The Fast Fourier Transform (FFT) algorithm can now be used to factor the Fourier matrix $\rm F_m$ and its Hermitian transpose. What is the cost?

Computing $\Lambda$ requires one DFT of length $m$. Using the FFT, that costs $O (m \log (m))$.
Computing $\rm F_m^* U$ requires $n$ DFTs of length $m$. Using the FFT, that costs $O (n \, m \log (m))$.
Multiplying $\rm (\Lambda + I_m)^{-1}$ and $\rm F_m^* U$ costs $m$ additions and $m \, n$ divisions.
Computing $\rm F_m \, (\Lambda + I_m)^{-1} \, F_m^* U$ requires $n$ DFTs of length $m$. Using the FFT, that costs $O (n \, m \log (m))$.
Dividing $\rm F_m \, (\Lambda + I_m)^{-1} \, F_m^* U$ by $m$ costs $m \, n$ divisions.

Thus, the total cost is $O (n \, m \log (m))$.

References

Robert M. Gray, Toeplitz and circulant matrices: a review.
Complex matrices; fast Fourier transform, MIT 18.06SC Linear Algebra, Fall 2011.