# An Efficient Algorithm for Overcomplete Sparsifying Transform Learning with Signal Denoising

Authors: Beiping Hou; Zhihui Zhu; Gang Li; Aihua Yu
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

This paper deals with the problem of overcomplete transform learning. An alternating minimization based procedure is proposed for solving the formulated sparsifying transform learning problem. A closed-form solution is derived for the minimization involved in transform update stage. Compared with existing ones, our proposed algorithm significantly reduces the computation complexity. Experiments and simulations are carried out with synthetic data and real images to demonstrate the superiority of the proposed approach in terms of the averaged representation and denoising errors, the percentage of successful and meaningful recovery of the analysis dictionary, and, more significantly, the computation efficiency.


---

## Body

## 1. Introduction

Signal representations play a role of cornerstones in signal processing. They have experienced the evolution from the Fourier transforms, which are traced back to the later part of the 19th century, to the wavelet expansions in the 1980s and then the sparse representations which have been one of the hotspots in signal processing community for the last two decades. Instead of the traditional expansions in terms of bases and frames, sparse redundant representations look for the best approximation of a signal vector with a linear combination of few atoms from an overcomplete set of well-designed vectors [1]. This topic is closely related to thesparsifying dictionary learning [2, 3] and the compressed sensing (CS) [4–6], an emerging area of research in which the sparsity of the signals to be recovered is a prerequisite.Letv
∈
R
N
×
1. The l
p-norm of v is defined as (1)
v
p
≜
∑
k
=
1
N
v
k
p
1
/
p
,
p
≥
1
.For convenience, v
0 is used to denote the number of nonzero elements in v though l
0 is not a norm as it does not satisfy the homogeneity.


### 1.1. Signal Models and Dictionary Learning

Lety
=
x
+
ξ
∈
R
N
×
1 be the signal under consideration, where x is the clean signal and ξ is an additive noise or modelling error. There are basically two models used for the clean signal x. The first one is calledsynthesis model which assumes that the clean signal x is given by (2)
x
=
∑
k
=
1
J
s
k
ψ
k
≜
Ψ
s
,where Ψ
∈
R
N
×
J is usually called asynthesis dictionary (throughout this paper, a lower-case letter, say q, is usually used to denote a column vector of the matrix Q. Furthermore, q
k and q
i
j denote the kth column vector and (
i
,
j
)th entry of matrix Q, resp.) and is saidovercomplete if N
<
J, and s
∈
R
J
×
1 is the coefficient vector of x in Ψ. The signal s is said κ-sparse in Ψ if s
0
≤
κ.The procedure for findings with x and Ψ given is referred to asanalysis. Clearly, (3)
s
=
Ω
xis a solution to this problem as long as Ψ
Ω
=
I
N, where I
N is the identity matrix of dimension N. It should be pointed out that (2) and (3) are equivalent when Ψ is square and nonsingular. When Ψ is overcomplete, such a claim does not hold as (2) is not an one-to-one mapping between x and s. Equation (3) is referred to asanalysis model [7–12] and Ω
∈
R
J
×
N is calledanalysis dictionary oranalysis operator. When J
>
N, which is assumed throughout this paper, the dictionary is saidredundant orovercomplete.Let{
y
n
} be a set of signal samples belonging to a certain class of signals. For a given sparsity level, say κ, sparsifying synthesis dictionary learning is to find Ψ such that y
n
=
Ψ
s
n
+
ξ
n, ∀
n with s
n
0
≤
κ. Methods for learning the dictionary Ψ from signal samples {
y
k
} have been one of the hottest topics in signal processing, among which the method of optimal directions (MOD) [13] and K-singular value decomposition (KSVD) based algorithms [14, 15] have been considered as state of the art in this area. The basic idea behind these methods is based on the alternating minimization strategy [16–18], in which the dictionary Ψ and the sparse coefficients s are updated alternatively using the orthogonal matching pursuit (OMP) techniques [19–21].Compared with the synthesis model-based sparsifying dictionary learning, its analysis model-based counterpart has received less attention. The analysis dictionary learning problem can be stated as(4)
Ω
,
X
≜
arg
⁡
m
i
n
Ω
~
,
X
~
Y
-
X
~
F
2
,
Ω
~
x
~
m
0
≤
J
-
l
≜
κ
,
∀
m
,where l, denoting the number of zero elements in s, is the cosparsity given and x
~
m denotes the mth column of X
~. (throughout this paper, we use “tilde” (say Q
~) to denote a (matrix/vector) variable, while Q is used for the optimal/true one in a problem). As X is the clean signal, (4) is also signal denoising.The first attempt to address such a problem was given in [7], where it is suggested to incrementally learn the dictionary Ω one row at a time. The performance of the proposed algorithm, however, depends heavily on a randomized initialization strategy. A different approach to the dictionary learning was proposed in [8–10]. As proposed in [10], (4) is investigated with the following formulation: (5)
m
i
n
Ω
~
,
X
~
Ω
~
X
~
1
s.t.
Ω
~
∈
C
,
Y
-
X
~
F
≤
ϵ
0
,where C is a specified set within which the optimal operator is searched and ϵ
0
≥
0 is the parameter corresponding to the noise level. Such a problem is alternatively addressed using a regularized version in [10].Another signal model, which is very much close to the analysis model, is thetransform model raised in [22]. Such a model is characterized with (6)
W
y
=
s
+
ϵ
,where W
∈
R
J
×
N is called atransform and s is sparse with s
0
≤
κ, while ϵ is the representation/modeling error in transform domain and is assumed to bevery small. It is interesting to note that, for an analysis signal y
=
x
+
ξ with Ω
x
=
s sparse, one has Ω
y
=
s
+
Ω
ξ which is always in the form specified in (6), while if y belongs to the transform model-based class in Ω, which means there exists a pair (
s
,
e
) with s being sparse such that Ω
y
=
s
+
e, Ω
x
ξ
=
s
e may have no solution for x, ξ at all given that the number of equations is bigger than that of unknowns due to N
<
J assumed.WhenW is available, the signal y is transformed to κ-sparse vector s: (7)
s
≜
arg
⁡
m
i
n
s
~
s
~
-
W
y
2
2
,
s.t.
s
~
0
≤
κ
.The solution can be achieved simply by thresholding the product W
y and retaining only the κ largest coefficients (magnitudewise) of this product. The estimate y
^ of the original signal y can be found from s and W by minimizing s
-
W
y
~
2
2 with respect to y
~. This yields (8)
y
^
=
W
T
W
-
1
W
T
sas long as W is full rank (under the assumption that J
≥
N), where T denotes the transpose operator.The sparsifying transform learning problem is formulated as(9)
W
,
S
≜
arg
⁡
m
i
n
W
~
,
S
~
S
~
-
W
~
Y
F
2
s
~
m
0
≤
κ
,
∀
m
.This problem is not well defined and has the problem of so-calledscale ambiguity. It is observed that a trivial solution for the above problem is S
=
0 and W
=
0. One approach to prevent such solutions is to add some penalty functions into the objective. In [22], the penalty ν
W
~
F
2
-
log
⁡
det
⁡
(
W
~
) was included in the cost function for learning a complete transform with dimension N
×
N. The term log
⁡
det
⁡
(
W
~
) is used to force the full rank of W and the term W
~
F
2 is applied to avoid the scaling issue. See (10)
L
R
B
W
~
,
S
~
≜
S
~
-
W
~
Y
F
2
+
λ
ν
W
~
F
2
-
log
⁡
det
⁡
W
~
,where λ and ν are two positive constants, while log
⁡
det
⁡
(
·
) is the natural logarithmic operator applied to the determinant of a matrix. The optimal transform design was then reformulated as (11)
W
,
S
≜
arg
⁡
m
i
n
W
~
,
S
~
L
R
B
W
~
,
S
~
,
s
~
m
0
≤
κ
,
∀
msubject to det
⁡
(
W
)
>
0. Such a problem was solved in [22, 23].


### 1.2. Contributions

In this paper, we will consider solving (9) with given (
Y
,
κ
). An efficient algorithm is proposed to solve the optimal sparsifying transform design problem. It should be pointed out that our proposed approach is more general than those in [22, 23] as it deals with overcomplete transform learning and is more efficient than [24] since the transform is updated using a closed-form solution rather than the one obtained with a gradient-based numerical procedure. All of these will be elaborated further in next section. Roughly speaking, (4) has been attacked using alternative optimization, in which with X
~ obtained, Ω is updated based on Ω
~
x
~
m
0
≤
J
-
l
≜
κ, ∀
m. This can be considered with (12)
Ω
≜
arg
⁡
m
i
n
Ω
~
,
S
~
S
~
-
Ω
~
X
~
F
2
,
s
~
m
0
≤
κ
,
∀
m
.As seen, this problem is exactly the same as (9) with Y replaced with X
~. So, the efficiency of the analysis KSVD proposed in [12] is expected to be improved using our proposed algorithm for analysis dictionary update. We will demonstrate the potential of this proposed algorithm in a series of experiments using both synthetic and real data and compare it with the approach in [12]. As to be seen, our approach yields a performance comparable to that by the latter in terms of denoising and recovery of the reference dictionary and more importantly a much improved computation efficiency.The outline of this paper is given as follows. In Section2, after briefing some related work, we formulate the problem to be investigated in the context of transform learning and how it is related to analysis dictionary learning. Section 3 is devoted to solving the optimal transform design. An efficient algorithm is derived for that purpose, in which a closed-form solution is obtained for updating the transform. The performance of this algorithm and its application to signal denoising are discussed in this section. In order to demonstrate the superiority of the proposed algorithm, experiments are carried out in Section 4, which show promising performance of the proposed approach. Some concluding remarks are given in Section 5.


## 1.1. Signal Models and Dictionary Learning

