text:
\section{Introduction}
The Plug-and-Play (PnP) ADMM is a variant of the alternating direction method of multiplier (ADMM) algorithm. Since its introduction in 2013 \cite{Venkatakrishnan_Bouman_Wohlberg_2013}, the algorithm has demonstrated extremely promising results in image restoration and signal recovery problems \cite{Dar_Bruckstein_Elad_2016, Metzler_Maleki_Baraniuk_2016,Rond_Giryes_Elad_2016,Kamilov_Mansour_Wohlberg_2017,Chang_Marchesini_2016,Ono_2017,Sreehari_Venkatakrishnan_Wohlberg_2016}. However, despite the enormous applications of PnP ADMM and several studies on its convergence \cite{Sreehari_Venkatakrishnan_Wohlberg_2016, Chan_Wang_Elgendy_2017, Teodoro_Bioucas_Figueiredo_2018}, it is generally unclear why the algorithm is performing so well. On the one hand, one can argue that the superior performance is attributed to the underlying image denoisers, which is evident from applications using convolutional neural networks \cite{Heide_Steinberger_Tsai_2014,Zhang2017_cvpr,Zhang2017_tip,Borgerding_Schniter_Rangan_2017}. Yet on the other hand, unless we can explicitly write down the optimization which describes the interaction of the denoisers, e.g., in the form of maximum-a-posteriori (MAP), it will be extremely difficult to quantify the solution even if we know the solution exists. The contribution of this paper is to present a set of formal performance analysis results for a family of denoisers known as the graph filters, so that we can understand the PnP ADMM behavior.

In order to state the problem more concretely, it would be useful to first review the PnP ADMM. The starting point of the PnP ADMM is the ADMM algorithm, which has become a standard tool for minimizing a sum of two separable functions of the form
\begin{equation}
\vxhat = \argmin{\vx} \;\; f(\vx) + \lambda g(\vx).
\label{eq:basic optimization}
\end{equation}
Here, $f$ represents the \emph{objective function} arising from the physical measurement models, e.g., blur, sampling, projection, or other linear transformations. The function $g$ is the \emph{regularization function}, representing the prior model of the underlying signals, e.g., $\ell_p$-norm, sparsity, or total variation. The parameter $\lambda$ is the regularization constant, and is assumed known a-priori and is fixed in this paper. Problems taking the form of \eref{eq:basic optimization} is broad, encompassing various linear inverse problems in imaging, deblurring, super-resolution, CT, MRI, and compressed sensing, to name a few.

ADMM solves \eref{eq:basic optimization} by converting the unconstrained optimization into a constrained problem
\begin{equation}
(\vxhat,\vvhat) = \argmin{\vx,\vv} \;\; f(\vx) + \lambda g(\vv), \;\; \subjectto \;\; \vx = \vv,
\label{eq:basic optimization 2}
\end{equation}
and considers the augmented Lagrangian function:
\begin{equation}
\calL(\vx,\vv,\vu) = f(\vx) + \lambda g(\vv) + \vu^T(\vx-\vv) + \frac{\rho}{2}\|\vx - \vv\|^2.
\label{eq:augmented lagrangian function}
\end{equation}
Then, the algorithm finds the solution by seeking a saddle point of $\calL$, which involves solving a sequence of subproblems in the form
\begin{align}
\vx^{(k+1)} &= \argmin{\vx\in \R^n} \;\; f(\vx) +\frac{\rho}{2} \|\vx - \vxtilde^{(k)}\|^2, \label{eq:ADMM2,x}\\
\vv^{(k+1)} &= \argmin{\vv\in \R^n} \;\; \lambda g(\vv) + \frac{\rho}{2}\|\vv - \vvtilde^{(k)}\|^2,\label{eq:ADMM2,v}\\
\vubar^{(k+1)} &= \vubar^{(k)} + (\vx^{(k+1)} - \vv^{(k+1)}),\label{eq:ADMM2,u}
\end{align}
where $\vubar^{(k)} \bydef (1/\rho)\vu^{(k)}$ is the scaled Lagrange multiplier, $\vxtilde^{(k)} \bydef \vv^{(k)}-\vubar^{(k)}$ and $\vvtilde^{(k)} \bydef \vx^{(k+1)}+\vubar^{(k)}$ are the intermediate variables. Under mild conditions, e.g., when both $f$ and $g$ are closed, proper and convex, and if a saddle point of $\calL$ exists, one can show that the iterates returned by \eref{eq:ADMM2,x}-\eref{eq:ADMM2,u} converges to the solution of \eref{eq:basic optimization}. Readers interested in knowing more about the theoretical properties of ADMM can consult tutorials such as \cite{Boyd_Parikh_Chu_Peleato_Eckstein_2011}.

The idea of PnP ADMM is to modify \eref{eq:ADMM2,v} by observing that it is essentially a denoising step, if we treat $\vvtilde^{(k)}$ as a ``noisy'' version of $\vv$ and $g(\vv)$ as a regularization for $\vv$. Based on this observation, we can replace \eref{eq:ADMM2,v} by a denoiser $\calD_\sigma: \R^n \rightarrow \R^n$ such that
\begin{equation}
\vv^{(k+1)} = \calD_\sigma\left( \vvtilde^{(k)} \right),
\label{eq:ADMM2,v rewrite 2}
\end{equation}
where $\sigma = \sqrt{\lambda/\rho}$ is the denoising strength or a hypothesized ``noise level''. The choice of $\calD_\sigma$ is broad. $\calD_\sigma$ can be a proximal map such as total variation denoising, or an off-the-shelf image denoisers such as BM3D, non-local means, and more recently neural networks.

Replacing the explicit optimization in \eref{eq:ADMM2,v} by an off-the-shelf denoiser in \eref{eq:ADMM2,v rewrite 2} leads to many interesting questions:
\begin{itemize}
\item \textbf{Existence of $g$}. Given an arbitrary denoiser $\calD_\sigma$, is it possible to have an equivalent $g$ such that \eref{eq:ADMM2,v} and \eref{eq:ADMM2,v rewrite 2} coincides? If not, under what conditions of $\calD_\sigma$ would this happen? Some recent studies provide an answer to this question \cite{Sreehari_Venkatakrishnan_Wohlberg_2016, Reehorst_Schniter_2018}, which will be elaborated in later sections of this paper. An alternative approach is the Regularization by Denoising (RED) \cite{Romano_Elad_Milanfar_2017}, which has advantages and limitations.
\item \textbf{Convergence}. Does PnP ADMM converge? Clearly the convergence does not hold for arbitrary denoisers. For global convergence, the most basic requirement is that $\calD_{\sigma}$ has symmetric gradient and is non-expansive \cite{Sreehari_Venkatakrishnan_Wohlberg_2016}. For fixed point convergence, $\calD_{\sigma}$ should be asymptotically invariant as $\sigma \rightarrow 0$ \cite{Chan_Wang_Elgendy_2017}.
\item \textbf{Performance}. If $g$ does exist, what is it and why is it an effective regularization function? If $g$ does not exists, then what does PnP ADMM solve, and why is it good? The focus of this paper is to address this performance issue. We will take two paths, one through the MAP formulation, and the other one through a concept called the consensus equilibrium.
\item \textbf{Generalization}. How to generalize the algorithm to multiple denoisers? The generalization can be done using the consensus equilibrium technique. However, can we explain the performance, and what would be the optimal combination weights?
\end{itemize}
The main assumption we make in this paper is that the denoisers are in the family of symmetric smoothing filters \cite{Milanfar_2013a, Milanfar_2013b, Chan_Zickler_Lu_2017} or graph filters \cite{Meyer_Shen_2012,Talebi_Zhu_Milanfar_2013,Taylor_Meyer_2012,Kheradmand_Milanfar_2014,Talebi_Milanfar_2014,Pang_Cheung_2017,Bai_Cheung_Liu_2018}. These filters are simple in the sense that they are (pseudo) linear and hence they can be expressed using matrices. However, they are also representative enough to cover a sufficiently large group of image denoisers that we encounter in the literature.

The remaining of the paper is arranged as follows. We begin with a brief review of the construction of graph filters and their properties in Section 2. The main results are presented in Section 3 and 4. Section 3 analyzes the problem from a MAP perspective, whereas Section 4 approaches the problem from the consensus equilibrium perspective. Further discussions of open problems are presented in Section 5.

\section{Graph Filters}
\subsection{Constructing a Graph Filter}
Graph filters is a class of image denoisers mostly studied in the field of graph signal processing. The application of graph filters is very broad, e.g., deblurring using graph Laplacian \cite{Kheradmand_Milanfar_2014}, boosting image denoisers \cite{Talebi_Zhu_Milanfar_2013}, image denoising using global similarity \cite{Talebi_Milanfar_2014}, JPEG compression using random walk \cite{Liu_Cheung_Wu_2017}, and  blind deconvolution using reweighted graph total variation \cite{Bai_Cheung_Liu_2018}. In terms of theory, there are studies of the graphs in continuous domain \cite{Pang_Cheung_2017}, sampling of graphs \cite{Anis_Gadde_Ortega_2016}, and filter banks of graphs \cite{Tremblay_Goncalves_Borgnat_2017}, to name a few.

Represented as matrices, graph filters take the form:
\begin{equation}
\vxhat = \calD_{\sigma}(\vy) = \mW \vy,
\label{eq: filtering output}
\end{equation}
where $\vy \in \R^n$ is the noisy image, and $\vxhat \in \R^n$ is the denoised image. The matrix $\mW \in \R^{n \times n}$ is the filter. There are multiple ways of constructing $\mW$. One of the most commonly used approaches is to define a non-local weight kernel $\mK \in \R^{n \times n}$:
\begin{equation}
[\mK]_{ij} = \exp\left\{-\frac{\|\vy_i - \vy_j\|^2}{2h^2}\right\},
\label{eq: K}
\end{equation}
where $\vy_i \in \R^d$ is a $d$-dimensional patch centered at pixel $i$, and $h$ is a parameter characterizing the denoising strength of the filter. Modifications to \eref{eq: K} are common, e.g., by incorporating spatial distance (which gives the non-local means filter), restricting to pixels instead of patches (which gives the bilateral filter), or extending the norm to a weighted norm (which gives the kernel regression filter). For more examples of the kernel matrices, we refer the readers to \cite{Milanfar_2013a}.

The $\mK$ matrix defined in \eref{eq: K} is symmetric but the row sum $\mK\vone$ is not 1. Thus, $\mK$ cannot be used directly as a denoising filter because it amplifies or attenuates signals. To normalize the matrix while preserving the symmetry, the most popular method is by means of the Sinkhorn-Knopp balancing algorithm \cite{Sinkhorn_Knopp_1967, Milanfar_2013b, Chan_Zickler_Lu_2017}, which iteratively normalizes the rows and columns of $\mK$ until convergence. When $\mK$ is symmetrized, the resulting matrix is called a symmetric smoothing filter, given by
\begin{equation}
\mW = \mD^{-1/2}\mK\mD^{-1/2},
\end{equation}
where the diagonal matrix $\mD$ is defined such that $\mW\vone = \vone$ and $\mW^T\vone =  \vone$. Since the columns and rows of $\mW$ sum to 1, $\mW$ is a doubly stochastic matrix.

