Summary
OP's problem has a polynomial-time algorithm via reduction to min-cost bipartite matching. (Lemma 1, below.)
Alternatively, one can strengthen OP's relaxation QP directly, by modeling the cost more carefully, to make QP it into a linear program (LP) with no integrality gap, so that any basic feasible solution to the LP gives an optimal solution to the original program (IQP). (Lemma 2, below.)
An approximate solution to IQP can be obtained by randomized rounding: 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|$.) (Lemmas 4 and 5, below.)
1. Reduction to min-cost matching
Lemma 1. The problem in the post has a polynomial-time algorithm via reduction to min-cost bipartite matching.
Proof of Lemma 1. Given an instance $G=(V,E)$ of the problem in the post, the reduction constructs an edge-weighted bipartite graph $G'=(U, W, E')$ as follows.
For each edge $e\in E$, add a vertex $\alpha_{e}$ to $U$.
For each vertex $v\in V$, add $d_G(v)$ vertices $\beta^v_1, \beta^v_2, \ldots, \beta^v_{d_v}$ to $W$, where $d_G(v)$ is the degree of $v$.
For each edge $e\in E$ and vertex $v\in e$, add edges $(\alpha_e, \beta^v_i)$ for $i\in[d_G(v)]$ to $E'$. Give each edge $(\alpha_e, \beta^v_i)$ cost $i^2-(i-1)^2$.
The problem is then to find a minimum-cost maximum matching $M$ in $G'$, that is, a matching of minimum cost among those having exactly one edge incident to each vertex $\alpha_e\in U$. This problem can be solved in polynomial time by standard algorithms for min-cost matching, or by reduction to min-cost flow. Clearly the reduction can be implemented in polynomial time.
To see that the reduction is correct, note that, given any orientation of the edges of $G$, there is a corresponding matching $M$: for each edge $e\in E$, match the vertex $\alpha_e$ to the vertex $\beta^v_i$, where $v\in e$ is the vertex that the oriented $e$ leaves, and $\beta^v_i$ is chosen greedily to be the first not-already-matched vertex in $\beta^v_1, \beta^v_2, \ldots, \beta^v_{d_G(v)}$. Then $d_{\text{out}}(v)$ of these vertices are matched, at total cost $\sum_{i=1}^{d_{\text{out}}(v)} i^2 - (i-1)^2$, which (as the sum telescopes) equals $d_{\text{out}}(v)^2$, so the total cost of the matching $M$ in $G'$ equals $\sum_{v\in V} d_{\text{out}}(v)^2$, as desired.
Conversely, let $M$ be any matching in $G'$ among those with an edge incident to each vertex $\alpha_e\in U$. Construct a corresponding orientation of $E$ as follows. For each edge $e\in E$, let $\beta^v_i$ be the vertex that $\alpha_e$ is matched to (so $v\in e$). Then orient $e$ out of $v$. The cost of $M$ is then at least $\sum_{v\in V} d_{\text{out}}(v)^2$, because in $M$ for each $v\in V$ there are $d_{\text{out}}(v)$ edges matched to vertices in $\{\beta^v_i : i\in d_G(v)\}$, so (by the convexity of the function $i\mapsto i^2$) the cost of these edges in $M$ is at least $d_{\text{out}}(v)^2$. So the cost $\sum_{v\in V} d_{\text{out}}(v)^2$ of the orientation is at most the cost of $M$, as desired.
It follows that, given a minimum-cost matching in $G'$, the corresponding orientation in $G$ minimizes $\sum_{v\in V} d_{\text{out}}(v)^2$, as desired. $~~~\Box$
2. Adjusting the quadratic program yields an LP with no integrality gap
Following the idea underlying Lemma 1, consider modifying OP's quadratic program by replacing the quadratic objective $\sum_{u\in V} y_u^2$ by $\sum_{u\in V} f(y_u)$, where $f(z)$ is the piecewise-linear function with integer breakpoints such that $f(z) = z^2$ for integer $z$. (Note that $z^2 < f(z)$ for non-integer $z$, yet $f$ is still convex, so this strengthens the relaxation.)
This modified program can in fact be formulated as a linear program by introducing additional non-negative variables $\{\Delta^u_i : u\in V, i\in [d_G(u)]\}$ to model the cost, as follows:
$$
\begin{align}
\text{minimize} & \sum_{u\in V} \sum_{i=1}^{d_G(u)} (i^2-(i-1)^2)\Delta^u_i \\
(\forall u\in V) & \sum_{i=1}^{d_G(u)} \Delta^u_i \ge y_u\\
& \vdots~~~~~(\text{the rest is the same as OP's program QP})
\end{align}
$$
Lemma 2. The optimal basic feasible solutions to the above LP are 0/1 solutions.
Proof idea. One can show that the LP equivalent to the (standard relaxation of) the matching problem in Lemma 1, in that the solutions of the two programs correspond by a correspondence that preserves cost, integrality, and the property of being a basic feasible solution. $~~~\Box$
3. Approximation ratio for randomized rounding of QP (Lemmas 3-5)
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 3. 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 4. 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 3,
$$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 3 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 5. 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$
