2138 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO.

6, JUNE 2013

Analysis Operator Learning and its Application to

Image Reconstruction
Simon Hawe, Martin Kleinsteuber, and Klaus Diepold

Abstract— Exploiting a priori known structural information sampling process. In many cases, reconstructing s by simply
lies at the core of many image reconstruction methods that can be inverting Equation (1) is ill-posed because either the exact
stated as inverse problems. The synthesis model, which assumes measurement process and hence A is unknown as in blind
that images can be decomposed into a linear combination of very
few atoms of some dictionary, is now a well established tool for image deconvolution, or the number of observations is much
the design of image reconstruction algorithms. An interesting smaller compared to the dimension of the signal, which is
alternative is the analysis model, where the signal is multiplied the case in Compressive Sensing or image inpainting. To
by an analysis operator and the outcome is assumed to be sparse. overcome the ill-posedness and to stabilize the solution, prior
This approach has only recently gained increasing interest. The knowledge or assumptions about the general statistics of
quality of reconstruction methods based on an analysis model
severely depends on the right choice of the suitable operator. images can be exploited.
In this paper, we present an algorithm for learning an analysis
operator from training images. Our method is based on  p -norm B. Synthesis Model and Dictionary Learning
minimization on the set of full rank matrices with normalized
columns. We carefully introduce the employed conjugate gradient One assumption that has proven to be successful in image
method on manifolds, and explain the underlying geometry of the reconstruction, cf. [6], is that natural images admit a sparse
constraints. Moreover, we compare our approach to state-of-the- representation x ∈ Rd over some dictionary D ∈ Rn×d
art methods for image denoising, inpainting, and single image with d ≥ n. A vector x is called sparse when most of
super-resolution. Our numerical results show competitive per-
formance of our general approach in all presented applications its coefficients are equal to zero or small in magnitude.
compared to the specialized state-of-the-art techniques. When s admits a sparse representation over D, it can be
expressed as a linear combination of only very few columns
Index Terms— Analysis operator learning, geometric
conjugate gradient, image reconstruction, inverse problems, of the dictionary {di }di=1 , called atoms, which reads as
oblique manifold.
s = Dx. (2)

I. I NTRODUCTION For d > n, the dictionary is said to be overcomplete or

redundant.
A. Problem Description
Now, using the knowledge that (2) allows a sparse solution,
Linear inverse problems are ubiquitous in the field of image an estimation of the original signal in (1) can be obtained from
processing. Prominent examples are image denoising [1], the measurements y by first solving
inpainting [2], super-resolution [3], or image reconstruction
from few indirect measurements as in Compressive Sensing x = arg min g(x) subject to ADx − y22 ≤ , (3)
x∈Rd
[4]. Basically, in all these problems the goal is to reconstruct
an unknown image s ∈ Rn as accurately as possible from a set and afterwards synthesizing the signal from the computed
of indirect and may be corrupted measurements y ∈ Rm with sparse coefficients via s = Dx . Therein, g : Rd → R is
n ≥ m, see [5] for a detailed introduction to inverse problems. a function that promotes or measures sparsity, and  ∈ R+ is
Formally, this measurement process can be written as an estimated upper bound on the noise power e22 . Common
choices for g include the  p -norm
y = As + e, (1) 
p
v p := |v i | p , (4)
where the vector e ∈ Rm
models sampling errors and noise,
i
and A ∈ Rm×n is the measurement matrix modeling the
with 0 < p ≤ 1 and differentiable approximations of (4).
Manuscript received April 23, 2012; revised December 14, 2012; accepted As the signal is synthesized from the sparse coefficients, the
January 22, 2013. Date of publication February 11, 2013; date of current
version March 29, 2013. This work was supported by the Cluster of Excellence reconstruction model (3) is called the synthesis reconstruction
CoTeSys – Cognition for Technical Systems, funded by the German Research model [7].
Foundation. The associate editor coordinating the review of this manuscript To find the minimizer of Problem (3), various algorithms
and approving it for publication was Prof. Joseph P. Havlicek.
The authors are with the Department of Electrical Engineering, Technische based on convex or non-convex optimization, greedy pursuit
Universität München, Munich 80290, Germany (e-mail: [Link]@[Link]; methods, or Bayesian frameworks exist that may employ dif-
kleinsteuber@[Link]; kldi@[Link]). ferent choices of g. For a broad overview of such algorithms,
Color versions of one or more of the figures in this paper are available
online at [Link] we refer the interested reader to [8]. What all these algorithms
Digital Object Identifier 10.1109/TIP.2013.2246175 have in common, is that their performance regarding the
1057-7149/$31.00 © 2013 IEEE
HAWE et al.: ANALYSIS OPERATOR LEARNING AND ITS APPLICATION TO IMAGE RECONSTRUCTION 2139

reconstruction quality severely depends on an appropriately the analyzed vector s describe the subspace containing the
chosen dictionary D. Ideally, one is seeking for a dictionary signal. To emphasize this difference, the term cosparsity has
where s can be represented most accurately with a coefficient been introduced in [19], which simply counts the number of
vector x that is as sparse as possible. Basically, dictionaries zero elements of s. As the sparsity in the synthesis model
can be assigned to two major classes: analytic dictionaries and depends on the chosen dictionary, the cosparsity of an analyzed
learned dictionaries. signal depends on the choice of the analysis operator .
Analytic dictionaries are built on mathematical models of Different analysis operators proposed in the literature
a general type of signal, e.g. natural images, they should include the fused Lasso [20], the translation invariant wavelet
represent. Popular examples include Wavelets [9], Bandlets transform [21], and probably best known the finite difference
[10], and Curvlets [11] among several others, or a con- operator which is closely related to the total-variation [22].
catenation of various such bases/dictionaries. They offer the They all have shown very good performance when used within
advantages of low computational complexity and of being the analysis model for solving diverse inverse problems in
universally applicable to a wide set of signals. However, this imaging. The question is: Can the performance of analysis
universality comes at the cost of not giving the optimally based signal reconstruction be improved when a learned
sparse representation for more specific classes of signals, e.g. analysis operator is applied instead of a predefined one, as it
face images. is the case for the synthesis model where learned dictionaries
It is now well known that signals belonging to a specific outperform analytic dictionaries? In [7], it has been discussed
class can be represented with fewer coefficients over a dictio- that the two models differ significantly, and the naïve way
nary that has been learned using a representative training set, of learning a dictionary and simply employing its transposed
than over analytic dictionaries. This is desirable for various or its pseudo-inverse as the learned analysis operator fails.
image reconstruction applications as it readily improves their Hence, different algorithms are required for analysis operator
performance and accuracy [12]–[14]. Basically, the goal is to learning.
find a dictionary over which a training set admits a maxi-
mally sparse representation. In contrast to analytic dictionaries,
D. Contributions
which can be applied globally to an entire image, learned
dictionaries are small dense matrices that have to be applied In this work, we introduce a new algorithm based on geo-
locally to small image patches. Hence, the training set consists metric optimization for learning a patch based analysis oper-
of small patches extracted from some example images. This ator from a set of training samples, which we name GOAL1
restriction to patches mainly arises from limited memory, and (GeOmetric Analysis operator Learning). The method relies
limited computational resources. on a minimization problem, which is carefully motivated in
Roughly speaking, starting from some initial dictionary the Section II-B. Therein, we also discuss the question of what is
learning algorithms iteratively update the atoms of the dictio- a suitable analysis operator for image reconstruction, and how
nary, such that the sparsity of the training set is increased. This to antagonize overfitting the operator to a subset of the training
procedure is often performed via block-coordinate relaxation, samples. An efficient geometric conjugate gradient method
which alternates between finding the sparsest representation on the so-called oblique manifold is proposed in Section III
of the training set while fixing the atoms, and optimizing for learning the analysis operator. Furthermore, in Section IV
the atoms that most accurately reproduce the training set we explain how to apply the local patch based analysis
using the previously determined sparse representation. Three operator to achieve global reconstruction results. Section V
conceptually different approaches for learning a dictionary sheds some light on the influence of the parameters required
became well established, which are probabilistic ones like by GOAL and how to select them, and compares our method
[15], clustering based ones such as K-SVD [16], and recent to other analysis operator learning techniques. The quality of
approaches which aim at learning dictionaries with specific the operator learned by GOAL on natural image patches is
matrix structures that allow fast computations like [17]. For further investigated in terms of image denoising, inpainting,
a comprehensive overview of dictionary learning techniques and single image super-resolution. The numerical results show
see [18]. the broad and effective applicability of our general approach.

