$\newcommand{\scp}[2]{\langle #1,#2\rangle}\newcommand{\id}{\mathrm{Id}}$ Let $f$ and $g$ be two proper, convex and lower semi-continuous functions (on a Hilbert space $X$ or $X=\mathbb{R}^n$) and let $g$ be continuously differentiable. Consequently, the subdifferential $\partial f$ and the gradient $\nabla g$ are monotone operators, i.e. $$ \scp{\nabla g(x)-\nabla g(y)}{x-y}\geq 0 $$ for all $x,y$ and for $u\in\partial f(x)$ and $v\in\partial f(y)$ $$ \scp{u-v}{x-y}\geq 0. $$ My question:

Is the operator $x\mapsto x - (\id + \gamma\partial f)^{-1}(x-\gamma\nabla g(x))$ monotone for $\gamma\geq 0$?

Note that I ask for all values $\gamma\geq 0$ specifically for large values.

An equivalent question is

Does for $\gamma\geq 0$ it holds that $$ \scp{(\id + \gamma\partial f)^{-1}(x-\gamma\nabla g(x)) - (\id + \gamma\partial f)^{-1}(y-\gamma\nabla g(y))}{x-y}\leq \|x-y\|^2 $$

Thoughts:

For $f=0$ and $g=0$ it's clear (for $g=0$ this follows since the "proximal operator" $x\mapsto (\id + \gamma\partial f)^{-1}(x)$ is non-expansive, for $f=0$, the monotonicity of $\gamma\nabla g$ works in the right direction and does the trick).

Also for small $\gamma$ (smaller that $2/L$ if $L$ is the Lipschitz constant of $\nabla g$, if is has one) the thing is clear as then the mapping $x\mapsto (\id + \gamma\partial f)^{-1}(x - \gamma\nabla g(x))$ is again non-expasive (and used as iteration in the so-called proximal gradient method). However, for large $\gamma$ I could not make use of any of these observation since using Cauchy-Schwarz for the inner product ruins the estimate then.

Moreover all examples I tried numerically (in various dimensions and for various functions $f$ and $g$) suggested that the claim holds.

Intuitively, all these together makes me think that the answer to my questions is yes, but I failed to prove it. Also all inequalities in Bauschke/Combettes "Convex analysis and Monotone Operator Theory in Hilbert spaces" I found were not helpful.