Lety
=
x
+
ξ
∈
R
N
×
1 be the signal under consideration, where x is the clean signal and ξ is an additive noise or modelling error. There are basically two models used for the clean signal x. The first one is calledsynthesis model which assumes that the clean signal x is given by (2)
x
=
∑
k
=
1
J
s
k
ψ
k
≜
Ψ
s
,where Ψ
∈
R
N
×
J is usually called asynthesis dictionary (throughout this paper, a lower-case letter, say q, is usually used to denote a column vector of the matrix Q. Furthermore, q
k and q
i
j denote the kth column vector and (
i
,
j
)th entry of matrix Q, resp.) and is saidovercomplete if N
<
J, and s
∈
R
J
×
1 is the coefficient vector of x in Ψ. The signal s is said κ-sparse in Ψ if s
0
≤
κ.The procedure for findings with x and Ψ given is referred to asanalysis. Clearly, (3)
s
=
Ω
xis a solution to this problem as long as Ψ
Ω
=
I
N, where I
N is the identity matrix of dimension N. It should be pointed out that (2) and (3) are equivalent when Ψ is square and nonsingular. When Ψ is overcomplete, such a claim does not hold as (2) is not an one-to-one mapping between x and s. Equation (3) is referred to asanalysis model [7–12] and Ω
∈
R
J
×
N is calledanalysis dictionary oranalysis operator. When J
>
N, which is assumed throughout this paper, the dictionary is saidredundant orovercomplete.Let{
y
n
} be a set of signal samples belonging to a certain class of signals. For a given sparsity level, say κ, sparsifying synthesis dictionary learning is to find Ψ such that y
n
=
Ψ
s
n
+
ξ
n, ∀
n with s
n
0
≤
κ. Methods for learning the dictionary Ψ from signal samples {
y
k
} have been one of the hottest topics in signal processing, among which the method of optimal directions (MOD) [13] and K-singular value decomposition (KSVD) based algorithms [14, 15] have been considered as state of the art in this area. The basic idea behind these methods is based on the alternating minimization strategy [16–18], in which the dictionary Ψ and the sparse coefficients s are updated alternatively using the orthogonal matching pursuit (OMP) techniques [19–21].Compared with the synthesis model-based sparsifying dictionary learning, its analysis model-based counterpart has received less attention. The analysis dictionary learning problem can be stated as(4)
Ω
,
X
≜
arg
⁡
m
i
n
Ω
~
,
X
~
Y
-
X
~
F
2
,
Ω
~
x
~
m
0
≤
J
-
l
≜
κ
,
∀
m
,where l, denoting the number of zero elements in s, is the cosparsity given and x
~
m denotes the mth column of X
~. (throughout this paper, we use “tilde” (say Q
~) to denote a (matrix/vector) variable, while Q is used for the optimal/true one in a problem). As X is the clean signal, (4) is also signal denoising.The first attempt to address such a problem was given in [7], where it is suggested to incrementally learn the dictionary Ω one row at a time. The performance of the proposed algorithm, however, depends heavily on a randomized initialization strategy. A different approach to the dictionary learning was proposed in [8–10]. As proposed in [10], (4) is investigated with the following formulation: (5)
m
i
n
Ω
~
,
X
~
Ω
~
X
~
1
s.t.
Ω
~
∈
C
,
Y
-
X
~
F
≤
ϵ
0
,where C is a specified set within which the optimal operator is searched and ϵ
0
≥
0 is the parameter corresponding to the noise level. Such a problem is alternatively addressed using a regularized version in [10].Another signal model, which is very much close to the analysis model, is thetransform model raised in [22]. Such a model is characterized with (6)
W
y
=
s
+
ϵ
,where W
∈
R
J
×
N is called atransform and s is sparse with s
0
≤
κ, while ϵ is the representation/modeling error in transform domain and is assumed to bevery small. It is interesting to note that, for an analysis signal y
=
x
+
ξ with Ω
x
=
s sparse, one has Ω
y
=
s
+
Ω
ξ which is always in the form specified in (6), while if y belongs to the transform model-based class in Ω, which means there exists a pair (
s
,
e
) with s being sparse such that Ω
y
=
s
+
e, Ω
x
ξ
=
s
e may have no solution for x, ξ at all given that the number of equations is bigger than that of unknowns due to N
<
J assumed.WhenW is available, the signal y is transformed to κ-sparse vector s: (7)
s
≜
arg
⁡
m
i
n
s
~
s
~
-
W
y
2
2
,
s.t.
s
~
0
≤
κ
.The solution can be achieved simply by thresholding the product W
y and retaining only the κ largest coefficients (magnitudewise) of this product. The estimate y
^ of the original signal y can be found from s and W by minimizing s
-
W
y
~
2
2 with respect to y
~. This yields (8)
y
^
=
W
T
W
-
1
W
T
sas long as W is full rank (under the assumption that J
≥
N), where T denotes the transpose operator.The sparsifying transform learning problem is formulated as(9)
W
,
S
≜
arg
⁡
m
i
n
W
~
,
S
~
S
~
-
W
~
Y
F
2
s
~
m
0
≤
κ
,
∀
m
.This problem is not well defined and has the problem of so-calledscale ambiguity. It is observed that a trivial solution for the above problem is S
=
0 and W
=
0. One approach to prevent such solutions is to add some penalty functions into the objective. In [22], the penalty ν
W
~
F
2
-
log
⁡
det
⁡
(
W
~
) was included in the cost function for learning a complete transform with dimension N
×
N. The term log
⁡
det
⁡
(
W
~
) is used to force the full rank of W and the term W
~
F
2 is applied to avoid the scaling issue. See (10)
L
R
B
W
~
,
S
~
≜
S
~
-
W
~
Y
F
2
+
λ
ν
W
~
F
2
-
log
⁡
det
⁡
W
~
,where λ and ν are two positive constants, while log
⁡
det
⁡
(
·
) is the natural logarithmic operator applied to the determinant of a matrix. The optimal transform design was then reformulated as (11)
W
,
S
≜
arg
⁡
m
i
n
W
~
,
S
~
L
R
B
W
~
,
S
~
,
s
~
m
0
≤
κ
,
∀
msubject to det
⁡
(
W
)
>
0. Such a problem was solved in [22, 23].


## 1.2. Contributions

In this paper, we will consider solving (9) with given (
Y
,
κ
). An efficient algorithm is proposed to solve the optimal sparsifying transform design problem. It should be pointed out that our proposed approach is more general than those in [22, 23] as it deals with overcomplete transform learning and is more efficient than [24] since the transform is updated using a closed-form solution rather than the one obtained with a gradient-based numerical procedure. All of these will be elaborated further in next section. Roughly speaking, (4) has been attacked using alternative optimization, in which with X
~ obtained, Ω is updated based on Ω
~
x
~
m
0
≤
J
-
l
≜
κ, ∀
m. This can be considered with (12)
Ω
≜
arg
⁡
m
i
n
Ω
~
,
S
~
S
~
-
Ω
~
X
~
F
2
,
s
~
m
0
≤
κ
,
∀
m
.As seen, this problem is exactly the same as (9) with Y replaced with X
~. So, the efficiency of the analysis KSVD proposed in [12] is expected to be improved using our proposed algorithm for analysis dictionary update. We will demonstrate the potential of this proposed algorithm in a series of experiments using both synthetic and real data and compare it with the approach in [12]. As to be seen, our approach yields a performance comparable to that by the latter in terms of denoising and recovery of the reference dictionary and more importantly a much improved computation efficiency.The outline of this paper is given as follows. In Section2, after briefing some related work, we formulate the problem to be investigated in the context of transform learning and how it is related to analysis dictionary learning. Section 3 is devoted to solving the optimal transform design. An efficient algorithm is derived for that purpose, in which a closed-form solution is obtained for updating the transform. The performance of this algorithm and its application to signal denoising are discussed in this section. In order to demonstrate the superiority of the proposed algorithm, experiments are carried out in Section 4, which show promising performance of the proposed approach. Some concluding remarks are given in Section 5.


## 2. Preliminaries and Problem Formulation

### 2.1. Related Works

The first attempt to address the analysis dictionary learning was given in [7], where it is suggested to incrementally learn the dictionary Ω one row at a time. The performance of the proposed algorithm, however, depends heavily on a randomized initialization strategy. A different approach to the dictionary learning was proposed in [8–10], where the goal of sparsifying the representations Ω
X
~ is formulated using an l
1-based penalty function rather than the l
0-based one. To avoid the trivial solutions, the dictionary is constrained to be uniform normalized tight frame; that is, Ω
T
Ω
=
I
N, the identity matrix of dimension N. Such a constraint clearly limits the searching space for the optimal dictionary.Recently, it was argued in [12] that, in the analysis model, one should emphasize the zeros of Ω
x. Based on this argument, (4) was reformulated in [12] as (13)
Ω
,
X
,
Λ
i
i
=
1
L
≜
arg
⁡
m
i
n
Ω
~
,
X
~
,
Λ
~
i
i
=
1
L
Y
-
X
~
F
2
,
Ω
~
Λ
~
i
x
~
i
=
0
,
R
a
n
k
Ω
~
Λ
~
i
=
N
-
r
,
i
=
1
,
…
,
L
,
ω
~
m
†
2
=
1
,
∀
m
,where q
m
† denotes the mth row vector of matrix Q, Λ
~
i is thecosupport of x
~
i, defined as the set of indices for the l rows of Ω
~ that are orthogonal to x
~
i, and Ω
~
Λ
~
i is the submatrix of Ω
~ that contains only the rows indexed by Λ
~
i, while r is the dimension of the subspace to which x
~
i belongs.An algorithm, called analysis KSVD, was proposed in [12] to address (13). The analysis KSVD, denoted as AKSVD, was developed in a similar way as the one adopted by its counterpart KSVD [14] for the synthesis model and as shown in [12], outperforming the one reported in [7]. It has been observed, however, that the AKSVD is very time consuming and therefore more efficient algorithms are needed.In order to restrict the solution set of the transform learning problem (9) to an admissible set, W
~
T
W
~
=
I
J was applied in [25], which means forcing the row of the transform W
~ to be a tight frame. However, only tight frame constraint is insufficient to prevent undesired solutions of (9) for learning an overcomplete transform W
~ (J
>
N), for which it is highly possible that practical algorithms result in a local minimal solution of the form W
=
W
1
T
0
T with W
1 being a tight frame. In [24], the overcomplete transform learning (9) was attacked with the following penalty function to be embedded to the cost function: (14)
ν
W
~
F
2
-
log
⁡
W
~
T
W
~
+
η
∑
n
≠
n
′
w
~
n
†
w
~
n
′
†
T
q
,where ν and η are weighting factors and q is usually set to 1 and 2. An iterative procedure using the alternating minimization strategy was proposed to solve the corresponding problem, where a gradient-based algorithm is utilized for updating the transform. Clearly, such a method is time consuming and may suffer from local minimum issue due to the use of gradient-based algorithm. Furthermore, the possibility of having same rows in W
~ can be reduced with the second term of the penalty but zero rows may occur if W
~ is not (row) normalized.


### 2.2. Problem Formulation