\subsection{Properties of Graph Filters}
There are a few interesting properties of $\mW$.
\begin{itemize}
\item $\mW$ can be considered as an undirected graph. $[\mW]_{ij}$ is the weight on the edge linking node $i$ and node $j$.
\item $0 \le \lambda_{n}(\mW) \le \ldots \le \lambda_{1}(\mW) = 1$. This follows from the fact that $\mW$ is doubly stochastic, and so the eigenvalues are all non-negative.
\item $\mI - \mW$ is the graph Laplacian, with the zero eigenvalue associated to the vector $\frac{1}{\sqrt{n}}\vone$.
\item The regularization defined through the graph Laplacian $\vx^T(\mI - \mW)\vx$ can be interpreted as
\begin{equation}
\vx^T(\mI - \mW)\vx = \sum_{i=1}^n \sum_{j\not=i}^n [\mK]_{ij}\left( \frac{x_i}{\sqrt{[\mD]_{ii}}} - \frac{x_j}{\sqrt{[\mD]_{jj}}}\right)^2,
\end{equation}
which is measuring the smoothness of the signal $\vx$ with the weights defined by the graph $\mW$.
\item The symmetric matrix $\mW = \mD^{-1/2}\mK\mD^{-1/2}$ has a better denoising performance than its non-symmetric counter part $\mW = \mDtilde^{-1}\mK$ where $\mDtilde = \mbox{diag}(\mK\vone)$. This is attributed to the implicit clustering of the pixels during the Sinkhorn-Knopp balancing algorithm \cite{Chan_Zickler_Lu_2017}, among other reasons such as reduced degree of freedom so that the drop in variance overrides the gain in bias \cite{Milanfar_2013b}.
\end{itemize}

In order to present the theoretical results of this paper, we make an assumption about $\mW$:
\begin{assumption}
\label{assumption: W invertible}
We assume that the matrix $\mW$ is invertible, i.e., there exists a $\mW^{-1}$ such that
\begin{equation}
\mW^{-1}\mW = \mI, \quad\mbox{and}\quad \mW\mW^{-1} = \mI.
\end{equation}
\end{assumption}

Assumption~\ref{assumption: W invertible} is very strong and could be unrealistic for practical $\mW$'s. However, for the derivation of the main results, at least in Section 3, we need $\mW^{-1}$ to simplify some of the technical details. In Section IV when we discuss the performance from the consensus equilibrium perspective, we will lift the assumption.

\vspace{2ex}
\noindent\textbf{Remark: Pseudo-Linearity}. An important remark we should pay attention to is the pseudo-linearity of $\mW$ in the actual PnP ADMM implementation. That is, the matrix $\mW$ is a function of the input noisy image $\vy$, or more precisely $\vxhat = \mW(\vy)\vy$ instead of \eref{eq: filtering output}. However, as mentioned in \cite{Romano2015}, the dependency of $\mW$ on $\vy$ is typically small if $\mW$ is estimated from a \emph{pre-filtered} version of $\vy$, for example applying a baseline denoising algorithm to $\vy$ and constructing a new weight matrix $\mW$ based on the baseline estimate. Practically this is also justified in applications such as \cite{Sreehari_Venkatakrishnan_Wohlberg_2016}, where $\mW$ is updated for the first tens of iterations and is kept fixed for the remaining iterations. Following the same line of arguments, we assume that $\mW$ is pre-defined and is independent of $\vy$. Of particular interest is the \emph{oracle} case where $\mW$ is estimated from the ground truth.

\section{Performance Analysis via MAP}
In this section we present the first set of analytic results. We call it the maximum-a-posteriori (MAP) analysis because we need to explicitly derive the MAP optimization.

\subsection{Existence of $g$}
Given a graph filter $\calD_{\sigma}: \R^n \rightarrow \R^n$ with $\calD_\sigma(\vx) = \mW\vx$, the first question we should ask is whether it is possible to obtain a regularization function $g$ defined in \eref{eq:ADMM2,v}. The rationale is that if we are able to do so, then we can characterize the MAP optimization in \eref{eq:basic optimization}, and hence understand the performance.

Consider the minimization
\begin{equation}
\vxhat = \argmin{\vx \in \R^n} \;\; g(\vx) + \frac{1}{2\sigma^2}\|\vx - \vxtilde\|^2,
\label{eq: rewrite G}
\end{equation}
the solution of which is characterized by
\begin{equation}
\vxhat \in \left( \calI + \sigma^2 \partial g \right)^{-1} \vxtilde,
\label{eq: resolvent}
\end{equation}
where $\partial g$ is the sub-differential of $g$, and the inverse operator is the resolvent of $g$ \cite{Boyd_Parikh_Chu_Peleato_Eckstein_2011}. We assume that $\calD_{\sigma}$ is an invertible matrix $\mW \in \R^{n \times n}$. Therefore, $g$ must be continuously differentiable and so the subgradient at every point has a unique gradient. This replaces $\partial g$ by $\nabla g$, and substituting $\vxhat = \calD_\sigma(\vxtilde)$ into \eref{eq: resolvent} yields
\begin{equation}
\calD_{\sigma}(\vxtilde) = \left( \calI + \sigma^2 \nabla g \right)^{-1}(\vxtilde), \quad \forall \vxtilde.
\label{eq: resolvent 3}
\end{equation}
Define $\mG \in \R^{n \times n}$ as the matrix representation of $\nabla g$ such that $\mG\vx = \nabla g(\vx)$. Since \eref{eq: resolvent 3} holds for all $\vxtilde$, this further implies that the mappings are identical:
\begin{equation}
\mW = \left( \mI + \sigma^2 \mG \right)^{-1}.
\end{equation}
Applying $\mI + \sigma^2 \mG$ to both sides yields $\left( \mI + \sigma^2 \mG \right) \mW = \mI$, and hence
\begin{equation}
\mG = \frac{1}{\sigma^2} \left(\mW^{-1} - \mI\right).
\label{eq: resolvent 4}
\end{equation}
The identity in \eref{eq: resolvent 4} offers a useful way of characterizing the existence of $g$.

\begin{theorem}[Sreehari et al. \cite{Sreehari_Venkatakrishnan_Wohlberg_2016} Theorem A.I, and Reehorst-Schniter \cite{Reehorst_Schniter_2018} Theorem 1.]
\label{thm: existence}
The mapping $g$ exists if and only if $\mG$ is symmetric.
\end{theorem}
\begin{proof}
The idea is to treat $g$ as a vector field which must be conservative. Therefore, the gradient $\nabla g$ must be symmetric using a classic result in vector calculus \cite{Kantorovitz_2016} (in particular Theorem 4.3.8 and 4.3.10).
\end{proof}

For graph filters, the symmetry of $\mG$ is guaranteed because $\mW$ is symmetric. In addition, we also have that $\mW$ is non-expansive, as the eigenvalues of $\mW$ are bounded above by 1. Therefore, if we use a graph filter and if $\mW^{-1}$ exists, then according to \cite{Sreehari_Venkatakrishnan_Wohlberg_2016} the corresponding minimization is a proximal map, and $g$ is convex, closed and proper. However, an arbitrary denoiser $\calD_{\sigma}$ does not necessarily have a symmetric Jacobian, and hence $g$ is not guaranteed to exist.

\subsection{Deriving $g$}
We can now state the first key result of this paper.
\begin{theorem}
\label{prop: s}
Let $\mW \in \R^{n \times n}$ be a graph filter with full rank so that the inverse of $\mW$ exists. If $g: \R^n \rightarrow \R$ is defined by
\begin{equation}
g(\vx) = \frac{1}{2} \vx^T \mG\vx = \frac{1}{2 \sigma^2} \vx^T (\mW^{-1} - \mI)\vx,
\label{eq: g expression}
\end{equation}
then the solution of the minimization of \eref{eq:ADMM2,v} is
\begin{equation}
\vxhat  \;\bydef\;  \argmin{\vx \in \R^n} \;\; g(\vx) + \frac{1}{2\sigma^2}\|\vx - \vxtilde\|^2 = \mW\vxtilde.
\label{eq: nlm v sub}
\end{equation}
\end{theorem}
\begin{proof}
The simplest proof is to substitute \eref{eq: g expression} into the optimization in \eref{eq: nlm v sub}. An alternative proof is to multiply both sides of \eref{eq: resolvent 4} by $\vx$ to obtain
\begin{equation}
\nabla g(\vx) = \mG\vx = \frac{1}{\sigma^2}(\mW^{-1} - \mI)\vx.
\end{equation}
Then integrating both sides with respect to $\vx$ yields $g(\vx) = \frac{1}{2} \vx^T \mG\vx$, which completes the proof.
\end{proof}

Can we replace $\mW^{-1}$ by its pseudo-inverse $\mW^+$ if $\mW$ is not invertible? The answer is no. If we replace $\mW^{-1}$ by $\mW^+$ in \eref{eq: g expression}, the optimization \eref{eq: rewrite G} becomes
\begin{equation}
\vxhat = \argmin{\vx} \;\; \frac{1}{2 \sigma^2} \vx^T (\mW^{+} - \mI)\vx + \frac{1}{2\sigma^2}\|\vx - \vxtilde\|^2.
\label{eq: pinv W optimization}
\end{equation}
The first order optimality condition of \eref{eq: pinv W optimization} is given by $\mW^+ \vxhat = \vxtilde$, which does not imply $\vxhat = \mW\vxtilde$ because $\mW^+\mW \not= \mI$ and $\mW\mW^+ \not= \mI$ unless $\mW$ has full rank. More specifically, $\vxhat = \mW\vxtilde$ requires $\vxhat$ to live in the range space of $\mW$. Since $\mW^+$ contains zero eigenvalues, there are multiple $\vxhat$ which can be mapped to the same $\vxtilde$ through $\mW^+ \vxhat = \vxtilde$. This is problematic, because the optimizer of \eref{eq: pinv W optimization} is no longer $\vxhat = \mW\vxtilde$, and hence it violates our assumption that the denoiser is a graph filter.

\subsection{The Role of $\rho$}
Knowing the explicit expression of $g$, we can now analyze its performance. But before we proceed to the discussion, we should clarify the role of the parameter $\rho$. By substituting \eref{eq: g expression} into \eref{eq:basic optimization 2}, and recalling $\sigma \bydef  \sqrt{\lambda/\rho}$, we can show that the original problem which PnP ADMM solves is
\begin{equation}
\minimize{\vx} \;\; f(\vx) + \frac{\rho}{2}\vx^T(\mW^{-1}-\mI)\vx.
\label{eq: solution PnP}
\end{equation}
There is an alternative optimization in the graph signal processing literature known as the graph Laplacian, which we will study extensively in this paper:
\begin{equation}
\minimize{\vx} \;\; f(\vx) + \frac{\lambda}{2}\vx^T(\mI - \mW)\vx.
\label{eq: solution Lps}
\end{equation}

In \eref{eq: solution PnP}, the original regularization parameter $\lambda$ is eliminated because $g$ contains its reciprocal. Thus, $\lambda$ is replaced by the ADMM internal parameter $\rho$. This internal parameter $\rho$ has no influence to the final solution in the classical ADMM because ADMM converges for all $\rho$ when $f$ and $g$ are closed, proper and convex. However, in PnP ADMM, the optimization cost is characterized by $\rho$. Thus, the solution would change as $\rho$ changes, a phenomenon reported in earlier papers, e.g., \cite{Chan_Wang_Elgendy_2017}.

