A Tourof Modern
Image Filtering
photo by david wagner[Peyman Milanfar]

R[New insights and methods,
both practical and theoretical]

ecent developments in computational imaging andrestoration have heralded the arrival and convergenceof several powerful methods for adaptive processing ofmultidimensional data. Examples include moving
least square (from graphics), the bilateral filter (BF)
and anisotropic diffusion (from computer vision), boosting, kernel, and spectral methods (from machine learning), nonlocal
means (NLM) and its variants (from signal processing), Bregman
iterations (from applied math), kernel regression, and iterative
scaling (from statistics). While these approaches found their inspirations in diverse fields of nascence, they are deeply connected.

Digital Object Identifier 10.1109/MSP.2011.2179329
Date of publication: 5 December 2012

In this article, I present a practical and accessible framework
to understand some of the basic underpinnings of these methods, with the intention of leading the reader to a broad understanding of how they interrelate. I also illustrate connections
between these techniques and more classical (empirical) Bayesian approaches.
The proposed framework is used to arrive at new insights
and methods, both practical and theoretical. In particular, several novel optimality properties of algorithms in wide use suchas block-matching and three-dimensional (3-D) filtering(BM3D), and methods for their iterative improvement (or nonexistence thereof) are discussed.
A general approach is laid out to enable the performance
analysis and subsequent improvement of many existing filtering

IEEE SIGNAL PROCESSING MAGAZINE [106] january 20131053-5888/13/$31.00©2013IEEE

algorithms. While much of the material discussed is applicableto the wider class of linear degradation models beyond noise(e.g., blur,) to keep matters focused, we consider the problem of
denoising here.

similar analysis to what follows can and should be carried out
for more general inverse problems such as deconvolution,interpolation, etc.

(2)
y = z + e.
Introduction
Multidimensional filtering is the most fundamental operation
It has been realized for some time now that effective restorain image and video processing, and low-level computer vision.
tion of signals will require methods which either model theIn particular, the most widely used canonical filtering operation
signal a priori (i.e., are Bayesian) or learn the underlying charis one that removes or attenuates the effect of noise. As such,
acteristics of the signal from the given data (i.e., learning,
the basic design and analysis of image filtering operations form
nonparametric, or empirical Bayes’ methods.) Most recently,
a very large part of the image processing literature; thethe latter category of approaches has become exceedingly popresulting techniques often
ular. Perhaps the most striking
quickly spreading to the wider
recent example is the popularity
Recent developments inrange of restoration and reconof patch-based methods [7]–[10].
computational imaging andstruction problems in imaging.
This new generation of algorestoration have heraldedOver the years, many approaches
rithms exploit both local andthe arrival and convergencehave been tried, but only recently
nonlocal redundancies or “selfof several powerful methods
in the last decade or so, a great
similarities” in the signals being
for adaptive processing ofleap forward in performance hastreated. Earlier on, the BF [8] wasmultidimensional data.
been realized. While largely unacdeveloped with very much theknowledged in our community,
same idea in mind, as were its
this phenomenal progress has been mostly thanks to the adopspiritually close predecessors: the Susan filter [11], normaltion and development of nonparametric point estimation proceized convolution [12], and the filters of Yaroslavsky [13]. The
dures adapted to the local structure of the given
common philosophy among these and related techniques ismultidimensional data. Viewed through the lens of the denoisthe notion of measuring and making use of affinities between
ing application, here we develop a general framework fora given data point (or more generally patch or region) of interunderstanding the basic science and engineering behind these
est, and others in the given measured signal y. These similaritechniques and their generalizations. Surely this is not the firstties are then used in a filtering context to give higher weights
article to attempt such an ambitious overview, and it will likely
to contributions from more similar data values, and to propnot be the last; but the aim here is to provide a self-contained
erly discount data points that are less similar. The pattern recpresentation that distills, generalizes, and puts into proper conognition literature has also been a source of parallel ideas. In
text many other excellent earlier works such as [1]–[5], and,particular, the celebrated mean-shift algorithm [14], [15] is in
more recently, [6]. It is fair to say that this article is, by necesprinciple an iterated version of point-wise regression as alsosity, not completely tutorial. Indeed it does contain several
described in [1] and [2]. In the machine learning community,
novel results; yet these are largely novel interpretations, forthe general regression problem has been carefully studied, andmalizations, or generalizations of ideas already known ordeep connections between regularization, least-squares regresempirically familiar to the community. Hence, I hope that thesion, and the support vector formalism have also been estabenterprising reader will find this article not only a good overlished [16]–[19].
view, but as should be the case with any useful presentation, aDespite the voluminous recent literature on techniquessource of new insights and food for thought.
based on these ideas, simply put, the key differences between
So to begin, let us consider the measurement model for thethe resulting practical filtering methods have been relatively
denoising problem
minor, but rather poorly understood. In particular, the underlying framework for each of these methods is distinct only tothe extent that the weights assigned to different data points is
y i = z i + e i, for i = 1, f, n,(1)
decided upon differently. To be more concrete and mathematically precise, let us consider the denoising problem (2) again.
where z i = z (x i) is the underlying latent signal of interest at a
The estimate of the signal z (x) at the position x is found using
position x i = 6 x i, 1, x i, 2@T, y i is the noisy measured signal (pixel
value), and e i is zero-mean, white noise with variance v 2 . Wea (nonparametric) point estimation framework; specifically,
the weighted least squares problem
make no other distributional assumptions for the noise. The
problem of interest then is to recover the complete set of samnp l e s o f z (x), w h i c h w e d e n o t e v e c t o r i a l l y a s
zt (x j) = arg min / 6 y i - z (x j)@2 K (x i, x j, y i, y j),(3)z (x j)i =1
z = 6z (x 1), z (x 2), f, z (x n)@T from the corresponding data set
y. To restate the problem more concisely, the complete meawhere the weight (or kernel) function K ($) is a a symmetricsurement model in vector notation is given by (surely afunction with respect to the indices i and j. K ($) is also aIEEE SIGNAL PROCESSING MAGAZINE [107] january 2013

Classical Regression Filters
positive valued and unimodal function that measures the
Naturally, the most naive way to measure the “distance”“similarity” between the samples y i and y j, at respective posibetween two pixels is to simply consider their spatial Euclidean
tions x i and x j . If the kernel function is restricted to be only adistance; specifically, using a Gaussian kernel,
function of the spatial locations x i and x j , then the resulting
formulation is what is known as (classical, or not data-adap- xi - x j 2
K (x i, x j, y i, y j) = exp eo.
tive) kernel regression in the nonparametric statistics literah 2x
ture [20], [21]. Perhaps more importantly, the key difference
between local and nonlocal patch-based methods lies essenSuch filters essentially lead to (possibly space-varying) Gausstially in the definition of the range of the sum in (3). Specifiian filters which are quite familiar from traditional image procally, indices covering a small spatial region around a pixel ofcessing [13], [20], [21], [25]. It is possible to adapt theinterest define local methods, and vice versa.
variance (or bandwidth) parameter h x to the local image staInterestingly, in the early 1980s, the essentially identical
tistics, and obtain a relatively modest improvement in perforconcept of moving least-squares emerged independently [22],mance. But the lack of stronger adaptivity to the underlying[23] in the graphics community. This idea has since been
structure of the signal of interest is a major drawback of thesewidely adopted in computerclassical approaches.
graphics [24] as a very effective
The humble weighted
tool for smoothing and interpolaThe Bilateral Filter
least-squares, with locallytion of data in three dimensions.
Another manifestation of the fordata-dependent weights
Surprisingly, despite the obvious
mulation in (3) is the BF [8], [13]connections between moving
where the spatial and photometbased on similarity kernels,
least-squares and the adaptive filric distances between two pixels
explains an enormous class ofters based on similarity, their
are taken into account in separarecently popularized and very
kinship has remained largely hidble fashion as follows:
effective filters.
den so far.
- xi - x j 2- (y i - y j ) 2
o exp eK (x i, x j, y i, y j) = exp eo2
Existing Algorithms
hx
h 2y

Over the years, the measure of similarity K (x i, x j, y i, y j) has- x i - x j 2 - (y i - y j ) 2
3= exp )+.
(4)
been defined in a number of different ways, leading to a
h 2x
h 2y
cacophony of filters, including some of the most well-known
recent approaches to image denoising [7]–[9]. Figure 1 gives a
As can be observed in the exponent on the right-hand side, andgraphical illustration of how different choices of similarity
in Figure 1, the similarity metric here is a weighted Euclidean
kernels lead to different classes of filters, some of which we
distance between the vectors (x i, y i) and (x j, y j) . This approach
discuss next.
has several advantages. Specifically, while the kernel is easy toconstruct and computationally simple
to calculate, it yields useful local adaptivity to the given data. In addition, it
Bilateral Filter
has only two control parameters
K (xi - x) . K (yi - y )(h x, h y), which make it very convenient
to use. However, as is well known,ythis filter does not provide effective
The Euclidean Distance
The Photometricperformance in low signal-to-noise: dx 2 + dy 2
Distance
scenarios [3].
dy = |yi - y|K (yi - y)

The Geodesic
Distance

Nonlocal Meansx

LARKdx = |xi - x|
The Spatial DistanceK (xi - x)
Classic Kernel Regression[Fig1] Similarity metrics and the resulting filters.

IEEE SIGNAL PROCESSING MAGAZINE [108] january 2013

NonLocal Means
The NLM algorithm [7], [26], [27], originally proposed in [28] and [29], has
stirred a great deal of interest in thecommunity in recent years. At its core,
however, it is a relatively simple generalization of the BF; specifically, the photometric term in the bilateral similarity
kernel, which is measured point-wise, issimply replaced with one that is patchwise. A second difference is that (at leastin theory) the geometric distancebetween the patches (corresponding to the first term in thebilateral similarity kernel), is essentially ignored, leading tostrong contribution from patches that may not be physicallynear the pixel of interest (hence, the name nonlocal). To summarize, the NLM kernel is
K (x i, x j, y i, y j) = exp e- xi - x j 2- yi - y j 2
oexp e
o (5)
2
hx
h 2y

with h x " 3, where y i and y j refer now to patches of pixels centered at pixels y i and y j, respectively. In practice, two implementation details should be observed. First, the patch-wisephotometric distance y i - y j 2 in the above is in fact measured
as (y i - y j) T G (y i - y j), where G is a fixed diagonal matrix containing Gaussian weights, which give higher importance to thecenter of the respective patches. Second, it is rather computationally impractical to compare all the patches y i to y j, so
although the NLM approach in Buades et al. [28] theoreticallyforces h x to be infinite, in practice typically the search is limited
to a reasonable spatial neighborhood of y j . Consequently, in
effect, the NLM filter too is more or less local; or said anotherway, h x is never infinite in practice. The method in Awate et al.
[29], on the other hand, proposes a Gaussian-distributed sample that comes closer to the exponential weighting on Euclidean
distances in (5).
Despite its popularity, the performance of the NLM filter
leaves much to be desired. The true potential of this filtering
scheme was demonstrated only later with the optimal spatial
adaptation (OSA) approach of Boulanger and Kervrann [26]. In
their approach, the photometric distance was refined toinclude estimates of the local noise variances within eachpatch. Specifically, they computed a local diagonal covariancematrix, and defined the locally adaptive photometric distanceas (y i - y j) T V -j 1 (y i - y j) in such a way as to minimize an estimate of the local mean squared error (MSE). Furthermore,they considered iterative application of the filter as discussed
in the section “Improving the Estimate by Iteration.”
Locally Adaptive Regression(Steering) Kernels
The key idea behind this measure of similarity, originally proposedin [9], is to robustly obtain the local structure of images by analyzing the photometric (pixel value) differences based on estimatedgradients, and to use this structure information to adapt the shapeand size of a canonical kernel. The locally adaptive regression kernel (LARK) is defined as follows:
K (x i, x j, y i, y j) = exp " - (x i - x j) T C i (x i - x j) ,,(6)

where the matrix C i = C (y i, y j) is estimated from the given
data asCi = / >jz 2xi,1(x j)z xi,1(x j) z x i,2(x j)z xi,1(x j) z xi,2(x j)z 2x i,2(x j)H.