Based on the previous discussions, the optimal sparsifying transform learning problem can be formulated as(15)
W
,
S
≜
arg
⁡
m
i
n
W
~
,
S
~
S
~
-
W
~
Y
F
2
,
s
~
m
0
≤
κ
,
∀
m
,
W
~
T
W
~
≻
0subject to (16)
w
~
k
†
0
=
1
,
∀
k
.To implement (15), we consider the following cost function: (17)
L
W
~
,
S
~
≜
S
~
-
W
~
Y
F
2
+
λ
F
W
~
,where F
(
W
~
) is a penalty term (i.e., a regularizer) defined as (18)
F
W
~
≜
ν
W
~
F
2
-
log
⁡
det
⁡
W
~
T
W
~
.Both λ and ν are positive parameters to be discussed later.The problem to be addressed in this paper is finally formulated as(19)
W
,
S
=
arg
⁡
m
i
n
W
~
,
S
~
L
W
~
,
S
~
,
s
~
m
0
≤
κ
,
∀
m
,
w
~
k
†
2
=
1
,
∀
k
.Comment 1.
(i) The additional regularizerF
(
W
~
) in (17) serves as a constraint on the space for searching the true transform W. The term W
~
F
2 is to ensure that W is not too big, while −log
⁡
det
⁡
(
W
~
T
W
~
) is to guarantee that W is not too small and hence the constraint of full rank is met. The combination of the two terms ensures that the optimal dictionary is searched within a well behaved space such that the undesired scaling problem can be avoided.
(ii) As to be seen in the next subsection, the minimum of this regularizer is a set of tight frames. Therefore, if (4) has a solution W close to a tight frame, it can be well estimated by (19).
(iii) Comparing (19) with (11), one realizes that (19) can be used for optimal overcomplete sparsifying transform learning. Besides, instead of −log
⁡
det
⁡
(
W
~
), −log
⁡
det
⁡
(
W
~
T
W
~
) is adopted in our approach. By doing so, the constraint det
⁡
(
W
~
)
>
0 is no longer needed.


### 2.3. Analysis Dictionaries and Transforms

Consider the analysis model-based dictionary learning (4). Define E
≜
Y
-
X
~ and E
≜
Ω
~
Y
-
S
~ as therepresentation error and thesparsification error, respectively, where Ω
~
X
~
=
S
~. Then (20)
Ω
~
Y
-
S
~
=
Ω
~
Y
-
X
~
.It is assumed in the sequel that the true dictionary Ω is full rank. Such an assumption ensures that the row vectors of Ω can span the entire signal space R
N
×
1.The following lemma specifies a class of true dictionariesΩ for which the equivalence between (15) and (4) holds.Lemma 1.
Equation (15) is equivalent to (4) if the true dictionary Ω is a c-tight frame (TF) (c
>
0); that is, Ω has an SVD of form (21)
Ω
=
c
U
I
N
0
T
V
T
,where both U and V are orthonormal matrices of proper dimension.Proof.
With (20) and (21), we have (22)
S
~
-
Ω
Y
F
2
=
Ω
Y
-
X
~
F
2
=
c
2
U
I
N
0
T
V
T
Y
-
X
~
F
2
=
c
2
I
N
0
T
V
T
Y
-
X
~
F
2
=
c
2
V
T
Y
-
X
~
F
2
=
c
2
Y
-
X
~
F
2
.Then, Lemma 1 follows immediately from the fact that the two constraints s
~
m
0
≤
κ and Ω
x
~
m
0
≤
κ are equivalent as S
~
=
Ω
X
~.Comment 2.
As seen from Lemma1, the equivalence between the two problems holds if the true dictionary Ω is a c-tight frame. This is in general not true. However, it is observed that for given Y the problem defined by (4) has in general more than one solution (
Ω
,
X
). In fact, define D
α
≜
diag
⁡
(
α
1
,
…
,
α
m
,
…
,
α
J
) and D
β
≜
diag
⁡
(
β
1
,
…
,
β
m
,
…
,
β
L
). It can be seen that if (
Ω
0
,
X
0
) is a solution, is (
Ω
,
X
): (23)
Ω
=
D
α
Ω
0
T
,
X
=
T
-
1
X
0
D
β
,where D
α is any given nonsingular diagonal matrix, while both T
∈
R
N
×
N and D
β are arbitrary except det
⁡
(
D
β
)
≠
0 such that Y
-
X
F
=
Y
-
X
0
F. This means that there exist a lot of degrees of freedom in the solution set of (4) and that by forcing finding the solution for which the dictionary Ω is the closest to a TF, (15) is then expected to yield a transform W that is a good estimate of a solution for the analysis dictionary learning problem (4).


### 2.4. Rule of Thumb for Choosingλ and ν

Now, we present a result regarding minimization ofF
(
W
~
). A similar result was given in [22].Lemma 2.
LetF
(
W
~
) be defined in (18). Then, (24)
F
W
~
≥
N
1
+
l
o
g
ν
≜
η
0with the equality achieved if and only if W
~ is a ν
-
1-tight frame.Proof.
LetW
~
=
U
Σ
N
0
T
V
T be the SVD of Ω
~, where Σ
N
=
diag
⁡
(
σ
1
,
…
,
σ
k
,
…
,
σ
N
) with σ
k
≥
σ
k
+
1
>
0, ∀
k. Clearly, F
(
Ω
~
)
=
∑
k
=
1
N
[
ν
σ
k
2
-
log
⁡
(
σ
k
2
)
]
≜
∑
k
=
1
N
f
(
σ
k
).
Note thatd
f
(
σ
k
)
/
d
σ
k
=
ν
2
σ
k
-
2
σ
k
/
σ
k
2 and σ
k
>
0. It follows from d
f
(
σ
k
)
/
d
σ
k
=
0 that the solution for minimizing f
(
σ
k
) is given by (25)
σ
k
∗
=
ν
-
1
⟹
f
σ
k
∗
=
1
+
log
⁡
ν
≤
f
σ
k
,
∀
kand hence (24) holds. Clearly, the equality is attainedif and only if
σ
k
∗
=
ν
-
1, ∀
k.Comment 3.
(i) Based on Lemma2, the cost function L
(
W
~
,
S
~
)
=
S
~
-
W
~
Y
F
2
+
λ
F
(
W
~
) is bounded from below with λ
η
0.
(ii) It follows from Lemmas1 and 2 that if the reference analysis operator happens to be a tight frame, then the solution of (19) should yield the same operator no matter what λ is taken.
(iii) As indicated in Comment2, the solution W can be expected to be close to a tight frame. This suggests that λ should take a value much larger than one since when λ is very large, minimizing L
(
W
~
,
S
~
) can be approximated with minimizing F
(
W
~
).
(iv) As to be seen from (19), we need to normalize the dictionary W
~ such that w
~
j
†
2
=
1, ∀
j. Under this constraint, the constant c in (21) turns out to be J
/
N. Since J
≥
N holds for overcomplete analysis dictionary, ν
<
1 should be ensured.
Before turning to next section, we should point out that though the proposed approach is not directly related to the overcomplete analysis dictionary learning (4) as (13) is in general, it can provide a good estimate of the true dictionary in a very efficient way, which is believed to be useful in improving the analysis KSVD algorithm.


## 2.1. Related Works

The first attempt to address the analysis dictionary learning was given in [7], where it is suggested to incrementally learn the dictionary Ω one row at a time. The performance of the proposed algorithm, however, depends heavily on a randomized initialization strategy. A different approach to the dictionary learning was proposed in [8–10], where the goal of sparsifying the representations Ω
X
~ is formulated using an l
1-based penalty function rather than the l
0-based one. To avoid the trivial solutions, the dictionary is constrained to be uniform normalized tight frame; that is, Ω
T
Ω
=
I
N, the identity matrix of dimension N. Such a constraint clearly limits the searching space for the optimal dictionary.Recently, it was argued in [12] that, in the analysis model, one should emphasize the zeros of Ω
x. Based on this argument, (4) was reformulated in [12] as (13)
Ω
,
X
,
Λ
i
i
=
1
L
≜
arg
⁡
m
i
n
Ω
~
,
X
~
,
Λ
~
i
i
=
1
L
Y
-
X
~
F
2
,
Ω
~
Λ
~
i
x
~
i
=
0
,
R
a
n
k
Ω
~
Λ
~
i
=
N
-
r
,
i
=
1
,
…
,
L
,
ω
~
m
†
2
=
1
,
∀
m
,where q
m
† denotes the mth row vector of matrix Q, Λ
~
i is thecosupport of x
~
i, defined as the set of indices for the l rows of Ω
~ that are orthogonal to x
~
i, and Ω
~
Λ
~
i is the submatrix of Ω
~ that contains only the rows indexed by Λ
~
i, while r is the dimension of the subspace to which x
~
i belongs.An algorithm, called analysis KSVD, was proposed in [12] to address (13). The analysis KSVD, denoted as AKSVD, was developed in a similar way as the one adopted by its counterpart KSVD [14] for the synthesis model and as shown in [12], outperforming the one reported in [7]. It has been observed, however, that the AKSVD is very time consuming and therefore more efficient algorithms are needed.In order to restrict the solution set of the transform learning problem (9) to an admissible set, W
~
T
W
~
=
I
J was applied in [25], which means forcing the row of the transform W
~ to be a tight frame. However, only tight frame constraint is insufficient to prevent undesired solutions of (9) for learning an overcomplete transform W
~ (J
>
N), for which it is highly possible that practical algorithms result in a local minimal solution of the form W
=
W
1
T
0
T with W
1 being a tight frame. In [24], the overcomplete transform learning (9) was attacked with the following penalty function to be embedded to the cost function: (14)
ν
W
~
F
2
-
log
⁡
W
~
T
W
~
+
η
∑
n
≠
n
′
w
~
n
†
w
~
n
′
†
T
q
,where ν and η are weighting factors and q is usually set to 1 and 2. An iterative procedure using the alternating minimization strategy was proposed to solve the corresponding problem, where a gradient-based algorithm is utilized for updating the transform. Clearly, such a method is time consuming and may suffer from local minimum issue due to the use of gradient-based algorithm. Furthermore, the possibility of having same rows in W
~ can be reduced with the second term of the penalty but zero rows may occur if W
~ is not (row) normalized.


## 2.2. Problem Formulation

Based on the previous discussions, the optimal sparsifying transform learning problem can be formulated as(15)
W
,
S
≜
arg
⁡
m
i
n
W
~
,
S
~
S
~
-
W
~
Y
F
2
,
s
~
m
0
≤
κ
,
∀
m
,
W
~
T
W
~
≻
0subject to (16)
w
~
k
†
0
=
1
,
∀
k
.To implement (15), we consider the following cost function: (17)
L
W
~
,
S
~
≜
S
~
-
W
~
Y
F
2
+
λ
F
W
~
,where F
(
W
~
) is a penalty term (i.e., a regularizer) defined as (18)
F
W
~
≜
ν
W
~
F
2
-
log
⁡
det
⁡
W
~
T
W
~
.Both λ and ν are positive parameters to be discussed later.The problem to be addressed in this paper is finally formulated as(19)
W
,
S
=
arg
⁡
m
i
n
W
~
,
S
~
L
W
~
,
S
~
,
s
~
m
0
≤
κ
,
∀
m
,
w
~
k
†
2
=
1
,
∀
k
.Comment 1.
(i) The additional regularizerF
(
W
~
) in (17) serves as a constraint on the space for searching the true transform W. The term W
~
F
2 is to ensure that W is not too big, while −log
⁡
det
⁡
(
W
~
T
W
~
) is to guarantee that W is not too small and hence the constraint of full rank is met. The combination of the two terms ensures that the optimal dictionary is searched within a well behaved space such that the undesired scaling problem can be avoided.
(ii) As to be seen in the next subsection, the minimum of this regularizer is a set of tight frames. Therefore, if (4) has a solution W close to a tight frame, it can be well estimated by (19).
(iii) Comparing (19) with (11), one realizes that (19) can be used for optimal overcomplete sparsifying transform learning. Besides, instead of −log
⁡
det
⁡
(
W
~
), −log
⁡
det
⁡
(
W
~
T
W
~
) is adopted in our approach. By doing so, the constraint det
⁡
(
W
~
)
>
0 is no longer needed.


