Solution (source code)

= Solution

For a differentiable <convex function> $F$ with $L$-<Lipschitz gradient> and a proper lower-semicontinuous <convex function> $H$, the <proximal gradient method> is
$$
\boxed{z^{k+1}=\operatorname{prox}_{\tau H}(z^k-\tau\nabla F(z^k)),\qquad0<\tau\le1/L.}
$$
The <proximal operator> is $\operatorname{prox}_{\tau H}(v)=\arg\min_z\{H(z)+\|z-v\|^2/(2\tau)\}$. To obtain an explicit iteration for <discrete isotropic total variation>, let $D=\nabla$ with the printed zero boundary differences and
$$
\mathcal P=\{p\in X\times X:|p_{ij}|_2\le1\text{ for every }i,j\}.
$$
The dual is minimization of $F(p)=\tfrac12\|g-D^*p\|_2^2$ over $\mathcal P$, and the primal minimizer is $u_*=g-D^*p_*$. Here $D^*=-\operatorname{div}$ must be the exact <adjoint operator> of these differences. For example,
$$
((D_x^+)^*p^1)_{ij}=\begin{cases}
-p^1_{1j},&i=1<N,\\
p^1_{i-1,j}-p^1_{ij},&1<i<N,\\
p^1_{N-1,j},&i=N>1,
\end{cases}
$$
with the analogous formula in $j$ for $(D_y^+)^*$; for $N=1$ both operators are zero. Equivalently use $p^1_{0j}=p^1_{Nj}=p^2_{i0}=p^2_{iN}=0$ in the backward-difference expression, ignoring unused edge components.

The <gradient> of the dual objective is $D(D^*p-g)$, with Lipschitz constant $\|D\|^2\le8$. Taking $H$ to be the <indicator functional of a constraint set> makes its <proximal operator> the <Euclidean projection onto a convex set>. Thus a valid explicit algorithm is
$$
v^k=p^k+\tau D(g-D^*p^k),\qquad
\boxed{p^{k+1}_{ij}=\frac{v^k_{ij}}{\max(1,|v^k_{ij}|_2)},\quad u^{k+1}=g-D^*p^{k+1},\quad0<\tau\le\tfrac18.}
$$
The bound follows by summing $(a-b)^2\le2a^2+2b^2$ over grid differences. Finite-dimensional <proximal gradient method> convergence gives a dual minimizer and convergence of the recovered $u^k$ to the unique primal minimizer. Projection is pointwise in the dual; applying pointwise <soft thresholding> directly to $Du$ would not compute the primal TV <proximal operator>.