Plug-and-play MMSE denoisers and stationarity

Since $0<{\rho }_{\sigma }(z-u)\le 1$, $p_{\sigma }(z)\le 1$. For fixed $z$, the measurable sets $E_k=\{u:{\rho }_{\sigma }(z-u)\ge 1\mathrm{/}k\}$, $k\ge 1$, cover ${\mathbb{R}}^n$. Their union is ${\mathbb{R}}^n$, so [F19] and [F2] give $1=\mu ({\mathbb{R}}^n)\le {\sum }_{k\ge 1}\mu (E_k)$; some $E_k$ therefore has positive $\mu$-measure, and $p_{\sigma }(z)\ge \mu (E_k)\mathrm{/}k>0$. For $j=1,2$, the function $u\mapsto \mathrm{\parallel }u{\mathrm{\parallel }}_2^j{\rho }_{\sigma }(z-u)$ is bounded because a polynomial factor is dominated at infinity by the Gaussian exponential, so its $\mu$-integral is finite. Put $w_{\sigma ,z}(u)={\rho }_{\sigma }(z-u)\mathrm{/}{p}_{\sigma }(z)$. This is a nonnegative Borel function. By [F14], $E\mapsto {\int }_E{w}_{\sigma ,z}d\mu$ is a measure; [F16] identifies it with ${\nu }_{\sigma ,z}$ and gives total mass $\int w_{\sigma ,z}d\mu =p_{\sigma }(z{)}^{-1}{p}_{\sigma }(z)=1$. Hence ${\nu }_{\sigma ,z}$ is a probability measure.

For a nonnegative simple Borel function $\phi ={\sum }_{j=1}^m{a}_j{\mathbf{1}}_{A_j}$ on disjoint Borel sets, $\phi w_{\sigma ,z}\in L^1(\mu )$ because $\phi$ is bounded and $w_{\sigma ,z}\in L^1(\mu )$ by step 1.1. By [F15] and [F17], $\int \phi d{\nu }_{\sigma ,z}={\sum }_{j=1}^m{a}_j{\nu }_{\sigma ,z}(A_j)=\int \phi w_{\sigma ,z}d\mu \mathrm{.}$ For arbitrary nonnegative Borel $\phi$, set ${\phi }_k(u)=2^{-k}\lfloor 2^k\min \{\phi (u),k\}\rfloor$. These are finite-valued Borel simple functions, increase pointwise, and converge to $\phi$. Applying [F18] under ${\nu }_{\sigma ,z}$ and $\mu$ yields $\int \phi d{\nu }_{\sigma ,z}={\lim }_k\int {\phi }_{k}d{\nu }_{\sigma ,z}={\lim }_k\int {\phi }_k{w}_{\sigma ,z}d\mu =\int \phi w_{\sigma ,z}d\mu =p_{\sigma }(z{)}^{-1}\int \phi {\rho }_{\sigma }(z-u)d\mu \mathrm{.}$ For real Borel $\phi$ with $\int \mathrm{\mid }\phi \mathrm{\mid }{\rho }_{\sigma }(z-u)d\mu <\mathrm{\infty }$, apply this identity to ${\phi }^{+}$ and ${\phi }^{-}$ and use [F17] to obtain the same formula for $\phi$. Taking $\phi (u)=\mathrm{\parallel }u{\mathrm{\parallel }}_2^2$ and using its weighted bound in step 1.1 proves the finite-second-moment claim; the finite first moment follows from $\mathrm{\parallel }u{\mathrm{\parallel }}_2\le 1+\mathrm{\parallel }u{\mathrm{\parallel }}_2^2$. Thus $D_{\sigma }(z)$ is well-defined.

