$\ell_0$ Sparsifying Transform Learning with Efficient Optimal Updates and Convergence Guarantees

0.41 MB·English

by  Saiprasad Ravishankar

1
ℓ
Sparsifying Transform Learning with Efficient
0
Optimal Updates and Convergence Guarantees
Saiprasad Ravishankar, Student Member, IEEE, and Yoram Bresler, Fellow, IEEE
Abstract—Manyapplicationsinsignal processingbenefitfrom The transform model is not only more general in its
5 the sparsity of signals in a certain transform domain or dictio- modeling capabilities than the analysis models, it is also
1 nary.Synthesissparsifyingdictionariesthataredirectlyadapted much more efficient and scalable than both the synthesis and
0 to data have been popularin applicationssuch as image denois- noisy signal analysis models. We briefly review the main
2 ing, inpainting, and medical image reconstruction. In this work,
we focus instead on the sparsifying transform model, and study distinctions between these sparse models (cf. [5] for a more
n the learning of well-conditioned square sparsifying transforms. detailed review, and for the relevant references) in this and
a The proposed algorithms alternate between a ℓ0 “norm”-based thefollowingparagraphs.Onekeydifferenceisin theprocess
J
sparsecodingstep,andanon-convextransformupdatestep.We
of finding a sparse representation for data given the model,
3 derive the exact analytical solution for each of these steps. The
or dictionary. For the transform model, given the signal y
1 proposed solution for the transform update step achieves the
global minimum in that step, and also provides speedups over and transform W, the transform sparse coding problem [5]
2
] iterative solutions involving conjugate gradients. We establish minimizes Wy x subject to x s, where s is a
L k − k2 k k0 ≤
thatouralternatingalgorithmsaregloballyconvergenttotheset given sparsity level. The solution xˆ is obtained exactly and
M of local minimizers of the non-convex transform learning prob-
cheaply by zeroing out all but the s coefficients of largest
lems. In practice, the algorithms are insensitive to initialization.
. magnitude in Wy 1. In contrast, for the synthesis or noisy
t We present results illustrating the promising performance and
a significantspeed-upsoftransformlearningoversynthesisK-SVD analysis models, the process of sparse coding is NP-hard
t
s in image denoising. (Non-deterministicPolynomial-timehard).While some of the
[ approximatealgorithmsthathave been proposedforsynthesis
Index Terms—Transform Model, Fast Algorithms, Image rep-
1 resentation, Sparse representation, Denoising, Dictionary learn- oranalysissparsecodingareguaranteedtoprovidethecorrect
v ing, Non-convex. solution under certain conditions, in applications, especially
9 thoseinvolvinglearningthemodelsfromdata,theseconditions
5
areoftenviolated.Moreover,thevarioussynthesisandanalysis
8 I. INTRODUCTION sparsecodingalgorithmstendtobecomputationallyexpensive
2
0 The sparsity of signals and images in a certain transform for large-scale problems.
. domain or dictionary has been widely exploited in numerous Recently,thedata-drivenadaptationofsparsemodelshasre-
1
0 applications in recent years. While transforms are a classical ceived much interest. The adaptation of synthesis dictionaries
5 tool in signal processing, alternative models have also been basedontrainingsignals[6]–[12]hasbeenshowntobeuseful
1 studied for sparse representation of data, most notably the in various applications [13]–[15]. The learning of analysis
v: popularsynthesismodel[1],[2],theanalysismodel[1]andits dictionaries, employing either the analysis model or its noisy
i morerealisticextension,thenoisysignalanalysismodel[3].In signal extension, has also received some recent attention [3],
X
this paper, we focus specifically on the sparsifying transform [16]–[18].
r
model [4], [5], which is a generalized analysis model, and Focusing instead on the transform model, we recently
a
suggests that a signal y Rn is approximately sparsifiable developed the following formulation [5] for the learning of
using a transform W R∈m×n, that is Wy = x+e where well-conditioned square sparsifying transforms. Given a ma-
x Rm is sparse in s∈ome sense, and e is a small residual. trix Y Rn×N, whose columns represent training signals,
∈ ∈
A distinguishing feature is that, unlike the synthesis or noisy our formulation for learning a square sparsifying transform
signal analysis models, where the residual is measured in the W Rn×n for Y is [5]
∈
signal domain, in the transform model the residual is in the
2 2
transform domain. (P0) min WY X +λ ξ W log detW
W,X k − kF k kF − | |
Copyright (c) 2014 IEEE. Personal use of this material is permitted. s.t. kXik0 ≤s ∀ i (cid:0) (cid:1)
However, permission to use this material for any other purposes must be
obtained [email protected]. where λ > 0, ξ > 0 are parameters, and X Rn×N is
ThisworkwasfirstsubmittedtotheIEEETransactionsonSignalProcess- ∈
a matrix, whose columns X are the sparse codes of the
ing on April 02, 2014. It was accepted in the current form on January 09, i
2
2015. corresponding training signals Yi. The term kWY −XkF in
ThisworkwassupportedinpartbytheNationalScienceFoundation(NSF) (P0) is called sparsification error, and denotes the deviation
undergrants CCF-1018660andCCF-1320953.
S. Ravishankar and Y. Bresler are with the Department of Electrical
and Computer Engineering and the Coordinated Science Laboratory, Uni- 1Moreover,givenW andsparsecodex,wecanalsorecoveraleastsquares
versity of Illinois, Urbana-Champaign, IL, 61801 USA e-mail: (ravisha3, estimate of the underlying signal as yˆ = W†x, where W† is the pseudo-
ybresler)@illinois.edu. inverseofW.
2
of the data in the transform domain from its sparse approxi- the positivity constraint by the absolute value in the for-
2
mation (i.e., the deviation of WY from the sparse matrix X). mulation in this paper. The W penalty in (P0) helps
Problem (P0) also has v(W) , log detW +ξ W 2 as remove a ‘scale ambiguity’ ikn thkeF solution, which occurs
− | | k kF
a regularizer in the objective to prevent trivial solutions. We when the data admit an exactly sparse representation[5]. The
2
have proposed an alternating algorithm algorithm for solving log detW and W penaltiestogether additionallyhelp
− | | k kF
(P0) [5], that alternates between updating X (sparse coding), control the condition number κ(W) of the learnt transform.
and W (transform update), with the other variable kept fixed. (Recall that the condition number of a matrix A Rn×n
∈
Because of the simplicity of sparse codingin the transform is defined as κ(A) = β1/βn, where β1 and βn denote the
model, the alternating algorithm for transform learning [5] largest and smallest singular values of A, respectively.) In
has a low computational cost. On the other hand, because, particular,badlyconditionedtransformstypicallyconveylittle
in the case of the synthesis or noisy signal analysis models, information and may degrade performance in applications
the learning formulations involve the NP-hard sparse coding, such as signal/image representation, and denoising [5]. Well-
such learning formulations are also NP hard. Moreover, even conditionedtransforms,ontheotherhand,havebeenshownto
when the ℓ0 sparse coding is approximated by a convex perform well in (sparse) image representation, and denoising
relaxation,thevarioussynthesisoranalysisdictionarylearning [5], [19].
problemsremainhighlynon-convex.Becausetheapproximate The condition number κ(W) can be upper bounded by a
algorithmsfor these problemsusually solve the sparse coding monotonically increasing function of v(W) (see Proposition
problem repeatedly in the iterative process of adapting the 1 of [5]). Hence, minimizing v(W) encourages reduction of
sparsemodel,thecostofsparsecodingismultipliedmanyfold. theconditionnumber.Theregularizerv(W)alsopenalizesbad
Hence,thesynthesisoranalysisdictionarylearningalgorithms scalings.GivenatransformW andascalarα R,v(αW)
∈ →
tendtobecomputationallyexpensiveinpracticeforlargescale as the scaling α 0 or α . For a fixed ξ, as λ
∞ → → ∞
problems. Finally, popular algorithmsfor synthesis dictionary is increased in (P0), the optimal transform(s) become well-
learningsuchasK-SVD [9],oralgorithmsforanalysisdictio- conditioned.Inthelimitλ ,theirconditionnumbertends
→∞
nary learning do not have convergenceguarantees. to 1, and their spectral norm (or, scaling) tends to 1/√2ξ.
In this follow-on work on transform learning, keeping the Specifically, for ξ = 0.5, as λ , the optimal transform
→ ∞
focus on the square transform learning formulation (P0), we tendstoanorthonormaltransform.Inpractice,thetransforms
make the following contributions. learntvia(P0)haveconditionnumbersverycloseto1evenfor
• We derive highly efficient closed-form solutions for the finiteλ[5].Thespecificchoiceofλdependsontheapplication
updatestepsinthealternatingminimizationprocedurefor and desired condition number.
(P0),thatfurtherenhancethecomputationalpropertiesof Inthispaper,toachieveinvarianceofthe learnedtransform
2
transform learning. to trivial scaling of the training data Y, we set λ=λ0kYkF
• We also consider the alternating minimization of an al- in (P0), where λ0 > 0 is a constant. Indeed, when the data
ternativeversionof(P0)in thispaper,thatisobtainedby Y are replaced with αY (α R, α = 0) in (P0), we
∈ 6
replacing the sparsity constraints with sparsity penalties. can set X = αX′. Then, the objective function becomes
• Importantly,we establishthatouriterativealgorithmsfor α2 kWY −X′k2F +λ0kYk2F v(W) , which is just a scaled
transform learning are globally convergent to the set of version of the objective in (P0) (for un-scaled Y). Hence,
(cid:0) (cid:1)
local minimizers of the non-convex transform learning its minimization over (W,X′) (with X′ constrained to have
problems. columns of sparsity s) yields the same solution(s) as (P0).
≤
Thus, the learnt transform for data αY is the same as for Y,
We organize the rest of this paper as follows. Section II
while the learnt sparse code for αY is α times that for Y.
briefly describes our transform learning formulations. In Sec-
We have shown [5] that the cost function in (P0) is lower
tion III, we derive efficient algorithms for transform learning,
anddiscussthealgorithms’computationalcost.InSectionIV, bounded by λv0, where v0 = n2 + n2 log(2ξ). The minimum
objective value in Problem (P0) equals this lower bound if
we present convergence guarantees for our algorithms. The
and only if there exists a pair (Wˆ,Xˆ) such that WˆY = Xˆ,
proof of convergenceis provided in the Appendix. Section V
with Xˆ Rn×N whose columns have sparsity s, and
presents experimental results demonstrating the convergence Wˆ Rn∈×n whose singular values are all equal t≤o 1/√2ξ
behavior, and the computational efficiency of the proposed (hen∈ce, the condition number κ(Wˆ) = 1). Thus, when an
scheme. We also show brief results for the image denoising
“error-free” transform model exists for the data, and the
application. In Section VI, we conclude.
underlying transform is unit conditioned, such a transform
modelisguaranteedtobeaglobalminimizerofProblem(P0)
II. LEARNING FORMULATIONS AND PROPERTIES (i.e., such a model is identifiable by solving (P0)). Therefore,
The transform learning Problem (P0) was introduced in it makes sense to solve (P0) to find such good models.
Section I. Here, we discuss some of its important properties. Another interesting property of Problem (P0) is that it ad-
The regularizer v(W) , log detW + ξ W 2 helps mitsanequivalenceclassofsolutions/minimizers.Becausethe
− | | k kF
prevent trivial solutions in (P0). The log detW penalty objectivein (P0) is unitarilyinvariant,then givena minimizer
eliminates degenerate solutions such as thos|e with|zero, or (W˜,X˜), the pair (ΘW˜,ΘX˜) is another equivalentminimizer
repeatedrows.WhileitissufficienttoconsiderthedetW >0 for all sparsity-preserving orthonormal matrices Θ, i.e., Θ
case [5], to simplify the algorithmic derivation we replace such that ΘX˜i 0 s i. For example, Θ can be a row
k k ≤ ∀
3
permutation matrix, or a diagonal 1 sign matrix. In the case of Problem (P1), we solve the following sparse
±
We note that a cost function similar to that in (P0), but coding problem
2
lacking the W penalty has been derived under certain
k kF N
assumptions in the very different setting of blind source min WY X 2 + η2 X (2)
separation [20]. However, the transforms learnt via Problem X k − kF i k ik0
i=1
(P0) perform poorly [5] in signal processing applications, X
2 AsolutionXˆ of(2)inthiscaseisobtainedasXˆ =Hˆ1(WY )
wwhheicnhthaesldeiasrcnuinssgedis,dhoenlpesexocvleurdcionmgethtehecrusccaialelkaWmbkiFgupietynaalntyd, ∀i,wherethe(hard-thresholding)operatorHˆη1(i·)isdηeifinedais
control the condition number. follows.
0 ,if b <η
In this work, we also consider an alternative version of Hˆ1(b) = | j| (3)
Problem (P0) by replacing the ℓ0 sparsity constraints with ℓ0 η j (bj ,if bj η
(cid:16) (cid:17) | |≥
penaltiesintheobjective(thisversionofthetransformlearning
Here, b Rn, and the subscriptj indexesvector entries. This
problem has been recently used for example in adaptive ∈
form of the solution to (2) has been mentioned in prior work
tomographicreconstruction[21],[22]).Inthiscase,we obtain
[21]. For completeness, we include a brief proof in Appendix
the following unconstrained(or, sparsity penalized) transform
A. When the condition (WY) = η occurs for some i
learning problem and j (where (WY) i|s the ejlie|ment iof WY on the jth
ji
N row and ith column), the corresponding optimal Xˆ in (2)
ji
(P1) Wm,iXn kWY −Xk2F +λv(W)+ ηi2kXik0 can be either (WY)ji or 0 (both of which correspond to the
i=1 minimumvalueofthecostin(2)).Thedefinitionin(3)breaks
X
whereη2, withη >0 i,denotetheweights(e.g.,η =η i the tie between these equally valid solutions by selecting the
i i ∀ i ∀ first. Thus, similar to Problem (1), the solution to (2) can be
for some η) for the sparsity penalties. The various aforemen-
computed exactly.
tioned properties for Problem (P0) can be easily extended to
2) Transform Update Step: The transform update step of
the case of the alternative Problem (P1).
(P0) or (P1) involvesthe followingunconstrainednon-convex
[5] minimization.
III. TRANSFORM LEARNING ALGORITHM
2 2
min WY X +λξ W λlog detW (4)
A. Algorithm W k − kF k kF − | |
NotethatalthoughNLCGworkswellforthetransformupdate
We have previously proposed [5] an alternating algorithm
step[5],convergencetotheglobalminimumofthenon-convex
forsolving(P0)thatalternatesbetweensolvingforX (sparse
transform update step has not been proved with NLCG. In-
coding step) and W (transform update step), with the other
stead,replacingNLCG,thefollowingpropositionprovidesthe
variable keptfixed. While the sparse codingstep has an exact
closed-formsolutionforProblem(4).Thesolutioniswrittenin
solution,thetransformupdatestepwasperformedusingitera-
terms of an appropriate singular value decomposition (SVD).
tive nonlinear conjugate gradients (NLCG). This alternating
algorithm for transform learning has a low computational Weuse()T todenotethematrixtransposeoperation,andM12
·
todenotethepositivedefinitesquarerootofapositivedefinite
costcomparedtosynthesis/analysisdictionarylearning.Inthe
following, we provide a further improvement: we show that matrix M. We let I denote the n n identity matrix.
×
both steps of transform learning (for either (P0) or (P1)) can Proposition 1: Given the training data Y Rn×N, sparse
in fact, be performed exactly and cheaply. code matrix X Rn×N, and λ > 0, ξ∈> 0, factorize
1) Sparse Coding Step: The sparse coding step in the YYT +λξI as L∈LT, with L Rn×n. Further, let L−1YXT
∈
alternating algorithm for (P0) is as follows [5] have a full SVD of QΣRT. Then, a global minimizer for the
transform update step (4) can be written as
2
min WY X s.t. X s i (1)
X k − kF k ik0 ≤ ∀ Wˆ =0.5R Σ+ Σ2+2λI 12 QTL−1 (5)
The above problem is to project WY onto the (non-convex)
set of matrices whose columns have sparsity s. Due to Thesolutionisunique(cid:16)ifand(cid:0)onlyifL−(cid:1)1Y(cid:17)XT isnon-singular.
the additivity of the objective, this correspond≤s to project- Furthermore, the solution is invariant to the choice of factor
ing each column of WY onto the set of sparse vectors L.
v{coexcm∈topruRtzned∈:bkRyxnkz,0etr≤hoeisno}gp,otwiumhtaicalhllpwrbouejtecctathileolntshzeˆcoso-enℓft0fiocbtiahelnelt.ssN-oℓof0wlb,aarfgollerisast asλtlroPgrWodoef:tYWTYhTe. W+obeλjteξhcIetinvWedeTfcuonmc−tpioo2nsetritn(hWe(4pY)oXsciatTinv)eb-+edetrfiren(-XiwterXimttTean)-
− (cid:8)| (cid:0) | (cid:1) (cid:9)
magnitude in z. If there is more than one choice for the s trix YYT +λξI as LLT (e.g., L can be the positive-definite
coefficientsoflargestmagnitudeinz(canoccurwhenmultiple square root, or the cholesky factor of YYT + λξI). The
entries in z have identical magnitude), then the optimal zˆ is objective function then simplifies as follows
not unique. We then choose zˆ= Hs(z), where Hs(z) is the tr WLLTWT 2WYXT +XXT λlog detW
projection, for which the indices of the s largest magnitude − − | |
elements (in z) are the lowest possible. Hence, an optimal Using(cid:0)a change of variables B = WL(cid:1), the multiplicativity
sparse code in (1) is computed as Xˆ =H (WY ) i. of the determinant implies log detB = log detW +
i s i
∀ | | | |
4
log detL. Problem (4) is then equivalent to Q , R or Q ,R , yielding equally valid alternative
k k k k
| | SVDs−ofL−1YX−T. However,becausezerosingularvaluesin
min tr BBT 2tr BL−1YXT λlog detB (6) (cid:0) (cid:1) (cid:0) (cid:1)
B − − | | the matrix Σ are mapped to non-zero singular values in the
Next, let(cid:0)B =(cid:1)UΓVT(cid:0), and L−1YX(cid:1) T = QΣRT be full matrixΓ(by(8)),wehavethatthefollowingtwomatricesare
equally valid solutions to (4).
SVDs (U,Γ,V,Q,Σ,R are all n n matrices), with γ and
i
σ denotingthe diagonalentriesof×Γ andΣ, respectively.The Wˆa = γ R QT +γ R QT L−1 (10)
i i6=k i i i k k k
unconstrained minimization (6) then becomes Wˆb =(cid:0)Pi6=kγiRiQTi −γkRkQTk(cid:1)L−1 (11)
n
mΓin "tr Γ2 −2 mU,aVx tr UΓVTQΣRT −λ log γi# wtrahnesrfeorγmk >is0n.oItt(cid:0)iusPnoiqbuveiouins tthhaist Wcˆaase6=. TWˆhebr,e(cid:1)if.oer.e,,thLe−o1pYtiXmaTl
i=1
(cid:0) (cid:1) (cid:8) (cid:0) (cid:1)(cid:9) X being non-singular is a necessary condition for uniqueness of
For the inner maximization, we use the result
the solution to (4).
max tr UΓVTQΣRT = tr(ΓΣ) [23], where the
U,V Next,we showthatthe non-singularityof L−1YXT isalso
maximum is attained by setting U = R and V = Q. The
(cid:0) (cid:1) asufficientconditionfortheaforementioneduniqueness.First,
remaining minimization with respect to Γ is then
ifthesingularvaluesofL−1YXT arenon-degenerate(distinct
n n n andnon-zero),thentheSVDofL−1YXT isuniqueuptojoint
min γ2 2 γ σ λ log γ (7)
{γi} Xi=1 i − Xi=1 i i− Xi=1 i tshcaatlinthgeosfoaluntyiopnaiWrˆ Q=i,Rin byγ±R1.QTThisLi−m1miesduiantieqluyeiminptlhieiss
Thisproblemisconvexinthenon-negativesingularvaluesγi, case. On the other(cid:0)hand, ifi(cid:1)=L1−1iYiXTi has some repeated but
and the solution is obtained by differentiating the cost in (7) (cid:0)P (cid:1)
stillnon-zerosingularvalues,thentheyaremappedtorepeated
with respect to the γ ’s and setting the derivative to 0. This
i (andnon-zero)singularvaluesin Γby(8). Letusassumethat
givesγ =0.5 σ σ2+2λ i. Since all the γ 0, the
i i± i ∀ i ≥ Σ has only one singular value that repeats (the proof easily
only feasible solution is
extends to the case of multiple repeated singular values) say
(cid:0) p (cid:1)
σ + σ2+2λ r times, and that these repeated values are arranged in the
γi = i 2i ∀i (8) bottom half of the matrix Σ (i.e., σn−r+1 = σn−r+2 = ... =
p σ =σˆ >0). Then, we have
Thus, a closed-form solution or global minimizer for the n
transform update step (4) is given as in (5). n−r
L−1YXT = σ Q RT +σˆ n Q RT (12)
The solution (5) is invariant to the specific choice of the i i i i=n−r+1 i i
matrix L. To show this, we will first show that if L1 Rn×n Xi=1 (cid:0)P (cid:1)
and L2 Rn×n satisfy YYT +λξI = L1LT1 = L2L∈T2, then Because the matrix defined by the first sum on the right (in
∈ (12)) corresponding to distinct singular values is unique, and
L2 = L1G, where G is an orthonormal matrix satisfying
GGT = I. A brief proof of the latter result is as follows. σˆ >0,sotooisthesecondmatrixdefinedbythesecondsum.
(Thisisalsoasimpleconsequenceofthefactthatalthoughthe
Since L1 and L2 are both n n full rank matrices (being
square roots of the positive de×finite matrix YYT +λξI), we singular vectors associated wtih repeated singular values are
have L2 =L1 L−11L2 =L1G, with G,L−11L2 a full rank notunique,thesubspacesspannedbythemare.)Thetransform
matrix. Moreov(cid:0)er, sinc(cid:1)e L1LT1 −L2LT2 =0, we have update solution (5) in this case is given as
L1 I−GGT LT1 =0 (9) Wˆ = in=−1rγiRiQTi +γn−r+1 ni=n−r+1RiQTi L(−113)
BecauseL1 hasfullran(cid:0)k,wemus(cid:1)tthereforehavethatGGT = Based o(cid:8)nPthe preceding argument(cid:0)s,Pit is clear that t(cid:1)h(cid:9)e right
I for (9) to hold. Therefore, G is an orthonormal matrix hand side of (13) is unique, irrespective of the particular
satisfying GGT =I. choice of (non-unique) Q and R. Thus, the transform update
Consider L1 and L2 as defined above. Now, if Q1 is the solution (5) is unique when L−1YXT has possibly repeated,
left singular matrix corresponding to L−11YXT, then Q2 = butnon-zerosingularvalues.Therefore,thenon-singularityof
GTQ1 isacorrespondingleftsingularmatrixforL−21YXT = L−1YXT is also a sufficient condition for the uniqueness of
GTL−11YXT. Therefore, replacing L1 by L2 in (5), making the solution to (4). (cid:4)
the substitutions L−1 = GTL−1, QT = QTG, and using the The transform update solution (5) is expressed in terms
2 1 2 1
orthonormalityofG,itisobviousthattheclosed-formsolution of the full SVD of L−1YXT, where L is for example, the
(5) involving L2 is identical to that involving L1. Cholesky factor, or alternatively, the Eigenvalue Decomposi-
Finally, we show that the solution (5) is unique if and tion (EVD) square root of YYT +λξI. Although in practice
only if L−1YXT is non-singular. First, the solution (5) can the SVD, or even the square root of non-negative scalars,
be written using the notations introduced above as Wˆ = are computed using iterative methods, we will assume in
n γ R QT L−1,whereR andQ aretheith columnsof the theoretical analysis in this paper, that the solution (5)
i=1 i i i i i
R and Q, respectively. We first show that the non-singularity is computed exactly. In practice, standard numerical methods
o(cid:0)Pf L−1YXT is(cid:1)a necessary condition for uniqueness of the are guaranteedto quicklyprovidemachineprecisionaccuracy
solution to (4). Now, if L−1YXT has rank < n, then a fortheSVD andothercomputations.Therefore,thetransform
singular vector pair Q ,R of L−1YXT corresponding to update solution (5) is computed to within machine precision
k k
a zero singular value (say σ = 0) can also be modified as accuracy in practice.
k
(cid:0) (cid:1)
5
Transform Learning Algorithms A1 and A2 and the solution is Xˆ = H (WY ) i. Moreover, the
i s i
Input: Y - training data, s - sparsity, λ - constant, ξ - transform update step for Problem (14) inv∀olves
constant, J0 - number of iterations.
Output: Wˆ - learned transform, Xˆ - learned sparse max tr WYXT s.t. WTW =I (15)
W
code matrix.
Initial Estimates: (W0,X0). DenotingthefullSVD(cid:0) ofYXT(cid:1)byUΣVT,whereU Rn×n,
Pre-Compute: L−1 = YYT +λξI −1/2. WΣˆ∈=RVn×UnT,.VTh∈isRsno×lunti,othneisoputnimiqaulesiofluantidonoinnlyPrifobthleem∈si(n1g5u)laisr
For k = 1: J0 Repeat(cid:0) (cid:1) values of YXT are non-zero.
1) Compute full SVD of L−1Y(Xˆk−1)T as QΣRT. For ξ = 0.5, Proposition 2 holds with the constraint
2) Wˆk =0.5R Σ+ Σ2+2λI 12 QTL−1. WTW =6 I in Problem (14) replaced by the constraint
3) XHˆˆik1(=WHkYs()W(cid:16)kiYfio)r(cid:0)∀AilgfoorrithAml(cid:1)gAor(cid:17)2it.hm A1, or Xˆik = W(14T)Wwi=th(t1h/e2ξm)oId.iTfiheedtrcaonnssftorarmintupWdaTteWsol=uti(o1n/f2oξr)PIroisbltehme
ηi i ∀ same as mentioned in Proposition 2, except for an additional
End
scaling of 1/√2ξ. The proof of Proposition 2 is provided in
Fig. 1. Algorithms A1 and A2 for solving Problems (P0) and (P1),
Appendix B.
respectively.Asuperscriptofkisusedtodenotetheiteratesinthealgorithms.
Although we begin with the transform update step in each iteration above, The orthonormaltransform case is special, in that Problem
onecouldalternatively startwiththesparsecodingstepaswell. (14) is also an orthonormal synthesis dictionary learning
problem, with WT denoting the synthesis dictionary. This
follows immediately, using the identity WY X =
F
ThealgorithmsA1andA2(correspondingtoProblems(P0) Y WTX , for orthonormal W. Henkce, Pr−oposiktion 2
F
and(P1),respectively)fortransformlearningareshowninFig. pkrov−ides an aklternating algorithm with optimal updates not
1. only for the orthonormal transform learning problem, but at
While Proposition 1 provides the closed-form solution to thesametimefortheorthonormaldictionarylearningproblem.
equation (4) for real-valued matrices, the solution can be
extended to the complex-valued case (useful in applications
C. Computational Cost
suchasmagneticresonanceimaging(MRI)[15])byreplacing
the ()T operation in Proposition 1 and its proof by ()H, the The proposed transform learning algorithms A1 and A2
Herm·itian transpose operation. The same proof appli·es, with alternate between the sparse coding and transform update
the trace maximization result for the real case replaced by steps. Each of these steps has a closed-form solution. We
max Re tr UΓVHQΣRH =tr(ΓΣ)forthecomplex now discuss their computational costs. We assume that the
U,V
case, where Re(A) denotes the real part of scalar A. matrices YYT +λI and L−1 (used in (5)) are pre-computed
(cid:8) (cid:0) (cid:1)(cid:9) (at the beginning of the algorithm) at total costs of O(Nn2)
and O(n3), respectively, for the entire algorithm.
B. The Orthonormal Transform Limit The computational cost of the sparse coding step in both
Algorithms A1 and A2 is dominated by the computation of
We have seen that for ξ = 0.5, as λ → ∞, the W theproductWY,andthereforescalesasO(Nn2).Incontrast,
minimizing (P0) tends to an orthonormal matrix. Here, we
theprojectionontothes-ℓ0ballinAlgorithmA1requiresonly
study the behavior of the actual sparse coding and trans-
O(nNlogn) operations,whenemployingsorting[5], andthe
form update steps of our algorithm as the parameter λ (or,
hardthresholding(asinequation(3))inAlgorithmA2requires
2
equivalently λ0, since λ = λ0kYkF) tends to infinity. The only O(nN) comparisons.
followingProposition2establishesthatasλ withξ held
→∞ For the transform update step, the computation of the
at 0.5, the sparse coding and transform update solutions for product YXT requires αNn2 multiply-add operations for an
(P0) approachthe correspondingsolutionsfor an orthonormal
X withs-sparsecolumns,ands=αn.Then,thecomputation
transform learning problem. Although we consider Problem
of L−1YXT, its SVD, and the closed-form transform update
(P0) here, a similar result also holds with respect to Problem
(5)requireO(n3)operations.Ontheotherhand,whenNLCG
(P1).
is employed for transform update, the cost (excluding the
Proposition 2: For ξ = 0.5, as λ , the sparse coding YXT pre-computation) scales as O(Jn3), where J is the
→ ∞
and transform update solutions in (P0) coincide with the number of NLCG iterations [5]. Thus, compared to NLCG,
corresponding solutions obtained by employing alternating the proposedupdate formula (5) allows for both an exact and
minimizationonthefollowingorthonormaltransformlearning potentially cheap (depending on J) solution to the transform
problem. update step.
Under the assumption that n N, the total cost per
min WY X 2 s.t. WTW =I, X s i (14) ≪
W,X k − kF k ik0 ≤ ∀ iteration (of sparse coding and transform update) of the
proposed algorithms scales as O(Nn2). This is much lower
Specifically, the sparse coding step for Problem (14) involves than the per-iteration cost of learning an n K overcomplete
×
(K > n) synthesis dictionary D using K-SVD [9], which
2
min WY X s.t. X s i
X k − kF k ik0 ≤ ∀ scales (assuming that the synthesis sparsity level s ∝ n)
6
as O(KNn2). Our transform learning schemes also hold R2. The local region defined by
a similar computational advantage over analysis dictionary ∆X <min φ (WY ): WY >s .
k k∞ i{ s i k ik0 }
learning schemes such as analysis K-SVD [5]. Furthermore, if we have WY s i, then ∆X can be
As illustrated in Section V-B, our algorithms converge in arbitrary. k ik0 ≤ ∀
fewiterationsinpractice.Therefore,theper-iterationcomputa-
The notation g∗ W0,X0 in Theorem 1 represents the
tionaladvantages(e.g.,overK-SVD)alsotypicallytranslateto
value to which the objective sequence g(Wk,Xk) con-
a net computationaladvantagein practice (e.g., in denoising).
verges, starting from(cid:0) an init(cid:1)ial (estimate) (W0,X0). Local
(cid:8) (cid:9)
region R2 in Theorem 1 is defined in terms of the scalar
IV. MAIN CONVERGENCE RESULTS
min φ (WY ): WY >s ,whichiscomputedbytaking
A. Result for Problem (P0) the cio{lusmns oif WkY wiitkh0spars}ity greater than s, and finding
Problem (P0) has the constraint X s i, which can thes-largestmagnitudeelementineachofthesecolumns,and
k ik0 ≤ ∀
instead be added as a penalty in the objective by using a choosing the smallest of those magnitudes. The intuition for
barrier function ψ(X) (which takes the value + when the this particular construction of the local region is provided in
∞
constraint is violated, and is zero otherwise). In this form, the proof of Lemma 9 in Appendix D.
Problem (P0) is unconstrained, and we denote its objective Theorem 1 indicates local convergence of our alternating
as g(W,X) = WY X 2 +λξ W 2 λlog detW + Algorithm A1. Assuming a particular initial (W0,X0), we
k − kF k kF − | |
ψ(X).Theunconstrainedminimizationprobleminvolvingthe have that every accumulation point (W,X) of the iterate
objective g(W,X) is exactly equivalent to the constrained sequence is a local optimum by equation (16), and satisfies
formulation (P0) in the sense that the minimum objective g(W,X) = g∗ W0,X0 . Thus, all accumulation points of
valuesaswellasthesetofminimizersofthetwoformulations the iterates (for a particular initial (W0,X0)) are equivalent
(cid:0) (cid:1)
are identical. To see this, note that whenever the constraint (in terms of their cost), or are equallygoodlocal minima.We
X s i is satisfied, the two objectives coincide. Other- thus have the following corollary.
k ik0 ≤ ∀
wise, the objective in the unconstrainedformulation takes the Corollary 1: For Algorithm A1, assuming a particular ini-
value+ andtherefore,itsminimumvalueisachievedwhere tial (W0,X0), the objective converges to a local minimum,
∞
the constraint X s i holds. This minimum value (and and the iterates converge to an equivalence class of local
k ik0 ≤ ∀
the corresponding set of minimizers) is therefore the same minimizers.
as that for the constrained formulation (P0). The proposed
The local optimality condition (16) holds for the algorithm
Algorithm A1 is an exact alternating algorithm for both the
irrespective of initialization. However, the local minimum
constrained and unconstrained formulations above. g∗ W0,X0 thattheobjectiveconvergestomaypossiblyde-
Problem(P0)istofindthebestpossibletransformmodelfor
pend on (i.e., vary with) initialization. Nonetheless, empirical
the given training data Y by minimizing the sparsification er- evi(cid:0)dence pre(cid:1)sented in Section V suggests that the proposed
ror,andcontrollingtheconditionnumber(avoidingtriviality).
transformlearning scheme is insensitive to initialization. This
We are interested to know whether the proposed alternating
leads us to conjecture that our algorithm could potentially
algorithm converges to a minimizer of (P0), or whether it
convergetotheglobalminimizer(s)ofthelearningproblemin
couldgetstuckinsaddlepoints,orsomenon-stationarypoints.
some(practical)scenarios.Fig.2providesasimpleillustration
Problem(P0)is non-convex,andthereforewell-knownresults
oftheconvergencebehaviorofouralgorithm.Wealsohavethe
on convergenceof alternatingminimization(e.g.,[24]) do not
followingcorollaryofTheorem1,where‘globallyconvergent’
applyhere.ThefollowingTheorem1providestheconvergence
refers to convergence from any initialization.
of our Algorithm A1 for (P0). We say that a sequence ak Corollary 2: Algorithm A1 is globally convergent to the
has an accumulation point a, if there is a subsequence that
set of local minimizers of the non-convex transform learning
(cid:8) (cid:9)
converges to a. For a vector h, we let φj(h) denote the objective g(W,X).
magnitude of the jth largest element (magnitude-wise) of h.
Note that our convergence result for the proposed non-
For some matrix B, B ,max B .
k k∞ i,j| ij| convex learning algorithm, is free of any extra conditions or
Theorem 1: Let Wk,Xk denote the iterate sequence
requirements. This is in clear distinction to algorithms such
generated by Algorithm A1 with training data Y and initial
as IHT [25], [26] that solve non-convexproblems,but require
(W0,X0). Then, th(cid:8)e objecti(cid:9)ve sequence g(Wk,Xk) is
extra stringent conditions (e.g., tight conditions on restricted
monotone decreasing, and converges to a finite value, say
isometry constants of certain matrices) for their convergence
g∗ =g∗ W0,X0 .Moreover,theiteratesequ(cid:8)enceisboun(cid:9)ded,
results to hold. Theorem 1 also holds for any choice of the
and each accumulation point (W,X) of the iterate sequence
(cid:0) (cid:1) parameter λ0 (or, equivalently λ) in (P0), that controls the
is a fixed point of the algorithm, and a local minimizer of
condition number.
the objective g in the following sense. For each accumulation
The optimality condition (16) in Theorem 1 holds true
point (W,X), ǫ′ =ǫ′(W)>0 such that
not only for local (small) perturbations in X, but also
∃
g(W +dW,X +∆X) g(W,X)=g∗ (16) for arbitrarily large perturbations of X in a half space.
≥ For a particular accumulation point (W,X), the condition
holds dW Rn×n satisfying dW ǫ′, and all ∆X tr (WY X)∆XT 0 in Theorem 1 defines a half-
Rn×N∀in the∈union of the followking rkeFgio≤ns. ∈ spaceofpe−rmissiblepertu≤rbationsinRn×N.Now,evenamong
(cid:8) (cid:9)
R1. The half-space tr (WY X)∆XT 0. the perturbations outside this half-space, i.e., ∆X satisfying
− ≤
(cid:8) (cid:9)
7
objectiveuinthefollowingsense.Foreachaccumulationpoint
(W,X), ǫ′ =ǫ′(W)>0 such that
∃
u(W +dW,X +∆X) u(W,X)=u∗ (17)
≥
holds dW Rn×n satisfying dW ǫ′, and all ∆X
Rn×N∀satisfy∈ing ∆X <mink ηk/F2 ≤. ∈
k k∞ i{ i }
The proofsof Theorems 1 and 2 are providedin Appendix
DandF.OwingtoTheorem2,resultsanalogoustoCorollaries
Fig.2. Possiblebehaviorofthealgorithmneartwohypotheticallocalminima
(marked with black dots) of the objective. The numbered iterate sequence 1 and 2 apply.
herehastwosubsequences(oneevennumbered,andoneoddnumbered)that Corollary 3: Corollaries 1 and 2 apply to Algorithm A2
converge to the two equally good (i.e., corresponding to the same value of
and the corresponding objective u(W,X) as well.
theobjective) localminima.
V. EXPERIMENTS
tr (WY X)∆XT >0 (and also outside the local region A. Framework
−
R2 in Theorem 1), we only need to be concerned about the In this section, we present results demonstrating the prop-
(cid:8) (cid:9)
perturbations that maintain the sparsity level, i.e., ∆X such erties of our proposed transform learning Algorithm A1 for
thatX+∆X hassparsity spercolumn.Foranyother∆X, (P0), and its usefulnessin applications.First, we illustrate the
≤
g(W+dW,X+∆X)=+ >g(W,X)trivially.Now,since convergence behavior of our alternating learning algorithm.
∞
X +∆X needs to have sparsity s per column, ∆X itself We consider various initializations for transform learning and
≤
canbe atmost2s sparsepercolumn.Therefore,thecondition investigate whether the proposed algorithm is sensitive to
tr (WY X)∆XT > 0 (corresponding to perturbations initializations. This study will provide some (limited) em-
−
that could violate (16)) essentially corresponds to a union of pirical understanding of local/global convergence behavior of
(cid:8) (cid:9)
lowdimensionalhalf-spaces(eachcorrespondingtoadifferent the algorithm. Then, we compare our proposed algorithm to
possible choice of supportof ∆X). In other words, the set of the NLCG-based transform learning algorithm [5] at various
“bad” perturbations is vanishingly small in Rn×N. patch sizes, in terms of image representation quality and
Note that Problem (P0) can be directly used for adaptive computationalcost of learning.Finally, we briefly discuss the
sparse representation (compression) of images [5], [27], in usefulness of the proposed scheme in image denoising.
whichcasetheconvergenceresultsherearedirectlyapplicable. All our implementations were coded in Matlab version
(P0) can also be used in applications such as blind denoising R2013a.All computationswere performedwith an Intel Core
[19], and blind compressedsensing [28]. The overallproblem i5 CPU at 2.5GHz and 4GB memory, employing a 64-bit
formulations [19], [28] in these applications are highly non- Windows 7 operating system.
convex (see Section V-D). However, the problems are solved Thedatainourexperimentsaregeneratedasthe2Dpatches
using alternating optimization [19], [28], and the transform ofnaturalimages.WeuseourtransformlearningProblem(P0)
learning Problem (P0) arises as a sub-problem. Therefore, by tolearnadaptivesparserepresentationsofsuchimagepatches.
using the proposed learning scheme, the transform learning The means of the patches are removed and we only sparsify
step of the alternating algorithms for denoising/compressed the mean-subtractedpatches which are stacked as columns of
sensing can be guaranteed (by Theorem 1) to converge 2. the training matrix Y (patches reshaped as vectors) in (P0).
The means are added back for image display. Mean removal
is typically adopted in image processing applications such as
B. Result for Penalized Problem (P1)
compression and denoising [19], [29]. Similar to prior work
When the sparsity constraints in (P0) are replaced with
[5], [19], the weight ξ =1 in all our experiments.
ℓ0 penalties in the objective with weights ηi2 (i.e., we solve We have previously introduced several metrics to evaluate
Problem(P1)),we obtainanunconstrainedtransformlearning
the quality of learnt transforms [5], [19]. The normalized
2 2
problemwithobjectiveu(W,X)= WY X +λξ W
c−oλnvloegrg|ednecteWgu|a+ranteNi=e1(sηii2mkiXlarikt0o.TIknhethoirsem−cas1e)ktfoFoor,AwlegkohraivthekmFa stkhpWearYcsiofil−cuamXtinoksn2FX/errko=Wr (YHNkS2F(EW,)wYfohr)eraeoftYrtahneissfomtrhametridWxatXaismddaeetnrfiionxte,edatnhades
P i s i
A2 that minimizes u(W,X). sparse codes [5]. The NSE measures the fraction of energy
Theorem 2: Let Wk,Xk denote the iterate sequence lost in sparse fitting in the transform domain, and is an
generated by Algorithm A2 with training data Y and initial interesting property to observe for the learnt transforms. A
(W0,X0). Then, th(cid:8)e objecti(cid:9)ve sequence u(Wk,Xk) is useful performance metric for learnt transforms in image
monotonedecreasing,andconvergestoafinitevalue,sayu∗ = representation is the recovery peak signal to noise ratio (or
u∗ W0,X0 . Moreover,the iterate sequence(cid:8)is bounded,(cid:9)and recovery PSNR), which was previously defined as 255√P/
each accumulation point (W,X) of the iterate sequence is Y W−1X in decibels (dB), where P is the number
a fi(cid:0)xed poin(cid:1)t of the algorithm, and a local minimizer of the of i−mage pixelFs and X is again the transform sparse code
(cid:13) (cid:13)
o(cid:13)f data Y [5(cid:13)]. The recovery PSNR measures the error in
2EvenwhendifferentcolumnsofX arerequiredtohavedifferentsparsity
recovering the patches Y (or equivalently the image, in the
levelsin(P0),ourlearningalgorithmandTheorem1canbetriviallymodified
toguarantee convergence. case of non-overlappingpatches) as W−1X from their sparse
8
codesX.TherecoveryPSNRservesasasimplesurrogatefor 10x 107 8x 107
the performance of the learnt transform in compression. Note
ittnhhfeaontrimfthtaehteiWopnr.ompoatsreidxatpoporowaocuhldwehraevteotboebuesetdrafnosrmciottmedpraessssiiodne, Objective Function 468 DKIRdLCaeTnnT dt IiInotnyimit tiIia naIlnilitiziziataaialtitizlioioaznnatiotionn Sparsification Error246 DKIRdLCaeTnnT dt IiInotnyimit tiIia naIlnilitiziziataaialtitizlioioaznnatiotionn
B. Convergence Behavior 120 0 101 102 100 0 101 102
Here, we study the convergence behavior of the proposed Iteration Number Iteration Number
(a) (b)
transform learning Algorithm A1. We extract the 8 8
(n = 64) non-overlapping (mean-subtracted) patches of×the 5
512 512 image Barbara [9]. Problem (P0) is solved to
learn×a square transform W that is adapted to this data. The mber4 DKLCTT IInniittiiaalliizzaattiioonn
data matrix Y in this case has N = 4096 training signals n Nu3 IRdaenndtiotym I nInitiiatializliazatiotionn
o
(patches represented as vectors). The parameters are set as diti
n
s = 11, λ0 = 3.1 10−3. The choice of λ0 here ensures Co2
×
well-conditioning of the learnt transform. Badly conditioned
transforms degrade performance in applications [5], [19]. 1 100 200 300 400 500 600
Iteration Number
Hence, we focus our investigation here only on the well- (c) (d)
conditioned scenario.
Fig. 3. Effect of different Initializations: (a) Objective function, (b)
We study the convergence behavior of Algorithm A1 for
Sparsification error, (c) Condition number, (d) Rows ofthe learnt transform
various initializations of W. Once W is initialized, the algo- shownaspatches forthecaseofDCTinitialization.
rithm iterates over the sparse coding and transform update
steps (this corresponds to a different ordering of the steps
in Fig. 1). We consider four different initializations (initial
The transforms learnt with different initializations in Fig. 3
transforms) for the algorithm. The first is the 64 64 2D
× alsoprovidesimilarrecoveryPSNRs(thatdifferbyhundredths
DCT matrix (obtained as the Kronecker productof two 8 8
× of a dB) for the Barbara image.
1D DCT matrices). Thesecondinitializationis the Karhunen-
Loe`ve Transform (KLT) (i.e., the inverse of PCA), obtained
here by inverting/transposing the left singular matrix of Y 3.
C. Image Representation
Thethirdandfourthinitializationsaretheidentitymatrix,and
a random matrix with i.i.d. Gaussian entries (zero mean and For the secondexperiment,we learn sparsifyingtransforms
standard deviation 0.2), respectively. fromthe√n √n(zeromean)non-overlappingpatchesofthe
×
Figure3showstheprogressofthealgorithmoveriterations image Barbara at various patch sizes n. We study the image
for the various initializations of W. The objective function representation performance of the proposed algorithm involv-
(Fig. 3(a)), sparsification error (Fig. 3(b)), and condition ing closed-form solutions for Problem (P0). We compare the
number(Fig.3(c)),allconvergequicklyforouralgorithm.The performance of our algorithm to the NLCG-based algorithm
sparsification error decreases over the iterations, as required. [5]thatsolvesaversion(withouttheabsolutevaluewithinthe
Importantly, the final values of the objective (similarly, the log-determinant)of (P0), and the fixed 2D DCT. The DCT is
sparsificationerror,andconditionnumber)arenearlyidentical a popular analytical transform that has been extensively used
for all the initializations. This indicates that our learning in compression standards such as JPEG. We set s=0.17 n
×
algorithm is reasonably robust, or insensitive to initialization. (roundedto nearestinteger),andλ0 is fixedtothe samevalue
Good initializations for W such as the DCT and KLT lead to asinSectionV-Bforsimplicity.TheNLCG-basedalgorithmis
fasterconvergenceoflearning.Thelearnttransformsalsohave executedwith 128NLCG iterationsforeachtransformupdate
identical Frobenius norms (5.14) for all the initializations. step, and a fixed step size of 10−8 [5].
Figure 3(d) shows the (well-conditioned) transform learnt Figure 4 plots the normalized sparsification error (Fig.
with the DCT initialization. Each row of the learnt W is 4(a)) and recovery PSNR (Fig. 4(b)) metrics for the learnt
displayed as an 8 8 patch, called the transform atom. The transforms, and for the patch-based 2D DCT, as a function
×
atoms here exhibit frequency and texture-like structures that of patch size. The runtimes of the various transform learning
sparsify the patches of Barbara. Similar to our prior work schemes (Fig. 4(c)) are also plotted.
[5], we observed that the transforms learnt with different The learnt transforms provide better sparsification and re-
initializations,althoughessentiallyequivalentinthesensethat covery than the analytical DCT at all patch sizes. The gap
they produce similar sparsification errors and are similarly in performance between the adapted transforms and the fixed
scaled and conditioned, appear somewhat different (i.e., they DCTalsoincreaseswithpatchsize(cf.[19]forasimilarresult
are not related by only row permutations and sign changes). and the reasoning). The learnt transforms in our experiments
areallwell-conditioned(conditionnumbers 1.2 1.6).Note
3We did not remove the means of the rows of Y here. However, we ≈ −
that the performancegap between the adapted transformsand
obtainalmostidenticalplotsinFig.3,whenthelearningalgorithmisinstead
initialized withtheKLTcomputedon(row)meancentered dataY. theDCTcanbeamplifiedfurtherateachpatchsize,byoptimal
9
malized Sparsification Error00000....0000.12468 CDlCoTsed Form Recovery PSNR (dB)3333334567 CICloAsed Form rWaIheClaspvhAForeiielnoemsoabelbwenslteystehaer,ovtriuivdenosedsnedcd[so3i(tmmct4hho]apie,mltaa[rFrpa3iarsl8pgeos]esotn.IrsrCfiitootoAhrnmm)sayaslonngftscohuereicimfsthhoisamrgadaseilisctnKetircFo-noSniamgVati.rvpDy4ea(rlb([eeb9)dau,]rtnwtpiosneelorgIfhw,CoawerAvmree).
Nor DCT
04 9 100 144 196 324 9 100 144 196 slightly better than the transform learning Algorithm A1 for
Patch Size (n) Patch Size (n)
the task of image representation. However, the learning and
(a) (b)
application of synthesis dictionaries also imposes a heavy
103 computationalburden(cf.[5]foracomparisonoftheruntimes
NLCG of synthesis K-SVD and NLCG-based transform learning).
me (sec)102 Closed Form Ifnodreaedco,amnpirmespsoiortnanatpapdlivcaantitoange(soifmoiularrtrtaoncsflaosrsmic-ablasaepdprsocahcehmees
Runti101 involving the DCT or Wavelets), is that the transform can be
applied as well as learnt very cheaply.
1040 9 100 144 196 While we adapted the transform to a specific image (i.e.,
Patch Size (n)
(c) image-specific transform) in Fig. 4, a transform adapted to a
varietyofimages(globaltransform)alsoperformswellintest
Fig. 4. Comparison of NLCG-based transform learning [5], Closed Form
images [27]. Both global and image-specific transforms may
transformlearningvia(P0),DCT,andICA[30]fordifferentpatchsizes:(a)
Normalizedsparsificationerror,(b)RecoveryPSNR,(c)Runtimeoftransform hold promise for compression.
learning. The plots for the NLCG and Closed Form methods overlap in (a)
and(b).Therefore,weonlyshowtheplotsfortheClosedFormmethodthere.
D. Image Denoising
Thegoalofdenoisingisto recoveranestimate ofanimage
x RP (2Dimagerepresentedasavector)fromitscorrupted
choice of λ0 (or, optimal choice of condition number 4). ∈
measurementy =x+h,wherehisthenoise.Weworkwithh
Theperformance(normalizedsparsificationerrorandrecov-
whose entriesare i.i.d.Gaussian with zero mean and variance
ery PSNR) of the NLCG-based algorithm [5] is identical to
σ2.Wehavepreviouslypresentedaformulation[19]forpatch-
that of the proposed Algorithm A1 for (P0) involving closed-
based image denoising using adaptive transforms as follows.
form solutions. However, the latter is much faster (by 2-11
times) than the NLCG-based algorithm. The actual speedups N
2 2
dsceapleensdwinithgernesepraelc,totnohNow/nJ. (thenumberofNLCGiterations) W,{xmi}in,{αi}Xi=1nkWxi−αik2+λiv(W)+τkRiy−xik2o
s.t. α s i (P2)
Inyetanothercomparison,we showinFig.4(b),therecov- k ik0 ≤ i ∀
ery PSNRs obtained by employing Independent Component Here, R Rn×P extracts the ith patch (N overlapping
i
Analysis (ICA – a method for blind source separation) [30]– ∈
patches assumed) of the image y as a vector R y. Vector
i
[34].SimilartopriorworkonICA-basedimagerepresentation x Rn denotes a denoised version of R y, and α Rn
i i i
[35], we learn an ICA model A (a basis here) using the ∈ ∈
is a sparse representation of x in a transform W, with an
i
FastICA algorithm[30], [36], to representthe training signals
apriori unknown sparsity s . The weight τ 1/σ [13], [19],
i
as Y =AZ, where the rows of Z correspond to independent ∝ 2
sources.NotethattheICAmodelenforcesdifferentproperties andλi issetbasedonthegivennoisydataRiy asλ0kRiyk2.
The net weighting on v(W) in (P2) is then λ= λ .
i i
(e.g., independence) than the transform model. Once the ICA
We have previously proposed a simple two-step iterative
modelislearnt(usingdefaultsettingsintheauthor’sMATLAB P
algorithm to solve (P2) [19], that also estimates the unknown
implementation [36]), the training signals are sparse coded in
s . Thealgorithmiteratesoveratransformlearningstep anda
i
the learnt ICA model A [35] using the orthogonal matching
variable sparsity update step (cf. [19] for a full description
pursuit algorithm [37], and the recovery PSNR (defined as
of these steps). We use the proposed alternating transform
in Section V-A, but with W−1X replaced by AZˆ, where Zˆ learningAlgorithmA1 (involvingclosed-formupdates) in the
is the sparse code in the ICA basis) is computed. We found
transform learning step. Once the denoised patches x are
i
that the A† obtained using the FastICA algorithm provides
found,the denoisedimagexis obtainedbyaveragingthex ’s
i
poornormalizedsparsificationerrors(i.e.,itisabadtransform
at their respective locations in the image [19].
model).Therefore,weonlyshowtherecoveryPSNRsforICA.
We now present brief results for our denoising framework
AsseeninFig.4(b),theproposedtransformlearningalgorithm
employing the proposed efficient closed-form solutions in
provides better recovery PSNRs than the ICA approach. This
transform learning. We work with the images Barbara, Cam-
illustrates the superiority of the transform model for sparse eraman, Couple 5, and Brain (same as the one in Fig. 1 of
[15]), and simulate i.i.d. Gaussian noise at 5 different noise
4The recovery PSNR depends on the trade-off between the sparsification levels (σ = 5, 10, 15, 20, 100) for each of the images.
errorandconditionnumber[5],[19].Fornaturalimages,therecoveryPSNR
We compare the denoising results and runtimes obtained by
using the learnt transform is typically better at λ values corresponding to
intermediate conditioning orwell-conditioning, ratherthanunitconditioning,
sinceunitconditioning istoorestrictive [5]. 5Thesethreewell-knownimageshavebeenusedinourpreviouswork[19].
10
Parameter Value Parameter Value Image σ NoisyPSNR K-SVD Transform
n 121 N′ 32000 5 34.15 38.09 38.28
λ0 0.031 τ 0.01/σ 10 28.14 34.42 34.55
C 1.04 s 12 Barbara 15 24.59 32.34 32.39
M′ 11 M 12 20 22.13 30.82 30.90
100 8.11 21.86 22.42
TABLEI
5 34.12 37.82 37.98
THEPARAMETERSETTINGSFOROURALGORITHM:n-NUMBEROF
10 28.14 33.72 33.87
PIXELSINAPATCH,λ0-WEIGHTIN(P2),C-SETSTHRESHOLDTHAT
Cameraman 15 24.60 31.50 31.65
DETERMINESSPARSITYLEVELSINTHEVARIABLESPARSITYUPDATESTEP
[19],M′-NUMBEROFITERATIONSOFTHETWO-STEPDENOISING 20 22.10 29.83 29.96
ALGORITHM[19],N′-TRAININGSIZEFORTHETRANSFORMLEARNING 100 8.14 21.75 22.01
STEP(THETRAININGPATCHESARECHOSENUNIFORMLYATRANDOM 5 34.14 42.14 42.74
FROMALLPATCHESINEACHDENOISINGITERATION)[19],M -NUMBER 10 28.12 38.54 38.78
OFITERATIONSINTRANSFORMLEARNINGSTEP,τ -WEIGHTIN(P2),s- Brain 15 24.62 36.27 36.43
INITIALSPARSITYLEVELFORPATCHES[19]. 20 22.09 34.70 34.71
100 8.13 24.73 24.83
5 34.16 37.29 37.35
10 28.11 33.48 33.67
our proposed algorithm with those obtained by the adaptive Couple 15 24.59 31.44 31.60
overcomplete synthesis K-SVD denoising scheme [13]. The 20 22.11 30.01 30.17
100 8.13 22.58 22.60
Matlab implementation of K-SVD denoising [13] available
from Michael Elad’s website [29] was used in our compar- TABLEII
PSNRVALUESINDECIBELSFORDENOISINGWITHADAPTIVE
isons, and we used the built-in parameter settings of that
TRANSFORMS,ALONGWITHTHECORRESPONDINGVALUESFOR64×256
implementation. OVERCOMPLETEK-SVD[13].THEPSNRVALUESOFTHENOISYIMAGES
We use 11 11 maximally overlapping image patches for (DENOTEDASNOISYPSNR)AREALSOSHOWN.
×
our transform-based scheme. The resulting 121 121 square
×
transform 6 has about the same number of free parameters
as the 64 256 overcomplete K-SVD dictionary [13], [29]. σ 5 10 15 20 100
The setting×s for the variousparameters(notoptimized)in our Average Speedup 9.82 8.26 4.94 3.45 2.16
transform-based denoising scheme are listed in Table I. At TABLEIII
σ = 100, we set the number of iterations of the two-step THEDENOISINGSPEEDUPSPROVIDEDBYOURTRANSFORM-BASED
SCHEME(INVOLVINGCLOSED-FORMSOLUTIONS)OVERK-SVD[13].
denoising algorithm [19] to M′ = 5 (lower than the value in THESPEEDUPSAREAVERAGEDOVERTHEFOURIMAGESATEACHNOISE
Table I), which also workswell, and providesslightly smaller LEVEL.
runtimes in denoising.
Table II lists the denoising PSNRs obtained by our
transform-based scheme, along with the PSNRs obtained by
K-SVD. The transform-based scheme provides better PSNRs cheaper than the K-SVD method.
than K-SVD for all the images and noise levels considered. Weobservethatthespeedupofthetransform-basedscheme
The average PSNR improvement (averaged over all rows of over K-SVD denoising decreases as σ increases in Table
Table II) provided by the transform-based scheme over K- III. This is mainly because the computational cost of the
SVD is 0.18 dB. When the NLCG-based transform learning transform-based scheme is dominated by matrix-vector mul-
[5] is used in our denoising algorithm, the denoising PSNRs tiplications (see [19] and Section III-C), and is invariant to
obtainedareverysimilarto theonesshownin TableIIforthe the sparsity level s. On the other hand, the cost of the
algorithm involving closed-form updates. However, the latter K-SVD denoising method is dominated by synthesis sparse
scheme is faster. coding,whichbecomescheaperasthesparsityleveldecreases.
We also show the average speedups provided by our Since sparsity levelsin K-SVD denoisingare set accordingto
transform-based denoising scheme 7 over K-SVD denoising an error threshold criterion (and the error threshold σ2)
∝
in Table III. For each image and noise level, the ratio of [13], [29], they decrease with increasing noise in the K-
the runtimes of K-SVD denoising and transform denoising SVD scheme. For these reasons, the speedupof the transform
(involving closed-form updates) is first computed, and these method over K-SVD is lower at higher noise levels in Table
speedups are averaged over the four images at each noise III.
level.Thetransform-basedschemeisabout10xfasterthanK- We would like to point out that the actual value of the
SVD denoising at lower noise levels. Even at very high noise speedup over K-SVD also depends on the patch size used
(σ =100),thetransform-basedschemeisstillcomputationally (by each method). For example, for larger images, a larger
patch size would be used to capture image informationbetter.
6Wehavepreviouslyshownreasonabledenoisingperformanceforadapted The sparsity level in the synthesis model typically scales as
(using NLCG-based transform learning [5]) 64×64 transforms [19]. The
a fraction of the patch size (i.e., s n). Therefore, the
denoisingperformanceusuallyimproveswhenthetransformsizeisincreased,
∝
butwithsomedegradation inruntime. actual speedup of transform-based denoising over K-SVD at
7Our MATLAB implementation is not currently optimized for efficiency. a particular noise level would increase with increasing patch
Therefore, the speedups here are computed by comparing our unoptimized
(andimage)size–aneffectthatisnotfullyexploredheredue
MATLABimplementation(fortransform-baseddenoising)tothecorrespond-
ingMATLABimplementation [29]ofK-SVDdenoising. to limitations of space.