IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

1

Understanding and Eliminating
the Large-kernel Effect in Blind Deconvolution

arXiv:1706.01797v3 [cs.CV] 19 Nov 2017

Li Si-Yao, Qian Yin*, and Ping Guo, Senior Member, IEEE

Abstract—Blind deconvolution consists of recovering a clear
version of an observed blurry image without specific knowledge
of the degradation kernel. The kernel size, however, is a required
hyper-parameter that defines the range of the support domain.
In this study, we experimentally and theoretically show how large
kernel sizes introduce noises to expected zeros in the kernel and
yield inferior results. We explain this effect by demonstrating
that sizeable kernels lower the squares cost in optimization. We
also prove that this effect persists with a probability of one for
noisy images. Using 1D simulation, we quantify the increment of
error of estimated kernel with its size. To eliminate this effect,
we propose a low-rank based penalty that reflects structural
information of the kernel. Compared to the generic `α , our
penalty can respond to even a small amount of random noise in
the kernel. Our regularization reduces the noise and efficiently
enhances the success rate of large kernel sizes. We also compare
our method to state-of-art approaches and test it using real-world
images.
Index Terms—image deblurring, blind deconvolution, kernel
size, regularization.

I. I NTRODUCTION

B

Lind deconvolution methods are to recover a latent clear
version from an observed blurry image without knowing
the degradation kernel. The model of image blur is
y = x ∗ k + n,

(1)

where y is the blurred image, x is the latent clear one, k is the
kernel (a.k.a. Point Spread Function (PSF)) and n represents
random Gaussian noise produced during photographing.
One popular approach is to maximize the posterior (MAP)
p(x, k|y) and reformulate the problem into a regularized least
squares optimization [1]:



x̂, k̂ = arg min kx ∗ k − yk2 + λg (x) + σh (k)
(2)
x,k

where g and h are prior function designed to prefer sharp
images and ideal kernels respectively. This two-variable optimization is implemented as two alternate sub-steps:


x̂(i+1) = arg min kx ∗ k̂ (i) − yk2 + λg (x)
(3)
x

This work is funded by National Natural Science Foundation of China
(NSFC) (61472043, 61375045) and Joint Research Fund in Astronomy under
cooperative agreement between the National Natural Science Foundation of
China (NSFC) and Chinese Academy of Sciences (CAS) (U1531242).
Li Si-Yao and Q. Yin are with the College of Information Science
and Technology, Beijing Normal University, Beijing, 100875 CN (e-mail:
lisiyao@mail.bnu.edu.cn, yinqian@bnu.edu.cn).
Ping Guo is with the School of System Science, Beijing Normal University,
Beijing, 100875, CN (email: pguo@ieee.org).
*Q. Yin is the corresponding author.

and


k̂ (i+1) = arg min kx̂(i+1) ∗ k − yk2 + σh (k) .
k

(4)

Although k is unknown, most implementations require
kernel size as a manually tuned hyper-parameter that limits the
support domain. The improper selection of kernel size yields
inferior results in MAP methods. When the preset kernel size
is too small, there will not be enough room for k̂ to expand
during iterations, which yields high errors of recovered images.
Contrastively, when kernel size is greater than ground truth,
MAP methods also yield inferior results. This phenomenon
was firstly mentioned by Fergus et al. [2]. Similarly, Cho and
Lee [3] showed that the residue cost of (2) increases with
kernel size when k is larger than the proper one. To show
this effct clearly we tested our implementation described in
Section 4 (but set µ and σ zero) using one-half to three times
the truth size iteratively. Figure 1(d-f) shows that the larger the
kernel size than ground truth the poorer the deblurred image.
At the same time, kernels become more noisy with higher
numerical errors. Figure 1(b) shows both the error ratio of
restored images and the Sum of Squared Difference (SSD) of
estimated kernels reach the lowest around the truth size and
increase afterwards.
It should be noted here that in this experiment, we ommited
regularization h, hence the k-step equals to a bare least squares
problem:
k̂ (i+1) = arg min kx̂(i+1) ∗ k − yk2 .
k

(5)

We also avoided using multi-scaling scheme and threshold in
this experiment. Parameters that performed well on the truth
size were kept identical for larger sizes during the experiment.
In this study, to explore the source of this effect, we transfer
convolution parts of (3) and (4) to affine transformations and
analyse their properties on kernel size. We show that for x
in sparse distributions of clear images, this large-kernel effect
remains with probabiliry 1. Further, we simulate and quantify
how the error increases with kernel size.
Ideal kernels are regarded to be sparse. A variety of generic
sparsity regularizations are used on k in the literature. The
most popular choice is `2 norm [4]–[7], taking advantage of
the conviniece of optimization. Some other implementations
used `1 norm [8]–[10] while Zuo et al. [11] used the hyperLaplacian prior (`α with 0 < α < 1) to maintain the sparsity
of the estimated kernel. However, previous works generally
treated regularization h as an accessory of the success of MAP
methods, e.g., Perrone and Favaro [12] omitted h directly. The
literature lacks a targeted discussion on the selection of kernel

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

2

6

0.02

3

0.015

0

13

23
Kernel size

47