Let $\alpha \in {\mathbb{N}}^n$ be a multi-index as in [F13]. At order zero, ${\mathrm{\partial }}^{\alpha }{\rho }_{\sigma }={\rho }_{\sigma }$. If ${\mathrm{\partial }}^{\alpha }{\rho }_{\sigma }(w)=P(w){\rho }_{\sigma }(w)$ for a polynomial $P$, then [F7] and [F8] give ${\mathrm{\partial }}_j(P(w){\rho }_{\sigma }(w))=({\mathrm{\partial }}_{j}P(w)-w_{j}P(w)\mathrm{/}{\sigma }^2){\rho }_{\sigma }(w)$, again a polynomial times ${\rho }_{\sigma }$. Induction on $\mathrm{\mid }\alpha \mathrm{\mid }$, in any coordinate order, proves the stated form for every finite derivative. On each compact set of $z$ values, both $\mathrm{\mid }{\mathrm{\partial }}^{\alpha }{\rho }_{\sigma }(z-u)\mathrm{\mid }$ and $\mathrm{\mid }{u}_i{\mathrm{\partial }}^{\alpha }{\rho }_{\sigma }(z-u)\mathrm{\mid }$ are bounded uniformly in $u$, since every polynomial times the Gaussian tends to zero at infinity. Apply [F3] coordinate by coordinate to $p_{\sigma }$ and to each numerator $N_i(z)=\int u_i{\rho }_{\sigma }(z-u)\mu (du)$; the uniform bounds supply an integrable majorant and [F4] justifies passage of each derivative limit through the integral. The same bounds for derivatives of order $\mathrm{\mid }\alpha \mathrm{\mid }+1$ give a uniform bound $M$ for the coordinate derivatives of each resulting derivative on any compact box $K$. For $z,w\in K$, change one coordinate at a time from $z$ to $w$. Each coordinate restriction is differentiable on an open interval containing its segment by [F3], hence continuous on the segment, so [F12] bounds its increment by $M\mathrm{\mid }{z}_j-w_j\mathrm{\mid }$. Summing gives a difference bound $M{\sum }_{j<n}\mathrm{\mid }{z}_j-w_j\mathrm{\mid }\le nM\mathrm{\parallel }z-w{\mathrm{\parallel }}_2$; thus each resulting derivative is continuous. Iterating over $\alpha$ makes $p_{\sigma }$ and every $N_i$ smooth. Since $p_{\sigma }>0$, the reciprocal $r(s)=1\mathrm{/}s$ is smooth on $(0,\mathrm{\infty })$: its first derivative is $-r(s{)}^2$ by the difference quotient, and repeated use of [F8] shows every higher derivative is a polynomial in $r$. Thus [F9] makes $D_{\sigma }=(N_i\mathrm{/}{p}_{\sigma }{)}_{i<n}$ smooth. Since $r=1\mathrm{/}s$ is smooth and [F6] gives $(\log {)}^{\mathrm{\prime }}=r$ on $(0,\mathrm{\infty })$, $\log$ is smooth there; then [F9] and $p_{\sigma }>0$ make $q_{\sigma }=-\log p_{\sigma }$ smooth.

For $a\in {\mathbb{R}}^n$, expand $u-a=(u-D_{\sigma }(z))+(D_{\sigma }(z)-a)$ and integrate. The cross term vanishes because $\int (u-D_{\sigma }(z)){\nu }_{\sigma ,z}(du)=0$, giving $\int \mathrm{\parallel }u-a{\mathrm{\parallel }}_2^2{\nu }_{\sigma ,z}(du)=\int \mathrm{\parallel }u-D_{\sigma }(z){\mathrm{\parallel }}_2^2{\nu }_{\sigma ,z}(du)+\mathrm{\parallel }a-D_{\sigma }(z){\mathrm{\parallel }}_2^2\mathrm{.}$ Thus $D_{\sigma }(z)$ is the unique posterior squared-error minimizer. Differentiating $N_i\mathrm{/}{p}_{\sigma }$ using [F3], [F5], [F7], [F8], and [F9] gives ${\mathrm{\partial }}_j{D}_{\sigma ,i}(z)=\frac{1}{{\sigma }^2}(\int u_i{u}_j{\nu }_{\sigma ,z}(du)-D_{\sigma ,i}(z)D_{\sigma ,j}(z)),$ because the terms containing $z_j$ cancel. This is exactly the stated symmetric covariance formula.