## 2.3. Analysis Dictionaries and Transforms

Consider the analysis model-based dictionary learning (4). Define E
≜
Y
-
X
~ and E
≜
Ω
~
Y
-
S
~ as therepresentation error and thesparsification error, respectively, where Ω
~
X
~
=
S
~. Then (20)
Ω
~
Y
-
S
~
=
Ω
~
Y
-
X
~
.It is assumed in the sequel that the true dictionary Ω is full rank. Such an assumption ensures that the row vectors of Ω can span the entire signal space R
N
×
1.The following lemma specifies a class of true dictionariesΩ for which the equivalence between (15) and (4) holds.Lemma 1.
Equation (15) is equivalent to (4) if the true dictionary Ω is a c-tight frame (TF) (c
>
0); that is, Ω has an SVD of form (21)
Ω
=
c
U
I
N
0
T
V
T
,where both U and V are orthonormal matrices of proper dimension.Proof.
With (20) and (21), we have (22)
S
~
-
Ω
Y
F
2
=
Ω
Y
-
X
~
F
2
=
c
2
U
I
N
0
T
V
T
Y
-
X
~
F
2
=
c
2
I
N
0
T
V
T
Y
-
X
~
F
2
=
c
2
V
T
Y
-
X
~
F
2
=
c
2
Y
-
X
~
F
2
.Then, Lemma 1 follows immediately from the fact that the two constraints s
~
m
0
≤
κ and Ω
x
~
m
0
≤
κ are equivalent as S
~
=
Ω
X
~.Comment 2.
As seen from Lemma1, the equivalence between the two problems holds if the true dictionary Ω is a c-tight frame. This is in general not true. However, it is observed that for given Y the problem defined by (4) has in general more than one solution (
Ω
,
X
). In fact, define D
α
≜
diag
⁡
(
α
1
,
…
,
α
m
,
…
,
α
J
) and D
β
≜
diag
⁡
(
β
1
,
…
,
β
m
,
…
,
β
L
). It can be seen that if (
Ω
0
,
X
0
) is a solution, is (
Ω
,
X
): (23)
Ω
=
D
α
Ω
0
T
,
X
=
T
-
1
X
0
D
β
,where D
α is any given nonsingular diagonal matrix, while both T
∈
R
N
×
N and D
β are arbitrary except det
⁡
(
D
β
)
≠
0 such that Y
-
X
F
=
Y
-
X
0
F. This means that there exist a lot of degrees of freedom in the solution set of (4) and that by forcing finding the solution for which the dictionary Ω is the closest to a TF, (15) is then expected to yield a transform W that is a good estimate of a solution for the analysis dictionary learning problem (4).


## 2.4. Rule of Thumb for Choosingλ and ν

Now, we present a result regarding minimization ofF
(
W
~
). A similar result was given in [22].Lemma 2.
LetF
(
W
~
) be defined in (18). Then, (24)
F
W
~
≥
N
1
+
l
o
g
ν
≜
η
0with the equality achieved if and only if W
~ is a ν
-
1-tight frame.Proof.
LetW
~
=
U
Σ
N
0
T
V
T be the SVD of Ω
~, where Σ
N
=
diag
⁡
(
σ
1
,
…
,
σ
k
,
…
,
σ
N
) with σ
k
≥
σ
k
+
1
>
0, ∀
k. Clearly, F
(
Ω
~
)
=
∑
k
=
1
N
[
ν
σ
k
2
-
log
⁡
(
σ
k
2
)
]
≜
∑
k
=
1
N
f
(
σ
k
).
Note thatd
f
(
σ
k
)
/
d
σ
k
=
ν
2
σ
k
-
2
σ
k
/
σ
k
2 and σ
k
>
0. It follows from d
f
(
σ
k
)
/
d
σ
k
=
0 that the solution for minimizing f
(
σ
k
) is given by (25)
σ
k
∗
=
ν
-
1
⟹
f
σ
k
∗
=
1
+
log
⁡
ν
≤
f
σ
k
,
∀
kand hence (24) holds. Clearly, the equality is attainedif and only if
σ
k
∗
=
ν
-
1, ∀
k.Comment 3.
(i) Based on Lemma2, the cost function L
(
W
~
,
S
~
)
=
S
~
-
W
~
Y
F
2
+
λ
F
(
W
~
) is bounded from below with λ
η
0.
(ii) It follows from Lemmas1 and 2 that if the reference analysis operator happens to be a tight frame, then the solution of (19) should yield the same operator no matter what λ is taken.
(iii) As indicated in Comment2, the solution W can be expected to be close to a tight frame. This suggests that λ should take a value much larger than one since when λ is very large, minimizing L
(
W
~
,
S
~
) can be approximated with minimizing F
(
W
~
).
(iv) As to be seen from (19), we need to normalize the dictionary W
~ such that w
~
j
†
2
=
1, ∀
j. Under this constraint, the constant c in (21) turns out to be J
/
N. Since J
≥
N holds for overcomplete analysis dictionary, ν
<
1 should be ensured.
Before turning to next section, we should point out that though the proposed approach is not directly related to the overcomplete analysis dictionary learning (4) as (13) is in general, it can provide a good estimate of the true dictionary in a very efficient way, which is believed to be useful in improving the analysis KSVD algorithm.


## 3. The Proposed Algorithm

This section is devoted to attacking the optimal sparsifying transform learning problem by deriving an algorithm to solve (19).


### 3.1. The Algorithm

Note that the problem defined in (19) is nonconvex. A practical approach to such problems is to adopt the alternating minimization strategy [16–18]. The basic idea is to update the dictionary W
~ and the (column) sparse matrix S
~ alternatively. The proposed algorithm is outlined below.A
l
g
P
r
o
p
o
s
e
dInitialization. Set the training data matrix Y, the number of iterations N
i
t
e, and the parameters λ
,
ν; initiate W with W
0. (there are several ways to initialize W. A simple way is to set W
=
r
a
n
d
n
(
J
,
N
). Other way is that for the mth row (m
=
1,2
,
…
,
J) we first randomly choose N
-
1 columns from the training data matrix Y to form a submatrix Y
-, then set w
m
†
=
e
-
e
Y
-
Y
-
‡, where e
=
r
a
n
d
n
(
1
,
N
) and Y
-
‡ denotes the pseudoinverse of Y
-. At last we normalize each row to unit l
2 norm. In the simulations to be presented, we choose the second approach to initialize W in order to have a fair comparison with [12]).Begin (
k
=
1,2
,
…
,
N
i
t
e
)Step 1.
UpdateS such that L
(
W
k
-
1
,
S
~
) is minimized under the constraint s
~
m
0
≤
κ, ∀
m. This is equivalent to (26)
S
k
≜
arg
⁡
m
i
n
S
~
S
~
-
W
k
-
1
Y
F
2
,
s
~
m
0
≤
κ
,
∀
m
.Step 2.
WithS
k just obtained, update W by solving (27)
W
~
k
≜
arg
⁡
m
i
n
W
~
L
W
~
,
S
k
.Step 3.
NormalizeW
~
k such that all the rows of W
k
=
T
s
c
W
~
k are of unit length. This can be achieved by setting the diagonal T
s
c with its (
j
,
j
)th entry equal to the inverse of the length for the jth row of W
~
k for j
=
1,2
,
…
,
J.End. Output W
^
=
W
N
i
t
e, the estimate of the true dictionary W.Now, let us elaborate the three steps involved in the proposed algorithm further.The first one involves updatingS with (26), whose solution S
k, as mentioned before, can be achieved directly with the product V
≜
W
k
-
1
Y in the following way.Begin (
m
=
1,2
,
…
,
L
)
(i)
Determineϱ
m which is equal to the absolute value of the κth largest coefficients (magnitudewise) in the mth column of V.
(ii)
Forj
=
1,2
,
…
,
J, set the (
j
,
m
)th element of S
k to (28)
v
j
m
,
v
j
m
≥
ϱ
m
,
0
,
v
j
m
<
ϱ
m
,
wherev
j
m denotes the (
j
,
m
)th element of matrix V.End. Output S
k. (In case where there are more than κ nonzero elements in S
k, just keep any set of the first κ largest ones and set the others to zero.)Now, let us consider (27) involved in Step 2. Define (29)
M
≜
S
~
Y
T
,
Z
≜
λ
ν
I
N
+
Y
Y
T
.It can be shown that the cost function L
(
W
~
,
S
~
) can be rewritten as (30)
L
W
~
,
S
~
=
S
~
F
2
-
2
t
r
M
W
~
T
+
t
r
W
~
Z
W
~
T
-
λ
log
⁡
det
⁡
W
~
T
W
~
≜
L
1
W
~
,where t
r
[
·
] denotes the trace operator of a matrix and L
1
(
W
~
) is used to emphasize the fact that S
~ is fixed in the minimization to be considered.Therefore, (27) is just a special case of the following minimization problem: (31)
W
∗
≜
arg
⁡
m
i
n
W
~
L
1
W
~
.Theorem 3.
LetM
∈
R
J
×
N, Z
∈
R
N
×
N, and S
~ three matrices independent of matrix W
~
∈
R
J
×
N with J
≥
N, and Z
=
U
z
Λ
2
U
z
T an SVD of the (symmetric) positive-definite Z and Z
s
q
r
t
≜
U
z
Λ
U
z
T. Furthermore, let M
~
≜
M
Z
s
q
r
t
-
1
=
U
0
[
Π
2
0
]
V
0
T be an SVD of M
~, where Π
=
diag
⁡
(
π
1
,
…
,
π
k
,
…
,
π
N
). A solution to (31) is then given by (32)
W
∗
=
U
0
Σ
2
0
V
0
T
Z
s
q
r
t
-
1
,where Σ
2
=
diag
⁡
(
σ
1
2
,
…
,
σ
k
2
,
…
,
σ
N
2
) with (33)
σ
k
2
=
π
k
2
+
π
k
4
+
4
λ
2
,
∀
k
.The detailed proof can be found in Appendix.As seen from (29), no matter what Y is taken the corresponding Z matrix in our problem is always (symmetric) positive-definite as long as λ
>
0 and ν
>
0. Therefore, (27) in Step 2 of the proposed algorithm can be solved using Theorem 3.Step3 is to normalize the transform such that W
k has all its rows of unit length. This is for avoiding zero rows in the transform W and thescaling ambiguity [22]. Such scaling has been popularly used in sparsifying dictionary design.Comment 4.
As noted, the proposed algorithm shares the same framework as the one used in [22]. It should be pointed out that besides the fact that our algorithm is more general in the sense that it can deal with analysis dictionaries (transforms) of any dimension rather those square nonsingular ones, the main difference between the two lies in the way of how to solve (27) in Step 2. In [22], a gradient-based algorithm was used to attack (27). Such an algorithm suffers in general from slow convergence and more seriously convergence to a local minimum of the problem. Therefore, our proposed algorithm is more efficient and faster convergence is then expected.