(b)

Kernel SSD

Error ratio

(a)

0.01

truth = 23
err = 1.0
(c)

size = 23
err = 1.6
(d)

size = 47
err = 5.9
(e)

size= 69
err = 69.9
(f)

Fig. 1. Large kernels produce inferior results. (a) The blurry image. (b) The numerical errors on kernel size. (c-f) Deblurred results: in the first row are
restored images; in the second row are corresponding estimated kernels; in the third row are the support domains (k > 0); in the last row are the domains of
k > max(k)/30, where adjacent positive pixels are colored identically and zeros are white.

regularization. In this study, we take a novel regularization
targeting to reduce the noise in k̂ and eliminate the unwanted
effect yielded by large kernels.
In addition to the regularization, operations after the kstep also enhance the quality of the estimated kernel. Cho
and Lee [6] removed any elements smaller than 1/20 of the
maximum of k̂ to purify the estimated kernel. Xu and Jia [4]
proposed an iterative support domain detector based on the
differences of elements of k̂. Beyond various refinements, two
constraints of k are widely recognized: (i) the sum of k equals
1 and (ii) the elements of k are no smaller than 0. With respect
to (i), Kundur et al. [13] indicates that the mean of a clear
image is preserved in the blurry; Levin et al. [14] emphasizes
that the normalization prevents scaling conflicts. For (ii), the
PSF is naturally greater than zero because the contribution of
light is positive. These constraints are realized as projections
after the k-step [8], [11]. Perrone and Favaro [12] showed that
projections after alternate Gradient Descent prevent the trivial
(delta) solution using 1D simulations.
II. A NALYSIS
Before analysing the source of larger-kernel effect, we firstly
introduce an interesting fact that we call Plumping Effect.
Theorem 1. (Plumping Effect) Let A = [v1 . . . vn ], where
vi ∈ Rm (m ≥ n + 1). Let B = [w1 A w2 ] where w1 , w2 ∈
Rm and rank(B) > rank(A). Given an m-D random vector

b whose elements are i.i.d. with the continuous probability
density function p, for u ∈ Rm

Pr inf{kBu − bk2 } < inf{kAu − bk2 } = 1.
Proof.
=
=


Pr inf{kBu − bk2 } ≥ inf{kAu − bk2 }

Pr b ∈ Rm \ span{w1 , w2 }
Z
dp(b)
Ω

where Ω = Rm \ span{w1 , w2 }. For w1 , w2 are linear
independent to {v1 , · · · , vn }, we have w1 6= 0 or w2 6= 0.
Hence, the Lebesgue measure of Ω is zero, and the probability
is zero.
Theorem 1 shows that padding linear independent colums to
a thin matrix yields lower residue cost; hence the least squares
solution will change.
The degradation (convolution) part in (1) is equivalent to
linear transformations:
x ∗ k = Tk x = Tx k.

(6)

where x, y ∈ RM ×N , k ∈ RL×K ; Tk and Tx are blocked
banded Toeplitz matrices [15], [16]; Tk ∈ RM N ×M N and
Tx ∈ RM N ×LK ; A represents the column-wise expanded
vector of matrix A; L and K are required to be odd.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

3

0.8

Dimensions of Tk will not change with kernel size. Especially, Tk remains the same when L and K increase by
wrapping a layer of zeros around k. The result of the x-step
keeps the same if k̂ expanded by zero padding, hence the xstep should not be blamed as the source of the large-kernel
effect. In contrast, when k is larger, Tx will become plumper
for the same x. In 1D cases, where N = K = 1, assume
L = 2l + 1, then

h

= J

(−l)

0

.

..

.

xM −1

x2
..
.
..
.
..
.

xM
0

xM −1
xM

..

.
···

x ··· J

(−1)

x1

0
x1

x x J

(1)

···
..
.
..
.
..
.

x2
..
.
..
.
..
.
xM −1

x ··· J

(l)

···
i
x .

0
.. 
. 


0 


x1 


x2 

.. 
. 

Pr (rank (Tx (L + 2)) > rank (Tx (L))) = 1,

Pr (rank (Tx (M )) = M ) .

0

23

1

23

(d)

0.4

0.6
0.3
0

1 13

(b)

0.2

0

35 47

0.8

ground truth
double size

1 13

0.4

ground truth
four times

0.6
0.3

35 47

(e)

0.2

0
−0.3
1

36

index

(c)

58

0
1

93

36

index

58

93

(f)

1

0.5

0

2347 93
Kernel size

253

(g)

(8)

which is the condition of the Plumping Effect. Firstly, we have
Pr (rank (Tx (L + 2)) > rank (Tx (L)))

(a)

0.8

−0.3

xM −l

1

Normalized kernels

..

x2
..
.
..
.
..
.