For each coordinate $j$, ${\mathrm{\partial }}_j{\rho }_{\sigma }(z-u)=(u_j-z_j){\rho }_{\sigma }(z-u)\mathrm{/}{\sigma }^2$. Differentiating under the integral by [F3] and using step 2.1 gives $\mathrm{\nabla }{p}_{\sigma }(z)=p_{\sigma }(z)(D_{\sigma }(z)-z)\mathrm{/}{\sigma }^2$. Since $p_{\sigma }(z)>0$, [F6] gives $\mathrm{\nabla }\log p_{\sigma }=\mathrm{\nabla }{p}_{\sigma }\mathrm{/}{p}_{\sigma }$, proving both Tweedie identities in the Statement.

For nonzero $v\in {\mathbb{R}}^n$, $v^{\mathsf{T}}J{D}_{\sigma }(z)v={\sigma }^{-2}\int ⟨v,u-D_{\sigma }(z){⟩}^2{\nu }_{\sigma ,z}(du)\mathrm{.}$ If this were zero, then for each $k\ge 1$ the set where $\mathrm{\mid }⟨v,u-D_{\sigma }(z)⟩\mathrm{\mid }\ge 1\mathrm{/}k$ would have ${\nu }_{\sigma ,z}$-measure zero, since the integral is at least $k^{-2}$ times that set's measure by [F16]; their union is the complement of the affine hyperplane $\{u:⟨v,u⟩=⟨v,D_{\sigma }(z)⟩\}$. Hence that hyperplane has full ${\nu }_{\sigma ,z}$-measure. The density ${\rho }_{\sigma }(z-u)\mathrm{/}{p}_{\sigma }(z)$ is strictly positive, so ${\nu }_{\sigma ,z}$ and $\mu$ have the same null sets: if ${\nu }_{\sigma ,z}(E)=0$, then for each $k$, [F16] and step 1.1 give $\mu (E\cap \{{\rho }_{\sigma }(z-u)\mathrm{/}{p}_{\sigma }(z)\ge 1\mathrm{/}k\})\le k{\nu }_{\sigma ,z}(E\cap \{{\rho }_{\sigma }(z-u)\mathrm{/}{p}_{\sigma }(z)\ge 1\mathrm{/}k\})=0$, and these sets cover $E$. Thus the hyperplane has full $\mu$-measure, contradicting the hypothesis. Therefore $J{D}_{\sigma }(z)$ is positive definite.

Fix $z$ and write ${\mathrm{\Psi }}_z(u)=\frac{1}{2}\mathrm{\parallel }{D}_{\sigma }(u)-z{\mathrm{\parallel }}_2^2-\frac{1}{2}\mathrm{\parallel }{D}_{\sigma }(u)-u{\mathrm{\parallel }}_2^2+{\sigma }^2{q}_{\sigma }(u)\mathrm{.}$ By the definition of $h_{\sigma ,\gamma }$, this equals $\frac{1}{2}\mathrm{\parallel }{D}_{\sigma }(u)-z{\mathrm{\parallel }}_2^2+\gamma h_{\sigma ,\gamma }(D_{\sigma }(u))$. Differentiate using [F7] and the symmetry of $J{D}_{\sigma }$ from step 4.1: $\mathrm{\nabla }{\mathrm{\Psi }}_z(u)=J{D}_{\sigma }(u)(D_{\sigma }(u)-z)-(J{D}_{\sigma }(u)-I)(D_{\sigma }(u)-u)+{\sigma }^2\mathrm{\nabla }{q}_{\sigma }(u)=J{D}_{\sigma }(u)(u-z),$ where the last equality uses $D_{\sigma }(u)-u=-{\sigma }^2\mathrm{\nabla }{q}_{\sigma }(u)$ from step 4.2.

If $z\mathrm{\ne }w$, put $v=w-z$ and $\psi (t)=⟨D_{\sigma }(z+tv),v⟩$ on $[0,1]$. By [F7] and step 5.1, ${\psi }^{\mathrm{\prime }}(t)=v^{\mathsf{T}}J{D}_{\sigma }(z+tv)v>0$. The mean-value theorem [F12] yields $⟨D_{\sigma }(w)-D_{\sigma }(z),w-z⟩>0$, so $D_{\sigma }(w)\mathrm{\ne }{D}_{\sigma }(z)$ and $D_{\sigma }$ is injective.