C. Analysis Model E. Notations

An alternative to the synthesis model (3) for reconstructing Matrices are written as capital calligraphic letters like X ,
a signal, is to solve column vectors are denoted by boldfaced small letters e.g. x,
s = arg min g(s) subject to As − y22 ≤ , (5) whereas scalars are either capital or small letters like n, N. By
s∈Rn v i we denote the i th element of the vector v, v i j denotes the
i th element in the j th column of a matrix V. The vector v:,i
which is known as the analysis model [7]. Therein,  ∈ Rk×n
denotes the i th column of V whereas vi,: denotes the transposed
with k ≥ n is called the analysis operator, and the analyzed
of the i th row of V. By Ei j , we denote a matrix whose i th entry
vector s ∈ Rk is assumed to be sparse, where sparsity is
in the j th column is equal to one, and all others are zero. Ik
again measured via an appropriate function g. In contrast to the
synthesis model, where a signal is fully described by the non- 1 Matlab implementation available at [Link]
zero elements of x, in the analysis model the zero elements of w00bhl/www/[Link]
2140 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO. 6, JUNE 2013

denotes the identity matrix of dimension (k × k), 0 denotes set of possible analysis operators is restricted to the set of
the zero-matrix of appropriate dimension, and ddiag(V) is the Uniform Normalized Tight Frames, i.e. matrices with uniform
 whose entries on the diagonal are those of V.
diagonal matrix row norm and orthonormal columns. The authors state that
By V2F = i, j v i2j we denote the squared Frobenius norm this algorithm has the limitation of requiring noiseless training
of a matrix V, tr(V) is the trace of V, and rk(V) denotes the samples whose analyzed vectors {si }i=1M are exactly cosparse.

rank. To overcome this restriction, the same authors propose an

extension of this algorithm that simultaneously learns the
II. A NALYSIS O PERATOR L EARNING analysis operator and denoises the training samples, cf. [26].
A. Prior Art This is achieved by alternating between updating the analysis
operator via the projected subgradient algorithm and denoising
The topic of analysis operator learning has only recently
the samples using an Augmented Lagrangian method. Therein,
started to be investigated, and only few prior work exists. In the
the authors state that their results for image denoising using
sequel, we shortly review analysis operator learning methods
the learned operator are only slightly worse compared to
that are applicable for image processing tasks. M employing the commonly used finite difference operator.
Given a set of M training samples si ∈ Rn i=1 ,
An interesting idea related to the analysis model, called
the goal of analysis operator learning is to find a matrix
Fields-of-Experts (FoE) has been proposed in [27]. The
 ∈ Rk×n with k ≥ n, which leads to a maximally cosparse
method relies on learning high-order Markov Random Field
representation si of each training sample. As mentioned in
image priors with potential functions extending over large
Subsection I-B, the training samples are distinctive vectorized
pixel neighborhoods, i.e. overlapping image patches. Moti-
image patches extracted from a set of example images. Let
vated by a probabilistic model, they use the student-t dis-
S = [s1 , . . . , s M ] ∈ Rn×M be a matrix where the training
tribution of several linear filter responses as the potential
samples constitute its columns, then the problem is to find
function, where the filters, which correspond to atoms from
 = arg min G(S) (6) an analysis operator point of view, have been learned from
 training patches. Compared to our work and the methods
where  is subject to some constraints, and G is some explained above, an underdetermined operator is used in their
function that measures the sparsity of the matrix S. In [23], experiments, i.e. k < n. Moreover, the learning process of
an algorithm is proposed in which the rows of the analysis FoE is computationally demanding, which limits the number
operator are found sequentially by identifying directions that and the size of trained filters. Note also, that these filters are
are orthogonal to a subset of the training samples. Starting learned independently of each other, while in contrast GOAL
from a randomly initialized vector ω ∈ Rn , a candidate row updates the analysis operator as a whole.
is found by first computing the inner product of ω with the
entire training set, followed by extracting the reduced training
set S R of samples whose inner product with ω is smaller than a B. Motivation of Our Approach
threshold. Thereafter, ω is updated to be the eigenvector corre- In the quest for designing an analysis operator learning
sponding to the smallest eigenvalue of S R S R . This procedure algorithm, the natural question arises: What is a good analysis
is iterated several times until a convergence criterion is met. If operator for our needs? Clearly, given a signal s that belongs
the determined candidate vector is sufficiently distinctive from to a certain signal class, the aim is to find an  such that s is
already found ones, it is added to  as a new row, otherwise it as sparse as possible. This motivates to minimize the expected
is discarded. This process is repeated until the desired number sparsity E[g(s)]. All approaches presented in Subsection II-
k of rows have been found. A can be explained in this way, i.e. for their sparsity measure
An adaption of the widely known K-SVD dictionary learn- g they aim at learning an  that minimizes the empirical mean
ing algorithm to the problem of analysis operator learning of the sparsity over all randomly drawn training samples. This,
is presented in [24]. As in the original K-SVD algorithm, however, does not necessarily mean to learn the optimal  if
G(S) = i si 0 is employed as the sparsifying function the purpose is to reconstruct several signals belonging to a
and the target cosparsity is required as an input to the diverse class, e.g. natural image patches. The reason for this
algorithm. The arising optimization problem is solved by is that even if the expected sparsity is low, it may happen
alternating between a sparse coding stage over each training with high probability that some realizations of this signal class
sample while fixing  using an ordinary analysis pursuit cannot be represented in a sparse way, i.e. that for a given
method, and updating the analysis operator using the optimized upper bound u, the probability Pr (g(s) ≥ u) exceeds a
training set. Then, each row of  is updated in a similar way tolerable value, cf. Fig. 1.
as described in the previous paragraph for the method of [23]. The algorithm presented here aims at minimizing the empir-
Interestingly, the operator learned on piecewise constant image ical expectation of a sparsifying function g(si ) for all
patches by [23] and [24] closely mimics the finite difference training samples si , while additionally keeping the empirical
operator.  variance moderate. In other words, we try to avoid that the
In [25], the authors use G(S) = i si 1 as the sparsity analyzed vectors of many similar training samples become
promoting function and suggest a constrained optimization very sparse and consequently prevent  from being adapted to
technique that utilizes a projected subgradient method for the remaining ones. For image processing, this is of particular
iteratively solving (6). To exclude the trivial solution, the interest if the training patches are chosen randomly from
HAWE et al.: ANALYSIS OPERATOR LEARNING AND ITS APPLICATION TO IMAGE RECONSTRUCTION 2141

1 Note, that this definition only yields a non-empty set if k ≥ n,

P r(g(Ωi s) ≤ x)
which is the interesting case in this work. Thus, from now on,
we assume k ≥ n. Remember that we require the rows of  to
have unit Euclidean norm. Hence, we restrict the transposed
0.5
of the learned analysis operator to be an element of OB(n, k).
i=2
Since OB(n, k) is open and dense in the set of matrices with
i=1 normalized columns, we need a penalty function that ensures
i=1
0 g3 g2 g1 the rank constraint (ii) and prevents iterates to approach the
u g1 g2 g1 u boundary of OB(n, k).
Lemma 1: The inequality 0 < det( 1k X X  ) ≤ ( n1 )n holds
Fig. 1. Illustration of two possible distributions Pr(g(i s) ≤ x) for two true for all X ∈ OB(n, k), where 1 < n ≤ k.
analysis operators. 1 : low expectation g 1 , high variance (dashed line); 2 :
moderate expectation g 2 , moderate variance. Although 1 yields a smaller
Proof: Due to the full rank condition on X , the product
expectation, there are more signals compared to 2 where the sparsity model X X  is positive definite, consequently the strict inequality
fails, i.e., Pr(g(1 s) ≥ u) > Pr(g(2 s) ≥ u) for a suitable upper bound u. 0 < det( 1k X X  ) applies. To see the second inequality of
Lemma 1, observe that
natural images, because there is a high probability of collecting
a large subset of very similar patches, e.g. homogeneous X 2F = tr(X X  ) = k, (12)
regions, that bias the learning process.
Concretely, we want to find an  that minimizes both the which implies tr( 1k X X  ) = 1. Since the trace of a matrix is
squared empirical mean equal to the sum of its eigenvalues, which are strictly positive
  2 in our case, it follows that the strict inequality 0 < λi < 1
g2 = M1 g(si ) (7) holds true for all eigenvalues λi of k1 X X  . From the well
i known relation between the arithmetic and the geometric mean
and the empirical variance we see
  