Specifically, z x i, (x j) are the estimated gradients of the underlying signal at point x i, computed from the given measurements
y j in a patch around the point of interest. In particular, the gradients used in the above expression can be estimated from thegiven noisy image by applying classical (i.e., nonadaptive)
locally linear kernel regression. Details of this estimation procedure are given in [30]. The reader may recognize the above
matrix as the well-studied “structure tensor” [31]. The advantage of the LARK descriptor is that it is exceedingly robust to
noise and perturbations of the data. The formulation is also theoretically well motivated since the quadratic exponent in (6)essentially encodes the local geodesic distance between the
points ( x i, y i) and ( x j, y j) on the graph of the function z (x, y)thought of as a two-dimensional (2-D) surface (a manifold)embedded in three dimensions. The geodesic distance was also
used in the context of the Beltrami-flow kernel in [32] and [33]in an analogous fashion.
*

Generalizations and Properties
The above discussion can be naturally generalized by definingTthe augmented data variable t i = 6 x Ti , y Ti @ and a general
Gaussian kernel as follows:
K (t i, t j) = exp " - (t i - t j) T Q (t i - t j) ,,Qx 0
E,
Q i, j = ;
0 Qy

(7)
(8)

where Q is symmetric positive definite (SPD).
Setting Q x = h12 I and Q y = 0, we have classical kernelxregression, whereas one obtains the BF framework whenQ x = h12 I and Q y = h12y diag[0, 0, f , 1, f , 0, 0]. The latterx
diagonal matrix picks out the center pixel in the elementwise difference of patches t i - t j . When Q x = 0 and Q y = h12y G,
we have the NLM filter and its variants. Finally, the LARK
kernel in (6) is obtained when Q x = C i and Q y = 0. More generally, the matrix Q can be selected so that it has nonzero
off-diagonal blocks. However, no practical algorithms withthis choice have been proposed so far. As detailed below, withan SPD Q, this general approach results in valid SPD kernels, a property that is used throughout the rest of our discussion. The definition of t given here is only one of many
possible choices. Our treatment in this article is equallyvalid when, for instance, t = T ^ x, y h is any feature derived
from a convenient linear or nonlinear transformation of theoriginal data.
The above concepts can be further extended using the powerful theory of reproducible kernels originated in functional
analysis (and later successfully adopted in machine learning[34], [35]) to present a significantly more general frameworkfor selecting the similarity functions. This will help us identify
the wider class of admissible similarity kernels more formally,
and to understand how to produce new kernels from ones
already defined [35]. Formally, a scalar-valued function K (t, s)over a compact region of its domain R n is called an admissible
kernel if

IEEE SIGNAL PROCESSING MAGAZINE [109] january 2013

■■ K is symmetric: K (t, s) = K (s, t)■■ K is positive definite. That is, fora n y c o l l e c t i o n o f p o i n t s t i,i = 1, f, n, the Gram matrix withelements K i, j = K (t i, t j) is positive
definite.
Such kernels satisfy some useful properties such as positivity, K (t, t) $ 0, and the Cauchy–Schwartz
inequality K (t, s) # K (t, t) K (s, s) .
With the above definitions in place, there are numerous
ways to construct new valid kernels from existing ones. We list
some of the most useful ones below, without proof [35]. Given
two valid kernels K 1 (t, s) and K 2 (t, s), the following constructions yield admissible kernels:
1) K (t, s) = aK 1 (t, s) + bK 2 (t, s) for any pair a, b $ 02) K (t, s) = K 1 (t, s) K 2 (t, s)3) K (t, s) = k (t) k (s), where k ($) is a scalar-valued function4) K (t, s) = p ^ K 1 (t, s) h, where p ($) is a polynomial with positive coefficients
5) K (t, s) = exp (K 1 (t, s)) .
Regardless of the choice of the kernel function, the weighted
least-square optimization problem (3) has a simple solution. In
matrix notation we can write
zt (x j) = arg min 6y - z (x j) 1 n@T K j 6y - z (x j) 1 n@,z (x j)(9)where 1 n = 61, f, 1@T, and K j = diag 6K (x 1, x j, y 1, y j), K (x 2, x j,
y 2, y j), f, K (x n, x j, y n, y j)@ . The closed-form solution to theabove is-1
zt (x j) = ^1 Tn K j 1 n h 1 Tn K j y(10)= c / K (x i, x j, y i, y j) m c / K (x i, x j, y i, y j) y i minonparametric approach in(3) can be further extended
to include a more generalexpansion of the signal z(x)in a desired basis. We briefly
discuss this case in “Generalization to Arbitrary Bases,”
but leave its full treatment for future research.
To summarize the discussion so far, we have presented ageneral framework that absorbs many existing algorithms as
special cases. This was done in several ways, including a generaldescription of the set of admissible similarity kernels, whichallows the construction of a wide variety of new kernels not considered before in the image processing literature. Next, we turn
our attention to the matrix formulation of the nonparametricfiltering approach. As we shall see, this provides a framework formore in-depth and intuitive understanding of the resulting filters, their subsequent improvement, and for their respective
asymptotic and numerical properties.
Before we end this section, it is worth saying a few words
about computational complexity. In general, patch-based methods are quite computationally intensive. Recently, many works
have aimed at both efficiently searching for similar patches, andmore cleverly computing the resulting weights. Among these,notable recent work appears in [36] and [37].

The performance of any
kernel-based denoising method
can be improved by some typeof iteration. The key is to use
the right iteration.

-1

i

K ij
=/
yi
K ij
/
i
i

The Martix Formulation and Its Properties
In this section, we analyze the filtering problems posed earlier
in the language of linear algebra and make several theoreticaland practical observations. In particular, we are able to not only
study the numerical/algebraic properties of the resulting filters,
but also to analyze some of their fundamental statisticalproperties.
To begin, recall the convenient vector form of the filters:zt (x j) = w Tj y,= / W ij y ii= w Tj y.

(11)

So in general, the estimate zt (x j) of the signal at position x j isgiven by a weighted sum of all the given data points y (x i),each contributing a weight commensurate with its similarity
as indicated by K ($), with the measurement y (x j) at the positionof interest. Furthermore, as should be apparent in (10), theweights sum to one. To control computational complexity, or todesign local versus nonlocal filters, we may choose to set theweight for some “sufficiently far-away” pixels to zero or a smallvalue, leading to a weighted sum involving a relatively small
number of data points in a properly defined vicinity of the sample of interest. This is essentially the only distinction between
locally adaptive processing methods (such as BL and LARK) andso-called nonlocal methods such as NLM. It is worth notingthat in the formulation above, despite the simple form of (10),in general we have a nonlinear estimator since the weights
W ij = W (x i, x j, y i, y j) depend on the noisy data. The(12)