At every $z$, the derivative $J{D}_{\sigma }(z)$ is invertible by step 5.1. Apply [F11] at each $z$: every point has an open neighborhood mapped diffeomorphically onto an open neighborhood. Their images cover $X_{\sigma }$, so $X_{\sigma }$ is open. Injectivity from step 6.1 makes the local inverse maps agree on overlaps, giving the global inverse $D_{\sigma }^{-1}:X_{\sigma }\to {\mathbb{R}}^n$, which is smooth. Since $q_{\sigma }$ is smooth by step 3.1, the formula in the Statement makes $h_{\sigma ,\gamma }$ smooth on $X_{\sigma }$.

For any $v\mathrm{\ne }0$, let $\chi (t)={\mathrm{\Psi }}_z(z+tv)$ on $[0,1]$. Step 5.2 gives ${\chi }^{\mathrm{\prime }}(t)=t{v}^{\mathsf{T}}J{D}_{\sigma }(z+tv)v>0$ for every $t\in (0,1]$ by step 5.1. By [F12], ${\mathrm{\Psi }}_z(z+v)>{\mathrm{\Psi }}_z(z)$. Thus $u=z$ is the unique global minimizer of ${\mathrm{\Psi }}_z$. Every $x\in X_{\sigma }$ has a unique form $D_{\sigma }(u)$ by step 7.1, while $h_{\sigma ,\gamma }=+\mathrm{\infty }$ outside $X_{\sigma }$; consequently $D_{\sigma }(z)$ is the unique global minimizer of $\frac{1}{2}\mathrm{\parallel }x-z{\mathrm{\parallel }}_2^2+\gamma h_{\sigma ,\gamma }(x)$ over all $x\in {\mathbb{R}}^n$.

For a point $s\in {\mathbb{R}}^n$, define the extended-real quadratic model ${\mu }_s(x)=g(s)+⟨\mathrm{\nabla }g(s),x-s⟩+\frac{1}{2\gamma }\mathrm{\parallel }x-s{\mathrm{\parallel }}_2^2+h_{\sigma ,\gamma }(x)$. By [F3] and $\gamma \le 1\mathrm{/}L$, its smooth part majorizes $g(x)$ and is tight at $x=s$. Adding $h_{\sigma ,\gamma }$ preserves the majorization; when $h_{\sigma ,\gamma }(x)=+\mathrm{\infty }$ both objective values are $+\mathrm{\infty }$. Thus ${\mu }_s(x)\ge f(x)$ for every $x$, and ${\mu }_s(s)=f(s)$ in extended-real order.

For $t\ge 2$ put $s=x_{t-1}\in X_{\sigma }$ and define $r_t(x)=g(s)+⟨\mathrm{\nabla }g(s),x-s⟩+\frac{1}{2\gamma }\mathrm{\parallel }x-s{\mathrm{\parallel }}_2^2-g(x)$. By [F3], $r_t(x)\ge 0$ for every $x$, $r_t(s)=0$, and $\mathrm{\nabla }{r}_t(s)=0$. Its gradient is $\mathrm{\nabla }{r}_t(x)=\mathrm{\nabla }g(s)+(x-s)\mathrm{/}\gamma -\mathrm{\nabla }g(x)$, so the triangle inequality and positive homogeneity [F6], the Lipschitz bound [F4], and the hypothesis on $\mathrm{\nabla }g$ give $\mathrm{\parallel }\mathrm{\nabla }{r}_t(x)-\mathrm{\nabla }{r}_t(y){\mathrm{\parallel }}_2\le (1\mathrm{/}\gamma +L)\mathrm{\parallel }x-y{\mathrm{\parallel }}_2\le (2\mathrm{/}\gamma )\mathrm{\parallel }x-y{\mathrm{\parallel }}_2$. On $X_{\sigma }$, ${\mu }_s-f=r_t$.