σg2 = M1
(g(si ) − g)2 (8)
n
λi ≤ n1 λi . (13)
i
Now, since the determinant ofa matrix is equal to the product
of the sparsity of the analyzed vectors. We achieve this by of its eigenvalues, and with λi = tr( 1k X X  ) = 1, we have
minimizing the sum of both, which is readily given by
 det( 1k X X  ) = λi ≤ ( n1 )n , (14)
g 2 + σg2 = M1 g(si )2 . (9)
i which completes the proof.
p Recalling that  ∈ OB(n, k) and considering Lemma 1,
Using g(·) = p · p ,
1
and introducing the factor 1
2 the function we can enforce the full rank constraint with the penalty
we employ reads as function

M
 
k
 
M
 p 2
J p (V) := |v j i | p 2
= h() := − n log(n)
1
log det( 1k  ). (15)
p v:,i  p ,(10)
1 1 1 1
2M p 2M
i=1 j =1 i=1 Regarding Condition (iii), the following result proves useful.
with 0 ≤ p ≤ 1 and V = S. Lemma 2: For a matrix X ∈ OB(n, k) with 1 < n ≤ k, the
 x | ≤ 1 applies, where equality holds true if
inequality |x:,i
Certainly, without additional prior assumptions on , the :, j
useless solution  = 0 is the global minimizer of Problem (6). and only if x:,i = ±x:, j .
To avoid the trivial solution and for other reasons explained Proof: By the definition of OB(n, k) the columns of
later in this section, we regularize the problem by imposing X are normalized, consequently Lemma 2 follows directly
the following three constraints on . from the Cauchy-Schwarz inequality.
1) The rows of  have unit Euclidean norm, i.e. ωi,: 2 = Thus, Condition (iii) can be enforced via the logarithmic
1 for i = 1, . . . , k. barrier function of the scalar products between all distinctive
2) The analysis operator  has full rank, i.e. rk() = n. rows of , i.e.
3) The analysis operator  does not have linear dependent 

r () := − log(1 − (ωi,: ω j,: )2 ). (16)
rows, i.e. ωi,: = ±ω j, for i = j .
1≤i< j ≤k
The rank condition (ii) on  is motivated by the fact that
different input samples s1 , s2 ∈ Rn with s1 = s2 should Finally, combining all the introduced constraints, our opti-
be mapped to different analyzed vectors s1 = s2 . With mization problem for learning the transposed analysis operator
Condition (iii) redundant transform coefficients in an analyzed reads as
vector are avoided.  = arg min J p (X  S) + κ h(X  ) + μ r (X  ). (17)
These constraints motivate the consideration of the set of X ∈OB(n,k)
full rank matrices with normalized columns, which admits a
Therein, the two weighting factors κ, μ ∈ R+ control the
manifold structure known as the oblique manifold [28]
influence of the two constraints on the final solution. The
OB(n, k) := {X ∈ Rn×k | rk(X ) = n, ddiag(X  X ) = Ik }. (11) following lemma clarifies the role of κ.
2142 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO. 6, JUNE 2013

Lemma 3: Let  be a minimum of h in the set of transposed

oblique matrices, i.e.
 ∈ arg min h(X  ), (18)
X ∈OB(n,k)
then the condition number of  is equal to one.
Proof: It is well known that equality of the arithmetic and
the geometric mean in Equation (13) holds true, if and only
if all eigenvalues λi of 1k X X  are equal, i.e. λ1 = . . . = λn .
Hence, if  ∈ arg min h(X  ), then det( 1k  ) = ( n1 )n ,
X ∈OB(n,k)
and consequently all singular values of  coincide. This
implies that the condition number of , which is defined as
the quotient of the largest to the smallest singular value, is
equal to one.
With other words, the minima of h are uniformly normalized
tight frames, cf. [25], [26]. From Lemma 3 we can conclude
that with larger κ the condition number of  approaches one.
Now, recall the inequality Fig. 2. This figure shows two points X and Y on a manifold M together
with their tangent spaces TX M and TY M. Furthermore, the Euclidean
σmin s1 − s2 2 ≤ (s1 − s2 )2 ≤ σmax s1 − s2 2 , (19) gradient ∇ f (X ) and its projection onto the tangent space TX M (∇ f (X )) are
depicted. The geodesic (X , H, t) in the direction of H ∈ TX M connecting
with σmin being the smallest and σmax being the largest the two points is shown. The dashed line typifies the role of a parallel transport
singular value of . From this it follows that an analysis of the gradient in TX M to TY M.
operator found with a large κ, i.e. obeying σmin ≈ σmax , carries
over distinctness of different signals to their analyzed versions.
The parameter μ regulates the redundancy between the rows properties, see e.g. [29] for a CG-type method on the oblique
of the analysis operator and consequently avoids redundant manifold.
coefficients in the analyzed vector s. To make this work self contained, we start by shortly
Lemma 4: The difference between any two entries of the reviewing the general concepts of optimization on matrix
analyzed vector s is bounded by manifolds. After that we present the concrete formulas and
implementation details for our optimization problem on the