### 3.2. Convergence

Consider the situation when the row normalization (i.e., Step3) is not executed in proposed A
l
g
P
r
o
p
o
s
e
d. As the solutions for the minimization involved in both Steps 1 and 2 are exact, the cost function L
(
W
~
,
S
~
) is ensured to be decreasing. Given that the cost function is bounded from below, the proposed A
l
g
P
r
o
p
o
s
e
d without Step 3 can ensure that the cost function L
(
W
~
,
S
~
) converges. As mentioned before, the obtained transform W may have zero rows during the iterative procedure. In order to avoid this phenomenon, Step 3 has to be applied when the transform is overcomplete and, in that case, the cost function is found empirically to converge as well.As argued in [22], the convergence of cost function does not indicate that of the iterates. There have been no theoretical results achieved regarding convergence of the iterates for the proposed algorithm. However, extensive experiments showed that the proposed algorithm does demonstrate good behavior of convergence.


### 3.3. Computation Efficiency

A thorough comparison between the algorithm proposed in [22], denoted as A
l
g
R
B, and some relevant algorithms in terms of computational cost can be found in [22]. Since our proposed algorithm shares the same framework with A
l
g
R
B and is derived as an alternative of A
l
g
A
K
S
V
D, the one recently reported in [12], we mainly focus on comparing our proposed A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D. For the sake of completeness, we present the computation issues of our proposed A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D as follows.Our algorithm is an iterative method. In each iteration, there are three steps. They are sparse coding, dictionary updating, and the dictionary normalizing. In the sparse coding step, we simply keep theκ largest entries of each column in the matrix W
Y by full sorting. This, as indicated in [26], requires O
(
J
L
log
⁡
N
) operations. Knowing that computing Ω
Y requires O
(
J
N
L
) operations, one concludes that the sparse coding step needs roughly O
(
J
L
(
N
+
log
⁡
N
)
) flops. In the dictionary update step, what we mainly do is just computing the SVD form of a matrix M. Thus, this step costs approximately O
(
N
J
2
+
N
3
) flops, while in the normalization step, it needs O
(
J
N
) flops. Thus the total cost per iteration of the proposed algorithm is the sum of the three, which is dominated by O
(
J
N
L
) under the assumption of N
≤
J
≪
L.As toA
l
g
A
K
S
V
D, there are two steps involved in each iteration, sparse coding and dictionary updating, and the computation cost is dominated by the sparse coding step. As given in [12], the sparse coding step requires O
(
L
N
3
J
) operations if the optimized backward greedy (OBG) algorithm is utilized. Therefore, the overall computation complexity scales as O
(
L
N
3
J
).Based on the computations above, one can see that our proposed algorithm is much more efficient thanA
l
g
A
K
S
V
D. This is confirmed by experiments to be given in the next section.


### 3.4. Signal Denoising

For a measurementy
i and a (row normalized) dictionary Ω
^ given, the clean signal x
i is estimated in A
l
g
A
K
S
V
D via solving the problem (34)
x
^
i
,
Λ
^
i
≜
arg
⁡
m
i
n
x
,
Λ
x
-
y
i
2
,
Ω
^
Λ
x
=
0
,
R
a
n
k
Ω
^
Λ
=
N
-
r
,where Λ
⊂
{
1,2
,
…
,
L
} is the set of row indices and r is the dimension of the subspace to which the signal belongs. The OBG is used to attack such a problem and it has to be run L times for each iteration in order to obtain the clear signal matrix X and the final estimate of the true dictionary Ω. That is why A
l
g
A
K
S
V
D is very time consuming.InA
l
g
P
r
o
p
o
s
e
d, instead of X the (column) sparse matrix S is used. As mentioned before, the latter can be updated efficiently using (28). Taking the transform obtained using A
l
g
P
r
o
p
o
s
e
d as the analysis dictionary, we can run the OBG to get the corresponding estimate x
^
i of x
i for i
=
1,2
,
…
,
L.Alternatively, following the same approach as in [22] we can consider the following: (35)
W
,
X
,
S
=
arg
⁡
m
i
n
W
~
,
X
~
,
S
~
J
W
~
,
X
~
,
S
~
,
s
~
m
0
≤
κ
,
∀
m
,where the objective function J
(
W
~
,
X
~
,
S
~
) is defined as (36)
J
W
~
,
X
~
,
S
~
≜
W
~
X
~
-
S
~
F
2
+
λ
F
W
~
+
μ
Y
-
X
~
F
2with F
(
W
~
) the same as defined before and μ
>
0 constant determining the importance of the representation error.This is simultaneous optimization of the dictionaryW, the clean signals X, and the sparse coefficients S. It is easy to understand that this problem can be addressed using alternating minimization and a similar algorithm to A
l
g
P
r
o
p
o
s
e
d can be derived. In fact, such an algorithm can be obtained from A
l
g
P
r
o
p
o
s
e
d with adding (i) X
0
=
Y to theinitialization and (ii) a new step between Steps 1 and 2 for updating the clean signals X: (37)
X
k
≜
arg
⁡
m
i
n
X
~
J
W
k
-
1
,
X
~
,
S
k
.The above is equivalent to minimizing W
k
-
1
X
~
-
S
k
F
2
+
μ
Y
-
X
~
F
2 with respect to X
~, a standard least squares problem.In this paper, we focus on the OBG based denoising approach as comparisons are given between proposedA
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D.


## 3.1. The Algorithm

Note that the problem defined in (19) is nonconvex. A practical approach to such problems is to adopt the alternating minimization strategy [16–18]. The basic idea is to update the dictionary W
~ and the (column) sparse matrix S
~ alternatively. The proposed algorithm is outlined below.A
l
g
P
r
o
p
o
s
e
dInitialization. Set the training data matrix Y, the number of iterations N
i
t
e, and the parameters λ
,
ν; initiate W with W
0. (there are several ways to initialize W. A simple way is to set W
=
r
a
n
d
n
(
J
,
N
). Other way is that for the mth row (m
=
1,2
,
…
,
J) we first randomly choose N
-
1 columns from the training data matrix Y to form a submatrix Y
-, then set w
m
†
=
e
-
e
Y
-
Y
-
‡, where e
=
r
a
n
d
n
(
1
,
N
) and Y
-
‡ denotes the pseudoinverse of Y
-. At last we normalize each row to unit l
2 norm. In the simulations to be presented, we choose the second approach to initialize W in order to have a fair comparison with [12]).Begin (
k
=
1,2
,
…
,
N
i
t
e
)Step 1.
UpdateS such that L
(
W
k
-
1
,
S
~
) is minimized under the constraint s
~
m
0
≤
κ, ∀
m. This is equivalent to (26)
S
k
≜
arg
⁡
m
i
n
S
~
S
~
-
W
k
-
1
Y
F
2
,
s
~
m
0
≤
κ
,
∀
m
.Step 2.
WithS
k just obtained, update W by solving (27)
W
~
k
≜
arg
⁡
m
i
n
W
~
L
W
~
,
S
k
.Step 3.
NormalizeW
~
k such that all the rows of W
k
=
T
s
c
W
~
k are of unit length. This can be achieved by setting the diagonal T
s
c with its (
j
,
j
)th entry equal to the inverse of the length for the jth row of W
~
k for j
=
1,2
,
…
,
J.End. Output W
^
=
W
N
i
t
e, the estimate of the true dictionary W.Now, let us elaborate the three steps involved in the proposed algorithm further.The first one involves updatingS with (26), whose solution S
k, as mentioned before, can be achieved directly with the product V
≜
W
k
-
1
Y in the following way.Begin (
m
=
1,2
,
…
,
L
)
(i)
Determineϱ
m which is equal to the absolute value of the κth largest coefficients (magnitudewise) in the mth column of V.
(ii)
Forj
=
1,2
,
…
,
J, set the (
j
,
m
)th element of S
k to (28)
v
j
m
,
v
j
m
≥
ϱ
m
,
0
,
v
j
m
<
ϱ
m
,
wherev
j
m denotes the (
j
,
m
)th element of matrix V.End. Output S
k. (In case where there are more than κ nonzero elements in S
k, just keep any set of the first κ largest ones and set the others to zero.)Now, let us consider (27) involved in Step 2. Define (29)
M
≜
S
~
Y
T
,
Z
≜
λ
ν
I
N
+
Y
Y
T
.It can be shown that the cost function L
(
W
~
,
S
~
) can be rewritten as (30)
L
W
~
,
S
~
=
S
~
F
2
-
2
t
r
M
W
~
T
+
t
r
W
~
Z
W
~
T
-
λ
log
⁡
det
⁡
W
~
T
W
~
≜
L
1
W
~
,where t
r
[
·
] denotes the trace operator of a matrix and L
1
(
W
~
) is used to emphasize the fact that S
~ is fixed in the minimization to be considered.Therefore, (27) is just a special case of the following minimization problem: (31)
W
∗
≜
arg
⁡
m
i
n
W
~
L
1
W
~
.Theorem 3.
LetM
∈
R
J
×
N, Z
∈
R
N
×
N, and S
~ three matrices independent of matrix W
~
∈
R
J
×
N with J
≥
N, and Z
=
U
z
Λ
2
U
z
T an SVD of the (symmetric) positive-definite Z and Z
s
q
r
t
≜
U
z
Λ
U
z
T. Furthermore, let M
~
≜
M
Z
s
q
r
t
-
1
=
U
0
[
Π
2
0
]
V
0
T be an SVD of M
~, where Π
=
diag
⁡
(
π
1
,
…
,
π
k
,
…
,
π
N
). A solution to (31) is then given by (32)
W
∗
=
U
0
Σ
2
0
V
0
T
Z
s
q
r
t
-
1
,where Σ
2
=
diag
⁡
(
σ
1
2
,
…
,
σ
k
2
,
…
,
σ
N
2
) with (33)
σ
k
2
=
π
k
2
+
π
k
4
+
4
λ
2
,
∀
k
.The detailed proof can be found in Appendix.As seen from (29), no matter what Y is taken the corresponding Z matrix in our problem is always (symmetric) positive-definite as long as λ
>
0 and ν
>
0. Therefore, (27) in Step 2 of the proposed algorithm can be solved using Theorem 3.Step3 is to normalize the transform such that W
k has all its rows of unit length. This is for avoiding zero rows in the transform W and thescaling ambiguity [22]. Such scaling has been popularly used in sparsifying dictionary design.Comment 4.
As noted, the proposed algorithm shares the same framework as the one used in [22]. It should be pointed out that besides the fact that our algorithm is more general in the sense that it can deal with analysis dictionaries (transforms) of any dimension rather those square nonsingular ones, the main difference between the two lies in the way of how to solve (27) in Step 2. In [22], a gradient-based algorithm was used to attack (27). Such an algorithm suffers in general from slow convergence and more seriously convergence to a local minimum of the problem. Therefore, our proposed algorithm is more efficient and faster convergence is then expected.