Completing the square gives ${\mu }_s(x)=\frac{1}{2\gamma }\mathrm{\parallel }x-(s-\gamma \mathrm{\nabla }g(s)){\mathrm{\parallel }}_2^2+h_{\sigma ,\gamma }(x)+g(s)-\frac{\gamma }{2}\mathrm{\parallel }\mathrm{\nabla }g(s){\mathrm{\parallel }}_2^2$. By [F1] and the exact update in [F2] with ${\tilde{D}}_{\sigma }=D_{\sigma }$, $x_t$ is the unique global minimizer of ${\mu }_{x_{t-1}}$. Hence $f(x_t)\le {\mu }_{x_{t-1}}(x_t)\le {\mu }_{x_{t-1}}(x_{t-1})=f(x_{t-1})$. Also $x_t\in X_{\sigma }$ for $t\ge 1$, where $h_{\sigma ,\gamma }$ is finite, so $f(x_t)$ is finite.

The sequence $(f(x_t){)}_{t\ge 1}$ is finite, nonincreasing by step 2.1, and bounded below by $\inf f\in \mathbb{R}$. It therefore converges: if $\mathrm{\ell }={\inf }_{t\ge 1}f(x_t)$, then for every $\epsilon >0$ some $N$ has $f(x_N)<\mathrm{\ell }+\epsilon$, and every $t\ge N$ has $\mathrm{\ell }\le f(x_t)\le f(x_N)$.

Because $x_t$ minimizes ${\mu }_s$, step 2.1 and the identity in step 1.2 yield $0\le r_t(x_t)={\mu }_s(x_t)-f(x_t)\le f(x_{t-1})-f(x_t)$. Summing from $t=2$ to $N$ gives ${\sum }_{t=2}^N{r}_t(x_t)\le f(x_1)-f(x_N)\le f(x_1)-\inf f$. The nonnegative summands therefore tend to zero.

Apply [F3] to $r_t$ with gradient Lipschitz constant $2\mathrm{/}\gamma$ and point $x_t$, taking $u_t=x_t-\frac{\gamma }{2}\mathrm{\nabla }{r}_t(x_t)$. Then $0\le r_t(u_t)\le r_t(x_t)-\frac{\gamma }{4}\mathrm{\parallel }\mathrm{\nabla }{r}_t(x_t){\mathrm{\parallel }}_2^2$. Step 3.2 implies $\mathrm{\parallel }\mathrm{\nabla }{r}_t(x_t){\mathrm{\parallel }}_2^2\le 4r_t(x_t)\mathrm{/}\gamma \to 0$.

For $t\ge 2$, $x_t$ is a global minimizer of ${\mu }_s$ and belongs to the open set $X_{\sigma }$, where ${\mu }_s$ is differentiable. For every direction $v$, the global minimum gives nonnegative directional difference quotients along $v$ and $-v$ for small positive increments; taking limits gives $⟨\mathrm{\nabla }{\mu }_s(x_t),v⟩=0$. Thus $\mathrm{\nabla }{\mu }_s(x_t)=0$. Since $r_t={\mu }_s-f$ on $X_{\sigma }$, $\mathrm{\nabla }{r}_t(x_t)=-\mathrm{\nabla }f(x_t)$; step 4.1 now proves $\mathrm{\parallel }\mathrm{\nabla }f(x_t){\mathrm{\parallel }}_2\to 0$.

Suppose an extended-real penalty $h$ that is finite somewhere and nowhere $-\mathrm{\infty }$ had $P(z)=-z$ as a global minimizer of $\frac{1}{2}\mathrm{\mid }x-z{\mathrm{\mid }}^2+\gamma h(x)$ for every $z$. Its objective has a finite competitor for every input, so each claimed global minimizer has finite value. For distinct $s,t\in \mathbb{R}$, compare the minimizer at input $s$ with competitor $P(t)$ and the minimizer at input $t$ with competitor $P(s)$, then add the inequalities. The finite penalty values cancel and expansion gives $(P(s)-P(t))(s-t)\ge 0$.

Substituting $P(s)-P(t)=-(s-t)$ into step 1.1 yields $-(s-t{)}^2\ge 0$, which is false for $s\mathrm{\ne }t$. Thus no such $h$ exists for this map, so the exact-MMSE construction and stationarity result cannot be transferred to an arbitrary denoiser without a separate proof.