We can also ask: For the best $\rho$ and $\lambda$, how would the solutions of \eref{eq: solution PnP} and \eref{eq: solution Lps} behave? \fref{fig: rho} shows an experimental result where $f(\vx) = \frac{1}{2}\|\vx - \vy\|^2$. We vary $\rho$ (and equivalently $\lambda$) to evaluate the mean squared error (MSE) for a fixed oracle $\mW$. As shown in \fref{fig: rho}, PnP ADMM and graph Laplacian have different operating regimes for $\rho$ and $\lambda$. PnP ADMM prefers smaller $\rho$, whereas graph Laplacian prefers larger $\lambda$. We will explain this observation in Section 4.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\linewidth]{rho}
\caption{Mean squared error of the solutions of \eref{eq: solution PnP} and \eref{eq: solution Lps} as $\rho$ (or $\lambda$) increases. PnP has a lower MSE than graph Laplacian for the best $\rho$ and $\lambda$.}
\label{fig: rho}
\end{figure}

\subsection{Comparison with Graph Laplacian}
We can now compare the performance of PnP ADMM and the conventional graph Laplacian regularization. In order to do so, we consider a denoising problem such that
$$
f(\vx) = \frac{1}{2}\|\vx - \vy\|^2,
$$
where $\vy = \vx + \veta$, where $\veta \sim \calN(\vzero, \sigma_{\eta}^2\mI)$. That is, $\vy$ is a noisy version of $\vx$. Such choice of $f(\vx)$ is also justified by the $\vv$-subproblem \eref{eq:ADMM2,v}, which is also a denoising problem.

We consider the following two optimizations:
\begin{align}
\vxhat_{L} &= \argmin{\vx} \;\; \frac{1}{2}\|\vx - \vy\|^2 + \frac{\alpha}{2} \vx^T (\mI - \mW) \vx, \label{eq: Lps analysis 1}\\
\vxhat_{P} &= \argmin{\vx} \;\; \frac{1}{2}\|\vx - \vy\|^2 + \frac{\alpha}{2} \vx^T (\mW^{-1} - \mI) \vx.  \label{eq: PnP analysis 1}
\end{align}
Here, the first optimization \eref{eq: Lps analysis 1} uses graph Laplacian as the regularization, and the second optimization \eref{eq: PnP analysis 1} uses the PnP regularization. We use a parameter $\alpha$ instead of $\rho$ or $\lambda$ to keep the notation simple. The solutions of the two optimizations are
\begin{align}
\vxhat_{L} &= \left[(1+\alpha)\mI - \alpha\mW \right]^{-1} \vy \bydef \mA_L^{-1}\vy, \label{eq: Lps analysis 2}\\
\vxhat_{P} &= \left[(1-\alpha)\mI + \alpha\mW^{-1} \right]^{-1} \vy \bydef \mA_P^{-1}\vy, \label{eq: PnP analysis 2}
\end{align}
where we defined matrices $\mA_L \bydef (1+\alpha)\mI - \alpha\mW$ and $\mA_P \bydef (1-\alpha)\mI + \alpha\mW $ to simplify notations.

Our goal is to analyze the mean squared error (MSE) of these two estimates. To this end, we first eigen-decompose $\mW$ as $\mW = \mU\mS\mU^T$, where $\mS = \mbox{diag}(s_1,\ldots,s_n)$ is the eigenvalue matrix. Define $\vb \bydef \mU^T\vx$ as the projection of $\vx$ by columns of $\mU$. Then, the MSE of $\vxhat_{L}$ and $\vxhat_{P}$ can be found as follows.

\begin{proposition}
Define $\vxhat_{L}$ and $\vxhat_{P}$ via \eref{eq: Lps analysis 2} and \eref{eq: PnP analysis 2}. Let $\mW = \mU\mS\mU^T$, where $\mU$ is the eigenvector matrix, and $\mS = \mbox{diag}(s_1,\ldots,s_n)$ is the eigenvalue matrix. Define $\vb \bydef \mU^T\vx$. Then, the MSE of $\vxhat_{L}$ and $\vxhat_{P}$ are respectively
\begin{align}
\mbox{MSE}_{L} &\bydef \E\|\vxhat_{L}-\vx\|^2 = \sum_{i=1}^n \frac{\alpha^2(1-s_i)^2 b_i^2 + \sigma_{\eta}^2}{(1+\alpha - \alpha s_i)^2},\\
\mbox{MSE}_{P} &\bydef \E\|\vxhat_{P}-\vx\|^2 = \sum_{i=1}^n \frac{\alpha^2(1-s_i)^2 b_i^2 + \sigma_{\eta}^2 s_i^2}{(\alpha + (1-\alpha) s_i)^2}.
\end{align}
\end{proposition}

\begin{proof}
The MSE of $\vxhat_{L}$ is
\begin{align}
\mbox{MSE}_{L}
&= \E\|\mA_L^{-1}(\vx+\veta)-\vx\|^2\\
&= \|(\mA_L^{-1}-\mI)\vx\|^2 + \sigma_\eta^2 \mbox{Tr} \left[ \left( \mA_{L}^{T} \mA_{L} \right)^{-1} \right] \notag \\
&= \left\|\mU\left(1 - \frac{1}{1+\alpha - \alpha s_i} \right) \mU^T \vx \right\|^2 + \cdots  \notag \\
&\quad\quad\quad + \sigma_\eta^2 \left\| \mU \left(\frac{1}{1+\alpha - \alpha s_i} \right) \mU^T\right\|_F^2 \notag \\
&= \sum_{i=1}^n \left(1 - \frac{1}{1+\alpha - \alpha s_i} \right)^2 b_i^2 + \sigma_{\eta}^2 \left(\frac{1}{1+\alpha - \alpha s_i} \right)^2 \notag \\
&= \sum_{i=1}^n \frac{\alpha^2(1-s_i)^2 b_i^2 + \sigma_{\eta}^2}{(1+\alpha - \alpha s_i)^2}.
\end{align}
Similarly, the MSE of $\vxhat_{P}$ is
\begin{align}
\mbox{MSE}_{P} &= \|(\mA_P^{-1}-\mI)\vx\|^2 + \sigma_\eta^2 \mbox{Tr} \left[  \left( \mA_{P}^{T} \mA_{P} \right)^{-1} \right] \notag  \\
&= \sum_{i=1}^n \frac{\alpha^2(1-s_i)^2 b_i^2 + \sigma_{\eta}^2 s_i^2}{(\alpha + (1-\alpha) s_i)^2}.
\end{align}
\end{proof}

Now, we can consider the $i$-th term of the sums:
\begin{align}
\mbox{MSE}_{L}^i &\bydef \frac{\alpha^2(1-s_i)^2 b_i^2 + \sigma_{\eta}^2}{(1+\alpha - \alpha s_i)^2},\\
\mbox{MSE}_{P}^i &\bydef \frac{\alpha^2(1-s_i)^2 b_i^2 + \sigma_{\eta}^2 s_i^2}{(\alpha + (1-\alpha) s_i)^2}.
\end{align}
It holds that for $s_i \rightarrow 0$,
\begin{align}
\lim_{s_i \rightarrow 0} \mbox{MSE}_{L}^i &= \frac{\alpha^2 b_i^2 + \sigma_\eta^2}{(1+\alpha)^2}, \label{eq:MSE L}\\
\lim_{s_i \rightarrow 0} \mbox{MSE}_{P}^i &= b_i^2. \label{eq:MSE P}
\end{align}
The limits in \eref{eq:MSE L} and \eref{eq:MSE P} provide important information about the performance. We make the following observations.

\vspace{2ex}
\noindent\textbf{1. Energy Concentration by $\mW$}. The graph filter $\mW$, when applied to a signal $\vx$, projects the signal onto the space spanned by its columns. The effectiveness of this projection is determined by the quality of $\mW$. In the oracle setting, $\mW$ ensures that most of the energy is captured by the first few eigenvectors. This is illustrated in the \fref{fig: projection}, where we consider a 1D signal $\vx$ and its noisy version $\vy$. The weight matrix $\mW$ is constructed from $\vx$ using the non-local means kernel defined in \eref{eq: K}. By projecting $\vx$ and $\vy$ using the eigenvector matrix $\mU$, we see that $\mU^T\vx$ has significantly fewer non-zeros than $\mU^T\vy$ for large eigen-indices. Since $\mU^T\vx$ is the vector $\vb$, we have the observation:

\vspace{2ex}
\noindent \emph{Observation 1}: For large $i$, $b_i$ is typically very small.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\linewidth]{projection}
\caption{Projecting a noisy signal $\vy$ and the ground truth signal $\vx$ by the principle eigenvectors. The two curves show $\mU^T\vy$ and $\mU^T\vx$, where $\mU$ is the eigenvector matrix of $\mW$.}
\label{fig: projection}
\end{figure}

\vspace{2ex}
\noindent\textbf{2. Eigenvalues at High Frequency}. Based on Observation 1, an immediate conclusion from \eref{eq:MSE L} and \eref{eq:MSE P} is that when $b_i \rightarrow 0$, the two limits are respectively
\begin{align*}
\lim_{b_i, s_i \rightarrow 0} \mbox{MSE}_{L}^i = \frac{\sigma_\eta^2}{(1+\alpha)^2}, \quad\mbox{and}\quad
\lim_{b_i, s_i \rightarrow 0} \mbox{MSE}_{P}^i = 0.
\end{align*}
In other words, the PnP regularization ensures that the MSE goes to zero for high frequency (i.e., large eigen-index $i$) whereas the graph Laplacian has residue MSE for high frequency. These can also be seen from the eigenvalues of $\mA_L^{-1}$ and $\mA_P^{-1}$:
\begin{align*}
\lambda_i(\mA_L^{-1}) = \frac{1}{1+\alpha-\alpha s_i}, \quad\mbox{and}\quad \lambda_i(\mA_P^{-1}) = \frac{s_i}{(1-\alpha)s_i+\alpha},
\end{align*}
where $\lambda_i(\cdot)$ denotes the $i$-th eigenvalue of the matrix. If we plot $\lambda_i(\mA_L^{-1})$ and $\lambda_i(\mA_P^{-1})$ as functions of the eigen-index $i$, we will obtain a result shown in \fref{fig: eigen value}. In this experiment we set $\sigma_\eta = 0.05$ and $\alpha = 0.2$. As $s_i \rightarrow 0$, we have
\begin{align*}
\lim_{s_i \rightarrow 0} \lambda_i(\mA_L^{-1}) = \frac{1}{1+\alpha}, \quad\mbox{and}\quad \lim_{s_i \rightarrow 0} \lambda_i(\mA_P^{-1}) = 0,
\end{align*}
which explain why $\mbox{MSE}_{L}^i$ reaches a steady state whereas $\mbox{MSE}_{P}^i$ reaches zero as $s_i \rightarrow 0$.

\vspace{2ex}
\noindent \emph{Observation 2}: As $i$ increases so that $s_i \rightarrow 0$ and $b_i \rightarrow 0$, the eigenvalue $\lambda_i(\mA_L^{-1})$ approaches a non-zero constant, whereas $\lambda_i(\mA_P^{-1})$ approaches 0. The same observation applies to MSEs.

\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{eigen_value}
\caption{Comparison of $\lambda_i(\mA_L^{-1})$ and $\lambda_i(\mA_P^{-1})$.}
\label{fig: eigen value}
\end{figure}