where w Tj = 6W (x 1, x j, y 1, y j), W (x 2, x j, y 2, y j), f , W (x n, x j, y n,y j)@ is a vector of weights for each j. Writing the above at once
for all j we have
R TV
Sw 1 WSw T2 W(13)
zt = S W y = W y.
ShW
Sw TnWT X
As such, the filters defined by the above process can be analyzed as the product of a (square, n # n) matrix of weights Wwith the vector of the given data y. First, a notational matter:
W is in general a function of the data, so strictly speaking, thenotation W (y) would be more descriptive. But as we will
describe later in more detail, the typical process for computing
these weights in practice involves first computing a preliminary denoised “pilot,” or “prefiltered” version of the image,from which the weights are then calculated. This preprocessing, done only for the purposes of computing the parameters

IEEE SIGNAL PROCESSING MAGAZINE [110] january 2013

Generalization to Arbitrary Bases
The nonparametric approach in (3) can be further extended
to include a more general model of the signal z (x) in someappropriate basis. Specifically, expanding the regressionfunction z (x) in a desired basis z l, we can formulate the following optimization problem:
n

N

i =1

l =0

2

zt (x j) = arg min / =y i - / b l (x j) { l (x i, x j)G K ( y i, y, x i, x j), (S1)b l (x j )where N is the model (or regression) order. For instance,with the basis set z l (x i, x j) = (x i - x j) l, we have the Taylor
series expansion, leading to the local polynomial approachesof classical kernel regression [20], [71], [9]. Alternatively, thebasis vectors could be learned from the given image using amethod such as K-SVD [72]. In matrix notation we havebt (x j) = arg min 6y - U j b (x j)@T K j 6y - U j b (x j)@ ,b (x j )

where

R
V
Sz 0 (x 1, x j) z 1 (x 1, x j) f z N (x 1, x j)W
Sz 0 (x 2, x j) z 1 (x 2, x j) f z N (x 2, x j)W
Uj = SWhhh
SW
Sz 0 (x n, x j) z 1 (x n, x j) f z N (x n, x j)W
T
X

(S2)

(S3)

is a matrix containing the basis vectors in its columns; andwhere b (x j) = 6b 0, b 1, f, b N@T (x j) is the coefficient vector of
the local signal representation in this basis. Meanwhile, K j isthe diagonal matrix of weights as defined in (9). To maintainthe ability to represent the “DC” pixel values, we can insist
that the matrix U j contain the vector 1 n = [1, 1, f, 1] T as oneof its columns. Without loss of generality, we assume this tobe the first column so that by definition z 0 (x i, x j) = 1 for all iand j.

of the filter, reduces the effective variance of the input noiseand stabilizes the calculation of the weights so that they are
not a strong function of the input noise level. As such, it ispractically more useful to think of the weights as depending
more directly on the underlying (unknown) image W (z) . As
such, we shall simply denote the weight matrix going forward
as W. A more rigorous justification for this is provided later in
the next section, and supported in “Stability of Filter Coefficients to Perturbations Due to Noise.” Further justification isprovided in [75].
Next, we highlight some important properties of thematrix W, which lend much insight to our analysis. Referring to(10), W can be written as the product
W = D -1 K,(14)

where 6K@ij = K (x i, x j, y i, y j) = K ij and D is a positive definite
diagonal matrix with the normalizing factors D jj =diag {R i K ij}
along its diagonals. We can write
W = D -1 K = D -1/2 D -1/2 KD -1/2 D 1/2 .
1 44 2 44 3
L(15)

Denoting the jth row of U j as z Tj = 6z 0 (x j, x j), z 1 (x j, x j), f ,z N (x j, x j)@ we have the closed-form solution

zt (x j) = z Tj bt (x j)=

U K j U j h U K j y = w y.
1 4 4 4 44 2T 4 4 4 44 3
T
zj ^

-1

T
j

T
j

T
j

(S4)
(S5)

wj

This again is a weighted combination of the pixels, though arather more complicated one than the earlier formulation.
Interestingly, the filter vector w j still has elements that sum
to one. This can be seen as follows. We havew Tj U j = z Tj ^U Tj K j U j h U Tj K j U j = z Tj .
-1

Writing this explicitly, we observe
RV
S1 z 1 (x 1, x j) f z N (x 1, x j)W
S1 z 1 (x 2, x j) f z N (x 2, x j)W
w Tj S
W =61, z 1 (x j, x j),f,z N (x j, x j)@ .
hhShWS1 z 1 (x n, x j) f z N (x n, x j)WT
X
Considering the inner product of w Tj with the first columnof U j, we have w Tj 1 n = 1 as claimed.
The use of a general basis as described above carries one disadvantage, since the coefficient vectors w j are influenced
now not only by the choice of the kernel, but also by the
choice of the basis, the elements of w j can no longer be guaranteed to be positive. While this can be useful in the context
of filtering as it affords additional degrees of freedom, it doessignificantly complicate the analysis of the resulting algorithms. Most notably, the resulting weight matrix W will nolonger have nonnegative elements, and hence the powerfulmachinery afforded by the Perron–Frobenius theory is no longer directly applicable. As such, we must resort to a broader
definition of a generalized stochastic matrix [73].

It is evident that since K and D are SPD, then so is L. Meanwhile, W and L are related through a diagonal similarity transformation, and therefore have the same eigenvalues. Hence,
interestingly, W has real, nonnegative eigenvalues (even though
it is not symmetric). An alternative, though not necessarily
more intuitive, description of the action of W is possible interms of graphical models. Specifically, if we consider the data
( x i, y i) to be nodes on a (finite, undirected, weighted) graphwith weights between nodes i and j given by the kernel values
K (x i, x j, y i, y j), the matrix L = L - I is precisely the “graphLaplacian” [38], [5], [39]. The name is not surprising. In the
canonical case, W is a low-pass filter (though possibly nonlinear
and space-varying.) Since W and L share the same spectrum, L
is also a low-pass filter. Therefore, L = L - I is generically ahigh-pass filter, justifying the Laplacian label from the point ofview of filtering. This is illustrated in Figure 2.
Next, we note that W is a positive row-stochastic matrix—
specifically, by definition, the sum of each of its rows is equal
to one, as should be clear from (10). Furthermore, the structure of the pixel data as samples of a function on a finite,
fixed, lattice means that W can also be interpreted [5], [39] as
the transition matrix of an aperiodic and irreducible Markov

IEEE SIGNAL PROCESSING MAGAZINE [111] january 2013

Stability of Filter Coefficients to Perturbations Due to NoiseThe statistical analysis of the nonparametric filters W is quite
the noise variance is sufficiently small, the fourth-ordercomplicated when W is considered random. Fortunately, howmoment is even smaller; hence, the change in the resulting
ever, the stability results in [47] allow us to consider thecoefficient matrix is bounded by a constant multiple of theresulting filters as approximately deterministic. More specifismall noise variance.
cally, we assume that the noise e corrupting the data has i.i.d.
Approximations to the moments of the perturbation
samples from a symmetric (but not necessarily Gaussian) disdW = P
W - W are also given in [47] and [78]tribution, with zero-mean and finite variance v 2 . The resultsin [47] imply that if the noise variance v 2 is small relative toE6d W@ . E 6dD@ D -2 K - D -1 E6dK@,
the clean data z (i.e., when signal-to-noise ratio is high), theE6d W 2@ . E6dD 2@ D -4 K 2 + D -2 E6dK 2@ -E6dD dK@ D -3 5 K,
W = W (y) computed from the noisy data
weight matrix P
y = z + e is near the latent weight matrix W (z) . That is, with
D - D, and dK = M
K - K, and & denotes an elementwhere dD = N
the number of samples n sufficiently large [47],wise product. A thorough description of the perturbation anal2
ysis for the class of Gaussian data-adaptive weights is given inP(S6)
W - W F # p c 1 v (e2) + c 2 v (e4),[47] and [78]. To summarize, for a sufficiently small noise varifor some constants c 1 and c 2, where v (e2) and v (e4) denote theance, the matrix of weights can be treated, with high confisecond and fourth order moments of e , respectively, anddence, as a deterministic quantity that depends only on the“ # p ” indicates that the inequality holds in probability. When
content of the underlying latent image z.

chain defined on the entries of the data vector y. Hence, theclassical Perron–Frobenius theory is at our disposal [40],[41], which completely characterizes the spectra of such
matrices. In particular [42], [43], we have that W has spectral
radius 0 # m (W) # 1. Indeed, the largest eigenvalue of W issimple (unrepeated) and is exactly m 1 = 1, with corresponding eigenvector v 1 = (1/ n ) [1,1, f ,1] T = (1/ n )1 n . Intuitively, this means that filtering by W preserves the average
grey level, and will leave a constant signal (i.e., a “flat”
image) unchanged. In fact, with its spectrum inside the unit
disk, W is an ergodic matrix [40], [41], which is to say thatits powers converge to a matrix of rank one, with identicalrows. More explicitly, using the eigen-decomposition of W we
can writen

W k = VS k V -1 = VS k U T = / m ik v i u Ti ,i=1

“Symmetric Approximation of W and its Properties” and citations therein, where we illustrate and justify both the approximation procedure and its fidelity in detail. The obtained
matrix is a close approximation of a given W which is bothSPD, and now doubly (i.e., row- and column-) stochastic.
Working with a symmetric (or rather “symmetrized”) W
makes possible much of the analysis that follows in theremainder of this article. Interestingly, there is an inherent
and practical advantage in using symmetrization. As we will
see in the experimental results, given any W, employing its
symmetrized version generally improves performance in the
MSE sense. We do not present a proof of this observation here,
but refer the interested reader for details in [75].
Therefore, from this point forward, we shall consider W to be
an SPD, and (doubly) stochastic matrix. The spectrum of W
determines the effect of the filter on the noisy signal y. We write
the eigen-decomposition

where we have defined U T = V -1, while denoting by u Ti the left
eigenvectors of W which are the columns of U. Therefore,lim W k = 1 n u T1 .

k"3

So u 1 summarizes the asymptotic effect of applying the filter W
many times.
As we made clear earlier, W = D -1 K is generically not asymmetric matrix, though it always has real, positive eigenvalues. As such, it should not come as a surprise that W is quite
close to an SPD matrix. As we illustrate in “Symmetric
Approximation of W and its Properties,” this is indeed thecase. As such, we shall find it useful to approximate W withsuch an SPD matrix. To effect this approximation requires a
bit of care, as the resulting matrix must still be row-stochastic.
Fortunately, this is not difficult to do. Sinkhorn’s algorithm[44], [45] for matrix scaling provides a provably optimal [46]
way to accomplish this task. We refer the reader to [75] andW = VSV T,(16)

where S = diag6m 1, f , m n@ contains the eigenvalues in decreasing order 0 # m n # g # m 1 = 1, and V is an orthogonal matrix.
As an example, we illustrate the spectrum of the LARK [9]KijK
W = D-1 K

ji

L = D1/2 WD-1/2 - I[Fig2] Construction of the graph Laplacian L from the kernel
weights K.

IEEE SIGNAL PROCESSING MAGAZINE [112] january 2013

Symmetric Approximation of W and its Properties
As elaborated in detail in [75], we approximate approximate
the matrix W = D -1 K with a doubly stochastic (symmetric) positive definite matrix. The algorithm we use to effect thisapproximation is due to Sinkhorn [44], [45]. He proved thatgiven an irreducible matrix with nonnegative elements (such asW), there exist diagonal matrices R = diag(r) and C = diag(c)such that

perturbation of A that is Hermitian with eigenvalues
W
m2 $ g $ W
m n, then
m1 $ W2
A F.
/ m i - Wm i 2 # 2 A - X
ni =1

W and normalSpecializing the result for W, defining T = W - Z
izing by n, we haveZ
W = R WC

1 / m -W
1 T 2.
mi 2 #
i
F
n i=1
n
n

is doubly stochastic. That is,Z
W1 n = 1 n and 1 Tn Z
W = 1 Tn.

(S7)

Furthermore, the vectors r and c are unique to within a scalar(i.e., a r, c/a.) Sinkhorn’s algorithm for obtaining r and c ineffect involves repeated normalization of the rows and columns (see Algorithm 1 below for details) so that they sum toone, and is provably convergent and optimal in the crossentropy sense [46]. To see that the resulting matrix Z
W is SPD,
we note a corollary of Sinkhorn’s result: When a symmetricmatrix is considered, the diagonal scalings are in fact identical[44], [45]. Applied to K, we see that the symmetric diagonalscaling W
K = MKM yields a symmetric doubly stochastic matrix.
But we have W = D -1 K, so

The right-hand side can be further bounded as described in [75].
Hence, for patches of reasonable size (say n $ 11), the boundon the right-hand side and the eigenvalues of the respective
matrices are quite close, as demonstrated in Figure S1. Bounds
may be established on the distance between the eigenvectors of
W and ZW as well [75]. Of course, since both matrices are row-stochastic, they share their first right eigenvector vt 1 = v 1, corresponding to a unit eigenvalue. The second eigenvectors arearguably more important, as they tell us about the dominant
structure of the region being filtered, and the corresponding
effect of the filters. Specifically, as indicated in [47], applying Theorem 5.2.8 from [43] gives a bound on the perturbation of thesecond eigenvectors

Z
W = RD -1 KC.

vt 2 - v 2 #Since the Z
W and the diagonal scalings are unique (to within ascalar), we must have RD -1 = C = M, and therefore Z
W=W
K. Thatis to say, applying Sinkhorn’s algorithm to either K or its row-sum
normalized version W yields the very same SPD doubly stochastic)
result, Z
W.
Algorithm 1: Scaling a matrix A to a nearby doubly stochastic
matrix X
A
Given a matrix A, let (n, n) = size(A) and initializer = ones (n, 1); for k = 1: iter;c = 1./ (A T r);
r = 1./ (A c);
end

4T F
o- 2 T F,

where o is the gap between the second and the third eigenvalues. It is worth noting that the bound essentially does notdepend on the dimension n of the data y.
As a side note, we mention the interesting fact that, whenapplied to a symmetric squared distance matrix, Sinkhorn’s algorithm has a rather nice geometric interpretation [76], [77]. Specifically, if P is the set of points that generated the symmetricdistance matrix A, then elements of its Sinkhorn scaled version
B = MAM correspond to the square distances between the corresponding points in a set Q, where the points in Q are obtainedby a stereographic projection of the points in P. The points in Q
are confined to a hypersphere of dimension d embedded in R d +1,
where d is the dimension of the subspace spanned by the points
in P.

C = diag(c); R = diag(r);X
A=R A C

1.2
Symmetric Filter

0.8
0.6
0.4

0.6
0.4
0.2
0

0

-0.2

100
200
300
Eigenvalue Index
(a)

400

Asymmetric Filter

0.8

0.2
0

Symmetric Filter

1

Asymmetric Filter
Eigenvalue

The convergence of this iterative algorithm is known to be linear [45], with
rate given by the subdominant eigenvalue m 2 of W. It is of interest to know
how much the diagonal scalings will
perturb the eigenvalues of W. This
will indicate how accurate our approximate analysis will be, which uses theW instead of W. Fortueigenvalues of Z
nately, Kahan [74], [43] provides a niceresult: Suppose that A is a non-Hermitian matrix with eigenvalues
A be a| m 1 | $ | m 2 | $ g $ | m n |. Let X

Eigenvalue

1

0

100
200
300
Eigenvalue Index
(b)

400

[Figs1] Eigenvalues of row-stochastic (asymmetric) filter matrices and the corresponding
doubly stochastic (symmetric) filter matrices symmetrized using Sinkhorn’s algorithm:(a) an NLM filter example and (b) a LARK filter example.

IEEE SIGNAL PROCESSING MAGAZINE [113] january 2013

Eigenvalue

Statistical Analysis of theFilters’ Performance1
We are now in a position to carry out ananalysis of the performance of filters
0.9
defined earlier. As we just explained, to
0.8
Edge 2carry out this analysis, we first note
Texture 1
that W is, strictly speaking, a function
0.7of the noisy data y. Indeed, the depen0.6
Corner 1
dence on the given data is what gives
Texture3
Flat 1
strength to these adaptive filters. On
0.5
the other hand, one may legitimately
0.4
worry that the effect of noise on theCorner 2
computation of these weights may be
0.3
Texture 2
Flat 2dramatic, resulting in too much sensi0.2
tivity for the resulting filters to be effec0.1
tive. This issue has indeed been
considered in the context of most
Edge 1
0
patch-based adaptive filters. The con50
100 150 200 250 300 350 400
sensus approach, in keeping with simiEigenvalue Index
lar thinking in the nonparametricstatistics literature, is to compute the[Fig3] Spectrum of the LARK filter on different types of patches.
weights from a “pilot” (or prefiltered)version zt 0 of the noisy data; that is,
W (zt 0) . The preprocessing step yielding
t
the
pilot
estimate
kernel in Figure 3 for different types of patches. As can be
z 0 can be carried out with any one of manyobserved, the spectrum of W decays rather quickly for “flat”
simpler (nonadaptive) filters, and is done only for the purpose
patches, indicating that the filter is very aggressive in denoisingof stably computing the parameters of the filter. The net effectthese types of patches. As we consider increasingly more comof this step, however, is significant: it reduces the effective
plex anisotropic patches containing edges, corners, and various
variance of the input noise and stabilizes the calculation of the
textures, the spectrum of the filter is (automatically) adjusted toweights so that they are not a strong function of the input
be less aggressive in particular directions. To take the specificnoise level any longer. In this setting, the noise affecting thecase of a patch containing a simple edge, the filter is adjusted in
calculation of the weights can, for all practical intents, beunsupervised fashion to perform strong denoising along theassumed small. As such, the stability results in [47] (summaedge, and to preserve and enhance the structure of the patch
rized in “Stability of Filter Coefficients to Perturbations Due
across the edge. This type of behavior is typical of the adaptive
to Noise”) give strong guarantees that the weights are consenonparametric filters we have discussed so far, but it is particuquently disturbed relatively little. That is, W (zt 0) . W (z) . For
larly striking and stable in the case of LARK kernels shown in
the purposes of performance analysis below, it is therefore reaFigure 3.
sonable (and practically verified in Figures 4–7) to consider
Of course, the matrix W is a function of the noisy data y andthe weights as depending more directly on the underlying
hence it is strictly speaking a stochastic variable. But another(unknown) image W (z), rather than W (y) . As such, we alsointeresting property of W is that when the similarity kernel K ($)
drop the dependence on its argument for the sake of simplifying the notation. But the reader is advised to always remember
is of the general Gaussian type as in (7), as described in “Symthat the weights are adapted to the (estimated) structure of
metric Approximation of W and its Properties” the resulting filthe underlying image.
ter coefficients W are quite stable to perturbations of the data byNow, let us compute an approximation to the bias and varinoise of modest size [47]. From a practical point of view, andance of the estimator zt = Wy. The bias in the estimate is
insofar as the computation of the matrix W is concerned, it isalways reasonable to assume that the noise variance is relatively
bias = E (zt ) - z = E (Wy) - z
small, because in practice we typically compute W on a “prefil. Wz - z = (W - I)z.
tered” version of the noisy image y anyway. This small noise
assumption is made only for the analysis of the filter coeffiRecalling that the matrix W has a unit eigenvalue correspondcients, and is not invoked in rest of the article. Going forward,
ing to a (constant) vector, we note that zt is an unbiased estimate
therefore, we consider W as essentially nonstochastic. Someif z is a constant image, but is biased for all other underlying
experiments in the section “Improving the Estimate by Iteralatent images. The squared magnitude of the bias is given bytion” and in [75] confirm the validity of this assumption. We
also elaborate a bit more on this point in the next section.
bias 2 = (W - I) z 2 .
(17)IEEE SIGNAL PROCESSING MAGAZINE [114] january 2013

80
70

6

60

4
6
8
Iteration Number
(a)

4
3

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)
2