(7)
During bind deconvolution iterations, for identical values
of x̂(i) , a larger L introduces a greater number of columns
onto both sizes of Tx̂(i) and results in different solutions.
To illustrate this point, we tested a 1D version of blind
deconvolution without kernel regularization and took different
values of L (truth and double and four times the truth size)
for the 50th k-step optimization after 49 truth-size iterations
(see Figure 2). Figure 2(a-c) show that the optimal solutions in
different sizes differ slightly on the main body that lies within
the ground truth size (colored in red), but greatly outside this
range (colored in green), where zeros are expected. Figure 2(df) compare ground truth to estimated kernels in (a-c) after
nonnegativity and sum-to-one projections. Larger sizes yield
a greater number of positive noises; hence, they lower the
weight of the main body after the projection and change the
outlook of kernel.
Even if x̂(i) successfully iterates to truth x, Theorem 1
implicates the large-kernel effect remains for random Gaussian
noise n. We now show

≥

0.2

Residue square cost ratio

0

0.3

−0.3

···

ground truth
truth size

0.6

Optimized kernels

Tx (L)

xl+1
 ..
 .


xM −1


=  xM


 0

 .
 ..

0.4

(9)

Kaltofen and Lobo [17] proved that for an M-by-M Toeplitz
matrix composed of finite filed of q elements,