\vspace{2ex}
\noindent\textbf{3. Trade offs in Bias and Variance}. MSE is the sum of bias and variance. So we can decompose the MSE and inspect
\begin{align*}
\mbox{Bias}_P^i &= \frac{\alpha^2(1-s_i)^2 b_i^2}{(\alpha + (1-\alpha) s_i)^2}, \quad\mbox{and}\quad
\mbox{Var}_P^i  = \frac{\sigma_{\eta}^2 s_i^2}{(\alpha + (1-\alpha) s_i)^2}\\
\mbox{Bias}_L^i &= \frac{\alpha^2(1-s_i)^2 b_i^2}{(1+\alpha - \alpha s_i)^2}, \quad\mbox{and}\quad
\mbox{Var}_L^i  = \frac{\sigma_{\eta}^2}{(1+\alpha - \alpha s_i)^2}.
\end{align*}
Since by construction $s_i$ is decreasing as $i$ grows, we can show that $\mbox{Var}_L^i \ge \mbox{Var}_P^i$.
\begin{proposition}
For any $0 \le s_i \le 1$, and $\alpha \le 1$, it holds that $\mbox{Var}_L^i \ge \mbox{Var}_P^i$ for all $i$.
\end{proposition}
\begin{proof}
To prove this result we just need to show the function
\begin{align*}
f(s_i)
&= s_i^2(1+\alpha-\alpha s_i)^2 - (\alpha + s_i - \alpha s_i)^2 \\
&= \alpha (s_i - 1)^2 (\alpha s_i^2 - 2 s_i - \alpha),
\end{align*}
is non-positive. In this expression, the first two terms are non-negative. The last term $\alpha s_i^2 - 2 s_i - \alpha$ is non-positive because $\alpha s_i^2 - 2 s_i - \alpha \le s_i - 2s_i - \alpha \le 0$, where the first inequality holds since $\alpha s_i \le 1$.
\end{proof}
The conclusion of this proposition is that the variance of PnP is always lower than the variance of graph Laplacian. However, the bias of PnP is larger than the bias of graph Laplacian, as we show in \fref{fig: bias var}. More specifically, we can prove the following result:
\begin{proposition}
For any $0 \le s_i \le 1$, and $\alpha < 1/2$, it holds that $\mbox{Bias}_L^i \le \mbox{Bias}_P^i$ for all $i$.
\end{proposition}
\begin{proof}
Since the numerators of $\mbox{Bias}_L^i$ and $\mbox{Bias}_P^i$ are identical, we only need to compare the denominator. Thus, we have
\begin{align*}
f(s_i) &= (\alpha + (1-\alpha) s_i)^2 - (1+\alpha - \alpha s_i)^2 \\
&= (1-s_i)[(2\alpha + 1) - (2\alpha-1)s_i] \le 0
\end{align*}
if $s_i \ge \frac{2\alpha+1}{2\alpha-1}$, which is true when $\alpha < 1/2$.
\end{proof}

The reason why the overall MSE is lower for PnP is attributed to the observation below.

\vspace{2ex}
\noindent \emph{Observation 3}: The drop in variance of PnP is more significant than the gain in bias, and hence PnP's overall MSE is lower than that of graph Laplacian.

\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{bias_var}
\caption{Bias and variance of PnP and graph Laplacian.}
\label{fig: bias var}
\end{figure}

\vspace{2ex}
\noindent\textbf{4. When is Graph Laplacian Better?} The analysis above suggests conditions under which PnP will perform worse than graph Laplacian. This happens when the vector $\vb = \mU^T\vx$ is large, because then the denominator of $\mbox{Bias}_P^i$ will amplify more than that of $\mbox{Bias}_L^i$. This also suggests that PnP is more sensitive to the correctness of $\mW$ than graph Laplacian. That is, when constructing $\mW$ from a pre-filtered signal, PnP demands higher quality of the pre-filtered signal.

To confirm this claim, we conduct an experiment by constructing a pre-filtered signal via $\vxhat_{\mathrm{pre}} = \vx + \vepsilon$, where $\vx \in \R^n$ is the true signal, and $\vepsilon \sim \calN(\vzero, \sigma_\epsilon^2\mI)$. That is, the pre-filtered signal is a noisy version of the clean signal, with perturbation specified by the noise variance $\sigma_\epsilon^2$. We then construct a graph filter $\mW$ using $\vxhat_{\mathrm{pre}} $. Technically speaking, $\mW$ is now $\mW(\vxhat_{\mathrm{pre}})$. Our goal is to see how $\mbox{MSE}_{L}$ and $\mbox{MSE}_{P}$ would change when $\sigma_\epsilon$ increases. \footnote{An alternative experiment would be to pick a specific baseline denoising algorithm and adjust its internal parameters. However, different baseline denoising algorithms have different characteristics and this could cause bias in the analysis. Thus, we choose to design a pre-filtered signal by adding iid noise.}

The result of this experiment is shown in \fref{fig: prefilter}, where we observe that PnP has a considerably lower MSE than graph Laplacian when $\sigma_\epsilon$ is small. But as $\sigma_\epsilon$ increases, $\mbox{MSE}_{P}$ quickly increases and eventually goes above $\mbox{MSE}_{L}$. This explains the relatively weak robustness of PnP compared to graph Laplacian.

\begin{figure}[h]
\centering
\includegraphics[width=1\linewidth]{prefilter}
\caption{The influence of pre-filtered signal quality to the final MSE. In this plot, $\sigma_\epsilon$ determines the amount of perturbation are there in the pre-filtered signal. As $\sigma_\epsilon$ increases, PnP shows a more sensitive behavior than graph Laplacian.}
\label{fig: prefilter}
\end{figure}

However, we should emphasize that such conclusion is only valid when $\mW$ is \emph{fixed}. In practical PnP ADMM, the weight matrix $\mW$ is updated in every ADMM iteration. If $\sigma_\epsilon$ decreases as PnP iterates, then $\mW$ could become a good filter.

\vspace{2ex}
\noindent \emph{Observation 4}: The performance of PnP could be worse than graph Laplacian when $\mW$ is estimated wrongly.

\vspace{2ex}
\noindent\textbf{5. Geometry: Truncating Eigenvectors of $\vx$}. If $\mW$ is oracle, is there any simple geometric interpretation, beside the bias-variance analysis, that can explain the performance?

By expressing the eigen-decomposition of the $\mW$ as
\begin{equation}
\mW =
\begin{bmatrix}
\mU_1 \\ \mU_2
\end{bmatrix}
\begin{bmatrix}
\mS_1 & \\
& \mS_2
\end{bmatrix}
\begin{bmatrix}
\mU_1^T & \mU_2^T
\end{bmatrix},
\label{eq: decompose W}
\end{equation}
the PnP ADMM regularization in \eref{eq: PnP analysis 1} can be written as
\begin{align}
&\vx^T (\mW^{-1} - \mI) \vx\label{eq: decompose regularization PnP} \\
&\quad\quad = \vx^T \left( \mU_1 (\mS_1^{-1} - \mI) \mU_1^T + \mU_2 (\mS_2^{-1} - \mI) \mU_2^T \right) \vx.  \notag
\end{align}
If the eigenvalues $\mS_2 = \vzero$, then the second term in \eref{eq: decompose regularization PnP} will dominate as $\mS_2^{-1}$ is infinity. In this case, we must have $\mU_2^T\vx = \vzero$ in order for the regularization to stay finite. Therefore, the overall optimization is now constrained:
\begin{align}
\vxhat_{P} &= \argmin{\vx} \;\; \frac{1}{2}\|\vx - \vy\|^2 + \frac{\alpha}{2} \vx^T \mU_1 (\mS_1^{-1} - \mI)\mU_1^T \vx \notag \\
&\quad\;\; \mbox{subject to} \;\; \mU_2^T \vx = \vzero.
\label{eq: decompose PnP optimization}
\end{align}

The constraint provides a simple geometric interpretation. While the graph Laplacian regularization searches in the entire column space, PnP ADMM only searches for a subspace spanned by eigenvectors of non-zero eigenvalues. As we see in Observation 1, this provides a very strong denoising power when $\mW$ is the oracle. Referring to the decomposition in \eref{eq: decompose W}, the signal is spanned only by $\mU_1$ whereas the noise is spanned by both $\mU_1$ and $\mU_2$. By restricting $\mU_2^T\vx = \vzero$, we ensure that any noise living in $\mU_2$ will be eliminated. See \fref{fig: eigen diagram} for illustration. As a result, the denoising is significantly more effective than graph Laplacian.

\begin{figure}[h]
\centering
\includegraphics[width=0.7\linewidth]{eigen_diagram}
\caption{PnP ADMM decomposes the optimization variable $\vx$ into two subspaces: $\mU_1$, the subspace encoding both signal and noise, and $\mU_2$, the subspace encoding only the noise. By restricting $\vx$ such that $\mU_2^T\vx = \vzero$, we force PnP ADMM to eliminate any noise living in $\mU_2$. }
\label{fig: eigen diagram}
\end{figure}

\noindent \emph{Observation 5}: Geometrically, PnP ADMM restricts the search space so that the pure noise components are forced to zero.

\subsection{Questions Not Answered}
To summarize, the fundamental reason why PnP works better is that the graph filter $\mW$ has effectively removed the unwanted noise from $\vy$ before proceeding to additional processing. Thus, a stronger denoiser leads to better PnP ADMM algorithm.

One major limitation of the above MAP-based analysis is that it tries to match $\calD_{\sigma}$ with a proximal map associated with a regularization function $g$. This is a particularly strong requirement, because many denoisers are not symmetric and hence there is no corresponding proximal map. It also requires $\mW^{-1}$ to exist, which is difficult to have in practice, e.g., a uniform average operator does not have an inverse.

The question we ask is that if we are willing to give up the proximal map requirement, would it be possible to obtain another interpretation that can explain the performance for denoisers that do not have $\mW^{-1}$? This leads to the next section of equilibrium-based analysis.

\section{Performance Analysis via Equilibrium}
The MAP optimization in \eref{eq:basic optimization} is not the only way to obtain $\vxhat$. In \cite{Buzzard_Chan_Bouman_2017}, Buzzard et al. generalize PnP ADMM by a concept called consensus equilibrium, which allows us to characterize the solution of PnP ADMM by analyzing the equilibrium state of a dynamical system. In this section we discuss the PnP ADMM performance from such equilibrium perspective.

\subsection{Equilibrium}
The starting point of the equilibrium analysis is the PnP ADMM algorithm listed in \eref{eq:ADMM2,x}-\eref{eq:ADMM2,u}. By defining a pair of functions $F$ and $G$ using \eref{eq:ADMM2,x} and \eref{eq:ADMM2,v} respectively, we can show that the PnP ADMM solution $\vxhat$ satisfies
\begin{align}
\vxhat &= F(\vxhat + \vuhat), \label{CE, eq1}\\
\vxhat &= G(\vxhat - \vuhat). \label{CE, eq2}
\end{align}
A solution $(\vxhat, \vuhat)$ satisfying \eref{CE, eq1} and \eref{CE, eq2} is called a consensus equilibrium solution. To make this more specific to the readers, consider the following denoising example.
\begin{proposition}
\label{prop: CE example}
Consider two operators $F$ and $G$:
\begin{align*}
F(\vx) &= \argmin{\vz\in \R^n} \;\;  \frac{1}{2}\|\vz-\vy\|^2 + \frac{\alpha}{2} \|\vz - \vx\|^2,\\
G(\vx) &= \mW\vx
\end{align*}
The consensus equilibrium solution $(\vxhat,\vuhat)$ satisfies
\begin{equation}
\left((1-\alpha)\mI + \alpha \mW^{-1} \right)\vxhat = \vy
\label{eq: CE example solution}
\end{equation}
and $\vuhat = (\mI - \mW^{-1})\vxhat$.
\end{proposition}