MSE

5
MSE

MSE

7
2.6
2.4
2.2
2
1.8
1.6
1.4
1.2
1

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

2
1

10

2

4
6
8
Iteration Number
(b)

50
40
30
20
10
0

10

2

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)
46810
Iteration Number(c)[Fig4] Parts (a)–(c) show plots of predicted MSE in diffusion based on ideally estimated symmetric LARK filters, and MSE of noisy
symmetric/asymmetric LARK filters through Monte-Carlo simulations. The input noise variance v 2 = 0.5, and for each simulation 100
independent noise realizations are implemented.

Writing the image z in the column space of the orthogonal
matrix V (which contains the eigenvectors as its columns) as
z = Vb, we can rewrite the squared bias magnitude asbias 2 = (W - I) z 2 = V (S - I) b 2n= (S - I) b 2 = / (m i - 1) 2 b i2 .

(18)

i =1

The ideal spectrum: BM3D
Let us pause for a moment and imagine that we can optimallydesign the matrix W from scratch in such a way as to minimize
the MSE expression above. That is, consider the set of eigenvalues for which (19) is minimized. Differentiating the expression
for MSE with respect to m i and setting equal to zero leads to a
familiar expression

We note that the last sum can be expressed over the indicesi = 2, f, n, because the first term in fact vanishes since m 1 = 1.
Next, we consider the variance of the estimate. We havecov (zt ) = cov (W y) . cov (W e) = v 2 W W T (nvar (zt ) = tr (cov (zt )) = v 2 / m i2 .
i=1Overall, the MSE is given byn

MSE = bias 2 + var (zt ) = / (m i - 1) 2 b 2i + v 2 m 2i .

(19)

m *i =b 2i
1=.
1 + snr -i 1
b 2i + v 2(20)This is, of course, 2nothing but a manifestation of the Wienerbfilter with snri = vi2 denoting the signal-to-noise ratio at eachcomponent i of the signal. This observation naturally raises
the question of whether any existing filters in fact get close tothis optimal performance [49]. Two filters that are designed
to approximately achieve this goal are BM3D [48] and patchwise locally optimal Wiener (PLOW) [50], which are currently
considered to be the state of the art in denoising. The BM3D

i =1

6.5

35

16

6

30

14

25

12

5
4.5

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

4
3.5

2

4
6
8
Iteration Number
(a)

10

20

MSE

MSE

MSE

5.5

15
10

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

5
0

2

4
6
8
Iteration Number
(b)

10

10
8
6
4
2
0

2

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)
46810
Iteration Number(c)[Fig5] Parts (a)–(c) show plots of predicted MSE in diffusion based on ideally estimated symmetric NLM filters, and MSE of noisy
symmetric/asymmetric NLM filters through Monte-Carlo simulations. The input noise variance v 2 = 0.5 and for each simulation 100
independent noise realizations are implemented.

IEEE SIGNAL PROCESSING MAGAZINE [115] january 2013

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

1.1

MSE

MSE

1
0.9
0.8
0.7
0.6
0.5

0

2

4
6
8
Iteration Number

10

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)
MSE

1.2

0

2

(a)

4
6
8
Iteration Number
(b)

10

4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

0

2

4
6
8
Iteration Number(c)10[Fig6] Parts (a)–(c) show plots of predicted MSE in twicing based on ideally estimated symmetric LARK filters, and MSE of noisy
symmetric/asymmetric LARK filters through Monte-Carlo simulations. The input noise variance v 2 = 0.5 and for each simulation 100
independent noise realizations are implemented.

patches for the aggregation process. In practice, of course,
algorithm can be briefly summarized by the following three
one does not have access to {b 2i }, so they are estimatedsteps:
1) Patches from an input image are classified according tofrom the given noisy image y, as is done also in [48]. The
their similarity, and are grouped into 3-D clusters
BM3D is not based on the design of a specific kernel K ($),
accordingly.
or even on a regression framework, as are other patch2) A so-called “3-D collaborative Wiener filtering” is implebased methods (such as NLM and LARK). In fact, there
mented to process each cluster.
need not exist any closed-form kernel that gives rise to the3) Filtered patches inside all the clusters are aggregated to
BM3D filter. Still, with the present framework, the diagoform an output.
nal Wiener shrinkage matrix
The core process of the BM3D
S = diag[m 1, f , m n] a n d t h e
One way to remedy the lackalgorithm is the collaborative
overall 3-D collaborative Wieof optimality of the choice ofWiener filtering that transner filtering process can now
kernel is to apply the resulting
forms a patch cluster through abe understood using a symmetfilters iteratively.
fixed orthonormal basis V (e.g.,ric, positive-definite filter
DCT) [48], where the coeffimatrix W with orthonormal
cients of each component are then scaled by some “shrinkeigen-decomposition VSV T . If the orthonormal basis V conage” factor {m i}, which are precisely chosen according totains a constant vector v il = 1n 1 n, one can easily make W(20). Following this shrinkage step, the data are transdoubly stochastic by setting its corresponding shrinkage
formed back to the spatial domain to generate denoised
factor m il = 1.

4.5

3

5
4
MSE

MSE

3.5
2.5
2
1.5
1
0

2

4

6

Iteration Number(a)8

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

10

3

1.4
1.2
1

2

0.8

1

0.6

0

0

2

4
6
8
Iteration Number(b)

Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

1.6

MSE

4

0.5

1.8

6
Ideally Predicted
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

10

0.4

0

2

4

6

8

10

Iteration Number(c)[Fig7] Parts (a)–(c) show plots of predicted MSE in diffusion based on ideally estimated symmetric NLM filters, and MSE of noisy
symmetric/asymmetric NLM filters through Monte-Carlo simulations. The input noise variance v 2 = 0.5 and for each simulation 100
independent noise realizations are implemented.

IEEE SIGNAL PROCESSING MAGAZINE [116] january 2013

To summarize this section, we have presented a widely
applicable matrix formulation zt = W y for the denoising problem, where W is henceforth considered symmetric, positivedefinite, and doubly stochastic. This formulation has thus far
allowed us to characterize the performance of the resulting filters in terms of their spectra and to obtain insights as to therelative statistical efficiency of existing filters such as BM3D.
One final observation worth making here is that the Birkhoff–
von Neumann theorem [41] notes that W is doubly stochastic
if and only if it is a convex combination of permutation matrin
ces. Specifically, W = / l = 1 a l Pl, where Pl are permutation
matrices; a l are nonnegative scalars with / a l = 1, and nr is at
most (n - 1) 2 + 1. This means that regardless of the choice of(an admissible) kernel, each pixel of the estimated image
zt = W y = / a l ^Pl y hn[51], whereas the second class includes the recently popularized Bregman iterations [52], which is closely related to , 2-boosting [53]. The latter is a generalization of twicing introduced by Tukey [54] more than 30 years ago. In fact, as we
will illustrate later, this latter class of iterations is also related
to biased anisotropic diffusion [55], otherwise known as
anisotropic reaction-diffusion.
Diffusion
So far, we have described a general formulation for denoising as
a spatially adaptive, data-dependent filtering procedure (3)amounting to
zt = W y.

-l=1

is always a convex combination of at most n= (n - 1) 2 + 1 (notnecessarily nearby) pixels of the noisy input.
Going forward, we first use the insights gained thus far to
study further improvements of the nonparametric approach.
Then, we describe the relationship between the nonparametricapproaches and more familiar parametric (empirical) Bayesian
methods.

Now consider applying the filter multiple times. That is, define
zt 0 = y, and the iteration
zt k = Wzt k - 1 = W k y.

The net effect of each application of W is essentially a step of
anisotropic diffusion [51], [5]. Rewriting the iteration in (21),
we have
zt k = Wzt k - 1,

Improving the Estimate by Iteration
As we described above, the optimal spectrum for the filter
matrix W is dictated by the Wiener condition (20), which
requires clairvoyant knowledge of the signal or at least careful
estimation of the SNR at each spatial frequency. In practice,
however, nearly all similarity-based methods we described so far(with the notable exception of the BM3D) design a matrix W
based on a choice of the kernel function, and without regard to
the ultimate consequences of this choice on the spectrum of theresulting filter. Hence, at least in the MSE sense, such choicesof W are always deficient. As former U.S. Secretary of Defense
Donald Rumsfeld might have said, we filter with the spectrum
that we have, not the spectrum that we would like to have! One
way to remedy the lack of optimality of the choice of kernel is to
apply the resulting filters iteratively, and that is the subject of
this section.
To put the problem in more practical terms, the estimate
derived from the application of a nonparametric denoisingfilter W will not only remove some noise, but invariably someof the underlying signal as well. A number of different
approaches have been proposed to reintroduce this “lost”
component back to the estimate. Interestingly, as we describe
later in the section “Relationship to (Empirical) Bayes’ andRegularization Approaches,” similar to the steepest descent(SD) methods employed in the solution of more classical
Bayesian optimization approaches, the iterative method
employed here involve multiple applications of a filter (orsequence of filters) on the data, and somehow pooling theresults. Next we discuss the two most well-known and frequently used approaches. The first is anisotropic diffusion

(21)

(22)

= zt k - 1 - zt k - 1 + Wzt k - 1,

(23)

= zt k - 1 + ^ W - I hzt k - 1,

(24)

zt k - zt k - 1 = ^ W - I h zt k - 1 .

(25)

which gives

Recall that W = D -1/2 LD 1/2 . Hence, we have
W - I = D -1/2 ^L - I hD 1/2 = D -1/2 L D 1/2,

(26)