## 3.2. Convergence

Consider the situation when the row normalization (i.e., Step3) is not executed in proposed A
l
g
P
r
o
p
o
s
e
d. As the solutions for the minimization involved in both Steps 1 and 2 are exact, the cost function L
(
W
~
,
S
~
) is ensured to be decreasing. Given that the cost function is bounded from below, the proposed A
l
g
P
r
o
p
o
s
e
d without Step 3 can ensure that the cost function L
(
W
~
,
S
~
) converges. As mentioned before, the obtained transform W may have zero rows during the iterative procedure. In order to avoid this phenomenon, Step 3 has to be applied when the transform is overcomplete and, in that case, the cost function is found empirically to converge as well.As argued in [22], the convergence of cost function does not indicate that of the iterates. There have been no theoretical results achieved regarding convergence of the iterates for the proposed algorithm. However, extensive experiments showed that the proposed algorithm does demonstrate good behavior of convergence.


## 3.3. Computation Efficiency

A thorough comparison between the algorithm proposed in [22], denoted as A
l
g
R
B, and some relevant algorithms in terms of computational cost can be found in [22]. Since our proposed algorithm shares the same framework with A
l
g
R
B and is derived as an alternative of A
l
g
A
K
S
V
D, the one recently reported in [12], we mainly focus on comparing our proposed A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D. For the sake of completeness, we present the computation issues of our proposed A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D as follows.Our algorithm is an iterative method. In each iteration, there are three steps. They are sparse coding, dictionary updating, and the dictionary normalizing. In the sparse coding step, we simply keep theκ largest entries of each column in the matrix W
Y by full sorting. This, as indicated in [26], requires O
(
J
L
log
⁡
N
) operations. Knowing that computing Ω
Y requires O
(
J
N
L
) operations, one concludes that the sparse coding step needs roughly O
(
J
L
(
N
+
log
⁡
N
)
) flops. In the dictionary update step, what we mainly do is just computing the SVD form of a matrix M. Thus, this step costs approximately O
(
N
J
2
+
N
3
) flops, while in the normalization step, it needs O
(
J
N
) flops. Thus the total cost per iteration of the proposed algorithm is the sum of the three, which is dominated by O
(
J
N
L
) under the assumption of N
≤
J
≪
L.As toA
l
g
A
K
S
V
D, there are two steps involved in each iteration, sparse coding and dictionary updating, and the computation cost is dominated by the sparse coding step. As given in [12], the sparse coding step requires O
(
L
N
3
J
) operations if the optimized backward greedy (OBG) algorithm is utilized. Therefore, the overall computation complexity scales as O
(
L
N
3
J
).Based on the computations above, one can see that our proposed algorithm is much more efficient thanA
l
g
A
K
S
V
D. This is confirmed by experiments to be given in the next section.


## 3.4. Signal Denoising

For a measurementy
i and a (row normalized) dictionary Ω
^ given, the clean signal x
i is estimated in A
l
g
A
K
S
V
D via solving the problem (34)
x
^
i
,
Λ
^
i
≜
arg
⁡
m
i
n
x
,
Λ
x
-
y
i
2
,
Ω
^
Λ
x
=
0
,
R
a
n
k
Ω
^
Λ
=
N
-
r
,where Λ
⊂
{
1,2
,
…
,
L
} is the set of row indices and r is the dimension of the subspace to which the signal belongs. The OBG is used to attack such a problem and it has to be run L times for each iteration in order to obtain the clear signal matrix X and the final estimate of the true dictionary Ω. That is why A
l
g
A
K
S
V
D is very time consuming.InA
l
g
P
r
o
p
o
s
e
d, instead of X the (column) sparse matrix S is used. As mentioned before, the latter can be updated efficiently using (28). Taking the transform obtained using A
l
g
P
r
o
p
o
s
e
d as the analysis dictionary, we can run the OBG to get the corresponding estimate x
^
i of x
i for i
=
1,2
,
…
,
L.Alternatively, following the same approach as in [22] we can consider the following: (35)
W
,
X
,
S
=
arg
⁡
m
i
n
W
~
,
X
~
,
S
~
J
W
~
,
X
~
,
S
~
,
s
~
m
0
≤
κ
,
∀
m
,where the objective function J
(
W
~
,
X
~
,
S
~
) is defined as (36)
J
W
~
,
X
~
,
S
~
≜
W
~
X
~
-
S
~
F
2
+
λ
F
W
~
+
μ
Y
-
X
~
F
2with F
(
W
~
) the same as defined before and μ
>
0 constant determining the importance of the representation error.This is simultaneous optimization of the dictionaryW, the clean signals X, and the sparse coefficients S. It is easy to understand that this problem can be addressed using alternating minimization and a similar algorithm to A
l
g
P
r
o
p
o
s
e
d can be derived. In fact, such an algorithm can be obtained from A
l
g
P
r
o
p
o
s
e
d with adding (i) X
0
=
Y to theinitialization and (ii) a new step between Steps 1 and 2 for updating the clean signals X: (37)
X
k
≜
arg
⁡
m
i
n
X
~
J
W
k
-
1
,
X
~
,
S
k
.The above is equivalent to minimizing W
k
-
1
X
~
-
S
k
F
2
+
μ
Y
-
X
~
F
2 with respect to X
~, a standard least squares problem.In this paper, we focus on the OBG based denoising approach as comparisons are given between proposedA
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D.


## 4. Experimental Results

In this section, we present some experiments to illustrate the performance of our proposed algorithmA
l
g
P
r
o
p
o
s
e
d and compare it with that by A
l
g
A
K
S
V
D given in [12]. For convenience, the dictionaries, obtained using these two algorithms, are denoted with Ω
P
r
o
p
o
s
e
d and Ω
A
K
S
V
D, respectively. All the experiments are performed on a computer with an Intel Core i7 CPU at 3.40 GHz, 16 GB memory, and 64-bit Windows 7 operating system.


### 4.1. Synthetic Experiments

