Maximum likelihood estimation complexity computation
Ordinary Least Squares problem
Your
$$\hat x = \arg\min_x \sum_{n=1}^{N_r} \left\lvert y_n - h_n x\right\rvert^2$$
is just a way of saying
$$\hat x = \arg\min_x \left\| y - Hx \right \|$$
and that is called an Ordinary Least Squares estimation problem.
Its objective function is of a quadratic form, so we don't have to actually iteratively optimize or anything. There's an analytic solution to the above equation:
$$ \hat x = \left(H^T H\right)^{-1}H^T y\text.$$
As you can see, this has absolutely nothing to do with the constellation order so far – and it doesn't have to. The constellation order only makes a difference to the decision complexity, and by the orthogonality principle, that's a step to be done after you've estimated the transmit signal.
Don't invert that matrix
Now, it's practically never desirable to calculate a matrix inverse, so the thing that you actually solve is
$$\left(H^T H\right)\hat x = H^T y\text.$$
How complex the solution to that is depends on the properties of your MIMO channel matrix!
The quadratic / full-rank MIMO case
In the (very special) case that $H$ is normal (that implies among other things quadratic, ie. as many receive as transmit antennas $N_r=N_t=:N$), you do that by applying a QR decomposition to $H^T H$:
\begin{align}
\left(H^T H\right)\hat x &amp;= H^T y\\
QR \hat x &amp;= H^T y\\
R\hat x &amp;= Q^T H^T y\tag{1}\\
&amp;= \nu y\tag{2}
\end{align}
Since $R$ is a triangular matrix, $(1)$ s the point where you just backsubstitute. Luckily, $R$, $Q^T$ and $H^T$ only need to be calculated once ever, so we can estimate many $\hat x$ from many $y$ as long as the channel doesn't change at fixed complexity.
Complexity
You do the usual paper trick of ignoring the fact that you need to acquire $H$ first. That can, and quite possibly will, dominate the complexity of this (so don't be a bad scientist – figure out how you get your $H$ and how compelx that is, if you're writing some kind of publication).
So, assuming you already have the $H$, the operations you are doing per channel realization (i.e. only once every time $H$ changes, which probably doesn't happen very often):

calculating one half of the hermitian matrix $H^TH$ (no need to calculate the whole thing, it's redundant).
a QR decomposition of a hermitian matrix ($H^TH$)
calculating $\nu = Q^TH^T$

The operation you do per receive vector (i.e. for every symbol time step, typically, so pretty often):

Backsubstituting in $R\hat x$ to finally get $\hat x$.

The complexities per channel realization are:

$\frac{N^3}2$
depends very much on your QR decomposition algorithm (which depends on matrix size, and need for numerical stability, and your platform), but typically no worse than $N^4$
$N^3$

So that's a total of $C_H = \mathcal O \left(\frac 32 N^3 + N^4\right)$ multiply-accumulates.
The complexity  per receive vector:

$\frac12 N^2$

So, $C_y = \mathcal O \left(\frac12 N^2\right)$.
If we introduce a coherence time (measured in symbol durations) $T_C$ for which $H$ doesn't change, you get an overall complexity of
\begin{align}
C_{\text{overall}} &amp;= \mathcal O \left(\frac {\frac 32 N^3 + N^4}{N T_C} + \frac1{2N} N^2\right)\\
 &amp;= \mathcal O \left(\frac {\frac 32 N^2 + N^3}{T_C} + \frac12 N\right)
\end{align}
per receive symbol (i.e. per element of $\hat x$).
Bringing the constellation order into this
As said above, the constellation order only makes a difference to the complexity of the decider. However, we know under ML rule that the likeliest transmit symbol vector is the one geometrically closest to the likeliest transmit symbol vector – and since all elements of that are to be assumed independent from each other, that minimum can be found individually. The complexity of the decision thus only depends on the (SISO) constellation decoder you pick; nothing has been said about that, and there's immensely different implementations of these, so I will not allow myself to give you a complexity estimate here – I'll just allow myself that I can't see the formula as given by the paper appearing anywhere in my head.
It is a very valid question to ask whether it makes sense to apply the ML detection principle to individual elements of the transmit / receive vector. Unless you're in the very massive MIMO case, it makes sense to assume that there's channel coding on these symbols, and that they are, in fact, not stochastically independent from each other. So, classically, a MLSE-based (typically: Viterbi) soft-input channel decoder would be what you use, or some kind of message-passing based channel decoder.