where L is the graph Laplacian operator mentioned earlier [39],[5]. Defining the normalized variable zr k = D 1/2 zt k, (25) embodies
a discrete version of anisotropic diffusion
zr k - zr k -1 = L zr k -1, )2zr (t)= d 2 zr (t) ( Diffusion Eqn) . (27)2twhere the left-hand side of the above is a discretization of the
2rz (t)derivative operator 2t .
Returning to the diffusion estimate (21), the bias in the estimate after k iterations isbias k = E (zt k) - z = E (W k y) - z = W k z - z = (W k - I) z.
Recalling that the matrix W has a unit eigenvalue corresponding to a (constant) vector, we note that since W k - I has thesame eigenvectors as W, the above sequence of estimators produce unbiased estimates of constant images.

IEEE SIGNAL PROCESSING MAGAZINE [117] january 2013

(a)

(b)

(c)

[Fig8] Example patches: (a) flat, (b) edge, and (c) texture.

As k grows, the variance tends to a constant value of v 2 . Overall,
the MSE is given by

The squared magnitude of the bias is given by2

bias k 2 = (W k - I) z .

(28)n

MSE k = bias k 2 + var (zt k) = / (m ik - 1) 2 b i2 + v 2 m i2k .

As before, writing the image z in the column space of V asz = Vb, we havenbias k 2 = (W k - I) z = V (S k - I) b = / (m ik -1) 2 b 2i . (29)
2

2

i =1

As k grows, so does the magnitude of the bias, ultimately converging to b 2 - b 21 . This behavior is consistent with what isobserved in practice; specifically, increasing iterations of diffusion produce more and more blurry (biased) results. Next, wederive the variance of the diffusion estimatorcov (zt k) = cov (W k y) = cov (W k e) = v 2 W k (W k) T,
which givesnvar (zt k) = tr (cov (zt k)) = v 2 / m i2k .

(30)

i =1

Diffusion experiments for denoising three types of patches(shown in Figure 8) are carried out to illustrate the evolutionof MSE. The variance of input (Gaussian) noise is set tov 2 = 25. Filters based on both LARK [9] and NLM [7] estimated from the latent noise-free patches are tested, and the
Sinkhorn algorithm (see “Symmetric Approximation of W and
its Properties”) is implemented to make the filter matrices
doubly stochastic. Predicted MSE k, var (zt k), and bias k 2
according to (30) through the diffusion iteration are shown in
Figure 9 (LARK) and Figure 10 (NLM), where we can see that
diffusion monotonically reduces the estimation variance andincreases the bias. In some cases [such as Figure 9(a) and (b)]diffusion further suppresses MSE and improves the estimation
performance. True MSEs for the standard (asymmetric) filters

i=1

MSE
Variance
Bias2

5
4

350
10

250

8

3

200

6

2

4

1

2
2

4
6
8
Iteration Number
(a)

10

MSE
Variance
Bias2

300

MSE
Variance
Bias2

150
100
50

2

4

6

8

10

0

Iteration Number
(b)

2

4

6

Iteration Number(c)[Fig9] Parts (a)–(c) show plots of predicted MSE k, var (zt k) and bias k 2 using (30) for LARK [9] filters in diffusion process.

IEEE SIGNAL PROCESSING MAGAZINE [118] january 2013

8

10

7

30

6

25

5

20

4

10

15

3

MSE
Variance
Bias2

2
1
0

15

2

4
6
8
Iteration Number(a)

MSE
Variance
Bias2

10
5

10

0

2

4
6
8
Iteration Number

MSE
Variance
Bias2

5

0

10

2

4
6
8
Iteration Number

(b)

10

(c)

[Fig10] Parts (a)–(c) show plots of predicted MSE k, var (zt k) and bias k 2 using (30) for NLM [7] filters in diffusion process.

and their symmetrized versions computed by Monte-Carlosimulations are also shown, where in each simulation 100
independent noise realizations are averaged. We see that theestimated MSEs of the symmetrized filters match quite wellwith the predicted ones (see Figures 11 and 12). The asymmetric filters, meanwhile generate slightly higher MSE, butbehave closely to the symmetrized ones especially in regions
around the optimal MSEs.
In the next set of Monte-Carlo simulations, we illustrate
the effect of W being computed from clean versus noisy
images (see Figures 4 and 5). Noise variance v 2 = 0.5, whichis relatively small, simulating the situation where “preprocessing” has been applied to suppress estimation variance. It
can be observed that the MSEs estimated from Monte-Carlosimulations are quite close to the ones predicted from theideal filters, which confirms the assumption in the section“The Martix Formulation and its Properties” that the filter
matrix W can be treated as deterministic under mostcircumstances.

5

10

4

(i)

(i)

MSE k = (m ik - 1) 2 b 2i + v 2 m 2i k .

468
Iteration Number
(a)

10

(i)

(i)

MSE k =
b

MSE k
v

2

= snri (m ik - 1) 2 + m 2i k,

2

350
300
250
200
150
100

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

50
2

(33)

where snri = vi2 is defined as the signal to noise ratio in the ith
mode.
(i)
One may be interested to know at which iteration MSE k isminimized. As we have seen experimentally, with an arbitrary

9.5

8

(32)

Equivalently, we can write

8.5
2

(31)

where MSE k denotes the MSE of the ith mode in the kth diffusion iteration, which is given by

9

3.5

(i)

i =1

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

10.5
MSE

MSE

4.5

3

11

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)n

MSE k = / MSE k ,

MSE

5.5

To further analyze the change of MSE in the diffusion process, let us consider the contribution of MSE in each component (or mode) separately

4
6
8
Iteration Number
(b)

10

0

2

4
6
8
Iteration Number(c)10[Fig11] Parts (a)–(c) show plots of predicted MSE and Monte-Carlo estimated MSE in diffusion process using LARK [9] filters. In theMonte-Carlo simulations, results for row-stochastic (asymmetric) filters and their symmetrized versions are both shown.

IEEE SIGNAL PROCESSING MAGAZINE [119] january 2013

6.5

35

6

30

16

25

14

5
4.5

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

4
3.5

2

4
6
8
Iteration Number
(a)

10

20

MSE

MSE

MSE

5.5

15
10

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

5
0

2

4
6
8
Iteration Number(b)Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

12
10
8
2

10

4
6
8
Iteration Number(c)10[Fig12] Parts (a)–(c) show plots of predicted MSE and Monte-Carlo estimated MSE in diffusion process using NLM [7] filters. In theMonte-Carlo simulations, result for row-stochastic (asymmetric) filters and their symmetrized versions are both shown.

filter W, there is no guarantee that even the very first iterationof diffusion will improve the estimate in the MSE sense. The
derivative of (33) with respect to k is(i)

2 MSE k= 2m ik log m i 6(snri + 1) m ik - snri@ .
2k

(34)

For i $ 2, since m i < 1, we note that the sign of the derivative
depends on the term inside the brackets, specifically,(snri + 1) m ik - snri . To guarantee that in the kth iteration the
MSE of the ith mode decreases (i.e., negative derivative), wemust have(snri + 1)qmi- snri 2 0 for any 0 # q 1 k.

(35)

In fact, the condition(snri + 1) m ik - snri 2 0(36)

guarantees that the derivative is always negative for q ! (0, k)qbecause given any scalar t = k ! (0, 1)(snri + 1) m ik - snri 2 0 & m tk
i 2`(snri + 1) m tk
i 2 snri `snri tj &snri + 1

snri t - 1
j 2 snri
snri + 1.

(37)