\begin{proof}
Since $G(\vx) = \mW\vx$, it holds that $G(\vxhat - \vuhat) = \mW(\vxhat - \vuhat)$, which implies that
\begin{align*}
\mW\vuhat = \mW\vxhat - \vxhat \quad \Rightarrow \vuhat = (\mI - \mW^{-1})\vxhat.
\end{align*}
The optimality of $F$ gives
\begin{equation*}
F(\vxhat + \vuhat) = (1+\alpha)^{-1}(\vy + \alpha(\vxhat + \vuhat))
\end{equation*}
Substituting $\vuhat = (\mI - \mW^{-1})\vxhat$ into $F(\vxhat + \vuhat) = \vxhat$ gives
\begin{align*}
(1+\alpha)^{-1} \left(\vy + \alpha(\vxhat + (\mI - \mW^{-1})\vxhat)\right) = \vxhat.
\end{align*}
Rearranging the terms yields the desired result.
\end{proof}

The interesting fact about Proposition~\ref{prop: CE example} is that, given $\vxhat$ we can easily check whether it is a consensus equilibrium point without going through the PnP ADMM algorithm. In other words, we can ``create'' an alternative optimization such that its solution is an equilibrium point, and we can interpret the meaning of the optimization to gain insights. To illustrate this point, let us continue with the example in Proposition~\ref{prop: CE example}.

\begin{proposition}
\label{prop: equilibrium}
Define the functions $\varphi$, $\psi_1$, $\psi_2$ and $\psi_3$:
\begin{align}
\varphi(\vx) &= \frac{1}{2}\|\vx - \vy\|^2 + \frac{\alpha}{2} \vx^T (\mW^{-1} - \mI) \vx, \label{eq: equilibrium phi}\\
\psi_1(\vx) &= \frac{1}{2}\|\vx - \mW\vy\|^2 - \frac{(1-\alpha)}{2}\vx^T (\mI - \mW)\vx, \label{eq: equilibrium psi}\\
\psi_2(\vx) &= \frac{1}{2}(\vx-\vy)^T\mW(\vx-\vy) + \frac{\alpha}{2} \vx^T(\mI - \mW)\vx,\\
\psi_3(\vx) &= \frac{1}{2}\left\| \vx - \left[(1-\alpha)\mI + \alpha \mW^{-1} \right]^{-1}\vy\right\|^2.
\end{align}
The minimizers of these four functions all satisfy the consensus equilibrium condition \eref{eq: CE example solution}.
\end{proposition}

\begin{proof}
Recalling \eref{eq: PnP analysis 2}, the minimizer of $\varphi$ is
\begin{align}
\vxhat = \left[(1-\alpha)\mI + \alpha \mW^{-1} \right]^{-1}\vy,
\label{eq: equilibrium 3}
\end{align}
which satisfies \eref{eq: CE example solution}. Rearranging the terms we obtain
\begin{align*}
\left[(1-\alpha)\mI + \alpha \mW^{-1} \right]^{-1}\vy
&= \left[(1-\alpha)\mW + \alpha \mI  \right]^{-1} \mW \vy \\
&= \left[\mI - (1-\alpha)(\mI - \mW) \right]^{-1} \mW \vy,
\end{align*}
which is the minimizer of $\psi_1$. The minimizers of $\psi_2$ and $\psi_3$ can be derived in a similar manner.
\end{proof}

Proposition~\ref{prop: equilibrium} shows four different optimization formulations that give the same consensus equilibrium solution. The first optimization in $\varphi$ is the MAP formulation we saw in Section 3. The other optimizations of $\psi_1$, $\psi_2$ and $\psi_3$ have some meaningful interpretations.

\vspace{2ex}
\noindent $\bullet$ \textbf{Interpretation of $\psi_1(\vx)$}. The optimization of $\psi_1$ is informative in the sense that it completely eliminates the inverse matrix $\mW^{-1}$, as opposed to the MAP in \eref{eq: PnP analysis 1}. The objective function $\frac{1}{2}\|\vx - \mW\vy\|^2$ says that instead of trying to minimize the residue between the optimization variable $\vx$ and the noisy observation $\vy$, we minimize $\vx$ with the filtered version $\mW\vy$. For oracle $\mW$'s where pure noise components are isolated at high frequencies (See Observation 5), $\mW\vy$ ensures that the pure noise components are eliminated. Therefore, when minimizing $\frac{1}{2}\|\vx - \mW\vy\|^2$, we are guaranteed to only search for the signal component.

The regularization term $-\frac{(1-\alpha)}{2}\vx^T (\mI - \mW)\vx$ is a concave function for any $0 \le \alpha \le 1$, because $\mI-\mW$ is positive semi-definite. The concave function \emph{injects} high frequency components that are lost by $\mW\vy$. In the extreme cases when $\alpha = 1$, we have
\begin{equation*}
\psi_1(\vx) = \frac{1}{2}\|\vx - \mW\vy\|^2,
\end{equation*}
and so the minimizer is just $\vxhat = \mW\vy$. When $\alpha = 0$, we have
\begin{align*}
\psi_1(\vx)
&= \frac{1}{2}\|\vx - \mW\vy\|^2 - \frac{1}{2}\vx^T (\mI - \mW)\vx \\
&= \frac{1}{2}\|\vx\|^2 - \vx^T\mW\vy + \frac{1}{2}\|\mW\vy\|^2 - \frac{1}{2}\|\vx\|^2 + \frac{1}{2}\vx^T\mW\vx\\
&= \frac{1}{2}(\vx - \vy)^T\mW(\vx - \vy) + \frac{1}{2}\vy^T\mW(\mW-\mI)\vy,
\end{align*}
and so the minimizer is $\vxhat = \vy$. In order to control the amount of high frequency injection, a small $\alpha$ is more preferred than a large $\alpha$. This explains why in \fref{fig: rho} a smaller $\alpha$ works better for PnP ADMM.

\vspace{2ex}
\noindent$\bullet$\textbf{Interpretation of $\psi_2(\vx)$}. The function $\psi_2$ offers an interpretation similar to $\psi_1$, e.g., the objective $(\vx-\vy)^T\mW(\vx-\vy)$ has no penalty over the subspace spanned by the eigenvectors of zero-eigenvalues, and so the regularization $\vx^T(\mI - \mW)\vx$ will force those dimensions to zero.

$\psi_2$ is also a special case of Kheradmand and Milanfar \cite{Kheradmand_Milanfar_2014}, who suggested
$$\psi(\vx) = \frac{1}{2}(\vx-\vy)^T\left(\mI + \beta(\mI - \mW)\right)(\vx-\vy) + \frac{\alpha}{2} \vx^T(\mI - \mW)\vx$$
for any $\beta \ge -1$. If we let $\beta = -1$, then we obtain $\psi_2$.

\vspace{2ex}
\noindent$\bullet$\textbf{Interpretation of $\psi_3(\vx)$}.
The last example $\psi_3$, while still giving the consensus equilibrium solution, is more or less useless (in terms of interpretation) because it is just a re-statement of the solution.

In summary, Proposition~\ref{prop: equilibrium} shows that the solution of the PnP ADMM is completely characterized by the consensus equilibrium. Solving the MAP optimization is just one of the many possible ways of achieving the equilibrium. If we depart from the MAP formulation, i.e., if we are willing to change the forward modeling $f(\vx)$ by integrating the denoiser into the forward model, then there is a wide range of optimization formulations that can offer meaningful interpretations while satisfying the equilibrium condition.

\subsection{General Linear Inverse Problems}
The equilibrium based approach also allows us to extend the analysis to general linear inverse problems which could be difficult to study by writing out the MSE in Section 3. The general linear inverse problem using the PnP ADMM with a graph filter requires two operators:
\begin{align*}
F(\vx) &= \argmin{\vz\in \R^n} \;\;  \frac{1}{2}\|\mA\vz-\vy\|^2 + \frac{\alpha}{2} \|\vz - \vx\|^2,\\
G(\vx) &= \mW\vx.
\end{align*}
of which the equilibrium solution $\vxhat$ satisfies
\begin{equation}
\left(\mA^T\mA + \alpha(\mW^{-1}-\mI)\right)\vxhat = \mA^T\vy.
\label{eq: CE linear condition}
\end{equation}
The condition in \eref{eq: CE linear condition} can be written as
\begin{align*}
\left(\mW^{\frac{1}{2}}\mA^T\mA\mW^{\frac{1}{2}} + \alpha(\mI-\mW)\right)\mW^{-\frac{1}{2}}\vxhat = \mW^{\frac{1}{2}}\mA^T\vy,
\end{align*}
which is also the equivalent to
\begin{equation}
\vxhat = \mW^{\frac{1}{2}} \cdot \left\{ \argmin{\vz \in \R^n } \;\; \frac{1}{2} \|\mA \mW^{\frac{1}{2}} \vz - \vy\|^2 + \frac{\alpha}{2} \vz^T(\mI - \mW)\vz \right\}.
\label{eq: CE linear}
\end{equation}

While \eref{eq: CE linear} may look unnatural at a first glance, in the sparse signal processing literature this is reminiscent to the MAP-synthesis minimization discussed by Elad \cite{Elad_Milanfar_Rubinstein_2007}.

How do we interpret \eref{eq: CE linear}? The inner minimization is the standard graph Laplacian minimization with $\mA\vz$ changed to $\mA \mW^{\frac{1}{2}} \vz$ in the objective. If $\mW$ has zero eigenvalues, then $\mW^{\frac{1}{2}}\vz$ ensures that the columns along those directions are untouched. Then by the virtue of $\vz^T(\mI - \mW)\vz$, the minimizer must be zero along those directions because $\mI-\mW$ is positive semi-definite. The outer step by multiplying the solution with $\mW^{\frac{1}{2}}$ ensures that the noise from the inner minimization is eliminated before creating the final solution.

\subsection{Performance of Multiple Priors}
Beyond linear inverse problem we can extend the analysis to PnP ADMM involving multiple priors. In this case, the MAP-based optimization is
\begin{equation}
\vxhat = \argmin{\vx \in \R^n} \;\; \varphi(\vx) \bydef \mu_0 f_0(\vx) + \sum_{i=1}^k \mu_i f_i(\vx),
\label{eq: CE main}
\end{equation}
where $f_0: \R^n \rightarrow \R $ is objective function associated with the forward model, and $\{f_i\}_{i=1}^k$ are the regularization functions associated with the prior models (if they exist).

Putting into the equilibrium framework, and if we restrict ourselves to graph filters, then we have
\begin{align*}
F_0(\vx) &= \argmin{\vz \in \R^n} \;\; f_0(\vz) + \frac{\alpha}{2}\|\vz - \vx\|^2,\\
F_i(\vx) &= \mW_i \vx, \quad\quad i = 1, \ldots, k.
\end{align*}
At equilibrium, the solution $(\vxhat, \{\vuhat_i\}_{i=0}^k)$ satisfies
\begin{align}
F_i(\vxhat + \vuhat_i) = \vxhat, \quad\mbox{and}\quad \sum_{i=0}^k \mu_i \vuhat_i = \vzero,
\label{eq: CE multi prior condition}
\end{align}
which is a generalization of the case $k=2$ in \eref{CE, eq1} and \eref{CE, eq2}. In \eref{eq: CE multi prior condition}, The mappings $\{F_i\}_{i=0}^k$ are called \emph{agents}, whose relative strengths are controlled by the parameters $\{\mu_i\}_{i=0}^k$. Without loss of generality we assume that $\sum_{i=1}^k \mu_i = 1$, and $\mu_0$ can take any positive value. The intuition behind \eref{eq: CE multi prior condition} is the concept of assembling individual experts to produce an overall better result. For example, $F_i$ could represent image denoisers trained by different image classes, e.g., human faces, buildings and plants, etc. Given an unknown scene, it is likely that one of the denoisers would perform better than the others. Therefore, an ensemble of weak experts can improve the robustness against model mismatch.