In order to have a fair comparison, the setup used here is similar to [12]. Let J
=
50, N
=
25, and L = 50,000. A J
×
N matrix W is generated with normally distributed entries and then Ω
=
D
α
W, where D
α is such a diagonal matrix that all the rows of Ω have unit l
2-norm. This Ω serves as the true overcomplete analysis dictionary. Then, a set of L vectors {
x
i
} of N
×
1 is produced with each living in a 4-dimensional subspace. Each vector x
i is obtained as follows: we first construct a submatrix Ω
(
i
)
∈
R
21
×
N by randomly choosing 21 rows from Ω and then compute the orthonormal basis {
ψ
k
(
i
)
} for the null space of Ω
(
i
), which can be obtained from the SVD of Ω
(
i
). With Ψ
(
i
)
≜
ψ
1
(
i
)
ψ
2
(
i
)
ψ
3
(
i
)
ψ
4
(
i
), calculate x
i
=
Ψ
(
i
)
s
i, where each element in s
i is randomly generated by Gaussian distribution of independent and identically distributed (i.i.d.) zero-mean and unit variance. At last, each vector x
i is normalized so that x
i
2
=
1 for all i
=
1,2
,
…
,
L. We can see that Ω
x
i has at least 21 zero components for all {
x
i
} (i
=
1,2
,
…
,
L) generated in such a way.LetI
κ
:
R
J
→
R
J denote a hard thresholding operator that keeps the largest (in magnitude) κ elements of a vector. We define the effective sparsity κ
~ of a length-J vector s as (38)
κ
~
=
m
i
n
κ
s.t.
I
κ
s
2
≥
0.99
s
2
.Figure 1 represents the histogram of the effective sparsity level of {
Ω
x
i
}
(
i
=
1,2
,
…
,
L
). As can be seen, the effective sparsity level is varying from 17 to 25and is most concentrated near 21.Figure 1
Histogram of the effective sparsity of{
Ω
x
i
}
(
i
=
1,2
,
…
,
L
) for J
=
50, N
=
25, and L
=
50000.The measurements{
y
i
} are generated with (39)
y
i
=
x
i
+
ϵ
i
,
∀
i
,where each entry of the additive noise vector ϵ
i is a white Gaussian noise with i.i.d. zero-mean and σ
2-variance. Define the signal-to-noise ratio (SNR) as ρ
s
n
r
≜
10
l
o
g
10
(
E
x
i
2
2
/
E
ϵ
i
2
2
), where E
[
·
] denotes the mathematical expectation operator. When σ
=
0.2
/
N
=
0.04, ρ
s
n
r
=
10
l
o
g
10
(
1
/
(
N
∗
σ
2
)
)
=
10
l
o
g
10
(
25
)
≈
14dB.We then estimate the analysis dictionaryΩ using A
l
g
A
K
S
V
D and A
l
g
P
r
o
p
o
s
e
d with X
~
=
Y (see (12)), respectively, for both noise-free case (i.e., σ
=
0) and noisy case with noise level ρ
s
n
r
=
14dB. In both cases, we fix κ
=
21 and N
i
t
e
=
200, while ν
=
0.01 and λ
=
300 for A
l
g
P
r
o
p
o
s
e
d. (As we said before, Ω
x
i has at least 21 zero components for all {
x
i
}. Thus, we can choose κ such that κ
≤
J
-
21
=
29. It is illustrated in Figure 1 that the effective sparsity is concentrated near 21, and therefore we choose κ
=
21.) Both A
l
g
P
r
o
p
o
s
e
d and AlgA
K
S
V
D are initialized with a dictionary Ω
0 in which each row is generated as follows: we first randomly choose N
-
1 columns from the training data matrix Y to form a submatrix Y
-; then set ω
m
†
=
e
-
e
Y
-
Y
-
‡, where e
=
r
a
n
d
n
(
1
,
N
) and Y
-
‡, as defined before, denotes the pseudoinverse of Y
-; at last we normalize each row to unit l
2 norm [12].Convergence Behavior of
A
l
g
P
r
o
p
o
s
e
d. The convergence behavior of our proposed A
l
g
P
r
o
p
o
s
e
d is shown in Figures 2 and 3 in terms of the objective function L
(
Ω
k
,
S
k
) and the relative change of iterates which is defined as Ω
k
-
Ω
k
-
1
F
/
Ω
k
-
1
F.Figure 2
Evolution of the objective functionL
(
Ω
k
,
S
k
) for J
=
50, N
=
25, and L
=
50000.Figure 3
Evolution of relative change of iteratesΩ
k
-
Ω
k
-
1
F
/
Ω
k
-
1
F for J
=
50, N
=
25, and L
=
50000.Comment 5.
(i) The results are self-explanatory. As expected, the objective functionL
(
Ω
k
,
S
k
) decreases smoothly in both noise-free case and noisy case. This is confirmed in Figure 2, from which it is seen that the latter case yields larger L
(
Ω
k
,
S
k
) than the noise-free case because of the larger sparsification error S
k
-
Ω
k
Y
F
2 in the latter case.
(ii) The convergence behavior of the iteratesΩ
k is depicted in Figure 3. One observes that the iterates of the proposed algorithm show a very good trend to convergence in both cases.Figures4 and 5 show the distribution of singular values of the two analysis dictionaries Ω
A
K
S
V
D and Ω
P
r
o
p
o
s
e
d.Figure 4
Distribution of singular values of reference dictionaryΩ, Ω
A
K
S
V
D, and Ω
P
r
o
p
o
s
e
d for the noise-free case.Figure 5
Distribution of singular values of reference dictionaryΩ, Ω
A
K
S
V
D, and Ω
P
r
o
p
o
s
e
d for the noisy case.As observed, the ratio between the maximal and the minimal singular values ofΩ
P
r
o
p
o
s
e
d is almost the same as that of Ω, smaller than the one of Ω
A
K
S
V
D. This implies that our proposed algorithm, as expected, yields a dictionary closer to a TF than Ω
A
K
S
V
D owing to the penalty term F
(
Ω
~
).Performance Evaluation of the Obtained Dictionaries. The averaged representation error is defined as Y
-
X
k
F
/
L
N for the dictionary Ω
k, where X
k is the estimate of the clean signals which are available for A
l
g
A
K
S
V
D and can be obtained using the OBG for A
l
g
P
r
o
p
o
s
e
d.The performance of a given dictionary for denoising is usually measured with the averaged denoising error defined asX
k
-
X
F
/
L
N, where X is the true (clean) signals.Figures6 and 7 show comparison of the two algorithms in terms of the averaged representation error and the averaged denoising error, respectively. (In Figure 7, the plots for the noise-free case are not presented as they are the same as those shown in Figure 6 due to the fact that Y
=
X.)Figure 6
Evolution of the averaged representation errorY
-
X
k
F
/
N
L versus iterations with J
=
50, N
=
25, and L
=
50000.Figure 7
Evolution of the average denoising errorX
k
-
X
F
/
L
N versus iterations with J
=
50, N
=
25, and L
=
50000.Now, let us look at how closeΩ
P
r
o
p
o
s
e
d and Ω
A
K
S
V
D are to the reference dictionary Ω.We say a rowω
m
† in the true analysis dictionary Ω is recovered in Ω
^ if (40)
m
i
n
n
1
-
ω
^
n
†
ω
m
†
T
<
0.01
.Denote J
k as the number of recovered rows in Ω
k. The percentage of the recovered rows denoted as ρ
p
r
r
(
k
) in Ω
k is calculated with ρ
p
r
r
(
k
)
≜
J
k
/
J. This is one of the performance indicators used in [12] for evaluating a learnt dictionary.With signals generated above and the same setup as the one used in [12] except N
i
t
e
=
200 instead of 100, we run each of the A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D and obtain a sequence of N
i
t
e analysis dictionaries {
Ω
k
}. Figure 8 shows the evolution of the percentage of the recovered rows for the two algorithms.Figure 8
Evolution of the percentage of the recovered rowsρ
p
r
r
(
k
) versus iterations with J
=
50, N
=
25, and L
=
50000.Comment 6.
(i) It is seen that our proposed algorithm surprisingly outperformsA
l
g
A
K
S
V
D in terms of averaged representation error, averaged denoising error, and percentage of recovered rows for both cases.
(ii) According to equation (6) of [12], the averaged denoising error in the oracle setup (i.e., the true dictionary and all the cosupports are known) is given by r
/
N
σ, which is equal to 0.4
σ for this example. As observed from Figure 7, the averaged denoising error for both algorithms seems to converge to a noise variance close to σ
=
0.04. This indicates these two algorithms just yield one of the local minimal points of the problems though very close to the reference one.Figures9–11, as another example, display the results for the case of N
=
20, J
=
40, and L
=
10000 and the vector x
i living in a 3-dimensional subspace, leading to κ
=
N
-
3
=
17. Here, we set λ
=
200 and ν
=
0.01 for our proposed algorithm. The results are self-explanatory.Figure 9
Evolution of the averaged representation errorY
-
X
k
F
/
N
L versus iterations with J
=
40, N
=
20, and L
=
10000.Figure 10
Evolution of the average denoising errorX
k
-
X
F
/
L
N versus iterations with J
=
40, N
=
20, and L
=
10000.Figure 11
Evolution of the percentage of the recovered rowsρ
p
r
r
(
k
) versus iterations with J
=
40, N
=
20, and L
=
10000.What is the most significant improvement of our proposed algorithm overA
l
g
A
K
S
V
D is the implementation efficiency. This can be seen from Figure 12 which yields the relationship between the time (in seconds) used and the number of training samples L for the setting: N
=
25, J
=
50, and N
i
t
e
=
200. As an example, one notes that, for L
=
50000, A
l
g
A
K
S
V
D requires about 1
0
5 seconds, while our proposed algorithm needs less than 3
×
1
0
2 seconds.Figure 12
Relationship between the time (in seconds) consumed and the number of the training samplesL.


### 4.2. Experiments on Image Denoising

Now, let us consider image denoising for five widely used images, “Lenna,” “Barbara,” “Boats,” “House,” and “Peppers.”For a given noisy image with noise levelσ
∈
{
5,10,15,20
} (which means each pixel is contaminated by a white Gaussian noise with i.i.d. zero-mean and σ
2-variance), we generate a training set consisting of L = 20,000 overlapping 8
×
8 patch examples. Then, we run A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D with λ
=
100, ν
=
0.01, r
=
7, κ
=
30, and N
i
t
e
=
20 to learn the analysis dictionaries of size J
×
N
=
128
×
64.We apply the learnt dictionaries for patch-based image denoising. Each noised patch is denoised by using error-based OBG for sparse coding, which is the same as that in [12]. The denoising results with respect to image PSNR in dB are given in Table 1.Table 1
Denoising performance (image PSNR in dB) evaluated with five different pictures in different noise levelσ for A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D.
σ
Algorithm
Lenna
Barbara
Boats
House
Peppers 256
5
A
l
g
A
K
S
V
D
37.32
36.88
36.60
38.04
36.89
A
l
g
P
r
o
p
o
s
e
d
37.54
37.11
36.65
38.11
36.98
10
A
l
g
A
K
S
V
D
33.35
32.00
32.07
33.89
32.11
A
l
g
P
r
o
p
o
s
e
d
33.46
32.85
32.07
33.88
32.29
15
A
l
g
A
K
S
V
D
31.05
29.22
29.62
31.35
29.53
A
l
g
P
r
o
p
o
s
e
d
31.30
29.73
29.91
31.60
29.94
20
A
l
g
A
K
S
V
D
29.45
27.30
27.92
29.58
27.79
A
l
g
P
r
o
p
o
s
e
d
29.49
28.11
28.54
29.64
28.18It is observed that denoising performance byA
l
g
P
r
o
p
o
s
e
d is in general better than that by A
l
g
A
K
S
V
D in terms of PSNR.Figure13 shows the visual effects of denoising on the image “Barbara” with a noise level of σ
=
10.Figure 13
“Barbara.” (a) The original; (b) noisy image with noiseσ
=
10 (PSNR = 28.13 dB); (c) denoised image using Ω
A
K
S
V
D (PSNR = 32.00 dB); (d) denoised image using Ω
P
r
o
p
o
s
e
d (PSNR = 32.85 dB).
(a)
(b)
(c)
(d)


## 4.1. Synthetic Experiments