Pr rank(TM ×M ) = M = 1 − 1/q.
(10)
In our case, clear images have statistical sparsity on derivative
fields [18], [19], and elements of x are modeled to be
continuous in hyper-Laplacian distributions [20]:
(
β exp (−γ|xj |α ) xj ∈ [−1, 1]
p(xj ) =
(11)
0
otherwise.

Fig. 2. The estimated kernels in 1D Blind deconvolution simulations. (ac) are optimized kernels of different sizes after the 50th iteration. (d-f) are
corresponding normalized kernels to left after nonegativity and sum-to-one
projections. (g) shows the residue squares cost of (5) of the optimized kernel
on its size. In this experiment, x is a 255×1 vector extracted from a real image
and the truth k is generated by column-wisely adding a 23 × 23 truth kernel
from Levin’s dataset [14]. The signal prior is ``1 . This figure is recommended
2
to view in color.

Then we have
Theorem 2.

Pr rank(Tx (M ) = M = 1.
Proof. See Appendix A.
Theorem 1 implicates that the larger the kernel size the

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

4

greater the error of k due to noise n. We now quantify this
increment.
Assume x̂ iterates to ground truth xgt during iterations, then,
for the estimated kernel, we have

During blind deconvolution iterations, implicit noise [9]
occurs when x̂ 6= xgt . Cho and Lee [3] indicated that x̂
should be regarded as a sparse approximation to xgt . Assume
x̂ = xgt + ∆x and n = 0, then
x̂ ∗ kgt =

k̂ = arg min kTxgt k − yk2
k

=

2

= arg min kTxgt k − Txgt kgt − nk
k

(12)

=kgt + arg min kTxgt k − nk2

and

k

= arg min k(Txgt + T∆x )k − Txgt kgt k2

Tx†gt n

Assume knk = 1, then

=(Txgt + T∆x )† Txgt kgt .
Then

SSD = k (Txgt + T∆x )† Txgt − I kgt k2 .

≤

kTx†gt nk

≤

sn (Tx†gt ),

(14)

1
1
≤ kTx†gt nk ≤
.
sn (Txgt )
s1 (Txgt )

(15)

Based on [21], [22], for arbitrary matrix A = [a1 , a2 , ..., ad ],
we have
height(a1 , a2 , ..., ad )
√
≤ s1 (A)
(16)
d
and
sn (A) ≤

√

To quantify how singular values of (Txgt + T∆x ) Txgt − I
changes with kernel size, we simulated 100 times, in each
of which we generated random x with length 254 and PDF
in (11) with M = 254, γ = 10 and α = 0.5, and generated random Gaussian vector ∆x where k∆xk = kxk/100.
Figure 4 (a) shows one example of generated x and ∆x.
Figure 4(b) shows means and standard deviations of s1 , sn
and s̄ of simulated (Tx + T∆x )† Tx − I on L, respectively. The
error of k̂ is expected to grow with kernel size even n = 0.
0.6

^x
x gt
Δx

0.3

d max kai k.

(17)

0

For arbitrary A = Tx (L), maxi kai k = kxk ≤ M . Hence,
√
1
L
√
≤ kTx† nk ≤
. (18)
height(J (−l) x, ..., x, ..., J (l) x)
ML

−0.3

i

Function height will decrease when L grows, because more
linear independent columns are introduced into A. To quantify
this decrement, we simulated 100 times, in each of which we
generated random x with length M = 254 and PDF as in (11)
with γ = 10 and α = 0.5. Figure 3(b) shows the decreasing
mean and standard deviation of simulated heights on L. Using
this simulated data, we plotted the lower and upper bound of
kTx† nk based on (18). We also extracted a line xex from a
clear image of Levin’s set [14], generated a normalized random
vector n and compared kTxex nk to theoretical boundaries.
Figure 3 (c) shows the error increases hyper-linearly.

height

0
−0.1
−0.2
1

index

(a)

254

||T x(L) +n||

1.8
1.2
0.6
0

25

L

(b)

253

50

upper bound
lower bound
experiment

−0.6
1

0.3

s1
s- n
s

0.2

√

0.1

(21)
†

where s1 and sn represent smallest and largest singular values.
Equavalently, we have

0.2

(20)

k

where A† represents Moore-Penrose pseudoinverse of matrix
A. Then
(13)
SSD = kk̂ − kgt k2 = kTx†gt nk2 .

s1 (Tx†gt )

(19)

xgt ∗ kgt + n∆x ,

k̂ = arg min kTx̂ k − yk2

k

=kgt +

xgt ∗ kgt + ∆x ∗ kgt

0.1

index

(a)

254

0
25

51

75
L

101

125

(b)

Fig. 4. (a) Simulated sparse signal examples. (b) Simulated singular values
of (Tx + T∆x )† Tx − I.

III. L OW- RANK REGULARIZATION
The large-kernel effect is yielded by noise in ultra-sized
kernels. Figure 1 and Figure 2 show that main bodies of
kernels are kept but noise take greater amounts when kernel
is larger. To constrain k̂ to be ideally clean, regularization h
is expected to distinguish noise from clean kernels efficiently.
To eliminate the noise in kernel, we take low-rank regularization on k such that (4) becomes


k̂ (i+1) = arg min kx̂(i+1) ∗ k − yk2 + σrank (k) . (22)
k

10
2
25

101

L

253

(c)

Fig. 3.
(a) The extracted row from a clear image. (b) The mean and
standard deviation of simulated heights. (c) Theoritical bounds using (b) and
experimental result of (a).

Because the direct rank optimization is an NP-hard problem,
continuous proximal functions are required. Fazel et al. [23]
proposed
log det ((X + δI))
(23)
as a heuristic proxy for X ∈ SN
+ where I is the N-by-N
identity matrix and δ is a small positive number. A 1D case
is shown in Figure 3.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

5

5

1
0

rank(x)
|x|
log(|x|+δ)

−3
−5

0
x

5

Fig. 5. 1D demonstration of the rank function, the nuclear norm and log det.

To allow this approximation to play a role in general
matrices, the low-rank object is substituted to (XX T )1/2 [24].
The regularization function then becomes
X
1
log(si + δ) (24)
h(X) = log det((XX T ) 2 + δI) =

are distinguished by high log det cost from real kernels (see
Figure 4(c)). Figure 4(d) shows that singular values of lowrank regularized kernel are distributed similarly as ground
truth, compared with the impure one.
Compared to previously used `α norms, low-rank regularization responds more efficiently to noise. To illustrate
this point, we generated a noisy kernel by adding a small
percentage () of Gaussian noise and 1 −  of the real
kernel. Figure 5 shows that the low-rank cost rapidly adjust
favorably to the noise but `α norms fail. That is because `α
only takes statistical information. An extrem example consists
of distrupting a truth kernel and randomly reorganizing its
elements, with `α cost unchanged. In contrast, rank (singular
values) corresponds to structural information of a matrix.
3

L

j
2.5

+

2
1
0.8

rank
logdet
nuclear norm

1.5
1
0.5
0

47×47

Sigular value

L

2
Cost ratio

where si is the i-th singular value of X.
The motivation to take low-rank regularization on k is based
on a common phenomenon of noise matrices. Figure 4(a-b)
shows a nonnegativity projected Gaussian noise matrix and
its singular values in decreasing order. For a noise matrix,
where light and dark alternate irregularly, the distribution
of singular values decays sharply at lower indices, then it
breaks and drag a relatively long and flat tail to the last. This
phenomenon is generic [15]. Based on this fact, noise matrices

L

0.1

0.2

ε

0.3

0.4

0.5

Fig. 7. Comparison on respond to noise. The cost ratio is calculated as
cost()−cost(0)
1+
. This figure is recommended to view in color.
|cost(0)|
¬ break aorund 2
2 ←

IV. I MPLEMENTATION
A. Updating k
0

(a)

1

Singular value ratio

logdet cost

−100

−400

−700

23

Index

(c)

47

It should be firstly noted that the log det function is nonconvex ( it is actually concave on S+ ). To solve the low-rank
regularized least squares problem 4 with h defined in (24),
we introduce the auxiliary variable Ψ = k and reformulate
the optimization into



1
min kx(i+1) ∗ Ψ − yk2 + σ log det (kk) 2 + δI
k,Ψ
(25)
s.t. Ψ = k

Index

(b)
impure kernel
clean kernel
low−rank regularizaed

0.5

0
1

long tail
¯↓
23
Index

47

(d)

Fig. 6. Singular values of clean kernels and noisy matrices. (a) shows the
support domain (black) of a 47 × 47 random Gaussian noise matrix after
nonnegativity projection. (b) is the distribution of singular values of (a). In
(c) we compare the log det cost of random Gaussian noise matrix (black),
the truth kernel from [14] after zero-padding (red), and the Gaussian PSF
with size/6 standard deviation (blue) on kernel size. In (d) we compare the
distribution of scaled singular values (the maximum to 1) of clean, impure
and regularized kernels. This figure is recommended to view in color.

Using the Lagrange method, (25) is solved by two alternate
sub-optimizations



(j+1)

= arg min kx̂(i+1) ∗ Ψ − yk2 + µkΨ − k̂ (j) k2
 Ψ̂
Ψ


1

(j+1)
(j+1)
2
 k̂
= arg min
kk − Ψ̂
k + σh (k)
k
2τ
(26)
where j is the iteration number while µ and τ are trade-off
parameters.
The Ψ-substep is convex and accomplished using the Conjugate Gradient method. For k-substep, the requirement for low
rank is limited; otherwise, the regularization may change the

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

6

Algorithm 1 Updating k with low-rank regularization
Input: x, y, µ, τ , OuterIterM ax, CGIterM ax,
innerIterM ax
Output: k̂
1: for j ← 0 to OuterIterM ax − 1 do
2:
if j = 0 then
3:
Ψ̂(j+1) ← minΨ kx ∗ Ψ − yk2 using CG with
maximum CGIterM ax iterations
4:
else

5:
µ(j) = µ ej /eOuterIterM ax

6:
Ψ̂(j+1) ← minΨ kx ∗ Ψ − yk2 + µ(j) kΨ − k (j) k2
using CG with maximum CGIterM ax iterations
7:
end if
8:
Initializing k (0) with all singular values equal to 1
9:
for t ← 0 to innerIterMax −1 do
10:
k (t+1) ← proxτ h (t) Ψ̂(j)
k
11:
end for
ax)
12:
k̂ (j+1) ← max{k (innerIterM
, 0}
P (j+1)
(j+1)
(j+1)
13:
k̂
← k̂
/ k̂
14: end for
15: k̂ ← k̂ (outerIterM ax)

main body of kernel—an extreme result is k̂ = 0. Thus, our
strategy is to lower the rank at Ψ̂ locally. Using the first-order
Taylor expansion of h at fixed matrix Z:
hZ (X) =h (Z)
+ tr



=h (Z) +

−1 
1
1
1 
ZZ T 2 + δI
XX T 2 − ZZ T 2
X si − ŝi
i

ŝi + δ

(27)
where ŝi is the i-th eigenvalue of Z, the k-substep is transformed into an iterative optimization


1
kk − Ψ̂(j+1) k2 + σhk(t) (k)
k (t+1) = arg min
(28)
k
2τ
where t is the inner iteration number. For convenience we set
σ as a flag (if σ = 0, the k-substep will be skipped) and only
tuned τ as the trade-off parameter.
Define the proximal mapping of function φ as follows:


1
2
ku − vk + φ(u) .
(29)
proxφ (v) = arg min
u
2

non-convex `1 /`2 regularized least squares problem is solved
by fixing `2 . Hence, the problem is transformed into a linear
optimization problem

min kx ∗ k − yk2 + λ0 kxk1
(31)
x

and is solved by iteratively minimizing a convex envelope (for
step η small enough) equivalent to


x(t+1) = proxηλkxk1 u(t) = sgn(u(t) )(|u(t) | − ηλ)+ (32)


where u(t) = x(t) − η ∇x kx ∗ k − yk2 x=x(t) .
It should be noted that we increased the fixed `2 norm
propotionally with sub-sampled image scale when taking an
multiscaling scheme.
C. Overall Implementation
The overall algorithms is shown in Algorithms 2.
Algorithm 2 Blind Deconvolution (single-scaling version)
Input: blurry image y, kernel size L, λ, η, τ , IterM ax
Output: clear image x and degradation kernel k
1: y ← [∇h y, ∇v y]
2: Initializing x ← y
3: Initializing k with an L × L zero matrix adding [0.5 0.5]
in the center
4: for t ← 1 to IterM ax do
5:
Updating x using Algorithm 3 in [8]
6:
Updating k using Algorithm 1
7: end for
8: x ← Non-blind deconvolution on k and y
We implement these algorithms on the derivative field of x
and y based on the fact that
∇ (x ◦ k) = ∇x ◦ k.

(33)

Kernels in our implementation are forced squares. We applied
the non-blind deconvolution method in [20].
V. E XPERIMENTS AND R ESULTS

Dong et al. [24] proved that one solution to the proximal
mapping of τ hZ is
proxτ hZ (X) = U (Σ − τ diag (w))+ V T

(30)

where U ΣV T is SVD of X, wi = 1/ (ŝi + δ) and (A)+ =
max {A, 0}. Local low-rank optimization is iteratively implemented via the given parameter τ (see Algorithm 1).
It should be noted here that µ is designed to exponentially
grow with j to allow more freedom of Ψ̂ for early iterations.
B. Updating x
As `α norms prefer the blurry images [12], [14], we take
g = `1 /`2 in this work, as proposed by Krishnan et al. [8]. The

size = 23, err = 1.55

size = 47, err = 1.56

size = 69, err = 2.14

Fig. 8. Deblurring results using low-rank regularization.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

7

Fig. 9. Tests on real-world images. From left to right: the blurry images, 31 × 31 kernel with low-rank, 61 × 61 without low-rank and 61 × 61 with low-rank.

1

1

0.8

0.8

0.6

Cho & Lee
Shan et al.
Perrone & Favaro
Krishnan et al.
ours
ours µ = σ = 0

0.4
0.2
0
0

1

2
3
Error ratio

(a)

4

Success rate

Success rate

In response to high error ratios of large kernels mentioned
in Figure 1 and to show the efficiency of our regularization,
we used the same parameters except µ and σ in Figure 1.
Figure 8 shows low-rank regularized kernels are much more
robust to kernel size. Noises in kernels are efficiently reduced
and qualities of restored images are enhanced.
We also tested our implementation on real-world images and
compared the results with and without using low-rank regularization (See Figure 9). Restored images without using lowrank regularization look muddy for strong noise in estimated
kernels. In contrast, regularized kernels are much cleaner, and
qualities of restored images match those restored by small
kernels.
It should be noted here that in experiments of Figure 8 and
Figure 9, we deliberately omitted multi-scaling scheme and
threshold to expose the effectiveness of low-rank regularization itself. Those operations may also reduce the large-kernel
effect by sub-sampling and purifying the estimated kernel.
Our regularization method performs robust on larger kernels

VI. C ONCLUSION

0.6

Cho & Lee
Shan et al.
Perrone & Favaro
Krishnan et al.
ours
ours µ = σ = 0

0.4
0.2

5

0
0

upon dataset from [14]. Figure 10 shows the success rates of
state-of-the-art methods versus our implementations with and
without (set µ and σ zero) low-rank regularization. Parameters
are fixed during the whole experiment, where σ = 1, µ = 1,
τ = 5 × 10−5 , OuterIterM ax = 20, CGIterM ax = 3 and
innerIterM ax = 10 and a 7-layer multi-scaling pyramid is
taken. Low-rank regularization works more effectively than the
regularization-free implementation and the state-of-art.
In addition, we test compared our implementation to several
state-of-arts on a real-world image to reveal the robustness of
our regularization on kernel sizes. Figure 11 shows that 185 ×
185 size yields strong noises in estimated kernels of previous
works (Cho & Lee, Krishnan et al.), and even changes main
bodies of kernels (Krishnan et al., Zuo et al.). Contrastively,
our implementation can keep the kernel stable for the larger
size.
Our implementation processed images of 255 × 255 pixels
with truth blur sizes in 85 sec on average on a Lenovo
ThinkCentre computer with Core i7 processor. Our MATLAB
codes and deblurring results will be distributed after the
present paper has been officially accepted for publication.

1

2
3
Error ratio

4

5

(b)

Fig. 10. Success rates. (a) Results on truth sizes. (b) Results on double sizes
of the truth.

In this paper, we demonstrate that kernel sizes larger than
ground truth produce increased noises. We attribute the largekernel effect to the Plumping Effect. We also take low-rank
regularization and reduce the large-kernel effect efficiently. In
practical implementations, even for noise-free y, the intermediate x̂(i) is unlikely to iterate to truth, hence some parts of
y will be treated as implicit noises, which may intensify the
effect even more than expected.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

blurred

8

Cho & Lee 85×85

Cho & Lee 185×185

Krishnan et al. 185×185

Krishnan et al. 85×85

Zuo et al. 185×185

Zuo et al. 85×85

ours 85×85

ours 185×185

Fig. 11. Tests on a real-world image. It’s recommanded to zoom into estimated kernels to compare the cleaness and main bodies of kernels.

det (Tx (M )), the item of largest degree is xM
m+1 with factor
1.

A PPENDIX
P ROOF TO T HEOREM 2
Assume M to be odd and M = 2m + 1 (m ∈ N+ ). Then,
Tx (M )

xm+1
 ..
 .


xM −1


=  xM


 0

 .
 ..
0

Lemma. Let X be a continuous r.v. in the finite support
domain [a, b]. Let P : R → R be a polynomial function
P (x) = xk + Q(x)

···

..

.

x2
..
.
..
.
..
.

..

.

xM −1

x2
..
.
..
.
..
.

xM
0

xM −1
xM

..

.
···

x1

0
x1
x2
..
.
..
.
..
.
xM −1

···
..
.
..
.
..
.

···

h
i
= J (−m) x · · · J (−1) x x J (1) x · · · J (m) x .

0
.. 
. 


0 


x1 


x2 

.. 
. 


where Q is a finite polynomial function with the largest degree
less than k. Generate a new r.v.
Y = P (X).
Then, for ∀y ∈ R, the Cumulative Distribution Function
(CDF) FY is continuous at y.
Proof.
Z
FY (y) = Pr(Y <= y) =

xm+1
(34)
(35)

For any M -by-M matrix A, rank(A) = M i.f.f. det(A) 6=
0. Thus,
Pr (rank (Tx (M )) = M ) = Pr (det (Tx (M )) 6= 0) . (36)
As we know, the explicit formula of determinant of a
Toeplitz matrix on its elements is unsolved in the current literature. Li [25] gives a concrete expression of det (Tx (M )) by
using LU factorization but fails to fit all situations (e.g. when
xm+1 = 0). However, it can be shown that det (Tx (M ))
equals a multivariate polynomial function without manipulating the whole expression. By using Laplace expansion on

fX (x)dx
Ξ(y)

where Ξ(y) = {x|P (x) <= y, x ∈ R}.
For ∀y ∈ R,
Ξ(y + ) = Ξ(y)
and
Ξ(y − ) = Ξ(y) \ Ω(y)
where Ω(y) = {x|P (x) − y = 0}.
Based on Beppo Levi’s Theorem,
Z
Z
lim+
fX (x)dx =
ξ→y

Ξ(ξ)

fX (x)dx.

Ξ(y)

Because P (x) 6≡ c (c is a constant), for ∀y ∈ R, zeros of
P (x) − y are finite, hence the Lebesgue measure of Ω(y) is
zero. We have
Z
Z
lim−
fX (x)dx =
fX (x)dx.
ξ→y

Ξ(ξ)

Ξ(y)

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

9

Thus

where
Y (ξ1 , . . . , ξm , ξm+2 , ξM )

Xm+1 · · ·
ξ2
..
 ..
 .
.


..
 ξM −1
.


.
.
..
..
= det  ξM


..
 0
. ξM −1

 .
..
 ..
.
ξM
0
···
0

FY (y + ) = FY (y − ) = F (y).

Theorem. Let X be a continuous r.v. with PDF
(
fX (x) =

β exp (−γ|x|α ) x ∈ [−1, 1]
0
otherwise.

For a sample of independent observations X1 , . . . , XM ,
generate a new r.v.

0

ξ2
..
.
..
.
..
.

ξ1

ξM −1
ξM

ξ2
..
.
..
.
..
.
ξM −1

···
..
.
..
.
..
.

0
..
.
0
ξ1
ξ2
..
.

···

Xm+1

M
=Xm+1
+ Qξ1 , ..., ξm , ξm+2 , ..., ξM (Xm+1 ) .

Qξ1 , ..., ξm , ξm+2 , ..., ξM is a polynomial function with the
largest degree less than M . Based on Lemma, we have

Z=


···

Xm+1
 ..
 .


XM −1


det  XM


 0

 .
 ..

..

.

X2
..
.
..
.
..
.

..

.

XM −1

X2
..
.
..
.
..
.

XM
0

XM −1
XM

..

.
···

0

X1

0
X1
X2
..
.
..
.
..
.
XM −1

···
..
.
..
.
..
.

0
..
.
0
X1
X2
..
.

···









.







Xm+1

Then,
Pr (Z = 0) = 0.
Proof. Based on the Law of Total Probability and Dominated
Convergence Theorem,
Pr (Z = 0)
= FZ (0) − FZ (0− )
Z∞
Z∞ Z∞
Z∞
=
···
···
FY (ξ1 ,...,ξm ,ξm+2 ,...,ξM ) (0)
−∞

−∞ −∞

−∞

fX (ξ1 ) · · · fX (ξm )fX (ξm+2 ) · · · fX (ξM )
dξ1 · · · dξm dξm+2 · · · dξM
Z∞
Z∞ Z∞
Z∞
lim−
···
···
FY (ξ1 ,...,ξm ,ξm+2 ,...,ξM ) (z)

−

ξ1

z→0

−∞

−∞ −∞

−∞

fX (ξ1 ) · · · fX (ξm )fX (ξm+2 ) · · · fX (ξM )
dξ1 · · · dξm dξm+2 · · · dξM

Z∞

Z∞ Z∞
···

=
−∞

−∞ −∞

Z∞
···

FY (ξ1 ,...,ξm ,ξm+2 ,...,ξM ) (0)

−∞


− FY (ξ1 ,...,ξm ,ξm+2 ,...,ξM ) (0− ) fX (ξ1 ) · · · fX (ξm )
fX (ξm+2 ) · · · fX (ξM )dξ1 · · · dξm dξm+2 · · · dξM

FY (ξ1 ,...,ξm ,ξm+2 ,...,ξM ) (0) − FY (ξ1 ,...,ξm ,ξm+2 ,...,ξM ) (0− ) = 0.
Hence,
Pr(Z = 0).

ACKNOWLEDGMENT
Li Si-Yao thank Dongwei Ren for his generous guidance
and constructive conversations.
We thank Junfeng Li, Yongjian Hu and Wenyi Zeng for
their suggestions and inspirations on the proof to Theorem 2.
R EFERENCES
[1] T. F. Chan and C.-K. Wong, “Total variation blind deconvolution,” IEEE
Trans. Image Process., vol. 7, no. 3, pp. 370–375, 1998.
[2] R. Fergus, B. Singh, A. Hertzmann, S. T. Roweis, and W. T. Freeman,
“Removing camera shake from a single photograph,” in ACM Trans.
Graph., vol. 25, no. 3, 2006, pp. 787–794.
[3] S. Cho and S. Lee, “Convergence analysis of map based blur kernel
estimation,” arXiv preprint arXiv:1611.07752, 2016.
[4] L. Xu and J. Jia, “Two-phase kernel estimation for robust motion
deblurring,” in Proc. ECCV, 2010, pp. 157–170.
[5] D. Gong, M. Tan, Y. Zhang, A. Van den Hengel, and Q. Shi, “Blind
image deconvolution by automatic gradient activation,” in Proc. CVPR,
2016, pp. 1827–1836.
[6] S. Cho and S. Lee, “Fast motion deblurring,” ACM Trans. Graph.,
vol. 28, no. 5, p. 145, 2009.
[7] L. Xu, S. Zheng, and J. Jia, “Unnatural l0 sparse representation for
natural image deblurring,” in Proc. CVPR, 2013, pp. 1107–1114.
[8] D. Krishnan, T. Tay, and R. Fergus, “Blind deconvolution using a
normalized sparsity measure,” in Proc. CVPR, 2011, pp. 233–240.
[9] Q. Shan, J. Jia, and A. Agarwala, “High-quality motion deblurring from
a single image,” ACM Trans. Graph., vol. 27, no. 3, p. 73, 2008.
[10] J. Pan, Z. Lin, Z. Su, and M.-H. Yang, “Robust kernel estimation with
outliers handling for image deblurring,” in Proc. CVPR, 2016, pp. 2800–
2808.
[11] W. Zuo, D. Ren, S. Gu, L. Lin, and L. Zhang, “Discriminative learning
of iteration-wise priors for blind deconvolution,” in Proc. CVPR, 2015,
pp. 3232–3240.
[12] D. Perrone and P. Favaro, “Total variation blind deconvolution: The devil
is in the details,” in Proc. CVPR, 2014, pp. 2909–2916.
[13] D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE
Signal Processing Mag., vol. 13, no. 3, pp. 43–64, 1996.
[14] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Understanding and
evaluating blind deconvolution algorithms,” in Proc. CVPR, 2009, pp.
1964–1971.
[15] H. C. Andrews and B. R. Hunt, Digital image restoration. Englewood
Cliffs, NJ: Prentice-Hall, 1977, ch. 5.2, pp. 102–103.

















IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. , NO. , DATE

[16] R. M. Gray, “Toeplitz and circulant matrices: A review,” Foundations
and Trends in Communication and Information Theory, vol. 2, no. 3,
pp. 155–239, 2006.
[17] E. Kaltofen and A. Lobo, “On rank properties of toeplitz matrices over
finite fields,” in Proc. Int. Symp. Symbolic and Algebraic Computation
(ISSAC), 1996, pp. 241–249.
[18] B. A. Olshausen and D. J. Field, “Emergence of simple-cell receptive
field properties by learning a sparse code for natural images,” Nature,
vol. 381, no. 6583, p. 607, 1996.
[19] Y. Weiss and W. T. Freeman, “What makes a good model of natural
images?” in Proc. CVPR, 2007, pp. 1–8.
[20] D. Krishnan and R. Fergus, “Fast image deconvolution using hyperlaplacian priors,” in Proc. NIPS, 2009, pp. 1033–1041.
[21] D. Spielman, “Lecture 2: The condition number,” MIT OpenCourseWare: Behavior of Algorithms, vol. 1, no. 2, p. 3, 2002.
[22] ——, “Lecture 3: The largest singular value of a matrix,” MIT OpenCourseWare: Behavior of Algorithms, vol. 1, no. 1, p. 1, 2002.
[23] M. Fazel, H. Hindi, and S. P. Boyd, “Log-det heuristic for matrix
rank minimization with applications to hankel and euclidean distance
matrices,” in Proc. American Control Conf. (ACC), vol. 3, 2003, pp.
2156–2162.
[24] W. Dong, G. Shi, X. Li, Y. Ma, and F. Huang, “Compressive sensing via
nonlocal low-rank regularization,” IEEE Trans. Image Process., vol. 23,
no. 8, pp. 3618–3632, 2014.
[25] H. Li, “On calculating the determinants of toeplitz matrices,” J. Appl.
Math. Bioinformatics, vol. 1, no. 1, p. 55, 2011.

Li Si-Yao received the B.Sc. degree in computer
science from Beijing Normal University, Beijing,
China, in 2016. He is currently pursuing the M.E. degree in computer science at Beijing Normal University. His research interest includes computer vision
and neural networks.

Qian Yin received her Ph.D in mathematic from
Beijing Normal University in 2006. She is currently
an Associate Professor in College of Information
Science and Technology, and the director of Image
Processing and Pattern Recognition Laboratory in
Beijing Normal University. Her current research
interest includes low-level vision and intelligent processing of astronomical data.

10

Ping Guo (SM’05) received his B.Sc. degree and
M.Sc. degree in optics from Peking University in
1980 and 1983 respectively, and the Ph.D degree in
computer science from The Chinese University of
Hong Kong in 2001.
From 1993 to 1994 he was with the Department of
Computer Science and Engineering at Wright State
University as a visiting faculty. From May to August
2000 he was with the National Laboratory of Pattern
Recognition at Chinese Academy of Sciences as a
research associate. From July to December 2005
he was with the Computer Science Department at The Chinese University
of Hong Kong as a research associate. He is currently a Professor with
the School of Systems Science at Beijing Normal University, and also a
Professor with the School of Computer Science and Technology at Beijing
Institute of Technology. He has authored over 260 papers in important
academic journals and conferences, including AAAI, the International Joint
Conference on Neural Networks, the IEEE TRANSACTIONS ON SYSTEMS,
MAN , AND CYBERNETICS : SYSTEMS and the IEEE TRANSACTIONS ON
NEURAL NETWORKS. His current research interests include computational
intelligence application on image processing, pattern recognition, software
reliability engineering, optical computing and spectral analysis.