To make the following discussion more concrete, let us focus on the denoising problem so that $f_0(\vx) = \frac{1}{2}\|\vx - \vy\|^2$ for some noisy signal $\vy$. Then, the agent $F_0$ is defined as
\begin{align*}
F_0(\vx) &= \argmin{\vz} \; \frac{1}{2}\|\vz - \vy\|^2 + \frac{\alpha}{2}\|\vz - \vx\|^2 = \frac{ \vy + \alpha \vx}{1 + \alpha}.
\end{align*}
The equilibrium condition can be determined by substituting $\vx = \vxhat + \vuhat_0$, leading to
\begin{equation*}
\vxhat = F_0(\vxhat + \vuhat_0) = \frac{\vy + (\alpha \vxhat + \vuhat_0)}{1 + \alpha}.
\end{equation*}
Rearranging the terms gives $\vuhat_0 = \frac{1}{\alpha}(\vxhat - \vy)$. Combining with other agents, i.e.,
$$
\vxhat = F_i(\vxhat + \vuhat_i) = \mW(\vxhat + \vuhat_i) \;\; \Rightarrow \;\; \vuhat_i = (\mW_i^{-1} - \mI)\vxhat,
$$
we can substitute into $\sum_{i=0}^k \mu_i \vuhat_i = \vzero$ to show that the consensus equilibrium solution satisfies
\begin{equation}
\left(\frac{\mu_0}{\alpha}\mI + \sum_{i=1}^k \mu_i \left(\mW_i^{-1}-\mI \right) \right)\vxhat = \frac{\mu_0}{\alpha}\vy, \quad i = 1,\ldots,k.
\label{eq: CE multi prior condition 2}
\end{equation}

If we want to write \eref{eq: CE multi prior condition 2} in terms of MAP, then a solution $\vxhat$ satisfying \eref{eq: CE multi prior condition 2} is the minimizer of
\begin{equation}
\vxhat = \argmin{\vx \in \R^n} \frac{\mu_0}{2}\|\vx - \vy\|^2 + \frac{\alpha}{2} \vx^T \left( \sum_{i=1}^k \mu_i\mW_i^{-1} - \mI \right) \vx.
\label{eq: CE multi prior condition 3}
\end{equation}
The difference of \eref{eq: CE multi prior condition 3} with the single-prior case, e.g., \eref{eq: equilibrium phi}, is that we replace $\mW^{-1}$ by a linear combination of multiple priors $\sum_{i=1}^k \mu_i\mW_i^{-1}$.

To visually see the benefit of the combination, we consider $k = 5$ filters $\{\mW_i\}_{i=1}^k$ where the filters $\mW_i$ are generated from kernels with different denoising strengths, see \eref{eq: K}. The combination weights are set as $\mu_i = 1/k$ for simplicity, and we set $\mu_0 = 1$ and $\alpha = 0.005$. The equilibrium condition of this combined filter gives a solution of the form
\begin{equation}
\vxhat = \left(\mu_0\mI + \alpha\left( \sum_{i=1}^k \mu_i \mW_i^{-1} - \mI  \right)\right)^{-1}\mu_0 \vy.
\label{eq: CE inversion matrix 1}
\end{equation}
For comparison, we also consider solutions obtained from a single filter. That is, we use only one $\mW_i$ with $\mu_i = 1$. This leads to a solution
\begin{equation}
\vxhat_i = \left(\mu_0\mI + \alpha( \mW_i^{-1} - \mI)\right)^{-1} \mu_0 \vy,
\label{eq: CE inversion matrix 2}
\end{equation}
for $i=1,\ldots,k$. The performance of the solutions are measured in terms of PSNR, and the PSNRs of this particular experiment are reported in \fref{fig: CE}. Here, the legend $\mW_i$ denotes the individual filters we use, and $\mW$ denotes the combined filter. In \fref{fig: CE} we also plot the eigenvalues of the respective inversion matrices. As one can see, the combined estimate staying in the middle $\mW_3$ and $\mW_4$, which are the two best performing matrices.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\linewidth]{CE}
\caption{Combining multiple priors. The curves show the eigenvalues of the inversion matrices defined in \eref{eq: CE inversion matrix 1} and \eref{eq: CE inversion matrix 2}. In this experiment we set $\rho = 0.005$, $\mu_0 = 1$, $\mu_i = 0.2$ for $i=1,\ldots,5$.}
\label{fig: CE}
\vspace{-2ex}
\end{figure}

\subsection{Choosing the Weights}
Is it possible to determine a set of weights $\{\mu_i\}_{i=1}^k$ so that they are optimal in some sense? In principle this is doable by minimizing the mean squared error between the $\vxhat$ defined in \eref{eq: CE inversion matrix 1} and the ground truth (if we have access to). However, such approach is computationally prohibited because the $\{\mu_i\}_{i=1}^k$ are located in the denominator of \eref{eq: CE inversion matrix 1}.

A slightly more ad hoc approach is to solve \eref{eq: CE inversion matrix 2} for every $\vxhat_i$, and define $\vxhat$ as a linear combination of them:
\begin{equation}
\vxhat = \sum_{i=1}^k \mu_i \vxhat_i = \sum_{i=1}^k \mu_i \mA_i^{-1} \vy,
\label{eq: CE multi prior condition 4}
\end{equation}
where $\mA_i \bydef \left(\mI + (\alpha/\mu_0) ( \mW_i^{-1} - \mI)\right)$. An interpretation of \eref{eq: CE multi prior condition 4} is the solution of an optimization
\begin{align}
&\min_{\vxhat, \{\vxhat_i\}_{i=1}^k} \quad \sum_{i=1}^k \mu_i  \left\| \vxhat - \vxhat_i \right\|^2 \notag \\
&\mbox{subject to}  \quad \left(\frac{\mu_0}{\alpha}\mI + \left(\mW_i^{-1}-\mI \right) \right)\vxhat_i = \frac{\mu_0}{\alpha}\vy,
\label{eq: CE multi prior condition 5}
\end{align}
which says that among the individual equilibrium solutions $\vxhat_i$, we are finding an average $\vxhat$ that is close to all $\vxhat_i$'s.

\begin{table*}[!]
\footnotesize{
\begin{tabular}{c|cc|cc|cc|cc|cc|cc|cc|cc}
\hline
          & \multicolumn{8}{c|}{Estimated $\mW$} & \multicolumn{8}{c}{Oracle $\mW$} \\
\hline
          & 80\%  &       & 60\%  &       & 40\%  &       & 20\%  &       & 80\%  &       & 60\%  &       & 40\%  &       & 20\%  &     \\
Img       & Lps   & PnP   & Lps   & PnP   & Lps   & PnP   & Lps   & PnP   & Lps   & PnP   & Lps   & PnP   & Lps   & PnP   & Lps   & PnP \\
\hline
Barb      & 26.86 & 28.97 & 26.30 & 27.88 & 25.75 & 26.80 & 23.64 & 24.36 & 50.47 & 53.50 & 48.72 & 52.41 & 47.14 & 50.58 & 44.26 & 49.11\\
Boat      & 26.64 & 28.54 & 25.93 & 27.31 & 25.16 & 25.96 & 23.04 & 23.46 & 50.16 & 51.92 & 48.84 & 51.89 & 47.10 & 50.01 & 43.64 & 48.01\\
Cmra      & 26.39 & 27.88 & 25.30 & 26.33 & 24.13 & 24.81 & 22.01 & 22.26 & 49.67 & 50.53 & 48.32 & 50.20 & 46.99 & 47.92 & 44.28 & 45.62\\
Coup      & 26.55 & 28.46 & 25.87 & 27.25 & 25.19 & 25.98 & 23.22 & 23.67 & 50.25 & 52.72 & 48.18 & 50.81 & 46.98 & 50.79 & 43.96 & 48.51\\
Hous      & 27.28 & 30.52 & 27.22 & 29.83 & 26.95 & 28.67 & 25.07 & 25.82 & 50.66 & 53.56 & 49.16 & 53.39 & 46.81 & 50.72 & 43.63 & 48.52\\
Lena      & 26.98 & 29.77 & 26.62 & 28.70 & 26.25 & 27.57 & 24.28 & 24.94 & 49.62 & 53.13 & 48.29 & 52.18 & 47.45 & 51.44 & 43.93 & 48.15\\
Man       & 26.81 & 29.10 & 26.39 & 28.12 & 25.80 & 26.85 & 23.95 & 24.58 & 50.19 & 53.27 & 49.05 & 52.06 & 47.23 & 51.17 & 44.08 & 49.25\\
Pepp      & 26.36 & 28.41 & 25.64 & 27.10 & 24.67 & 25.51 & 22.66 & 23.21 & 49.69 & 51.12 & 49.20 & 51.26 & 46.61 & 50.64 & 41.98 & 46.77\\
\hline
\end{tabular}}
\caption{Inpainting Experiment. The input images contains of uniformly random sampled pixels at a rate of 80\%, 60\%, 40\%, and 20\%. The reconstruction methods are graph Laplacian (Lps) and PnP ADMM (PnP). The regularization parameter $\rho$ and $\lambda$ are optimized for every image and every configuration. }
\label{tab: inpainting}
\end{table*}

The weighted average in \eref{eq: CE multi prior condition 4} provides a simple way of determining $\{\mu_i\}_{i=1}^k$. Assuming that we have access to the ground truth $\vx$, the weights $\{\mu_i\}_{i=1}^k$ can be found by minimizing
\begin{equation}
\min_{\{\mu_i\}_{i=1}^k} \; \left\|\vx - \sum_{i=1}^k \mu_i \vxhat_i \right\|^2 \;\; \mbox{subject to} \;\; \sum_{i=1}^k \mu_i = 1.
\label{eq: CE multi prior condition 6}
\end{equation}
Letting $\vmu = [\mu_1,\ldots,\mu_k]^T$, and defining the matrix $\mSigma \in \R^{k \times k}$ such that $[\mSigma]_{ij} = (\vx-\vxhat_i)^T(\vx-\vxhat_j)$, then \eref{eq: CE multi prior condition 6} can be written as
\begin{equation}
\min_{\vmu} \;\; \vmu^T\mSigma\vmu \;\;\;\; \mbox{subject to} \;\;\;\; \vmu^T\vone = 1.
\label{eq: CE multi prior condition 7}
\end{equation}
Closed-form solution of \eref{eq: CE multi prior condition 7} exists, and is given by
\begin{equation}
\vmu^* = \frac{\mSigma^{-1}\vone}{\vone^T \mSigma^{-1}\vone}.
\end{equation}
Readers interested in the details can refer to \cite{Choi_Elgendy_Chan_2018_TIP} by Choi et al., where they discuss the solution's geometry and its extension to cases where ground truth is not available.

\section{Further Discussions}
We close the paper by discussing a few issues that might be of interest to readers.

\vspace{2ex}
\noindent $\bullet$ \textbf{Experiment on Images}. Thus far the numerical results are all based on 1D signals. Readers who do image processing may wonder how well can the theory be translated to images. To articulate this problem, we consider an inpainting problem which has a forward model as
\begin{equation}
\vy = \mS\vx + \veta.
\end{equation}
In this equation, $\mS \in \R^{m \times n}$ is a sampling matrix which either selects or not selects a pixel from the input image $\vx$. The noise vector $\veta$ is assumed i.i.d. Gaussian with standard deviation $\sigma_{\eta} = 0.05$. Therefore, the objective function is $f(\vx) = \frac{1}{2}\|\mS\vx - \vy\|^2$.