|ωi,: s − ωj,: s| ≤  ω ) s .
2(1 − ωi,: j,: 2 (20) oblique manifold. For an in-depth introduction on optimization
Proof: From the Cauchy-Schwarz inequality we get on matrix manifolds, we refer the interested reader to [30].

|ωi,: s − ωj,: s| = |(ωi,: − ω j,: ) s| ≤ ωi,: − ω j,: 2 s2 .(21)
A. Optimization on Matrix Manifolds
Since by definition ωi,: 2 = ω j,: 2 = 1, it follows that
Let M be a smooth Riemannian submanifold of Rn×k with
ωi,: − ω j,: 2 =  ω ).
2(1 − ωi,: j,: the standard Frobenius inner product Q, P := tr(Q P),
The above lemma implies, that if the i th entry of the ana- and let f : Rn×k → R be a differentiable cost function. We
lyzed vector is significantly larger than 0 then a large absolute consider the problem of finding
ω
value of ωi,: th
j,: prevents the j entry to be small. To achieve
large cosparsity, this is an unwanted effect that our approach arg min f (X ). (22)
X ∈M
avoids via the log-barrier function r in (16). It is worth
mentioning that the same effect is achieved by minimizing the The concepts presented in this subsection are visualized in
analysis operator’s mutual coherence max |ωi,:  ω | and that Fig. 2 to alleviate the understanding.
j,:
i = j To every point X ∈ M one can assign a tangent space
our experiments suggest that enlarging μ leads to minimizing TX M. The tangent space at X is a real vector space containing
the mutual coherence. all possible directions that tangentially pass through X . An
In the next section, we explain how the manifold structure element ∈ TX M is called a tangent vector at X . Each
of OB(n, k) can be exploited to efficiently learn the analysis tangent space is associated with an inner product inherited
operator. from the surrounding Rn×k , which allows to measure distances
and angles on M.
III. A NALYSIS O PERATOR L EARNING A LGORITHM The Riemannian gradient of f at X is an element of the
Knowing that the feasible set of solutions to Problem (17) tangent space TX M that points in the direction of steepest
is restricted to a smooth manifold allows us to formulate ascent of the cost function on the manifold. As we require
a geometric conjugate gradient (CG-) method to learn the M to be a submanifold of Rn×k and since by assumption
analysis operator. Geometric CG-methods have been proven f is defined on the whole Rn×k , the Riemannian gradient
efficient in various applications, due to the combination of G(X ) is simply the orthogonal projection of the (standard)
moderate computational complexity and good convergence gradient ∇ f (X ) onto the tangent space TX M. In formulas,
HAWE et al.: ANALYSIS OPERATOR LEARNING AND ITS APPLICATION TO IMAGE RECONSTRUCTION 2143

this reads as current optimal point until some user-specified convergence

criterion is met, or a maximum number of iterations has been
G(X ) := TX M (∇ f (X )). (23)
reached.
A geodesic is a smooth curve (X , , t) emanating from
X in the direction of ∈ TX M, which locally describes the B. Geometric Conjugate Gradient for Analysis Operator
shortest path between two points on M. Intuitively, it can be Learning
interpreted as the equivalent of a straight line in the manifold In this subsection we derive all ingredients to implement
setting. the geometric conjugate gradient method as described in the
Conventional line search methods search for the next iterate previous subsection for the task of learning the analysis oper-
along a straight line. This is generalized to the manifold setting ator. Results regarding the geometry of OB(n, k) are derived
as follows. Given a current optimal point X (i) and a search e.g. in [30]. To enhance legibility, and since the dimensions
direction H(i) ∈ TX (i) M at the i th iteration, the step size α (i) n and k are fixed throughout the rest of the paper, the oblique
which leads to sufficient decrease of f can be determined by manifold is further on denoted by OB.
finding the minimizer of The tangent space at X ∈ OB is given by
α (i) = arg min f ( (X (i) , H(i) , t)). (24) TX OB = { ∈ Rn×k | ddiag(X  ) = 0}. (31)
t ≥0

Once α (i) has been determined, the new iterate is computed The orthogonal projection of a matrix Q ∈ Rn×k onto the
by tangent space TX OB is

X (i+1) = (X (i) , H(i) , α (i) ). (25) TX OB (Q) = Q − X ddiag(X  Q). (32)

Now, one straightforward approach to minimize f is to Regarding geodesics, note that in general a geodesic is
alternate Equations (23), (24), and (25) using H(i) = −G (i) , the solution of a second order ordinary differential equation,
with the short hand notation G (i) := G(X (i) ), which corre- meaning that for arbitrary manifolds, its computation as well
sponds to the steepest descent on a Riemannian manifold. as computing the parallel transport is not feasible. Fortunately,
However, as in standard optimization, steepest descent only as the oblique manifold is a Riemannian submanifold of a
has a linear rate of convergence. Therefore, we employ a product of k unit spheres S n−1 , the formulas for parallel
conjugate gradient method on a manifold, as it offers a transport and the exponential mapping allow an efficient
superlinear rate of convergence, while still being applicable implementation.
to large scale optimization problems with low computational Let x ∈ S n−1 be a point on a sphere and h ∈ Tx S n−1 be a
complexity. tangent vector at x, then the geodesic in the direction of h is
In CG-methods, the updated search direction H(i+1) ∈ a great circle
TX (i+1) M is a linear combination of the gradient G (i+1) ∈ x, if h2 = 0
TX (i+1) M and the previous search direction H(i) ∈ TX (i) M. γ (x, h, t) = sin(t h2 ) (33)
x cos(th2 ) + h h2 , otherwise.
Since adding vectors that belong to different tangent spaces
is not defined, we need to map H(i) from TX (i) M to The associated parallel transport of a tangent vector ξ ∈
TX (i+1) M. This is done by the so-called parallel transport Tx S n−1 along the great circle γ (x, h, t) reads as
T ( , X (i) , H(i) , α (i) ), which transports a tangent vector ∈
ξ h
TX (i) M along the geodesic (X (i) , H(i) , t) to the tangent space τ (ξ , x, h, t) = ξ − xh2 sin(th2 )
TX (i+1) M. Now, using the shorthand notation h22

T (i+1) := T ( , X (i) , H(i) , α (i) ), (26) +h(1 − cos(th2 )) . (34)

the new search direction is computed by As OB is a submanifold of the product of unit spheres,
H (i+1)
= −G (i+1) (i)
+ β TH(i) , (i+1)
(27) the geodesic through X ∈ OB in the direction of H ∈
TX OB is simply the combination of the great circles emerging
where β (i) ∈ R is calculated by some update formula adopted by concatenating each column of X with the corresponding
to the manifold setting. Most popular are the update formulas column of H, i.e.
by Fletcher-Reeves (FR), Hestenes-Stiefel (HS), and Dai-Yuan 
(i+1) (X , H, t) = γ (x:,1 , h:,1 , t), . . . , γ (x:,k , h:,k , t) . (35)
(DY). With Y (i+1) = G (i+1) − TG (i) , they read as
(i+1) (i+1)
(i)
βFR = G G (i) ,,G , (28) Accordingly, the parallel transport of ∈ TX OB along the
G (i)
(i+1) (i+1)
geodesic (X , H, t) is given by
(i)
βHS = G (i+1),Y (i+1) , (29)
T ,Y T ( , X , H, t) =
H(i) 
(i) G (i+1) ,G (i+1) (36)
βDY = . (30) τ (ξ:,1 , x:,1 , h:,1 , t), . . . , τ (ξ:,k , x:,k , h:,k , t) .
T (i+1)
(i) ,Y
(i+1)
H

Now, a solution to Problem (22) is computed by alternating Now, to use the geometric CG-method for learning the
between finding the search direction on M and updating the analysis operator, we require a differentiable cost function
2144 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO. 6, JUNE 2013

Algorithm 1 Backtracking Line Search on Oblique Manifold Algorithm 2 Geometric Analysis Operator Learning (GOAL)
Input:
(i)
t0 > 0, 0 < c1 < 1, 0 < c2 < 0.5, X (i) , G (i) , H(i) Input: Initial analysis operator init, training data S, parame-
(i)
Set: t ← t0 ters p, ν, κ, μ
while f ( (X (i) , H(i) , t)) > f (X (i) ) + tc2 G (i) , H(i) do Set: i ← 0, X (0) ←  init, H
(0) ← −G (0)

t ← c1 t repeat
end while α (i) ← arg min f ( (X (i) , H(i) , t)), cf. Algorithm 1 in
t ≥0
Output: α (i) ← t conjunction with Equation (41)
X (i+1) ← (X (i) , H(i) , α (i) ), cf. Equation (35)
G (i+1) ← TX (i+1) M (∇ f (X (i+1) )), cf. Equations (32) and
f . Since, the cost function presented in Problem (17) is (42)
not differentiable due to the non-smoothness of the ( p, q)-  (i) (i) 
β (i) ← max 0, min(βDY , βHS ) , cf. Equations (29), (30)
pseudo-norm (10), we exchange Function (10) with a smooth
approximation, which is given by H(i+1) ← −G (i+1) + β (i) TH(i+1)
(i) , cf. Equations (26), (36)
 k 2 i ←i +1

M 
J p,ν (V) := 2M
1 1 p
(v i2j + ν) 2 , (37) until X (i) −X (i−1)  F < 10−4 ∨ i = maximum # iterations
p
j =1 i=1
Output:  ← X (i)

with ν ∈ R+ being the smoothing parameter. The smaller ν

is, the more closely the approximation resembles the original
c2 = 10−2 . As an initial guess for the step size at the first
function. Again, taking V = S and with the shorthand
CG-iteration i = 0, we choose
notation z i j := (S)i j , the gradient of the applied sparsity
(0)
promoting function (37) reads as t0 = G (0) −1
F , (44)

1  1  2
M k
∂ p as proposed in [32]. In the subsequent iterations, the backtrack-
∂ J p,ν (S) = M p (z i j + ν) 2
j =1 i=1 ing line search is initialized by the previous step size divided
(i) (i−1)

k  p
 by the line search parameter, i.e. t0 = α c1 . Our complete
× z i j (z i2j + ν) 2 −1 Ei j S. (38) approach for learning the analysis operator is summarized in
i=1 Algorithm 2. Note, that under the conditions that the Fletcher-
The gradient of the rank penalty term (15) is Reeves update formula is used and some mild conditions on
the step-size selection, the convergence of Algorithm 2 to a
∂
∂ h() = − kn log(n)
2
( 1k  )−1 (39) critical point, i.e. lim inf i→∞ G (i)  = 0, is guaranteed by a
and the gradient of the logarithmic barrier function (16) is result provided in [33].
⎡ ⎤
 2ω ω
⎣ i,: j,: IV. A NALYSIS O PERATOR -BASED I MAGE
∂
∂ r () = (E + E j i )⎦ . (40)
 ω )2 i j