In order to have a fair comparison, the setup used here is similar to [12]. Let J
=
50, N
=
25, and L = 50,000. A J
×
N matrix W is generated with normally distributed entries and then Ω
=
D
α
W, where D
α is such a diagonal matrix that all the rows of Ω have unit l
2-norm. This Ω serves as the true overcomplete analysis dictionary. Then, a set of L vectors {
x
i
} of N
×
1 is produced with each living in a 4-dimensional subspace. Each vector x
i is obtained as follows: we first construct a submatrix Ω
(
i
)
∈
R
21
×
N by randomly choosing 21 rows from Ω and then compute the orthonormal basis {
ψ
k
(
i
)
} for the null space of Ω
(
i
), which can be obtained from the SVD of Ω
(
i
). With Ψ
(
i
)
≜
ψ
1
(
i
)
ψ
2
(
i
)
ψ
3
(
i
)
ψ
4
(
i
), calculate x
i
=
Ψ
(
i
)
s
i, where each element in s
i is randomly generated by Gaussian distribution of independent and identically distributed (i.i.d.) zero-mean and unit variance. At last, each vector x
i is normalized so that x
i
2
=
1 for all i
=
1,2
,
…
,
L. We can see that Ω
x
i has at least 21 zero components for all {
x
i
} (i
=
1,2
,
…
,
L) generated in such a way.LetI
κ
:
R
J
→
R
J denote a hard thresholding operator that keeps the largest (in magnitude) κ elements of a vector. We define the effective sparsity κ
~ of a length-J vector s as (38)
κ
~
=
m
i
n
κ
s.t.
I
κ
s
2
≥
0.99
s
2
.Figure 1 represents the histogram of the effective sparsity level of {
Ω
x
i
}
(
i
=
1,2
,
…
,
L
). As can be seen, the effective sparsity level is varying from 17 to 25and is most concentrated near 21.Figure 1
Histogram of the effective sparsity of{
Ω
x
i
}
(
i
=
1,2
,
…
,
L
) for J
=
50, N
=
25, and L
=
50000.The measurements{
y
i
} are generated with (39)
y
i
=
x
i
+
ϵ
i
,
∀
i
,where each entry of the additive noise vector ϵ
i is a white Gaussian noise with i.i.d. zero-mean and σ
2-variance. Define the signal-to-noise ratio (SNR) as ρ
s
n
r
≜
10
l
o
g
10
(
E
x
i
2
2
/
E
ϵ
i
2
2
), where E
[
·
] denotes the mathematical expectation operator. When σ
=
0.2
/
N
=
0.04, ρ
s
n
r
=
10
l
o
g
10
(
1
/
(
N
∗
σ
2
)
)
=
10
l
o
g
10
(
25
)
≈
14dB.We then estimate the analysis dictionaryΩ using A
l
g
A
K
S
V
D and A
l
g
P
r
o
p
o
s
e
d with X
~
=
Y (see (12)), respectively, for both noise-free case (i.e., σ
=
0) and noisy case with noise level ρ
s
n
r
=
14dB. In both cases, we fix κ
=
21 and N
i
t
e
=
200, while ν
=
0.01 and λ
=
300 for A
l
g
P
r
o
p
o
s
e
d. (As we said before, Ω
x
i has at least 21 zero components for all {
x
i
}. Thus, we can choose κ such that κ
≤
J
-
21
=
29. It is illustrated in Figure 1 that the effective sparsity is concentrated near 21, and therefore we choose κ
=
21.) Both A
l
g
P
r
o
p
o
s
e
d and AlgA
K
S
V
D are initialized with a dictionary Ω
0 in which each row is generated as follows: we first randomly choose N
-
1 columns from the training data matrix Y to form a submatrix Y
-; then set ω
m
†
=
e
-
e
Y
-
Y
-
‡, where e
=
r
a
n
d
n
(
1
,
N
) and Y
-
‡, as defined before, denotes the pseudoinverse of Y
-; at last we normalize each row to unit l
2 norm [12].Convergence Behavior of
A
l
g
P
r
o
p
o
s
e
d. The convergence behavior of our proposed A
l
g
P
r
o
p
o
s
e
d is shown in Figures 2 and 3 in terms of the objective function L
(
Ω
k
,
S
k
) and the relative change of iterates which is defined as Ω
k
-
Ω
k
-
1
F
/
Ω
k
-
1
F.Figure 2
Evolution of the objective functionL
(
Ω
k
,
S
k
) for J
=
50, N
=
25, and L
=
50000.Figure 3
Evolution of relative change of iteratesΩ
k
-
Ω
k
-
1
F
/
Ω
k
-
1
F for J
=
50, N
=
25, and L
=
50000.Comment 5.
(i) The results are self-explanatory. As expected, the objective functionL
(
Ω
k
,
S
k
) decreases smoothly in both noise-free case and noisy case. This is confirmed in Figure 2, from which it is seen that the latter case yields larger L
(
Ω
k
,
S
k
) than the noise-free case because of the larger sparsification error S
k
-
Ω
k
Y
F
2 in the latter case.
(ii) The convergence behavior of the iteratesΩ
k is depicted in Figure 3. One observes that the iterates of the proposed algorithm show a very good trend to convergence in both cases.Figures4 and 5 show the distribution of singular values of the two analysis dictionaries Ω
A
K
S
V
D and Ω
P
r
o
p
o
s
e
d.Figure 4
Distribution of singular values of reference dictionaryΩ, Ω
A
K
S
V
D, and Ω
P
r
o
p
o
s
e
d for the noise-free case.Figure 5
Distribution of singular values of reference dictionaryΩ, Ω
A
K
S
V
D, and Ω
P
r
o
p
o
s
e
d for the noisy case.As observed, the ratio between the maximal and the minimal singular values ofΩ
P
r
o
p
o
s
e
d is almost the same as that of Ω, smaller than the one of Ω
A
K
S
V
D. This implies that our proposed algorithm, as expected, yields a dictionary closer to a TF than Ω
A
K
S
V
D owing to the penalty term F
(
Ω
~
).Performance Evaluation of the Obtained Dictionaries. The averaged representation error is defined as Y
-
X
k
F
/
L
N for the dictionary Ω
k, where X
k is the estimate of the clean signals which are available for A
l
g
A
K
S
V
D and can be obtained using the OBG for A
l
g
P
r
o
p
o
s
e
d.The performance of a given dictionary for denoising is usually measured with the averaged denoising error defined asX
k
-
X
F
/
L
N, where X is the true (clean) signals.Figures6 and 7 show comparison of the two algorithms in terms of the averaged representation error and the averaged denoising error, respectively. (In Figure 7, the plots for the noise-free case are not presented as they are the same as those shown in Figure 6 due to the fact that Y
=
X.)Figure 6
Evolution of the averaged representation errorY
-
X
k
F
/
N
L versus iterations with J
=
50, N
=
25, and L
=
50000.Figure 7
Evolution of the average denoising errorX
k
-
X
F
/
L
N versus iterations with J
=
50, N
=
25, and L
=
50000.Now, let us look at how closeΩ
P
r
o
p
o
s
e
d and Ω
A
K
S
V
D are to the reference dictionary Ω.We say a rowω
m
† in the true analysis dictionary Ω is recovered in Ω
^ if (40)
m
i
n
n
1
-
ω
^
n
†
ω
m
†
T
<
0.01
.Denote J
k as the number of recovered rows in Ω
k. The percentage of the recovered rows denoted as ρ
p
r
r
(
k
) in Ω
k is calculated with ρ
p
r
r
(
k
)
≜
J
k
/
J. This is one of the performance indicators used in [12] for evaluating a learnt dictionary.With signals generated above and the same setup as the one used in [12] except N
i
t
e
=
200 instead of 100, we run each of the A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D and obtain a sequence of N
i
t
e analysis dictionaries {
Ω
k
}. Figure 8 shows the evolution of the percentage of the recovered rows for the two algorithms.Figure 8
Evolution of the percentage of the recovered rowsρ
p
r
r
(
k
) versus iterations with J
=
50, N
=
25, and L
=
50000.Comment 6.
(i) It is seen that our proposed algorithm surprisingly outperformsA
l
g
A
K
S
V
D in terms of averaged representation error, averaged denoising error, and percentage of recovered rows for both cases.
(ii) According to equation (6) of [12], the averaged denoising error in the oracle setup (i.e., the true dictionary and all the cosupports are known) is given by r
/
N
σ, which is equal to 0.4
σ for this example. As observed from Figure 7, the averaged denoising error for both algorithms seems to converge to a noise variance close to σ
=
0.04. This indicates these two algorithms just yield one of the local minimal points of the problems though very close to the reference one.Figures9–11, as another example, display the results for the case of N
=
20, J
=
40, and L
=
10000 and the vector x
i living in a 3-dimensional subspace, leading to κ
=
N
-
3
=
17. Here, we set λ
=
200 and ν
=
0.01 for our proposed algorithm. The results are self-explanatory.Figure 9
Evolution of the averaged representation errorY
-
X
k
F
/
N
L versus iterations with J
=
40, N
=
20, and L
=
10000.Figure 10
Evolution of the average denoising errorX
k
-
X
F
/
L
N versus iterations with J
=
40, N
=
20, and L
=
10000.Figure 11
Evolution of the percentage of the recovered rowsρ
p
r
r
(
k
) versus iterations with J
=
40, N
=
20, and L
=
10000.What is the most significant improvement of our proposed algorithm overA
l
g
A
K
S
V
D is the implementation efficiency. This can be seen from Figure 12 which yields the relationship between the time (in seconds) used and the number of training samples L for the setting: N
=
25, J
=
50, and N
i
t
e
=
200. As an example, one notes that, for L
=
50000, A
l
g
A
K
S
V
D requires about 1
0
5 seconds, while our proposed algorithm needs less than 3
×
1
0
2 seconds.Figure 12
Relationship between the time (in seconds) consumed and the number of the training samplesL.


## 4.2. Experiments on Image Denoising

Now, let us consider image denoising for five widely used images, “Lenna,” “Barbara,” “Boats,” “House,” and “Peppers.”For a given noisy image with noise levelσ
∈
{
5,10,15,20
} (which means each pixel is contaminated by a white Gaussian noise with i.i.d. zero-mean and σ
2-variance), we generate a training set consisting of L = 20,000 overlapping 8
×
8 patch examples. Then, we run A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D with λ
=
100, ν
=
0.01, r
=
7, κ
=
30, and N
i
t
e
=
20 to learn the analysis dictionaries of size J
×
N
=
128
×
64.We apply the learnt dictionaries for patch-based image denoising. Each noised patch is denoised by using error-based OBG for sparse coding, which is the same as that in [12]. The denoising results with respect to image PSNR in dB are given in Table 1.Table 1
Denoising performance (image PSNR in dB) evaluated with five different pictures in different noise levelσ for A
l
g
P
r
o
p
o
s
e
d and A
l
g
A
K
S
V
D.
σ
Algorithm
Lenna
Barbara
Boats
House
Peppers 256
5
A
l
g
A
K
S
V
D
37.32
36.88
36.60
38.04
36.89
A
l
g
P
r
o
p
o
s
e
d
37.54
37.11
36.65
38.11
36.98
10
A
l
g
A
K
S
V
D
33.35
32.00
32.07
33.89
32.11
A
l
g
P
r
o
p
o
s
e
d
33.46
32.85
32.07
33.88
32.29
15
A
l
g
A
K
S
V
D
31.05
29.22
29.62
31.35
29.53
A
l
g
P
r
o
p
o
s
e
d
31.30
29.73
29.91
31.60
29.94
20
A
l
g
A
K
S
V
D
29.45
27.30
27.92
29.58
27.79
A
l
g
P
r
o
p
o
s
e
d
29.49
28.11
28.54
29.64
28.18It is observed that denoising performance byA
l
g
P
r
o
p
o
s
e
d is in general better than that by A
l
g
A
K
S
V
D in terms of PSNR.Figure13 shows the visual effects of denoising on the image “Barbara” with a noise level of σ
=
10.Figure 13
“Barbara.” (a) The original; (b) noisy image with noiseσ
=
10 (PSNR = 28.13 dB); (c) denoised image using Ω
A
K
S
V
D (PSNR = 32.00 dB); (d) denoised image using Ω
P
r
o
p
o
s
e
d (PSNR = 32.85 dB).
(a)
(b)
(c)
(d)


## 5. Conclusions

In this paper, we have investigated the problem of overcomplete transform learning for analysis model. An alternating minimization based iterative procedure has been proposed to solve the formulated optimal transform design problem. A closed-form solution has been derived for updating the transform. The superiority of the proposed approach in terms of the averaged representation and denoising errors, the percentage of the recovered rows, and, more significantly, the computation efficiency has been demonstrated with a series of experiments using both synthetic and real data.As observed, both the AKSVD and our proposed algorithm demonstrate good convergence to a local minimum. Deeper studies along those lines [27, 28] are needed for theoretical guarantees of convergence to a global minimum and developing high performance algorithms for sparsifying dictionary and transform design. Furthermore, it is expected that the performance of the analysis KSVD algorithms in [12] can be improved with the proposed approach for updating dictionary embedded. Investigations along these lines are ongoing.



---

*Source: 6329740-2016-07-19.xml*