The result of this experiment is shown in Table~\ref{tab: inpainting}, where we compare graph Laplacian with PnP at different sampling ratios (20\%, 40\%, 60\%, 80\%). The $\mW$ matrix is the non-local means weight matrix symmetrized by the Sinkhorn-Knopp algorithm. For the estimated case, $\mW$ is constructed from a pre-estimated image, where the pre-estimation is done using the Sheperd interpolation algorithm \cite{Shepard_1968}. Tested on 8 ``standard'' images, the PSNR values are recorded and compared. The regularization parameters $\rho$ (or $\lambda$) are optimized for every image and for every testing condition. Therefore, the PSNRs shown in Table~\ref{tab: inpainting} are the best possible results of the methods. For both oracle and estimated $\mW$, PnP ADMM is consistently better than graph Laplacian, which is coherent to the theoretical arguments.

\vspace{2ex}
\noindent $\bullet$ \textbf{Complexity Consideration}. When we compare \eref{eq: Lps analysis 1} and \eref{eq: Lps analysis 2}, we tend to think that \eref{eq: Lps analysis 1} is computationally less intensive because it is a standard least squares which can be solved using conjugate gradient (or other linear solvers). The PnP solution \eref{eq: Lps analysis 2}, on the other hand, requires spectral decomposition of $\mW^{-1}$ so that we can obtain the solution via, e.g., by solving an optimization \eref{eq: decompose PnP optimization}. However, a complexity analysis based on this line of argument is not totally correct because in order to obtain the PnP solution we actually use the ADMM algorithm \eref{eq:ADMM2,x}-\eref{eq:ADMM2,u}. Thus, the real complexity bottleneck is at the inversion, i.e., \eref{eq:ADMM2,x}, not the graph filter $\mW$. For specific problems such as deblurring where $f(\vx)$ has structures, e.g., toeplitz matrix, the inversion could be done in $\calO(n \log n)$ via Fourier transform.

\vspace{2ex}
\noindent $\bullet$ \textbf{Statistical Interpretations of $\mW$ and $\mW^{-1}$}. If we regard $\mW$ as some kind of covariance matrix, then the inverse $\mW^{-1}$ is the precision matrix. The $(i,j)$-th element of the covariance matrix measures the \emph{correlation} between node $i$ and node $j$ of the graph, whereas the $(i,j)$-th element of the precision matrix measures the \emph{partial correlation} between node $i$ and node $j$ condition on all the other $n-2$ nodes. Partial correlation is useful in handling situations where a high correlation between $i$ and $j$ is caused by some some third party node $k$. Thus by conditioning on $k$ we can remove the dependency of node $k$, and leaving just the partial correlation between $i$ and $j$.

As an example, in \fref{fig: W and invW} we show a ground truth signal $\vx$ used to generate a filter matrix $\mW$. We then compute the pseudo-inverse $\mW^{+}$ by keeping on the 50 leading eigenvectors, and plot its magnitude in the log-scale. (Note: We cannot compute the actual inverse $\mW^{-1}$ when $\mW$ is not invertible.) As we can see, while $\mW$ shows a banded diagonal behavior due to a spatial constraint we put, the inverse matrix $\mW^+$ shows a very strong signal dependent behavior.

\begin{figure}[h]
\centering
\begin{tabular}{cc}
\includegraphics[width=0.35\linewidth]{W}&
\includegraphics[width=0.35\linewidth]{invW}
\end{tabular}
\includegraphics[width=0.8\linewidth]{z0}
\caption{[Top Left] $\mW$. [Top Right] $\mW^{+}$, the pseudo-inverse of $\mW$. We plot $\mW$ and $\mW^{+}$ in log absolute scale so that we can visualize them. [Bottom] The ground truth signal $\vx$ used to generate $\mW$.}
\label{fig: W and invW}
\end{figure}

\vspace{2ex}
\noindent $\bullet$ \textbf{Connection to RED}. A very similar idea to PnP ADMM is the regularization by denoising (RED) \cite{Romano_Elad_Milanfar_2017}. In RED, the idea is to define the regularization by
\begin{equation}
g(\vx) = \frac{1}{2}\vx^T(\vx - \calD_{\sigma}(\vx)).
\end{equation}
One advantage of such formulation is that we are able to interpret the regularization as the denoiser $\calD_{\sigma}$ is explicitly expressed in $g$. However, as recently discussed by Reehorst and Schniter \cite{Reehorst_Schniter_2018}, while RED offers an interpretable regularization, the existence of $g$ is problematic if $\calD_{\sigma}$ is not symmetric.

If we restrict $\calD_{\sigma}$ to symmetric smoothing filters so that $\calD_{\sigma}(\vx) = \mW\vx$, then the RED regularization is nothing but
\begin{equation*}
g(\vx) = \frac{1}{2}\vx^T(\vx - \mW \vx) = \frac{1}{2}\vx^T (\mI - \mW) \vx,
\end{equation*}
which is exactly the graph Laplacian regularization. As we have shown in the analysis above, under the oracle setting, this does not have the noise rejection capability as PnP has.

\vspace{2ex}
\noindent $\bullet$ \textbf{Dynamics}. In this paper we assume that $\mW$ is defined through some pre-filtered signal, yet in practice $\mW$ is updated with the ADMM iterations. If $\mW$ is updated dynamically, obtaining a closed-form expression of the consensus equilibrium point would become very difficult. This is an open question, and more studies are needed to fully understand how to characterize the solution of this type.

\vspace{2ex}
\noindent $\bullet$ \textbf{Duality}. There is an interesting duality relationship between the graph Laplacian $(\mI-\mW)$ and the PnP prior $(\mW^{-1}-\mI)$. In fact, one can show that
\begin{equation}
(\mI-\mW)^{-1} - (\mW^{-1}-\mI)^{-1} = \mI.
\end{equation}
This result is reminiscent to one of the recent theoretical studies by Feizi et al \cite{Feizi_Suh_Xia_2017} on generative adversarial network (GAN). The analogy is that if $(\mI-\mW)$ represents a discriminator for the real data, then $(\mI-\mW^{-1})$ represents a discriminator for the generated data. It would be interesting to derive similar duality interpretation from a graph perspective.

\vspace{2ex}
\noindent $\bullet$ \textbf{Beyond Image Restoration}. Thus far we have exclusively discussed imaging applications. If we consider $\mW$ as a forward diffusion, then $\mW^{-1}$ is the reversal diffusion. It would be interesting to investigate applications in general graph domain.

\section{Conclusion}
We presented a formal performance analysis for the PnP ADMM algorithm. The analysis is based on the class of graph filters. Our findings suggest a few reasons why PnP ADMM performs better than the conventional graph Laplacian regularization. From the direct mean squared error (MSE) perspective, the inverse $\mW^{-1}$ associated with PnP ADMM ensures that the pure noise components are strongly attenuated during the ADMM algorithm, whereas they are less attenuated by the conventional graph Laplacian. As a result, the MSE of the PnP ADMM estimate has a much lower variance than the conventional graph Laplacian. PnP ADMM shows gain in bias, but the bias is typically smaller than the drop in variance unless the filter $\mW$ is poorly estimated.

To further understand the problem, this paper also introduces a performance analysis approach from an equilibrium perspective. Instead of deriving the solution from a MAP optimization, the equilibrium approach allows us to derive the solution by inspecting the equilibrium condition of the ADMM subproblems. This significantly alleviates the difficulty of assuming $\mW^{-1}$ to exist, tackling general linear inverse problems, and problems with multiple priors. We presented a number of equivalent optimization formulations and explained the behavior of the equilibrium solution.

As more sophisticated image restoration modules, e.g., deep neural networks, are being integrated into traditional model based recovery frameworks, we hypothesize that equilibrium will play an important role in analyzing the performance of these methods.

\section{Acknowledgement}
The author thanks Prof. Charles Bouman and Prof. Greg Buzzard of Purdue University for many inspiring discussions, Dr. Brendt Wolhberg of Los Alamos National Lab for offering valuable feedbacks on this paper, and Prof. Gene Cheung of York University for sharing thoughts on graph signal processing. This work is supported, in part, by the National Science Foundation under grants CCF-1763896 and CCF-1718007.

\balance
\bibliographystyle{IEEEbib}

redpajama_set_name: RedPajamaArXiv


The ozone bed sterilizer KZ-X-DL1 made in Guang Dong, China has just been released into various markets around the globe over the past few months. On June 3, 2009, Guangzhou Kangzhen Medical Equipment Co. Ltd announced the increased availability of this portable machine. The manufacturer also took this opportunity to point out some of the equipment's benefits.
Although the bed and mattress must be encased in a specially sealed plastic bag for the system to work properly, sheets and other items remain in place in the room where they are used by the patient. Once hooked up, the equipment runs through an automatic cycle that creates its own ozone for sterilization. Through a set of vacuum and pump phases, deeply embedded pathogens are directly exposed to the ozone and killed.
This procedure can be used in addition to normal steam pressing of linens for sanitization purposes. In this way, ozone exposure creates an additional level of protection against cross-contamination in hospital settings. This method is superior to UV light as an additional step for sterilization of bedding because light only sanitizes the surface. Getting a pillow clean requires the sterilant to reach all the nooks and crannies where germs reside.
Hospitals are well known for their tendency to promote the spread of nosocomial diseases. Immune systems that are already stressed tend to be more susceptible to these secondary infections. The potential of being in close contact with other patients and the items that they may have touched magnifies the risks involved in any stay. The challenge medical personnel face is to create a healthy environment within a setting where hundreds if not thousands of people with different illnesses are routinely exposed to each other. Infection control and decontamination of bedding is a constant concern for hospitals all over the world.
Poor sanitization practices when dealing with patient linens is often the 'Achilles heel' of infection control programs. Lately, China in particular has been proactive in limiting the chances of cross infections. During the past few months they have sometimes been criticized for quarantining visitors from countries that have an ongoing rampant epidemic of the H1N1 flu.
Chinese health officials have set up special hotels that continue to be used strictly for this purpose. Everything in those rooms must be thoroughly decontaminated for the protection of quarantine patients and future tourists. This new product from Guangzhou Kangzhen would no doubt come in handy in this situation.
Actually, it wouldn't be a bad idea to use ozone sterilization as a standard practice in any type of high traffic setting. Hospitals are simply the logical place to start in the U.S. Perhaps hotels around the world will be next.
About Guangzhou Kangzhen Medical Equipment Co.
The company was founded in 2003. Since then it has built a reputation in China for excellent research, development, and application of advanced technology. Their goal is to bring to market innovative products using state of the art medical technology that streamline or improve patient care.
If you are interested in purchasing sterilization or autoclaving equipment, take a moment to browse our free report. It can help you discover which features will meet your needs and how to maximize yoru return on investment.