1≤i< j ≤k
1 − (ω i,: j,: R ECONSTRUCTION
In this section we explain how the analysis operator
Combining Equations (38), (39), and (40), the gradient of
 ∈ Rk×n is utilized for reconstructing an unknown image
the cost function
s ∈ R N from some measurements y ∈ Rm following the analy-
f (X ) := J p,ν (X  S) + κ h(X  ) + μ r (X  ) (41) sis approach (5). Here, the vector s ∈ R N denotes a vectorized
image of dimension N = wh, with w being the width and h
which is used for learning the analysis operator reads as being the height of the image, respectively, obtained by
∇ f (X ) = ∂∂X J p,ν (X  S) + κ ∂∂X h(X  ) + μ ∂∂X r (X  ).(42) stacking the columns of the image above each other. In the
following, we will loosely speak of s as the image.
Regarding the CG-update parameter β (i) , we employ a Remember, that the size of  is very small compared to
hybridization of the Hestenes-Stiefel Formula (29) and the Dai the size of the image, and it has to be applied locally to small
Yuan formula (30) image patches rather than globally to the entire image. Arti-
(i)  (i) (i)  facts that arise from naïve patch-wise reconstruction are com-
βhyb = max 0, min(βDY , βHS ) , (43)
monly reduced by considering overlapping patches. Thereby,
which has been suggested in [31]. As explained therein, each patch is reconstructed individually and the entire image
formula (43) combines the good numerical performance of is formed by averaging over the overlapping regions in a
HS with the desirable global convergence properties of DY. final step. However, this method misses global support during
Finally, to compute the step size α (i) , we use an adaption the reconstruction process, hence, it leads to poor inpainting
of the well-known backtracking line search to the geodesic results and is not applicable for e.g. Compressive Sensing
(X (i) , H(i) , t). In that, an initial step size t0(i) is iteratively tasks. To overcome these drawbacks, we use a method related
decreased by a constant factor c1 < 1 until the Armijo to the patch based synthesis approach from [12] and the
condition is met, see Algorithm 1 for the entire procedure. FoE algorithm [27], which provides global support from local
In our implementation we empirically chose c1 = 0.9 and information. Instead of optimizing over each patch individually
HAWE et al.: ANALYSIS OPERATOR LEARNING AND ITS APPLICATION TO IMAGE RECONSTRUCTION 2145

and combining them in a final step, we optimize over the problem is to solve
entire image demanding that a pixel is reconstructed such that
s = arg min 12 As − y22 + b(s) + λg( F s). (49)
the average sparsity of all patches it belongs to is minimized. s∈R N
When all possible patch positions are taken into account, this
procedure is entirely partitioning-invariant. Therein, λ ∈ R+ balances between the sparsity of the solu-
√ √ For legibility,
√ we
tion’s analysis coefficients and the solution’s fidelity to the
assume square patches i.e. of size ( n × n) with n being
a positive integer. measurements. The measurement matrix A ∈ Rm×N and the
Formally, let r ⊆ {1, . . . , h} and c ⊆ {1, . . . , w} denote measurements y ∈ Rm are application dependent.
sets of indices √with ri+1 − ri = dv , ci+1 − ci = dh and
1 ≤ dv , dh ≤ n. Therein, dv , dh determine the degree of V. E VALUATION AND E XPERIMENTS
overlap between two adjacent patches in vertical, and hori- The first part of this section aims at answering the ques-
zontal direction, respectively. We consider all image patches tion of what is a good analysis operator for solving image
whose center is an element of the cartesian product set r × c. reconstruction problems and relates the quality of an analysis
Hence, with | · | denoting the cardinality of a set, the total operator with its mutual coherence and its condition number.
number of patches being considered is equal to |r||c|. Now, This, in turn allows to select the optimal weighting parameters
let Prc be a binary (n × N) matrix that extracts the patch κ and μ for GOAL. Using this parameters, we learn one
centered at position (r, c). With this notation, we formulate general analysis operator  by GOAL, and compare its
the (global) sparsity promoting function as image denoising performance with other analysis approaches.
In the second part, we employ this  unaltered for solving

k
p two classical image reconstruction tasks of image inpainting
(( Prc s)2i + ν) 2 , (45)
and single image super-resolution, and compare our results
r∈r c∈c i=1
with the currently best analysis approach FoE [27], and state-
which measures the overall approximated  p -pseudo-norm of of-the-art methods specifically designed for each respective
the considered analyzed image patches. We compactly rewrite application.
Equation (45) as
K 
 p A. Global Parameters and Image Reconstruction
F 2
g( s) := ( F s)2i + ν , (46) To quantify the reconstruction quality, as usual, weuse
N
i=1 the peak signal-to-noise ratio PSNR = 10 log(2552 N/ i=1

(si − si ) ). Moreover, we measure the quality using the Mean
2
with K = k|r||c| and
⎡ ⎤ Structural SIMilarity Index (MSSIM) [35], with the same set
 Pr1 c1 of parameters as originally suggested in [35]. Compared to
⎢  Pr1 c2 ⎥ PSNR, the MSSIM better reflects a human observer’s visual
⎢ ⎥
 F := ⎢ .. ⎥ ∈ R K ×N (47) impression of quality. It ranges between zero and one, with
⎣ . ⎦
one meaning perfect image reconstruction.
 Pr|r| c|c|
Throughout all experiments, we fixed the size of the image
being the global analysis operator that expands the patch patches to (8 × 8), i.e. n = 64. This is in accordance to
based one to the entire image. We treat image boundary the patch-sizes mostly used in the literature, and yields a
effects by employing constant √ padding, i.e. replicating the good trade-off between reconstruction quality and numerical
values at the image boundaries  2n  times, where · denotes burden. Images are reconstructed by solving the minimization
rounding to the smaller integer. Certainly, for image processing problem (49) via the conjugate gradient method proposed in
applications  F is too large for being applied in terms of [34]. Considering the pixel intensity bounds, we used bl = 0
matrix vector multiplication. Fortunately, applying  F and and bu = 255, which is the common intensity range in
its transposed can be implemented efficiently using sliding 8-bit grayscale image formats. The sparsity promoting function
window techniques, and the matrix vector notation is solely (46) with p = 0.4 and ν = 10−6 is used for both learning the
used for legibility. analysis operator by GOAL, and reconstructing the images.
According to [34], we exploit the fact that the range of pixel Our patch based reconstruction algorithm as explained in
intensities is limited by a lower bound bl and an upper bound Section IV achieves the best results for the maximum possible
bu . We enforce this bounding constraint by minimizing the overlap dh = dv = 1. The Lagrange multiplier λ and the
N measurements matrix A depend on the application, and are
differentiable function b(s) := b(si ), where b is a penalty briefly discussed in the respective subsections.
i=1
term given as
⎧ B. Analysis Operator Evaluation and Parameter Selection
⎨ |s − bu |2 if s ≥ bu
b(s) = |s − bl |2 if s ≤ bl (48) For evaluating the quality of an analysis operator and
⎩ for selecting appropriate parameters for GOAL, we choose
0 otherwise.
image denoising as the baseline experiment. The images to be
Finally, combining the two constraints (46) and (48) with reconstructed have artificially been corrupted by additive white
the data fidelity term, the analysis based image reconstruction Gaussian noise (AWGN) of varying standard deviation σnoise .
2146 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO. 6, JUNE 2013

1.6
0.9

Mutual Coherence
0.8 1.4

0.7 1.2
Fig. 3. Five training images used for learning the analysis operator.
0.6
1
1.2 1.4 1.6 1.8 2 2.2 1 16 32 48 64
Condition Number

This baseline experiment is further used to compare GOAL (a) (b)

with other analysis operator learning methods. We would like Fig. 4. (a) Performance of 400 analysis operators learned by GOAL in
to emphasize that the choice of image denoising as a baseline relation to their mutual coherence and their condition number. Color ranges
experiment is not crucial neither for selecting the learning from dark blue (worst) to dark red (best). The green dot corresponds to the
best performing operator  . (b) Singular values of  .
parameters, nor for ranking the learning approaches. In fact,
any other reconstruction task as discussed below leads to the
same parameters and the same ranking of the different learning
algorithms.
For image denoising, the measurement matrix A in Equation
(49) is the identity matrix. As it is common in the denoising
literature, we assume the noise level σnoise to be known and
adjust λ accordingly. From our experiments, we found that λ =
σnoise
16 is a good choice. We terminate our algorithm after 6−30
iterations depending on the noise level, i.e. the higher the noise
level is the more iterations are required. To find an optimal
analysis operator, we learned several operators with varying
values for μ, κ, and k and fixed all other parameters according
to Subsection V-A. Then, we evaluated their performance for Fig. 5. Learned atoms of the analysis operator  . Each of the 128 atoms
the baseline task, which consists of denoising the five test is represented as an 8 × 8 square, where black corresponds to the smallest
negative entry, gray is a zero entry, and white corresponds to the largest
images, each corrupted with the five noise levels as given positive entry.
in Table I. As the final performance measure we use the
average PSNR of the 25 achieved results. The training set
consisted of M = 200 000 image patches, each normalized C. Comparison With Related Approaches
to unit Euclidean norm, that have randomly been extracted
from the five training images shown in Fig. 3. Certainly, The purpose of this subsection is to rank our approach
these images are not considered within any of the performance among other analysis operator learning methods, and to com-
evaluations. Each time, we initialized GOAL with a random pare its performance with state-of-the-art denoising algo-
matrix having normalized rows. Tests with other initializations rithms. Concretely, we compare the denoising performance
like an overcomplete DCT did not influence the final operator. using  learned by GOAL with total-variation (TV) [36]
The results showed, that our approach clearly benefits from which is the currently best known analysis operator, with the
over-completeness. The larger we choose k, the better the oper- recently proposed method AOL [25], and with the currently
ator performs with saturation starting at k = 2n. Therefore, best performing analysis operator FoE [27]. Note that we used
we fixed the number of atoms for all further experiments to the same training set and dimensions for learning the operator
k = 2n. Regarding κ and μ, note that by Lemma 3 and 4 by AOL as for GOAL. For FoE we used the same setup as
these parameters influence the condition number and the originally suggested by the authors. Concerning the required
mutual coherence of the learned operator. Towards answer- computation time for learning an analysis operator, for this
ing the question of what is a good condition number and setting GOAL needs about 10-minutes on an Intel Core i7 3.2
mutual coherence for an analysis operator, Fig. 4(a) shows GHz quad-core with 8 GB RAM. In contrast, AOL is approxi-
the relative denoising performance of 400 operators learned mately ten times slower, and FoE is the computationally most
by GOAL in relation to their mutual coherence and condition. expensive method requiring several hours. All three methods
We would like to mention that according to our experiments, are implemented in unoptimized Matlab code.
this relation is mostly independent from the degree of over- The achieved results for the five test images and the five
completeness. It is also interesting to notice that the best noise levels are given in Table I. Our approach achieves the
learned analysis operator is not a uniformly tight frame. The best results among the analysis methods both regarding PSNR,
concrete values, which led to the best performing analysis and MSSIM. For a visual assessment, Fig. 6 exemplarily shows
operator  ∈ R128×64 in our experiments are κ = 9000 some denoising results achieved by the four analysis operators.
and μ = 0.01. Its singular values are shown in Fig. 4(b) To judge the analysis methods’ denoising performance
and its atoms are visualized in Fig. 5. This operator  globally, we additionally give the results achieved by current
remains unaltered throughout all following experiments in state-of-the-art methods BM3D [37] and K-SVD Denoising
Subsections V-C–V-E. [12], which are specifically designed for the purpose of image
HAWE et al.: ANALYSIS OPERATOR LEARNING AND ITS APPLICATION TO IMAGE RECONSTRUCTION 2147

(a) (b) (a) (b)

(c) (d) (c) (d)

Fig. 6. Images exemplarily showing the typical artifacts created by the four Fig. 7. Results for reconstructing the “Lena” image from 10% of all pixels
compared analysis operators for image denoising (“Man” image, σnoise = 20). using  learned by GOAL, FoE, and [38]. (a) Masked 90% missing pixels.
For a better visualization, a close up is provided for each image. (a) GOAL, (b) Inpainted image GOAL, PSNR 28.57dB and MSSIM 0.840. (c) Inpainted
PSNR 30.44 dB and MSSIM 0.831. (b) AOL [25], PSNR 27.33 dB and image FoE, PSNR 28.06dB and MSSIM 0.822. (d) Inpainted image [38],
MSSIM 0.720. (c) TV [36], PSNR 29.63 dB and MSSIM 0.795. (d) FoE PSNR 27.63dB and MSSIM 0.804.
[27], PSNR 29.75 dB and MSSIM 0.801.

denoising. In most of the cases our method performs slightly In Fig. 7, we show exemplary results for reconstructing the
better than the K-SVD approach, especially for higher noise “lena” image from 10% of all pixels using GOAL, FoE, and
levels, and besides of the “barabara” image it is at most the recently proposed synthesis based method [38]. Table II
≈ 0.5dB worse than BM3D. This effect is due to the very gives a comparison of further images and further number of
special structure of the “barbara” image that rarely occurs in missing pixels. It can be seen that our methods performs best
natural images, which are smoothed by the learned operator. independent of the configuration.

D. Image Inpainting E. Single Image Super-Resolution

In image inpainting as originally proposed in [2], the goal In single image super-resolution (SR), the goal is to recon-
is to fill up a set of damaged or disturbing pixels such that struct a high resolution image s ∈ R N from an observed low
the resulting image is visually appealing. This is necessary resolution image y ∈ Rm . In that, y is assumed to be a blurred
for the restoration of damaged photographs, for removing and downsampled version of s. Mathematically, this process
disturbances caused by e.g. defective hardware, or for deleting can be formulated as y = DBs + e where D ∈ Rm×N is a
unwanted objects. Typically, the positions of the pixels to decimation operator and B ∈ R N×N is a blur operator. Hence,
be filled up are given a priori. In our formulation, when the measurement matrix is given by A = DB. In the ideal case,
N − m pixels must be inpainted, this leads to a binary m × N the exact blur kernel is known or an estimate is given. Here,
dimensional measurements matrix A, where each row contains we consider the more realistic case of an unknown blur kernel.
exactly one entry equal to one. Its position corresponds to a Therefore, to apply our approach for magnifying an image by a
pixel with known intensity. Hence, A reflects the available factor of d in both vertical and horizontal dimension, we model
image information. Regarding λ, it can be used in a way the blur via a Gaussian kernel of dimension (2d −1)×(2d −1)
that our method simultaneously inpaints missing pixels and and with standard deviation σblur = d3 .
denoises the remaining ones. For our experiments, we artificially created a low resolution
As an example for image inpainting, we disturbed some image by downsampling a ground-truth image by a factor
ground-truth images artificially by removing N − m pixels of d using bicubic interpolation. Then, we employed bicubic
randomly distributed over the entire image as exemplary interpolation, FoE, the method from [39], and GOAL to
shown in Fig. 7(a). In that way, the reconstruction quality can magnify this low resolution image by the same factor d. This
be judged both visually and quantitatively. We assumed the upsampled version is then compared with the original image in
data to be free of noise, and empirically selected λ = 10−2 . terms of PSNR and MSSIM. In Table III, we present the results
2148 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO. 6, JUNE 2013

TABLE I
A CHIEVED PSNR IN D ECIBELS (dB) AND MSSIM FOR D ENOISING F IVE T EST I MAGES C ORRUPTED BY F IVE N OISE L EVELS . E ACH C ELL C ONTAINS THE
A CHIEVED R ESULTS FOR THE R ESPECTIVE I MAGE W ITH S IX D IFFERENT A LGORITHMS , W HICH A RE : T OP L EFT G OAL , T OP R IGHT AOL [26],
M IDDLE L EFT TV [36], M IDDLE R IGHT FoE [27], B OTTOM L EFT K-SVD D ENOISING [12], AND B OTTOM R IGHT BM3D [37]

Lena Barbara Man Boat Couple

σnoise / PSNR PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM
5 / 34.15 38.65 36.51 0.945 0.924 37.96 35.95 0.962 0.944 37.77 35.91 0.954 0.932 37.09 35.77 0.938 0.926 37.43 35.55 0.951 0.932
37.65 38.19 0.936 0.938 35.56 37.25 0.948 0.958 36.79 37.45 0.944 0.949 36.17 36.33 0.925 0.917 36.26 37.06 0.940 0.944
38.48 38.45 0.944 0.942 38.12 38.27 0.964 0.964 37.51 37.79 0.952 0.954 37.14 37.25 0.939 0.938 37.24 37.14 0.950 0.951
10 / 28.13 35.58 32.20 0.910 0.856 33.98 31.27 0.930 0.883 33.88 31.33 0.907 0.851 33.72 31.24 0.883 0.842 33.75 30.87 0.903 0.844
34.24 35.12 0.890 0.901 30.84 32.91 0.886 0.923 32.90 33.44 0.884 0.893 32.54 33.23 0.863 0.868 32.32 33.37 0.878 0.889
35.52 35.79 0.910 0.915 34.56 34.96 0.936 0.942 33.64 33.97 0.901 0.907 33.68 33.91 0.883 0.887 33.62 33.86 0.901 0.909
20 / 22.11 32.63 28.50 0.869 0.772 30.17 27.26 0.880 0.791 30.44 27.33 0.831 0.720 30.62 27.16 0.819 0.711 30.39 26.95 0.833 0.727
31.09 31.97 0.827 0.856 26.79 28.39 0.773 0.849 29.63 29.75 0.795 0.801 29.30 29.96 0.778 0.793 28.87 29.77 0.783 0.807
32.39 32.98 0.861 0.875 30.87 31.78 0.881 0.905 30.17 30.59 0.814 0.833 30.44 30.89 0.805 0.825 30.08 30.68 0.817 0.847
25 / 20.17 31.65 27.47 0.854 0.742 29.05 26.08 0.856 0.750 29.43 26.28 0.801 0.677 29.61 26.08 0.792 0.671 29.32 25.81 0.802 0.679
30.05 30.87 0.796 0.836 25.73 27.05 0.724 0.813 28.66 28.62 0.759 0.761 28.32 28.87 0.744 0.758 27.87 28.57 0.746 0.767
31.33 32.02 0.842 0.859 29.59 30.72 0.850 0.887 29.14 29.62 0.780 0.804 29.36 29.92 0.772 0.801 28.92 29.65 0.780 0.820
30 / 18.59 30.86 26.50 0.839 0.717 27.93 24.95 0.818 0.706 28.64 25.30 0.774 0.638 28.80 25.07 0.769 0.630 28.46 24.79 0.780 0.633
29.40 30.00 0.786 0.823 24.91 25.97 0.690 0.787 27.95 27.85 0.736 0.740 27.56 28.01 0.720 0.737 27.09 27.70 0.715 0.743
30.44 31.22 0.823 0.843 28.56 29.82 0.821 0.868 28.30 28.87 0.750 0.780 28.48 29.13 0.744 0.779 27.95 28.81 0.746 0.795

TABLE II
R ESULTS A CHIEVED FOR I NPAINTING F IVE T EST I MAGES W ITH VARYING N UMBER OF M ISSING P IXELS U SING T HREE D IFFERENT M ETHODS . I N
E ACH C ELL , THE PSNR IN dB AND THE MSSIM A RE G IVEN FOR G OAL (T OP ), FoE [27](M IDDLE ), AND M ETHOD [38] (B OTTOM )

% of Missing Pixels Lena Barbara Boat Man House

PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM
90% 28.57 0.840 22.61 0.696 25.61 0.743 26.35 0.755 28.35 0.828
28.06 0.822 22.45 0.682 25.14 0.719 26.23 0.747 28.18 0.828
27.63 0.804 22.49 0.658 24.80 0.683 25.56 0.715 26.62 0.784
80% 31.82 0.895 24.90 0.814 28.55 0.833 28.93 0.847 32.00 0.887
31.09 0.880 23.48 0.762 27.76 0.804 28.51 0.836 31.36 0.880
30.95 0.878 24.72 0.780 27.80 0.804 28.24 0.821 30.20 0.874
50% 37.75 0.956 34.51 0.965 34.47 0.936 34.12 0.947 38.89 0.961
36.70 0.947 28.64 0.919 33.17 0.907 33.49 0.940 37.72 0.957
36.75 0.943 33.21 0.953 33.77 0.918 33.27 0.934 38.11 0.960
20% 43.53 0.985 42.12 0.991 41.04 0.982 40.15 0.985 45.43 0.990
42.29 0.981 36.03 0.981 38.45 0.963 39.15 0.982 44.21 0.989
40.77 0.965 40.63 0.983 39.45 0.966 39.06 0.977 42.95 0.978

TABLE III
R ESULTS IN T ERMS OF PSNR AND MSSIM FOR U PSAMPLING THE S EVEN T EST I MAGES BY A FACTOR OF d = 3 U SING F IVE D IFFERENT A LGORITHMS
G OAL , F O E [27], M ETHOD [39], AND B ICUBIC I NTERPOLATION

Method Face August Barbara Lena Man Boat Couple

PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM PSNR MSSIM
GOAL 32.37 0.801 23.28 0.791 24.42 0.731 32.36 0.889 29.48 0.837 28.25 0.800 27.79 0.786
FoE 32.19 0.797 22.95 0.782 24.30 0.727 31.82 0.885 29.17 0.832 28.00 0.797 27.64 0.782
Method [39] 32.16 0.795 22.90 0.771 24.25 0.719 32.00 0.881 29.29 0.829 28.04 0.793 27.56 0.778
Bicubic 31.57 0.771 22.07 0.724 24.13 0.703 30.81 0.863 28.39 0.796 27.18 0.759 26.92 0.743

for upsampling the respective images by d = 3. The presented upsampling and denoising. In Fig. 8 we present the result
results show that our method outperforms the current state-of- for simultaneously denoising and upsampling a low resolu-
the-art. We want to emphasize that the blur kernel used for tion version of the image “august” by a factor of d = 3,
downsampling is different from the blur kernel used in our which has been corrupted by AWGN with σnoise = 8. As
upsampling procedure. it can be seen, our method produces the best results both
Note that many single image super-resolution algorithms visually and quantitatively, especially regarding the MSSIM.
rely on clean noise free input data, whereas the general Due to high texture this image is hard to upscale even when
analysis approach as formulated in Equation (49) naturally no noise is present, see the second column of Table III.
handles noisy data, and is able to perform simultaneous Results obtained for other images confirm this good perfor-
HAWE et al.: ANALYSIS OPERATOR LEARNING AND ITS APPLICATION TO IMAGE RECONSTRUCTION 2149

[3] W. T. Freeman, T. R. Jones, and E. C. Pasztor, “Example-based super-

resolution,” IEEE Comput. Graph. Appl., vol. 22, no. 2, pp. 56–65, Apr.
2002.
[4] E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles:
Exact signal reconstruction from highly incomplete frequency informa-
tion,” IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, Feb. 2006.
[5] A. Kirsch, An Introduction to the Mathematical Theory of Inverse
Problems. New York, USA: Springer-Verlag, 1991.
[6] M. Elad, M. A. T. Figueiredo, and Y. Ma, “On the role of sparse and
(a) (b) redundant representations in image processing,” Proc. IEEE, vol. 98,
no. 6, pp. 972–982, Jun. 2010.
[7] M. Elad, P. Milanfar, and R. Rubinstein, “Analysis versus synthesis in
signal priors,” Inverse Problems, vol. 23, no. 3, pp. 947–968, Jun. 2007.
[8] J. A. Tropp and S. J. Wright, “Computational methods for sparse solution
of linear inverse problems,” Proc. IEEE, vol. 98, no. 6, pp. 948–958,
Jun. 2010.
[9] S. Mallat, “A theory for multiresolution signal decomposition: The
wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 11,
no. 7, pp. 674–693, Jul. 1989.
(c) (d)
[10] E. Le Pennec and S. Mallat, “Sparse geometric image representations
with bandelets,” IEEE Trans. Imag. Process., vol. 14, no. 4, pp. 423–438,
Apr. 2005.
[11] J.-L. Starck, E. J. Candès, and D. L. Donoho, “The curvelet transform
for image denoising,” IEEE Trans. Imag. Process., vol. 11, no. 6, pp.
670–684, Jun. 2002.
[12] M. Elad and M. Aharon, “Image denoising via sparse and redundant
representations over learned dictionaries,” IEEE Trans. Imag. Process.,
vol. 15, no. 12, pp. 3736–3745, Dec. 2006.
(e) (f) [13] J. Mairal, F. Bach, J. Ponce, and G. Sapiro, “Online learning for matrix
factorization and sparse coding,” J. Mach. Learn. Res., vol. 11, no. 1,
Fig. 8. Single image super-resolution results of four algorithms on noisy pp. 19–60, Mar. 2010.
data for magnifying a low resolution image by a factor of three together [14] M. Zhou, H. Chen, J. Paisley, L. Ren, L. Li, Z. Xing, D. Dunson,
with the corresponding PSNR and MSSIM. The low resolution image has G. Sapiro, and L. Carin, “Nonparametric Bayesian dictionary learning
been corrupted by AWGN with σnoise = 8. (a) Original image “august.” for analysis of noisy and incomplete images,” IEEE Trans. Imag.
(b) Noisy low resolution image. (c) Bicubic Interpolation, PSNR 21.63 dB Process., vol. 21, no. 1, pp. 130–144, Jan. 2012.
and MSSIM 0.653. (d) Method [39], PSNR 22.07 dB and MSSIM 0.663. [15] K. Kreutz-Delgado, J. F. Murray, B. D. Rao, K. Engan, T. W. Lee, and
(e) FoE [27], PSNR 22.17 dB and MSSIM 0.711. (f) GOAL, PSNR 22.45 dB T. J. Sejnowski, “Dictionary learning algorithms for sparse representa-
and MSSIM 0.726. tion,” Neural Comput., vol. 15, no. 2, pp. 349–396, Feb. 2003.
[16] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for
mance of GOAL but are not presented here due to space designing overcomplete dictionaries for sparse representation,” IEEE
limitation. Trans. Signal Process., vol. 54, no. 11, pp. 4311–4322, Nov. 2006.
VI. C ONCLUSION [17] R. Rubinstein, M. Zibulevsky, and M. Elad, “Double sparsity: Learning
sparse dictionaries for sparse signal approximation,” IEEE Trans. Signal
This paper deals with the topic of learning an analysis Process., vol. 58, no. 3, pp. 1553–1564, Mar. 2010.
operator from example image patches, and how to apply it for [18] I. Tošić and P. Frossard, “Dictionary learning,” IEEE Signal Process.
Mag., vol. 28, no. 2, pp. 27–38, Mar. 2011.
solving inverse problems in imaging. To learn the operator, we [19] S. Nam, M. Davies, M. Elad, and R. Gribonval, “Cosparse analysis
motivate an  p -minimization on the set of full-rank matrices modeling-uniqueness and algorithms,” in Proc. IEEE Int. Conf. Acoust.
with normalized columns. A geometric conjugate gradient Speech Signal Process., May 2011, pp. 5804–5807.
[20] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, and K. Knight, “Sparsity
method on the oblique manifold is suggested to solve the and smoothness via the fused lasso,” J. Royal Stat. Soc. Ser. B, vol. 67,
arising optimization task. Furthermore, we give a partition- no. 1, pp. 91–108, Feb. 2005.
ing invariant method for employing the local patch based [21] I. W. Selesnick and M. A. T. Figueiredo, “Signal restoration with
overcomplete wavelet transforms: Comparison of analysis and synthesis
analysis operator such that globally consistent reconstruction priors,” in Proc. Int. Soc. Opt. Photon. Wavelets XIII, Sep. 2009, p.
results are achieved. For the famous tasks of image denois- 74460D.
ing, image inpainting, and single image super-resolution, we [22] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based
provide promising results that are competitive with and even noise removal algorithms,” Physica D, vol. 60, nos. 1–4, pp. 259–268,
Nov. 1992.
outperform current state-of-the-art techniques. Similar as for [23] B. Ophir, M. Elad, N. Bertin, and M. D. Plumbley, “Sequential minimal
the synthesis signal reconstruction model with dictionaries, eigenvalues–An approach to analysis dictionary learning,” in Proc. Eur.
we expect that depending on the application at hand, the Signal Conf., 2011, pp. 1465–1469.
[24] R. Rubinstein, T. Faktor, and M. Elad, “K-SVD dictionary-learning for
performance of the analysis approach can be further increased the analysis sparse model,” in Proc. IEEE Int. Conf. Acoust. Speech,
by learning the particular operator with regard to the specific Signal, 2012, pp. 5405–5408.
problem, or employing a specialized training set. [25] M. Yaghoobi, S. Nam, R. Gribonval, and M. E. Davies, “Analysis
operator learning for overcomplete cosparse representations,” in Proc.
Eur. Signal Conf., 2011, pp. 1470–1474.
R EFERENCES [26] M. Yaghoobi, S. Nam, R. Gribonval, and M. E. Davies, “Noise aware
[1] J. Portilla, V. Strela, M. Wainwright, and E. Simoncelli, “Image denois- analysis operator learning for approximately cosparse signals,” in Proc.
ing using scale mixtures of Gaussians in the wavelet domain,” IEEE IEEE Int. Conf. Acoust. Speech, Signal, Mar. 2012, pp. 5409–5412.
Trans. Imag. Process., vol. 12, no. 11, pp. 1338–1351, Nov. 2003. [27] S. Roth and M. Black, “Fields of experts,” Int. J. Comput. Vis., vol. 82,
[2] M. Bertalmìo, G. Sapiro, V. Caselles, and C. Ballester, “Image inpaint- no. 2, pp. 205–229, Apr. 2009.
ing,” in Proc. Assoc. Comput. Mach. Spec. Interest Group Comput. [28] N. T. Trendafilov, “A continuous-time approach to the oblique Procrustes
Graph., 2000, pp. 417–424. problem,” Behaviormetrika, vol. 26, no. 2, pp. 167–181, 1999.
2150 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO. 6, JUNE 2013

[29] M. Kleinsteuber and H. Shen, “Blind source separation with compres- Martin Kleinsteuber received the Ph.D. degree
sively sensed linear mixtures,” IEEE Signal Process. Lett., vol. 19, no. 2, in mathematics from the University of Würzburg,
pp. 107–110, Feb. 2012. Würzburg, Germany, in 2006.
[30] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Since 2009, he has been an Assistant Profes-
Matrix Manifolds, Princeton, NJ, USA: Princeton University Press, Sep. sor of geometric optimization and machine learn-
2008. ing with the Department of Electrical Engineering
[31] Y. Dai and Y. Yuan, “An efficient hybrid conjugate gradient method and Information Technology, Technische Univer-
for unconstrained optimization,” Ann. Oper. Res., vol. 103, nos. 1–4, sität München, Munich, Germany. He was a Post-
pp. 33–47, Mar. 2001. Doctoral Fellow with National ICT Australia Ltd.,
[32] J. C. Gilbert and J. Nocedal, “Global convergence properties of con- the Australian National University, Canberra, Aus-
jugate gradient methods for optimization,” Soc. Indian Autom. Manuf. tralia, and the University of Würzburg. His current
J. Optim., vol. 2, no. 1, pp. 21–42, 1992. research interests include latent variable analysis and feature extraction, partic-
[33] W. Ring and B. Wirth, “Optimization methods on Riemannian manifolds ularly robust dimensionality reduction techniques and sparse representations,
and their application to shape space,” Soc. Indian Autom. Manuf. together with their applications to audio- and image processing.
J. Optimi., vol. 22, no. 2, pp. 596–627, 2012. Dr. Kleinsteuber was a recipient of the SIAM Student Paper Prize in 2004
[34] S. Hawe, M. Kleinsteuber, and K. Diepold, “Cartoon-like image recon- and the Robert-Sauer-Award of the Bavarian Academy of Science in 2008 for
struction via constrained  p -minimization,” in Proc. IEEE Int. Conf. his contribution to research on Jacobi-type methods on Lie algebras.
Acoust. Speech, Signal, Mar. 2012, pp. 717–720.
[35] Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality
assessment: From error visibility to structural similarity,” IEEE Trans.
Image Process., vol. 13, no. 4, pp. 600–612, Apr. 2004.
[36] J. Dahl, P. C. Hansen, S. Jensen, and T. L. Jensen, “Algorithms and
software for total variation image reconstruction via first-order methods,”
Numer. Algorithms, vol. 53, no. 1, pp. 67–92, Jan. 2010.
[37] 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.
[38] M. Zhou, H. Chen, J. Paisley, R. Lu, L. Li, Z. Xing, D. Dunson, Klaus Diepold was born, raised, and educated
G. Sapiro, and L. Carin, “Nonparametric Bayesian dictionary learning in Munich, Germany. He received the Dipl.-Ing.
for analysis of noisy and incomplete images,” IEEE Trans. Image and the Dr.-Ing. degrees in electrical engineering
Process., vol. 21, no. 1, pp. 130–144, Jan. 2012. from the Technische Universität München (TUM),
[39] J. Yang, J. Wright, T. Huang, and Y. Ma, “Image super-resolution via Munich, Germany, in 1987 and 1992, respectively.
sparse representation,” IEEE Trans. Image Process., vol. 19, no. 11, From 1993 to 2002, he was with the Research
pp. 2861–2873, Nov. 2010. Laboratory, DynaPel Systems, Inc., Munich, Ger-
many where he was lead research on video signal
processing and the television and video compression
industry. In 2002, he joined the Department of
Electrical Engineering and Information Technology,
Simon Hawe received the [Link]. degree in elec- TUM, as a Full Professor. His current research interests include computational
trical engineering from the Technische Universität aspects of motion picture technology, video processing and compression, com-
München, Munich, Germany, in 2008, where he puter graphics, computer vision, 3-D audio, machine learning, and cognition
is currently pursuing the Ph.D. degree in learning for technical systems.
sparse data models via geometric optimization with Dr. Diepold is a member of COST IC1003 Qualitnet, a European initiative
the Electrical Engineering Department. dedicated to quality of experience. He was the Principal Investigator in a
His current research interests include sparse meth- cooperative research center involved in research on high-fidelity telepresence
ods in machine learning, geometric optimization, and teleaction. He is a member on the Executive Board of the excellence
inverse problems, computer vision, and image cluster, Cognition for Technical Systems and the Scientific Director of the
processing. Centre for Digital Technology and Management.