Solution (source code)

= Solution

Compute the <convex conjugate> of $f_z$ at a general pair $(v,y)$. With $w=z+u-x$, the independent variables become $(x,w)$ and $u=w+x-z$, so
$$
\begin{aligned}
f_z^*(v,y)
&=\sup_{x,w}\{\langle v+y,x\rangle+\langle y,w\rangle
-\langle y,z\rangle-k(x)-h(w)\}\\
&=k^*(v+y)+h^*(y)-\langle y,z\rangle.
\end{aligned}
$$
Therefore the <Lagrange dual function> is
$$
\boxed{\psi_z(y)=\langle y,z\rangle-k^*(y)-h^*(y).}
$$
Using <strong duality> from the preceding solution gives
$$
\boxed{F(z)=\sup_y[\langle y,z\rangle-k^*(y)-h^*(y)],
\qquad F=(k^*+h^*)^*.}
$$
Equivalently the <conjugate of an infimal convolution> is $F^*=k^*+h^*$, and the continuous convex $F$ equals its <biconjugate>.

For practical <subgradient> computation, minimize the known convex dual objective
$$
k^*(y)+h^*(y)-\langle z,y\rangle.
$$
Every optimizer, and only an optimizer, belongs to $\partial F(z)$ by <Fenchel–Young inequality>. The full characterization is
$$
\boxed{\partial F(z)
=\arg\max_y\psi_z(y)
=\{y:z\in\partial(k^*+h^*)(y)\}.}
$$
This is <infimal-convolution dual subgradients>; the set is nonempty because $F$ is finite convex everywhere. If both conjugates are differentiable at the optimizer, solve $\nabla k^*(y)+\nabla h^*(y)=z$. For nonsmooth conjugates, use the displayed aggregate <subdifferential> or a convex minimization algorithm. Replacing it by $\partial k^*(y)+\partial h^*(y)$ requires the usual <subdifferential sum rule> qualification; it is not automatically justified solely by knowing the two conjugates.