image produced by the ith mode of the denoising filter. The
above condition implies that if the filter entropy exceeds orequals the channel capacity of the ith mode, then the kth iteration of diffusion will in fact produce an improvement in the corresponding mode. One may also be interested in identifying anapproximate value of k for which the overall MSE k is minimized. This depends on the signal energy distribution ({b 2i })over all the modes, which is generally not known, but for which
we may have some prior in mind. But for a given mode i, theminimum MSE (i) is achieved in the k*th iteration, wherek * = log `

(38)

(39)

Of course, in practice we can seldom find a single iteration
number k* minimizing MSE (i) for all the modes simultaneously; but in general, the take-away lesson is that optimizing agiven filter through diffusion can [under proper conditions(38)] make its corresponding eigenvalues closer to those of an
ideal Wiener filter that minimizes MSE.
Another interesting question is this: given a particular W,
how much further MSE reduction can be achieved by implementing diffusion? For example, can we use diffusion to
improve the 3-D collaborative Wiener filter in BM3D [48]? Let
us assume that this filter has already approximately achieved
the ideal Wiener filter condition in each mode, specifically

The condition (36) has an interesting interpretation. Rewriting(36), we havelog (1 + snri) 1 log c 1 m,
1 - mli
144 2
443
ei

snri
j log m i .
snri + 1

m *i .

snri.
snri + 1(40)

Then we can see that(snri + 1) (m *i ) k - snri 1 0 for any k 2 1,

(41)

k

where mli = m i denotes the ith eigenvalue of W k . The left-handside of the inequality is Shannon’s channel capacity [56] of theith mode for the problem y = z + e. We name the right-handside expression e i the entropy of the ith mode of the filter W k.
The larger e i, the more variability can be expected in the output

which means for all the modes, diffusion will definitely worsen(increase) the MSE. In other words, multi-iteration diffusion is
not useful for filers where the Wiener condition has been
(approximately) achieved.

IEEE SIGNAL PROCESSING MAGAZINE [120] january 2013

Extensions to More General Restoration Problems
Many of the ideas discussed in the context of denoising here canbe extended to more general image restoration problems. As asimple example, we briefly note the case where the measurement model includes a linear distortion such as blur, expressed asfollows:
y = Az + e.

(S8)

One approach would be to employ the framework of the residual iterations defined earlier. More specifically, define
zt k + 1 = zt k + A T W (y - Azt k),(S9)where we observe that the term W (y - Azt k) is filtering (i.e.,denoising) the residuals from the measurement model in a spatially adaptive, nonparametric manner using W. The operator A T
is the adjoint to the blur, and applied to these filtered residuals,

Twicing, , 2-Boosting, Reaction-Diffusion,
and Bregman Iterations
An alternative to repeated applications of the filter W is to consider the residual signals, defined as the difference between theestimated signal and the measured signal. The use of the residuals in improving estimates has a rather long history, dating atleast back to the work of John Tukey, who, in his landmark book[54], termed the idea twicing. More recently, the idea has been
proposed in the applied mathematics community for application
to image restoration problems under the rubric of Bregman
iterations [52]. And in the machine learning and statistics literature, [53] recently proposed the idea of , 2-boosting, which may
be viewed as an ( , 2 ) regression counterpart of the more familiar
notion of boosting in classification, originally proposed byFreund and Schapire [57]. All of the aforementioned ways oflooking at the use of residuals are independently quite rich andworth exploring on their own. Indeed, for a wide-angle view ofthe boosting approach to regression and regularization, and its
connections to functional gradient descent interpretations, we
can refer the interested reader to the nice treatment in [58].
The general idea of Bregman divergence, and its use in solving
large-scale optimization problems is also a very deep area ofwork, for which we can refer the interested reader to the book[59] in which both the foundations and some applications tosignal processing are discussed. Also, in [60], the ideas ofBregram divergence and boosting are connected. More general
inverse problems such as deconvolution can also be addressedin this framework. We give a very brief hint of this direction in“Extensions to More General Restoration Problems.”
In whatever guise, the basic notion—to paraphrase Tukey—
is to use the filtered residuals to enhance the estimate by adding
some “roughness” to it. Put another way, if the residuals contain some of the underlying signal, filtering them should
recover this leftover signal at least in part.
Formally, the residuals are defined as the difference between
the estimated signal and the measured signal: rk = y - zt k -1,
where here we define the initialization zt 0 = Wy. Note that thisacts as a “backprojection” operation. The asymptotic behavior ofthis iterative process is interesting. Specifically, as k " 3,
A T Wy = A T WAzt 3

(S10)

-1
zt 3 = ^A T WA h A T Wy.

(S11)or equivalently

Without the adaptive weights in W, this asymptotic solution is thefamiliar (global) least squares solution, which is unstable due to thedirect inverse. The inclusion of the kernel weights, and stopping theprocess after a finite number of iterations can result in a stabilized
solution with smaller MSE. The analysis of the stability and statistical characteristics of such estimates, along with the choice of the
optimal stopping point will require in depth understanding of thespectrum of the operator A T W, much in the same way as we didwith the spectrum of W in the denoising framework.

initialization is a once-filtered version of the noisy data using
some initial filter W, and different from (z 0 = y ) used in the diffusion iterations. With this definition, we write the iterated
estimates as
zt k = zt k - 1 + W rk = zt k - 1 + W (y - zt k - 1) .

(42)

These estimates also trade off bias against variance with increasing iterations, though in a fundamentally different way than the
diffusion estimator we discussed earlier. Indeed, we immediately
notice that this iteration has a very different asymptotic behavior than diffusion. Specifically, as k " 3, the estimate here
tends back toward the original data y, whereas the diffusion
tends to a constant!
For an intuitive interpretation, note that for k = 1, the(residual) iteration above gives
zt 1 = zt 0 + W (y - zt 0) = Wy + W (y - Wy) = ^2W - W 2 h y.
This first iterate zt 1 is in fact precisely the original “twicing” estimate of Tukey [54]. More recently, this particular filter has beensuggested (i.e., rediscovered) in other contexts. In [39], forinstance, Coifman et al. suggested it as an ad hoc way to enhancethe effectiveness of diffusion. The intuitive justification given [5]was that ^2W - W 2 h = ^2I - W h W can be considered as a twostep process; specifically blurring, or diffusion (as embodied by
W), followed by an additional step of inverse diffusion or sharpening (as embodied by 2I - W.) As can be seen, the ad hoc suggestion has a clear interpretation here. Furthermore, we see thatthe estimate zt 1 can also be thought of as a kind of nonlinear,
adaptive unsharp masking process, further clarifying its effect.
In [61], this procedure and its relationship to the Bregman iterations were studied in detail. For the sake of completeness, we
also mention in passing that in the nonstochastic setting andwhere W is rank one, the residual-based iterations (42) are
known as matching pursuit [62]. This approach has found many
applications and extensions for fitting overcomplete dictionaries.

IEEE SIGNAL PROCESSING MAGAZINE [121] january 2013

behavior of the estimates, we rewrite the iterative process in(42) more explicitly in terms of the data y. We have [40] and [53]To further understand the iteration (42) in comparison to
the earlier diffusion framework, consider again a very similar
algebraic manipulation of (42)zt k = zt k - 1 + W (y - zt k - 1),(43)= ^I - W h zt k - 1 + Wy,(44)

zt k - zt k - 1 = ^I - W hzt k - 1 + (zt 0 - zt k - 1) .

(45)

zt k = / W (I - W) j y = ^I - ^I - W hk + 1 h y.
k

Clearly the above iteration does not monotonically blur thedata; but a rather more interesting behavior for the biasvariance tradeoff emerges. Analogous to our earlier derivation,
the bias in the estimate after k iterations isbias k = E (zt k) - z = ^I - ^I - W hk + 1 h E (y) - z = - ^I - W hk + 1 z.

Redefining the graph Laplacian [note that this definition of thegraph Laplacian differs by a negative sign from the earlier onein (26)]The squared magnitude of the bias is given bybias k 2 = ^I - W hk +1 z

I - W = D -1/2 ^I - L hD 1/2 = D -1/2 L D 1/2,
with a change of variable zr k = D 1/2 zt k, we can write

2zr (t)= d 2 zr (t) + (zr (0) - zr (t)) .
t 2 4 44 3 1 4 4 2 4 4 3
1 4244
reaction term(46)

var (zt k) = tr (cov (zt k)) = v 2 / ^1 - (1 - m i) k + 1 h .
n

The above is a modified version of the diffusion equation studied
earlier, where the modification takes the form of a “forcing”
function, which is in fact the residual between the estimate attime t and a prefiltered version of the original data. This additional term reacts against the strict smoothing effect of thestandard diffusion, and results in a nontrivial steady-state. Not
coincidentally, the notion of a “biased” anisotropic diffusion [55]is essentially identical to the above formulation. More recently,
Farbman et al. [63] proposed a variational formulation which
also leads precisely to the present setup. This will become moreclear as we discuss the connection to (empirical) Bayesian
approaches in the section “Relationship to (Empirical) Bayes’
and Regularization Approaches.”
The residual-based iterations also have fundamentally different statistical properties than diffusion. To study the statistical

As k grows, the variance tends to a constant value of v 2. Overall,
the MSE is
MSE k = bias k 2 + var (zt k)= / (1 - m i) 2k +2 b 2i + v 2 ^1 - (1 - m i) k +1h .
nExperiments for denoising the patches given in Figure 8 using
the “twicing” approach are carried out, where the same doubly
stochastic LARK and NLM filters used for the diffusion tests in
Figures 9 and 10 are employed. Plots of predicted MSE k, var (zt k)
and bias k 2 with respect to iteration are illustrated in
Figure 13 (LARK) and Figure 14 (NLM). It can be observed that

35
30
25

15

20

10

15

10

0

0

2

4
6
8
Iteration Number
(a)

10

MSE
Variance
Bias2

5
0

(49)

2

i =1

20

5

2

i=1

25

MSE
Variance
Bias2

(48)

i =1

which gives

diffusion

15n= / (1 - m i) 2k +2 b i2 .

cov (zt k) = cov 6^I - ^I - W hk + 1h y@T= v 2 ^I - ^I - W hk + 1h^I - ^I - W hk + 1h ,reaction term

1

20

2

The behavior of the (decaying) bias in this setting is in stark
contrast to the (increasing) bias of the diffusion process.
The variance of the estimator and the corresponding MSE are

zr k - zr k -1 = L zr k -1 + (zr 0 - zr k -1),
1 4 4 44 2 4 4 44 3 1 44 2 44 3
diffusion

(47)

j=0

0

2

4
6
8
Iteration Number

10

MSE
Variance
Bias2

10
5
0

(b)

0

2

4
6
8
Iteration Number(c)[Fig13] Parts (a)–(c) show plots of predicted MSE k, var (zt k) and bias k 2 using (49) for LARK [9] filters in the twicing process.

IEEE SIGNAL PROCESSING MAGAZINE [122] january 2013

10

4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

MSE
Variance
Bias2

6

25

5

20

4

15

3

1
0

2

4
6
8
Iteration Number

10

0

10

MSE
Variance
Bias2

2

0

2

4

6

8

MSE
Variance
Bias2

5
0

10

0

2

4
6
8
Iteration Number

Iteration Number
(b)

(a)

10

(c)

[Fig14] Parts (a)–(c) show plots of predicted MSE k, var (zt k) and bias k 2 using (49) for NLM [7] filters in twicing process.

in contrast to diffusion, twicing monotonically reduces the estimation bias and increases the variance. So twicing may in fact
reduce the MSE in cases where diffusion fails to do so [suchas Figure 9(c)]. True MSEs, for the asymmetric filters and theirsymmetrized versions, are estimated through Monte-Carlo

20

24

40

22

35

20

10

0

2

4
6
8
Iteration Number

16
14

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

12
10

10

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

30

18

MSE

MSE

MSE

15

5

simulations, and they all match very well with the predicted
ones (see Figures 15 and 16).
Similar to the diffusion analysis before, we plot MSE resulting from the use of W filters directly estimated from noisy
patches through Monte-Carlo simulations with noise variance

0

2

(a)

4
6
8
Iteration Number
(b)

10

25
20
15
10

0

2

4
6
8
Iteration Number
(c)

10

6

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

25

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

5.5

20

5

MSE

4.5
4
3.5
3
2.5
2
1.5
1
0.5

MSE

MSE

[Fig15] Parts (a)–(c) show plots of predicted MSE and Monte-Carlo estimated MSE in twicing process using LARK [9] filters. In theMonte-Carlo simulations, row-stochastic (asymmetric) filters and their symmetrized versions are all tested.

4.5

15

Predicted (Sym.)
Monte-Carlo (Sym.)
Monte-Carlo (Asy.)

4
0

2

4
6
8
Iteration Number
(a)

10

3.5

0

2

4
6
8
Iteration Number
(b)

10

10

0

2

4
6
8
Iteration Number(c)[Fig16] Parts (a)–(c) show plots of predicted MSE and Monte-Carlo estimated MSE in twicing process using NLM [7] filters. In theMonte-Carlo simulations, row-stochastic (asymmetric) filters and their symmetrized versions are all tested.

IEEE SIGNAL PROCESSING MAGAZINE [123] january 201310

v 2 = 0.5, and compare them withthe predicted MSEs from the ideally estimated filters. Again, theMSEs estimated from MonteCarlo simulations are quite close
to the predicted data (see Figures
6 and 7).
Here again, the contribution to
MSE k of the ith mode can be written as

MSE k = (1 - m i) 2k +2 b 2i + v 2 ^1 - (1 - m i) k +1 h .
2

(i)

(50)

Proceeding in a fashion similar to the diffusion case, we can(i)analyze the derivative of MSE k with respect to k to see whether
the iterations improve the estimate. We have(i)

2 MSE k= 2 (1 - m i) k + 1 log (1 - m i)
2k
# 6(snri + 1) (1 - m i) k + 1 - 1@ .

(51)

Reasoning along parallel lines to the earlier discussion, the following condition guarantees that the kth iteration improves theestimate:(snri + 1) (1 - m i) k +1 2 1.

(52)

Rewriting the above, we havelog (1 + snri) 2 log cknowledge (or estimation) of thelatent signal energy {b 2i } (or theSNR,) and this is an ongoing challenge. Research on estimating the
SNR of image or video data without a reference “ground truth” isof broad interest in the community as of late [64]. In particular,
this problem points to potentially important connections withongoing work in no-reference image quality assessment [65]as well.

specifying the kernel functionor the corresponding weights
is essentially equivalent
to estimating a particular
type of empirical prior fromthe given data.

1 m,
1 - mli(53)

where now mli = 1 - (1 - m i) k + 1 is the ith eigenvalue ofI - (I - W) k +1, which is the equivalent filter matrix for the kth
twicing iteration. This is the reverse of the condition we derived
earlier; specifically, here we observe that if the entropy does notexceed the channel capacity, then iterating with the residuals
will indeed produce an improvement. This makes reasonable
sense because if the filter removes too much detail from thedenoised estimate, this “lost” detail is contained in the residuals, which the iterations will attempt to return to the estimate.
(i)
Again, the minimum MSE k is achieved at the k*th iteration,wherek * = - log ^1 + snri h / log (1 - mli ) - 1,(54)

which is when the i th eigenvalue of the matrix
I - (I - W) k * + 1 becomes snri / (snri + 1), the Wiener condition.
In most cases we cannot optimize all the modes with a uniform number of iterations, but under proper conditions (53)twicing can make the eigenvalues closer to those of the ideal
Wiener filter.
From the above analysis, we can see that it is possible toimprove the performance of many existing denoising algorithms in the MSE sense by implementing iterative filtering(see Figure 17). However, to choose either diffusion or twicing
and to determine the optimal iteration number require prior

Relationship to (Empirical) Bayes’
and Regularization Approaches
It is interesting to contrast the adaptive nonparametric framework described so far with the more classical Bayesian estimation approaches. Of course, Bayesian thinking requires that a“prior” be specified independently of the measured data, whichcharacterizes the aggregate behavior of the underlying latentimage of interest. In the context we have described so far, no
specific prior information was used. Instead, data-dependant
weights were specified in a least-squares framework, yielding
locally adaptive filters. On the surface then, it may seem that
there is no direct relationship between the two approaches.
However, as we will illustrate below, specifying the kernel function or the corresponding weights is essentially equivalent toestimating a particular type of “prior” from the given data. (Onecringes to call this estimate a “prior” any longer, but we will doso for lack of a better word.) This empirical Bayesian approach
has in fact a long and useful history [66] in the statistics literature. Furthermore, another class of modern patch-based filtering methods based on the empirical Bayesian framework wererecently pioneered in [67] and extended in [68] and [69] basedon the notion of “adaptive” (i.e., data-dependent) priors. In thesame spirit, we illustrate below that an empirical version of the
popular maximum a posteriori method has a data-dependent
counterpart in the nonparametric techniques we describedso far.
A canonical Bayesian approach is to specify a prior for the
unknown z, or to equivalently use a regularization functional. The so-called maximum a posteriori (MAP) estimate
is then given by the solution of the following optimization
problem:
MAP: zt = arg min 1 y - z 2 + m R (z),
z
2
2

(55)

where the first term on the right-hand side is the data-fidelity(log-likelihood) term, and the second (regularization, or logprior) term essentially enforces a soft constraint on the global
structure and smoothness of the signal. In the MAP approach,the regularization term is a functional R (z), which is typically(but not always) convex, to yield a unique minimum for theoverall cost. Particularly popular recent examples include
R (z) = dz , and R (z) = z 1 . The above approach implicitly
contains a global (Gaussian) model of the noise (captured by theIEEE SIGNAL PROCESSING MAGAZINE [124] january 2013

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

[Fig17] Denoising example using iterative LARK filtering: Parts (a)–(c) show input noisy patches; noise variance v 2 = 25. (d)–(f) show
output patches by filtering once. Parts (g)–(i) show MSE optimized outputs: part (g) is the sixth diffusion output; (h) is the fourth
diffusion output; and (i) is the third twicing output.

quadratic data-fidelity term) and an explicit model of the signal(captured by the regularization term.)
Regardless of what choice we make for the regularization
functional, this approach is global in the sense that the explicit
assumptions about the noise and the signal constrain the
degrees of freedom in the solution, hence limiting the global
behavior of the estimate. This often results in well-behaved (butnot always desirable) solutions. This is precisely the main

motivation for using (local/nonlocal) adaptive nonparametrictechniques in place of global regularization methods. Indeed,though the effect of regularization is not explicitly present
in the nonparametric framework, its work is implicitly done by
the design of the kernel function K ($), which affords us local
control, and therefore more freedom and often better adaptation to the given data, resulting in more powerful techniqueswith broader applicability.

IEEE SIGNAL PROCESSING MAGAZINE [125] january 2013

if and only if W is symmetric.) We observe therefore that the
Another important advantage of the local nonparaequivalent regularization term is a quadratic function of themetric methods is that the resulting algorithms inherently
data and residuals. In particular, the function R diffusion (zt k) can
operate in relatively small windows, and can hence be efficiently implemented and often easily parallelized. With global
be thought of as an estimated log-prior at each iteration as it
parametric and variational methods, the solutions are often
depends on the data y, the last estimate zt k, and weights W,
implemented using large-scale global optimization techwhich are computed from the data. More specifically, theniques, which can be hard (sometimes even impossible) to
implied empirical Bayes’ “prior” would be“kernelize,” for algorithmic efficiency [33].
Though we do not discuss them here, of course hybrid
pt (z) = c exp 6- R (zt )@,methods are also possible, where the first term in (55) is
replaced by the corresponding
where c is a normalization
weighted term (3), with regularconstant. The empirical (MAP)
Another important advantage
ization also applied explicitly
Bayesian interpretation of theof the local nonparametric[70]. As we shall see below, theresidual-based method follows
methods is that the resulting
iterative methods for improving
similarly, but is somewhat suralgorithms inherently operate in
the nonparametric estimatesprisingly more straightforward, as
relatively small windows.
(such as diffusion, twicing, etc.)
we illustrate below.
mirror the behavior of SD iterations in the empirical Bayesian case, but with a fundamental
Empirical Bayes’ Interpretationdifference. Specifically, the corresponding implied regularizaof the Residual Iterations
tion functionals deployed by the nonparametric methods are
We arrive at the empirical Bayes’ interpretation of the residual
significantly different, adapting directly to the given data.
iterations (42) in a similar manner. In particular, comparing
To connect MAP to the nonparametric setting, we proceed in(42) to the SD MAP iterations in (56) we have
similar lines as Elad [3] by considering the simplest iterative
approach: the SD method. The SD iteration for MAP, with fixed(61)
zt k + 1 = zt k - n 6(zt k - y) + mdR (zt k)@
step size µ is
(62)
zt k +1 = zt k + W (y - zt k) .
t
t
t
6
@
(56)
MAP:
z k + 1 = z k - n (z k - y) + mdR (z k) .
Again, equating the right-hand sides, we get
Empirical Bayes’ Interpretation of Diffusion
Let us compare the diffusion iterations (25) to the SD MAP iterations above in (56). Specifically, we equate the right-handsides ofzt k + 1 = zt k - n 6(zt k - y) + mdR (zt k)@(57)

zt k + 1 = zt k + ^ W - I h zt k,(58)

Solving for dR (zt k) we obtaindR (zt k) = 1 ^ W - nI h^ y - zt kh .
nm

(63)

With W symmetric, we have (to within an additive constant)which gives- n 6(zt k - y) + mdR (zt k)@ = ^ W - I hzt k .
Solving for dR (zt k) we obtaindR (zt k) = 1 ^ W - (1 - n) I h^ y - zt kh - 1 ^I - W h y . (59)nm

- n 6(zt k - y) + mdR (z k)@ = W (y - zt k) .

nm

When W is symmetric, integrating both sides we have (to within
an additive constant)
R diff. (z) = 1 ^ y - z hT ^(1 - n) I - W h^ y - z h + 1 y T ^I - W h z.
2nm
nm
		(60)(Since the left-hand side of (59) is the gradient of a scalar-valued
function, the right-hand side must be curl free. This is possible

R residual (z) = 1 ^ y - z hT ^ nI - W h^ y - z h .
2nm

(64)

Therefore, the (implicit) regularization term is an adaptivefunction of the data and residuals. Once again, the empiricallyestimated prior is of the same form pt (z) = c exp 6-R (zt )@, where
R (zt ) is computed directly on the residuals y - zt of the filtering
process at every iteration. It is especially interesting to note that
this empirical prior is a quadratic form in the residuals alone
where, not coincidentally, the matrix nI - W is closely related
to the Laplacian operator defined earlier. This hints directly atthe edge-preserving behavior of this type of iteration; a property
not shared by the diffusion process, which monotonically blurs
the data.
Conclusions
It has been said that in both literature and the arts, there are
only seven basic plots (comedy, tragedy, etc.)—that all stories

IEEE SIGNAL PROCESSING MAGAZINE [126] january 2013

AUthorare combinations or variations on
deep connections between
Peyman Milanfar (milanfar@ucsc.
these basic themes. Perhaps it istechniques used commonlyedu) received the B.S. degree intaking the analogy too far to say
in computer vision, image
electrical engineering and mathethat an exact parallel exists in our
processing, graphics, andmatics from the University of Califield as well. But it is fair to argue
machine learning exist.
fornia, Berkeley, and the M.S. andthat the basic tools of our trade in
Ph.D. degrees in electrical engirecent years have revolved aroundneering and computer science from the Massachusetts Institutea small number of key concepts as well, which I tried to highof Technology, respectively. Until 1999, he was a senior research
light in this article.
engineer at SRI International, Menlo Park, California. He is a■■ The most successful modern approaches in image andprofessor of electrical engineering and was an associate dean forvideo processing are nonparametric. We have drifted away
research at the Baskin School of Engineering at the University
from model-based methods, which have dominated signalof California, Santa Cruz. He is currently a visiting scientist atprocessing for decades.
Google[x]. His research interests include statistical signal,■■ In one-dimensional signal processing, there is a long
image and video processing, computational photography, andhistory of design and analysis of adaptive filtering algocomputer vision. He won the 2000 U.S. National Science Founrithms. A corresponding line of thought has only
dation CAREER Award and the 2010 IEEE Signal Processing
recently become the dominant paradigm in processing
Society (SPS) Best Paper Award. He is a member of the SPS’s
higher-dimensional data. Indeed, adaptivity to data is aAwards Board as well as the Image, Video, and Multidimensional
central theme of all the algorithms and techniques disSignal Processing Technical Committee. He is an IEEE Fellow.
cussed here, which represent a snapshot of the state of
the art.
■■ The traditional approach that every graduate student ofReferences[1] D. Barash, “A fundamental relationship between bilateral filtering, adaptive
our field learns is to carefully design a filter, apply it to thesmoothing and the nonlinear diffusion equation,” IEEE Trans. Pattern Anal.
Mach. Intell., vol. 24, no. 6, pp. 844–847, June 2002.
given data, and call it a day. In contrast, many of the most[2] D. Barash and D. Comaniciu, “A common framework for nonlinear diffusuccessful recent approaches involve repeated applications ofsion, adaptive smoothing, bilateral filtering and mean shift,” Image Vis. Comput.,a filter or sequence of filters to data, and aggregating thevol. 22, no. 1, pp. 73–81, 2004.
[3] M. Elad, “On the origin of the bilateral filter and ways to improve it,” IEEE
results. Such ideas have been around for some time in statisTrans. Image Process., vol. 11, no. 10, pp. 1141–1150, Oct. 2002.
tics, machine learning, and elsewhere, but we have just[4] O. Scherzer and J. Weickert, “Relations between regularization and diffusion
begun to make careful use of sophisticated iteration andfiltering,” J. Math. Imaging Vision, vol. 12, no. 1, pp. 43–63, Feb. 2000.
boosting mechanisms.
[5] A. Singer, Y. Shkolinsky, and B. Nadler, “Diffusion interpretation of nonlocal
neighborhood filters for signal denoising,” SIAM J. Imaging Sci., vol. 2, no. 1,
As I highlighted here, deep connections between techniques
pp. 118–139, 2009.
used commonly in computer vision, image processing, graph[6] L. Pizarro, P. Mrazek, S. Didas, S. Grewenig, and J. Weickert, “Generalised
nonlocal image smoothing,” Int. J. Comput. Vision, vol. 90, no. 1, pp. 62–87, 2010.
ics, and machine learning exist. The practitioners of these fields[7] A. Buades, B. Coll, and J. M. Morel, “A review of image denoising algorithms,
have been using each other’s techniques either implicitly orwith a new one,” Multiscale Modeling Simulation (SIAM Interdisciplinary J.),
vol. 4, no. 2, pp. 490–530, 2005.
explicitly for a while. The pace of this convergence has quick[8] C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,”
ened, and this is not a coincidence. It has come about throughin Proc. IEEE Int. Conf. Comput. Vision, Bombay, India, pp. 836–846,
many years of scientific iteration in what my late friend and colJan. 1998.
[9] H. Takeda, S. Farsiu, and P. Milanfar, “Kernel regression for image processleague Gene Golub called the “serendipity of science,” an aptiing and reconstruction,” IEEE Trans. Image Process., vol. 16, no. 2, pp. 349–366,
tude we have all developed for making desirable discoveries by
Feb. 2007.
happy accident.
[10] K. Dabov, A. Foi, V. Katkovnic, and K. Egiazarian, “Image denoising by sparse
Acknowledgments
I thank the editorial board of IEEE Signal Processing Magazine, and Prof. Antonio Ortega in particular, for their supportand valuable feedback. I am grateful to all my students in theMultidimensional Signal Processing research group at the
University of California, Santa Cruz for their insights andcomments on earlier versions of this manuscript. In particular, I acknowledge the diligent help of Xiang Zhu for many of
the simulations contained here. I am also thankful for feedback and encouragement of my friends and colleagues Prof.
Michael Elad of the Technion and Prof. Bernd Girod of Stanford. This work was supported in part by the U.S. Air Force
grant FA9550-07-1-0365 and the National Science Foundation
grant CCF-1016018.

3D transform-domain collaborative filtering,” IEEE Trans. Image Process., vol. 16,
no. 8, pp. 2080–2095, Aug. 2007.

[11] S. M. Smith and J. M. Brady, “Susan: A new approach to low level image processing,” Int. J. Comput. Vision, vol. 23, no. 1, pp. 45–78, 1997.
[12] H. Knutsson and C. F. Westin, “Normalized and differential convolution:
Methods for interpolation and filtering of incomplete and uncertain data,” inProc. IEEE Comput. Society Conf. Comput. Vision Pattern Regocnit. (CVPR),
June 1993, pp. 515–523.
[13] L. Yaroslavsky, Digital Picture Processing. Berlin, Germany: SpringerVerlag, 1987.
[14] K. Fukunaga and L. D. Hostetler, “The estimation of the gradient of a density
function, with applications in pattern recognition,” IEEE Trans. Inform. Theory ,vol. 21, no. 1, pp. 32–40, 1975.
[15] D. Comaniciu and P. Meer, “Mean shift: A robust approach toward feature space
analysis,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 24, no. 5, pp. 603–619,
May 2002.
[16] A. Smola and B. Scholkopf, “A tutorial on support vector regression,” Statistics Comput., vol. 14, no. 3, pp. 199–222, 2004.
[17] A. Smola, B. Scholkopf, and K. Muller, “The connection between regularization operators and support vector kernels,” Neural Netw., vol. 11, no. 4, pp. 637–
649, 1998.

IEEE SIGNAL PROCESSING MAGAZINE [127] january 2013

[18] R. Rosipal and L. J. Trejo, “Kernel partial least squares regression in reproducing kernel Hilbert space,” J. Mach. Learning Res., vol. 2, pp. 97–123, Dec.
2001.
[19] F. Girosi, “An equivalence between sparse approximation and support vector
machines,” Neural Comput., vol. 10, no. 6, pp. 1455–1480, 1998.
[20] M. P. Wand and M. C. Jones, Kernel Smoothing (Monographs on Statistics
and Applied Probability). London: Chapman & Hall, 1995.
[21] W. Hardle, Applied Nonparametric Regression. Cambridge, U.K.: Cambridge
Univ. Press, 1990.
[22] P. Lancaster and K. Salkauskasr, “Surfaces generated by moving least squares
methods,” Math. Comput., vol. 87, no. 155, pp. 141–158, 1981.
[23] D. Levin, “The approximation power of moving least-squares,” Math. Comput., vol. 67, no. 224, pp. 1517–1531, 1998.

[49] P. Chatterjee and P. Milanfar, “Is denoising dead?” IEEE Trans. Image
Process., vol. 19, no. 4, pp. 895–911, Apr. 2010.
[50] P. Chatterjee and P. Milanfar, “Patch-based near-optimal denoising,” IEEE
Trans. Image Process., vol. 21, no. 4, pp. 1635–1649, Apr. 2011.
[51] P. Perona and J. Malik, “Scale-space and edge detection using anistropic
diffusion,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 12, no. 9, pp. 629–639,
July 1990.
[52] S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, “An iterative regularization method for total variation-based image restoration,” Multiscale Modeling
Simulation, vol. 4, no. 2, pp. 460–489, 2005.
[53] P. Buhlmann and B. Yu, “Boosting with the l2 loss: Regression and classification,” J. Amer. Statist. Assoc., vol. 98, no. 462, pp. 324–339, 2003.
[54] J. W. Tukey, Exploratory Data Analysis. Reading, MA: Addison Wesley, 1977.

[24] M. Alexa, J. Behr, D. Cohen-or, S. Fleishman, and D. Levin, “Computing andrendering point set surfaces,” IEEE Trans. Visualization Comput. Graph., vol. 9,
no. 1, pp. 3–15, 2003.

[55] N. Nordstrom, “Biased anisotropic diffusion: A unified regularization anddiffusion approach to edge detection,” Image Vision Comput., vol. 8, no. 4, pp.
318–327, 1990.

[25] R. Gonzalez and J. Woods, Digital Image Processing, 3rd ed. Englewood
Cliffs, NJ: Prentice-Hall, 2008.

[56] D. P. Palomar and S. Verdu, “Gradient of mutual information in linear vector
Gaussian channels,” IEEE Trans. Inform. Theory, vol. 52, no. 1, pp. 141–154,
2006.

[26] C. Kervrann and J. Boulanger, “Optimal spatial adaptation for patch-based
image denoising,” IEEE Trans. Image Process., vol. 15, no. 10, pp. 2866–2878,
Oct. 2006.
[27] S. P. Awate and R. T. Whitaker, “Unsupervised, information-theoretic, adaptive
image filtering for image restoration,” IEEE Trans. Pattern Anal. Mach. Intell,vol. 28, no. 3, pp. 364–376, Mar. 2006.
[28] A. Buades, B. Coll, and J. Morel, “A non-local algorithm for image denoising,” in Proc. CVPR, June 2005, pp. 60–65.
[29] S. Awate and R. Whitaker, “Higher-order image statistics for unsupervised,
information-theoretic, adaptive, image filtering,” in Proc. CVPR , June 2005, pp.
44–51.
[30] H. Takeda, S. Farsiu, and P. Milanfar, “Kernel regression for image processing and reconstruction,” IEEE Trans. Image Process., vol. 16, no. 2, pp. 349–366,
Feb. 2007.
[31] H. Knutsson and C.-F. Westin, “Normalized and differential convolution:
Methods for interpolation and filtering of incomplete and uncertain data,” inProc. IEEE Comput. Soc. Conf. Comput. Vision Pattern Recognit., June 16–19
1993, pp. 515–523.
[32] N. Sochen, R. Kimmel, and A. Bruckstein, “Diffusions and confusions in signal
and image processing,” J. Math. Imaging Vision, vol. 14, no. 3, pp. 195–209, 2001.
[33] A. Spira, R. Kimmel, and N. Sochen, “A short time Beltrami kernel forsmoothing images and manifolds,” IEEE Trans. Image Process., vol. 16, no. 6,
pp. 1628–1636, 2007.
[34] J. Ham, D. Lee, S. Mika, and B. Schölkopf, “A kernel view of the dimensionality reduction of manifolds,” in Proc. 21st Int. Conf. Mach. Learning, 2004.
[35] T. Hofmann, B. Scholkopf, and A. J. Smola, “Kernel methods in machinelearning,” The Ann. Statistics, vol. 36, no. 3, pp. 1171–1220, 2008.
[36] N. Dowson and O. Salvado, “Hashed nonlocal means for rapid image filtering,”
IEEE Trans. Pattern Anal. Mach. Intell., vol. 33, no. 3, pp. 485–499, 2011.
[37] J. Baek and D. Jacobs, “Accelerating spatially varying Gaussian filters,” ACM
Trans. Graph., vol. 29, no. 6, Dec. 2010.
[38] M. Belkin and P. Niyogi, “Towards a theoretical foundation for Laplacianbased manifold methods,” in Proc. 18th Conf. Learning Theory (COLT), 2005,
pp. 486–500.
[39] A. Szlam, M. Maggioni, and R. Coifman, “Regularization on graphs withfunction-adapted processes,” J. Mach. Learning Res., vol. 9, pp. 1711–1739,
Aug. 2008.
[40] E. Seneta, Non-Negative Matrices and Markov Chains (Springer Series in
Statistics). Berlin, Germany: Springer, 1981.
[41] R. A. Horn and C. R. Johnson. Matrix Analysis. Cambridge, U.K.: Cambridge
Univ. Press, 1991.
[42] H. J. Landau and A. M. Odlyzko, “Bounds for eigenvalues of certain stochastic
matrices,” Linear Algebra Its Applicat., vol. 38, pp. 5–15, June 1981.
[43] W. J. Stewart, Introduction to the Numerical Solution of Markov Chains.
Princeton, NJ: Princeton Univ. Press, 1994.
[44] P. Knight, “The Sinkhorn-Knopp algorithm: Convergence and applications,”
SIAM J. Matrix Anal. Appl., vol. 30, no. 1, pp. 261–275, 2008.
[45] R. Sinkhorn, “A relationship between arbitrary positive matrices and doubly
stochastic matrices,” Annals Math. Statist., vol. 35, no. 2, pp. 876–879, 1964.
[46] J. Darroch and D. Ratcliff, “Generalized iterative scaling for log-linear models,” Ann. Math. Statist, vol. 43, no. 5, pp. 1470–1480, 1972.
[47] L. Huang, D. Yan, M. I. Jordan, and N. Taft, “Spectral clustering with perturbed data,” in Proc. Neural Information Process. Syst. (NIPS), 2008.
[48] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by
sparse 3-D transform-domain collaborative filtering,” IEEE Trans. Image
Process., vol. 16, no. 8, pp. 2080–2095, Aug. 2007.

[57] Y. Freund and R. E. Schapire, “A short introduction to boosting,” J. Japan.
Soc. Artif. Intell., vol. 14, no. 5, pp. 771–780, Sept. 1999.
[58] P. Buhlmann and T. Hothorn, “Boosting algorithms: regularization, prediction and model fitting,” Statist. Sci., vol. 22, no. 4, pp. 477–505, 2007.
[59] Y. Censor and S. Zenios, Parallel Optimization (Numerical Mathematics
and Scientific Computation). New York: Oxford Univ. Press, 1997.
[60] M. Collins, R. Schapire, and Y. Singer, “Logistic regression, Adaboost andBregman distances,” Mach. Learning, vol. 48, no. 1–3, 2002.
[61] M. Charest and P. Milanfar, “On iterative regularization and its application,”
IEEE Trans. Circuits Syst. Video Technol., vol. 18, no. 3, pp. 406–411, Mar. 2008.
[62] S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Trans. Signal Process., vol. 41, no. 12, pp. 3397–3415, Dec. 1993.
[63] Z. Farbman, R. Fattal, D. Lischinski, and R. Szeliski, “Edge-preserving decompositions for multi-scale tone and detail manipulation,” in Proc. ACM SIGGRAPH
ACM Trans. Graph., vol. 27, no. 3, Aug. 2008.
[64] F. Luisier, T. Blu, and M. Unser, “A new SURE approach to image denoising: Interscale orthonormal wavelet thresholding,” IEEE Trans. Image Process.,
vol. 16, no. 3, pp. 593–606, Mar. 2007.
[65] Z. Wang and A. Bovik, “Mean squared error: Love it or leave it? A new look at
signal fidelity measures,” IEEE Signal Process. Mag., vol. 26, no. 1, pp. 98–117,
Jan. 2009.
[66] B. Carlin and T. Louis, Bayes and Empirical Bayes Methods for Data
Analysis (Texts in Statistical Science). London: Chapman & Hall, 2000.
[67] T. Weissman, E. Ordentlich, G. Seroussi, S. Verdu, and M. Weinberger, “Universal discrete denoising: Known channel,” IEEE Trans. Info. Theory, vol. 51, no.
1, pp. 5–28, 2005.
[68] S. Awate and R. Whitaker, “Feature-preserving MRI denoising: A nonparametric empirical Bayes approach,” IEEE Trans. Med. Imaging, vol. 26, no. 9, pp.
1242–1255, 2007.
[69] K. Sivaramakrishnan and T. Weissman, “A context quantization approach to universal denoising,” IEEE Trans. Signal Process., vol. 57, no. 6, pp. 2110–2129, 2009.
[70] A. Bruckstein, “On globally optimal local modeling: From moving leastsquare to overparametrization,” in Proc. SIAM Conf. Imaging Sci., Chicago, IL.
2000.
[71] L. Wasserman, All of Nonparametric Statistics. Berlin, Germany: Springer,
2005.
[72] M. Aharon, M. Elad, and A. Bruckstein, “The K-SVD: An algorithm for designing of overcomplete dictionaries for sparse representation,” IEEE Trans.
Signal Process., vol. 54, no. 11, pp. 4311–4322, Nov. 2006.
[73] R. Khoury, “Closest matrices in the space of generalized doubly stochastic
matrices,” J. Math. Anal. Applicat., vol. 222, no. AY985970, pp. 562–568, 1998.
[74] W. Kahan, “Spectra of nearly Hermitian matrices,” in Proc. Amer. Math.
Soc., vol. 48, pp. 11–17, 1975.
[75] P. Milanfar, “Symmetrizing smoothing filters,” SIAM J. Imaging Sci., to be
published. DOI: 10.1137/120875843[76] H. Jegou, C. Schmid, H. Harzallah, and J. Verbeek, “Accurate image search using the contextual dissimilarity measure,” IEEE Trans. Pattern Anal. Mach. Intell.,vol. 32, no. 1, pp. 2–11, Jan. 2010.
[77] C. Johnson, R. Masson, and M. Trosset, “On the diagonal scaling of Euclidean
distance matrices to doubly stochastic matrices,” Linear Algebra Its Applicat., vol.
397, no. 1, pp. 253–264, 2005.
[78] D. Yan, L. Huang, and M. I. Jordan, “Fast approximate spectral clustering,”
in Proc. 15th ACM Conf. Knowledge Discovery Data Mining (SIGKDD), 2009,
pp. 907–916.

[SP]

IEEE SIGNAL PROCESSING MAGAZINE [128] january 2013