Is that edge orientation optimization problem NP-hard? Is the following optimization problem NP-hard?

Problem. For a given undirected graph $G=(V,E)$, find an orientation of the edges that minimizes the objective value $\sum_\limits{u\in V} ~\left( d_{out}(u) \right)^2$, where $d_{out}$ is the number of out-neighbors.

Motivation: The problem is related to listing triangles in graphs. A motivation can be found here where the optimization problem is similar, but where the edge orientation needs to induce a DAG.
One can observe that the DAG constraint is not necessary: the actual constraint in order to have a correct algorithm (i.e. an algorithm that will list all triangles) is that the edge orientation should not induce a cycle of length 3.
I think that the constraint &quot;no induced cycle of length 3&quot; is too complicated and decided to formulate the above problem as a stepping stone.
Another formulation: The problem can be written as the following Integer Quadratic Programming.
\begin{align}
 \text{Minimize} &amp;\qquad \sum_{u\in V} (y_u)^2 &amp;\qquad  (IQP)\\
 \text{s.t.} &amp;\qquad \forall u \in V, \forall v \in N_u, x_{u,v}  \in \{0,1\}  \\
 &amp;\qquad \forall (u,v) \in E, x_{u,v}+x_{v,u}=1\\
 &amp;\qquad \forall u \in V, y_u=\sum_{v \in N_u} x_{u,v}\\
\end{align}
with $N_u$ the set of neighbors of node $u$ in the input undirected graph.
Is (IQP) NP-hard?
It can be relaxed into the following Quadratic Programming:
\begin{align}
 \text{Minimize} &amp;\qquad \sum_{u\in V} (y_u)^2 &amp;\qquad  (QP)\\
 \text{s.t.} &amp;\qquad \forall u \in V, \forall v \in N_u, x_{u,v}  \in [0,1]  \\
 &amp;\qquad \forall (u,v) \in E, x_{u,v}+x_{v,u}=1\\
 &amp;\qquad \forall u \in V, y_u=\sum_{v \in N_u} x_{u,v}\\
\end{align}
for which an optimal solution can be obtained in polynomial time using a quadratic programming solver as the objective is convex. Can an exact or approximate solution to (IQP) be obtained out of the (QP) solution? Using randomized rounding maybe?
The greedy algorithm is not optimal: The algorithm &quot;fix some orientation then repeatedly flip the orientation of an edge that decrease the score&quot; is not optimal. The following is a counter example found by alt-tab.
We answer OP's last question: can an approximate solution to IQP be obtained by randomized rounding?
We show that the natural randomized-rounding scheme gives a 2-approximation, and a $(1+1/\overline d)$-approximation in graphs with average degree $\overline d$. (Note that any connected graph has average degree at least $2-2/|V|$.)
For the record, here's the rounding scheme. Given a fractional solution $x$ for QP, define 0/1 solution $x'$ as follows: for each edge $(u,w)$ independently: take $x'_{uw} = 1$ and $x'_{wu}=0$ with probability $x_{uw}$; otherwise take $x'_{uw} = 0$ and $x'_{wu} = 1$.
Let $f(x)$ be the objective function of QP.
Lemma 1. For any feasible $x$,
$$E[f(x')] = f(x) + 2\sum_{(u,w)\in E} x_{uw}x_{wu} \le f(x) + |E|/2.$$
Proof. This is a standard calculation.  Note that
$$f(x) = \sum_{u\in V} \Big(
\sum_{w\in N(u)} x^2_{uw} + 2\sum_{w'\in N(u), w'\ne w} x_{uw}x_{uw'}\Big)$$
and by linearity of expectation, the independence of $x'_{uw}$ and $x'_{uw'}$, and $E[x'^2_{uw}] = x_{uw}$,
$$E[f(x')] = \sum_{u\in V} \Big(
\sum_{w\in N(u)} x_{uw} + 2\sum_{w'\in N(u), w'\ne w} x_{uw}x_{uw'}\Big).$$
Using this and $x_{uw} = 1 - x_{wu}$,
$$
E[f(x')] - f(x) 
= \sum_{(u,w)\in E} x_{uw} - x^2_{uw} + x_{wu} - x^2_{wu}
= 2\sum_{(u,w)\in E} x_{uw} x_{wu}.$$
Lemma 2. For any feasible $x$, $E[f(x')] \le 2 f(x)$. That is, the rounding scheme gives a 2-approximation.
Proof.
By inspection of the expansion of $f(x)$ in the proof of Lemma 1,
$$f(x) \ge \sum_{(u,w)\in E} x^2_{uw} + x^2_{wu}.$$
Applying $p^2 + (1-p)^2 \ge 2p(1-p)$ for $p\in[0,1]$ (with $p=x_{uw}$)
to each term gives
$$f(x) \ge 2\sum_{(u,w)\in E} x_{uw}x_{wu}.$$
This and Lemma 1 imply $E[f(x')] \le f(x) + f(x) = 2 f(x)$. $~~~\Box$
The bound is tight in that for any graph with maximum degree 1,
the fractional solution $x_{uw} = x_{wu} = 1/2$ obtains value $|E|/2$,
but any integer solution has value $|E|$.
(Also note that $|E|$ is a lower bound on OPT,
using $\sum_{u\in V} d^2_{out}(v) \ge \sum_{u\in V} d_{out}(v) = |E|$.)
Lemma 3. In any graph with average degree $\overline d = 2|E|/|V|$, the rounding scheme gives at least a $(1+1/\overline d)$-approximation.
Proof.  Fix an optimal solution. Let $d'_v$ denote the out-degree of vertex $v$ in the optimal solution,
so $d'$ is a vector in $\mathbb R_+^{|V|}$.
Let $\mathbf 1$ denote the all-ones vector in $\mathbb R^{|V|}$.
By Cauchy-Schwarz,
$$\text{OPT} = d' \cdot d' 
\ge (d' \cdot \mathbf 1)^2/(\mathbf 1\cdot \mathbf 1) 
= \textstyle \big(\sum_v d'_v\big)^2/|V| 
= |E|^2/|V| = |E|\overline d/2.$$
By this and Lemma 1
$$E[f(x')] \le f(x) + |E|/2 \le f(x) + \text{OPT} / \overline d \le \text{OPT}(1+1/\overline d).$$
$~~~\Box$