Below is FINRA's press release of the $200M restitution order the SEC and FINRA levied on Regions Morgan Keegan (RMK) related to bond funds it misrepresented to investors. Big for investors. We had blogged about this in April 2010 when the SEC charged Morgan Keegan & Company, Morgan Asset Management and two employees with fraudulently overstating the value of securities backed by subprime mortgages. These were the first federal government allegations related to the Regions Morgan Keegan bond funds which lost significant value in 2008. Many law firms around the U.S., including our firm, have been retained by investors who lost money in the Regions Morgan Keegan funds. Big fight will be to get fines into evidence during the pending arbitrations.
WASHINGTON – The Financial Industry Regulatory Authority (FINRA), the Securities and Exchange Commission (SEC) and five state regulators from Alabama, Kentucky, Mississippi, South Carolina and Tennessee announced today that each has settled enforcement proceedings against Morgan Keegan & Company, Inc. Morgan Keegan will pay restitution of $200 million for customers who invested in seven affiliated bond funds, including the Regions Morgan Keegan Select Intermediate Bond Fund (Intermediate Fund). Morgan Keegan's affiliate, Morgan Asset Management, managed the funds.
FINRA's settlement includes findings that Morgan Keegan failed to establish, maintain and enforce an adequate supervisory system, including written supervisory procedures reasonably designed to achieve compliance with NASD rules. Morgan Keegan's supervisory system and written procedures were not reasonably designed to ensure that its sales literature disclosed certain information as to risk and did not contain exaggerated claims. As a result, Morgan Keegan failed to adequately describe the nature, holdings and certain risks of the Intermediate Fund. In addition, beginning in 2007 when the particular risks associated with the Intermediate Fund's holdings began to impact negatively the holdings in the fund, Morgan Keegan failed to take steps reasonably designed to revise its advertising materials to inform customers of the specific risks of investing in the fund under the current market conditions.
This investigation was conducted by Gino Ercolino, David Fenimore, Gregory Firehock, Theresa Ridder and Richard Santiago of the Enforcement Department, with the assistance of FINRA's Advertising Regulation Department.
Investors can obtain more information about, and the disciplinary record of, any FINRA-registered broker or brokerage firm by using FINRA's BrokerCheck. FINRA makes BrokerCheck available at no charge. In 2010, members of the public used this service to conduct 17.2 million reviews of broker or firm records. Investors can access BrokerCheck at www.finra.org/brokercheck or by calling (800) 289-9999. Investors may find copies of this complaint, as well as other disciplinary documents, in FINRA's Disciplinary Actions Online database.


Client referrals are the best way Ken Nissen Hairstylist can connect with new customers. Current clients are invited to refer two friends and, when they book an appointment, you receive a $20 credit. Or ask your friends to mention you when scheduling their first appointment.


text:
Sunny Hostin '94 J.D. to host 'Truth About Murder' series
Author: Notre Dame Law School
Sunny Hostin's latest project is very different from what she does as co-host of ABC's daytime show "The View," reports USA Today.
Hostin, who received her J.D. from Notre Dame Law School in 1994, will host a new series called "Truth About Murder" on Investigation Discovery. The six-episode series, which starts October 22, relies on Hostin's legal knowledge from her career as a prosecutor and her own experience with crime as a witness to her uncle's stabbing when she was a child.
Read the USA Today feature about Hostin and her new show.

redpajama_set_name: RedPajamaCommonCrawl


Automatically release of units at a specified degree of list.
Different configuration of the list angle and different activation sequences of outputs are available.
Can also release units manually by local or remote push button.
Vessels with increased risk of capsize due to working conditions such as fishing vessels and workboats.
With current float-free release systems on the market, safety equipment is released first when it reaches a depth of 1,5 to 4 metres. When a ship or a vessel capsizes without sinking, there is therefore a risk that liferafts and EPIRBs are trapped under the vessel, or never released at all if not submerged to a sufficient water depth.
With List Angle Detection (LAD) liferafts and Epirbs can be automatically released at a specified degree of list when a vessel capsizes. The released safety units reach the surface before the ship flips around, significantly reducing the risk for it to be trapped or entangled in constructions on deck.
ERRS LAD can be controlled from several remote release positions or automatically by a HPS (Hydrostatic Pressure Switch) or a water sensor. Output number two is activated automatically 2 seconds after output number one has been activated as a standard configuration. For activation of relay outputs the unit has to be configured for this function and an ERRS Addon modul (HM-0464) is needed.
ERRS LAD automatically performs a system check once every hour. The system check monitors battery, emergency power voltage, ERU circuits and the wiring to an external activation switch. The battery condition is checked every 24 hours. If the system check detects a deviation an alarm message will appear on the LED indicator and the alarm output is activated. It is possible to initiate a manual system check which is useful at raft service.
ERRS is powered by the ship's emergency power supply. Each control panel is also equipped with a back-up battery to ensure operation in all emergency conditions.


Cathy has worked in the arena since the very start – a staggering 18 years and 3 days! (Yes, she still loves it!) She explains that she started 6 months before the official opening of the arena, situated in temporary buildings right across the water from the arena.
For the first 6 months, she explains it was all VERY different. There was no computerised system for ticketing so when the day of an on-sale rolled around they had books of tickets behind the counter. To make sure they knew that customer were sitting in the correct seats they would be ticking the seats off using a picture of the floor plan. (Phew, we're glad those days are over!) Cathy explains the arena and Ticketmaster system has now become fully digital, so thankfully none of that happens any longer!
Cathy explains there are many more words she could use, however she settles on meticulous (could you imagine getting the seating of a show wrong? We wouldn't want the stress!), enjoyable and challenging, as every show is completely different!
What's a normal day for Box Office Manager?
Every day is different for Cathy. As Box Office Manager, she oversees all ticketing sales for every show that comes to the arena! To start off, she'll get working on the executive sales reports through the Ticketmaster system which we run off. These will be analysed and sent to the respective show promoters.
Then comes the fun part of setting up shows from scratch! (That's usually around 10,000 seats per show!) Again, using the TM's system she'll implement all pricing structures in each open area, cordoning off certain stage holds that are required due to the production of the show. (Ever wondered why some shows will have 'production holds' released a few days before a show? This is why). She then will do the same for the Heineken and West Lounge, which are our executive lounges!
Throughout the day Cathy will be liaising with show promoters to gauge how their shows are selling and if there is anything they'd like to change around the shows. She'll then liaise with Ticketmaster Ireland's offices in Belfast and Dublin reviewing ticketing strategy and digital developments for the Box Office.
Another part of her job is on an event day. Box Office is one of the busiest departments with people picking up tickets, sometimes losing tickets (a BIG disaster when you want to come to see your favourite show) and facilitating VIP Guest lists. Then at the end of the event, she'll go backstage to help settle the show accounts with the tour accountant!
Since she's been here for the long haul –18 years and one day – she has lots! Firstly, she explains it is A LOT busier than 18 years ago. The pace of events coming to Belfast has increased and this is something she really encourages artists to do! With that comes the diversity of events coming to the arena. Now we can see a Giants game, a fully standing Kylie show and a fully seated André Rieu show all in 3 days!
The favourite part of her role is to make sure the audience is happy! From setting up the show to then seeing it through when the main act takes the stage, she explains the excitement is what keeps her job so enjoyable! She explains that on the morning of an on-sale she is greeted by the excitement of customers as they buy tickets. Then on an event day, it's even more intense as everyone anticipates the main act coming on stage! When everything runs smoothly with ticketing on the day of the show that makes Cathy's life easy!
Favourite show to work on?
Usually, no-one can give me an answer to this question but Cathy has one straight away! She still remembers Bruce Springsteen like it was yesterday. He only played two shows, one in London and one in Belfast that year.
Cathy explains this was another diverse show as they were selling tickets BEHIND the stage. She wasn't sure how this was going to work until the production went underway and the stage went into place. She still explains it as "the best experience" for the 64 people seated there, as it was as if they were part of the show, rather than in the crowd!
Cathy names Pink as one of the most visually and unusual shows that has come to the arena. From the acrobatics to her singing, she was able to be so diverse and give the audience an amazing show. However, she goes back to The Boss for her ultimate show. She still describes as "the best show she's ever seen of ALL TIME" (that's how excited she is talking about it). The show in question was Bruce Springsteen Seeger Sessions. She explains that at one point there must have been a good "60 musicians on the stage, I'm not exaggerating" (she was exaggerating, I googled it and it was 20).
The dream for Cathy would be Elvis. She still watches clips of him online and said she would have done anything to see him play live (this may show her age)!
Cathy has no qualms about this question.
Keep moving forward…and in this new age of ticketing and Box Office, it's very apt!


text:
Angry President Demands Obedience
On September 12, 2021 By AZ ConservativeIn Uncategorized
By: John Semmens
An angry President Biden took to the airwaves to announce new marching orders for America. Despite the fact that far less than 1% of the population has been felled by COVID, 70% of the adult population having been vaccinated, and an unknown, but probably large percentage of the unvaccinated having natural immunity to the virus, the President unveiled a nation-wide mandate to force more people to take the vaccine.
"I've urged everyone to get vaccinated, but a distinct minority have refused to take my advice," Biden complained. "My patience is wearing thin. If people won't make the right decision voluntarily it is my obligation to take the decision out of their hands—as I had to do numerous times when my young sons Beau and Hunter got uppity with me and defiantly disobeyed my good advice. Accordingly, I am ordering all federal employees to take the vaccination. I am ordering all businesses with of 100 or more employees to make getting vaccinated a requirement for keeping their jobs. I am also recommending that all state and local governments impose a vaccination requirement on their employees."
The President dismissed the contention that he lacks the constitutional authority to unilaterally impose this requirement, saying "polls show that a majority of Americans favor universal vaccination. Enforcing the will of the majority on a minority is not tyranny. It is democracy."
Numerous Republican governors spoke out against the vaccine mandate. Gov. Greg Gianforte (R-Mt) asserted "President Biden's vaccination mandate is unlawful and un-American. We are committed to protecting Montanans' freedoms and liberties against this gross federal overreach." Gov. Kim Reynolds (R-Iowa) said "President Biden is taking dangerous and unprecedented steps to insert the federal government even further into our lives while dismissing the ability of Iowans and Americans to make healthcare decisions for themselves."
Undaunted by this outspoken opposition to his decree, Biden boasted "If these governors won't help us beat the pandemic, I will use my power as president to get them out of the way. I can start by withholding federal highway payments. If necessary I can send troops like President Eisenhower did to enforce school desegregation in Arkansas. If they still don't get the message, I can have them arrested like President Lincoln did to critics of his handling of the Civil War. We'll see if any of them have the nerve to face off against their very determined president."
Dr. Anthony Fauci called Biden's vaccination mandate plan "disappointingly moderate. If I were devising a mandate I would not have allowed for any opt outs. That includes getting regular COVID tests, having had a previous infection, existing medical conditions that might make vaccination more risky, and religious or moral objections. Our situation is too perilous for us to compromise. Freedom must give way to discipline. In my opinion, everyone must be vaccinated whether they like it or not."
In a seemingly discordant opinion, Vice-President Kamala Harris was quoted saying "when people are able to make choices without government interference for themselves in terms of their well being and the well-being of their family in consultation with whomever they may choose, we are a stronger society." Harris hastened to clarify that "I was talking about a woman's right to abort her unwanted baby, which is an entirely different thing than mandating everyone be vaccinated." The most striking difference is that the fatality rate for babies from abortion is nearly 100%, while the fatality rate for not being vaccinated is less than 1%.
In related news, Rutgers University barred a student from attending virtual classes because he has not been vaccinated. Rutgers spokeswoman Dory Devlin called the university's rule "fair. As President Biden has made clear, the disloyal who refuse to get vaccinated endanger us all. Even though this student will not directly infect anyone via virtual classes he is still a living embodiment of the kind of selfishness that is disunifying America. For this he must pay a price—either get vaccinated of face exclusion from all normal activities."
San Fran to Pay Criminals to Not Shoot People
Pfizer Announces "Covicure"