Proof that the difference between the input and output of the proximal operator belongs to the subdifferential The proximal operator $\text{prox}$ is defined as follows, for a function $f$, at point $x$:
$$\text{prox}_f(x) = \underset{u}{\text{argmin}} (f(u) + \frac{1}{2} ||x-u||^2)$$
The subdifferential $\partial f(x)$ of a function $f$ at a point $x$ is the set of subgradients of $f$ at point $x$, where a subgradient is any vector such that the line generated by it is always under the function.
It is claimed that $\forall x \in \text{dom} (f), x - \text{prox}_f(x) \in \partial f(x)$, yet I can't figure out why.
 A: EDIT: the following is wrong, it uses that $(f+g)^*=f^*+g^*$ which is horrendously wrong! I cannot delete an accepted answer, unfortunately. Is there any way to un-accept it?

You might be looking for: $\text{prox}_f(x) - x \in \partial f^*(x)$, where * denotes convex conjugate.
A slightly hacky way to see this: since the convex conjugate of $\frac{\|\cdot\|^2}{2}$ is itself and $\partial \frac{\|\cdot\|^2}{2}(x) = x$,
$$\text{prox}_f(x) = \underset{u}{\text{argmax}} \langle x, u \rangle - \left( f + \frac{\|\cdot\|^2}{2} \right) \\
\in \partial \left( f + \frac{\|\cdot\|^2}{2} \right)^*(x)
= \partial f^*(x) + x$$
by using that the Fenchel inequality is saturated at the subdifferential.
(This also used the subdifferential sum rule, which applies since the domain of $\frac{\|\cdot\|^2}{2}$ is the whole space)
Alternatively, starting from $x - u \in \partial f(u)$ where $u = \text{prox}_f(x)$,
$$x \in \partial \left( f + \frac{\|\cdot\|^2}{2} \right)(u) \\
u \in \partial \left( f + \frac{\|\cdot\|^2}{2} \right)^*(x) 
= \partial f^*(x) + x$$
by using that $u \in \partial g(x) \iff x \in \partial g^*(u)$
