Frame-based Sparse Analysis and Synthesis Signal Representations and Parseval K-SVD

Frames are the foundation of the linear operators used in the decomposition and reconstruction of signals, such as the discrete Fourier transform, Gabor, wavelets, and curvelet transforms. The emergence of sparse representation models has shifted of …

Authors: Wen-Liang Hwang, Ping-Tzan Huang, Tai-Lang Jong

Frame-based Sparse Analysis and Synthesis Signal Representations and   Parseval K-SVD
1 Frame-based Sparse Analysis and Synthesis Signal Representations and P arse v al K-SVD W en-Liang Hwang, Ping-Tzan Huang, and T ai-Lang Jong Abstract Frames are the f ounda tio n of the linear o p erators used in th e decomp o sition and recon - struction of signals, such as the discrete F ourier transfo rm, Gabor , wa velets, and cur velet transform s. Th e emergence o f sparse rep resentation models has shifted of the emphasis in f rame theory tow ard sparse l 1 -minimization problems. In this pap er , we a pply f rame theory to the sparse re presentation of signals in which a synthesis dictionar y is used fo r a fr ame and an an alysis dic tio nary is used fo r a dual f rame. W e sought to for mulate a n ovel dual frame design in which the sparse vector o btained thr o ugh the d e composition of any signal is also th e sparse solution re presenting signals based on a rec onstruction frame. Our findings demon strate th at this type of dual frame cann ot be constructed for over -comp le te frames, th ereby pre cluding the use o f any linear ana lysis op e rator in dri ving the sparse synthesis coefficient fo r signal represen tatio n. No netheless, the best approxim a tio n to the sparse synthesis solution can b e derived from the analysis coefficient using the canon ical dual frame. In this study , we developed a novel dictionary lea r ning algorithm ( called Parsev al K-SVD) to learn a tight-fram e d ictionary . W e then leveraged the ana ly sis and synthesis perspectives of signal representatio n with frames to d eriv e optimizatio n for mulations for problem s pertaining to im age recovery . Our p reliminary , r esults dem onstrate that the images recovered using this ap proach are corr elated to the frame bou nds of d ictionaries, thereby demonstra tin g the impo rtance of u sing different d ictionaries for different ap plications. W en-Liang Hwang is at Institute of Information Science, Academia Sinica, T aipei 115 29, T aiwan. Ping-Tzan Huang and T ai-L ang Jong are at Department of Electrical Engineering, National Tsi ng Hua Univ ersity , Hsinchu , T aiwan. DRAFT 2 I . I N T R O D U C T I O N Signal representation is a fundamental aspect of signal processi ng, serving as the basis of sig n al decomposition (analysis), processing, and rec onstruction (s y nthesis) [1]. A s ignal is first mapped as a vector of coefficients in the transform d o main using a decomposition operator (analysis operator). The coef ficients are then modified and mapped to the sign al space using a re construction operator (synt hesis operator). Th e fundamental prerequisite i n the design of analys i s and synthesis operators i s the perfect reconstruction of the signal [2], [3]. Frame theory refers to the branch of mathematics tasked with the design of linear decom p osition and reconstruction operators capable of perfectly reconstructi n g any signal i n a vector s pace [4]–[7]. A frame comprises of a set of linearly in dependent vectors spanning a vector space. Any signal in the vector space can be repre sented as a linear combinatio n of elements i n the frame. Analysis coef ficients (also referred t o as decomposit ion coef ficients or frame coef ficients) are deriv ed by applying the inp u t produ ct between the signal and elements in a dual frame. A dual frame ass ociated wit h the frame us ed for t his type of signal representation is not unique. Milestones i n the de velopment of fr ame theory include the const ruction of the canoni cal dual frame via the frame operator , the constructions of discrete Fourier transform, Gabor , wav elets, and curvelet transforms [6], [8]–[10], and the exploration of frame redundancy in various signal processing app lications [3], [11], [12]. In cases where the signal space is a finite-dimensi o nal discrete vector space, th e frame and canonical dual frame h ave a connecti on t o pseudo-in verse, singular value decompo sition in matrix linear algebra. The emergence of the sparse representatio n m odel has shi fted the emphasis i n frame theory toward sparse l 1 -minimizatio n probl em s [13]–[15], as evidenced by t h e plethora of algorithms aim ed at eluci d ating the model and making practical use of it. Suf ficient conditions for a unique k -sparse solution to the sparse m odel can be derived through the analysis of mutual i ncoherence [15 ], [16], the null space property [17], the restricted isometry property (RIP) [18], and spark [19], ( ψ ) , which is defined as th e sm allest possible number k such that th ere exists a subset o f k columns of ψ that are linearly dependent. Throughout the paper , we assume that spark ( ψ ) > 2 k (1) to ensure t hat the solut i on to the l 1 -synthesis problem i s uniqu e. W e also assume that DRAFT 3 ψ (denoting a fr ame) and φ (denoting a dual frame) are bot h in R m × n with m ≥ n and the ranks of ψ and φ are n . In t h is paper , we introd u ce a novel method by which to establ i sh a connection between sparse representation problems and frame theory . If we regard the synthesis dictionary for sparse representation as a frame, then the analysis di ctionary can be regarded as its dual frame. This novel perspective leads to the following conclusions: (1) It is impossib l e to construct a dual frame for any over -complete frame with the aim of obtaining t he minimizer of the ℓ 1 -synthesis-based prob lem; and, (2) The canonical dual frame is best linear decomposition operator by whi ch to obtain t he approximation o f the ℓ 1 -minimizer . This provides t heoretical sup p ort for the trend tow ard the use of non- linear d ecom position operators to derive th e ℓ 1 -minimizer in a straightforward manner; i.e., without relying on iterati ve algo rithms [20], [21]. The relationshi p between sparse representatio n problems and frame theory shi fts our perspectiv e with regard to dictionary l earning p rob lems. T o the best o f our knowledge, no existin g d i ctionary (or frame) learning algorithm has addressed the issue of frame bounds. Frame bounds correspond to the lar gest and smallest eigenv al u es of the frame operator ψ ψ ⊤ of frame ψ . The ratio of t he bou nds (i.e., the condition num ber) determines the degree o f redundancy between frame coef ficients [22], which means that it i s fundamentally ass ociated with numerical stability and the performance of various signal processing problems. In this paper , we demonstrate that the properti es of frame theory can be appli ed to the development of a learning algorithm to learn a Parse v al d i ctionary from a s et o f observations. Th i s is a frame who se numerical properties are cl o sest to an orth o normal matrix, as the ratio of the frame bounds is 1 . Finally , we describe the preliminary application of the resulting learned dictio nary to probl ems associated with the rest o rati on of images (denoising , image compression, and filli ng-the-missing - pixels). Our objective is to illu s trate ho w the frame bound s of dictionaries aff ects the performance of optimization problems that leverage the synthesis and analysis perspectiv es of fr ame coeffi cients of images. The remainder of th e paper i s or ganized as follows. Section II presents a revie w of studies on frame t heory as well as analysis and synthesis sparsity . Our main theoret- ical cont ri butions are presented in Section III and Section IV. Section V outlines our approach to training a P arse val ti ght frame from observations. Section VI presents our experiment results. Concludi n g remarks are drawn in Section VII. DRAFT 4 I I . B AC K G RO U N D A N D R E L A T E D W O R K S A. F rames A frame comprises a set of linearly independent vectors in a vector space spanned by the vectors. Frames are the cornerstone o f signal p rocessing in the formulat i on of perfect reconstruction pairs of linear o p erators used to decompose signals into the transform domain and reconstruct the original signals from the t ransform coef ficients. Frame studies d esi gning a set of vectors ψ i ∈ V span the vector space by representing any f ∈ V as f = m X i =1 h f , ˜ c i i ψ i , (2) where { ˜ c i } and {h f , ˜ c i i} are called the dual frame and the frame coef ficients of ψ , respectiv ely [4], [6]. If frame ψ is over -complete, then there are different choices for ˜ g i 6 = ˜ c i in which f = m X i =1 h f , ˜ g i i ψ i . (3) The frame operator S of frame f i with i = 1 , · · · , m is defined as follo ws: S f = ψ ⊤ ψ f = m X i h f , ψ i i ψ i , (4) where ψ ⊤ maps signal f in V to coeffic ients in C m . Ensuring that S has an in verse requires that the frame condition is met for an y f ∈ V : A k f k 2 ≤ h f , S f i = h f , ψ ψ ⊤ f i = h ψ ⊤ f , ψ ⊤ f i ≤ B k f k 2 , (5) where ∞ > B ≥ A > 0 are frame boun ds. If A = B , then t h e frame is tight; and if A = B = 1 , then the frame is Parse v al. The canonical dual frame i s defined as sequence S − 1 f 1 , · · · , S − 1 f m , where S − 1 is the i nv erse of t h e frame operator . The canonical dual frame s at i sfies the frame condition for any f ∈ V with bounds A − 1 and B − 1 : B − 1 k f k 2 ≤ h f , S − 1 f i ≤ A − 1 k f k 2 . (6) A mil est one of the frame th eory indicates that t h e frame coefficients deri ved from the canonical dual frame is the solution of ℓ 2 -synthesis problem    min u k u k 2 f = P m i =1 u i ψ i for all f ∈ V , (7) with u = [ u 1 · · · u m ] ⊤ . This sequence of coef ficients has the smallest l 2 -norm of all frame coef ficients of ψ representing any signal f ∈ V . DRAFT 5 B. Analysis and S ynt hesis Sparse Models The sy nthesis-based sp arse model comprises a dictionary ψ o f size n × m with m ≥ n in w h i ch it is assumed that a s i gnal i s a l i near com bination o f fe wer than k atom s in the dictionary ( x = ψ u and k u k 0 ≤ k ). Signals of interest lie within a union o f k -subspaces of the space sp ann ed by all atoms in the di ctionary . The para llel analysis-based sparse model (also called the co-sparse analysis model) is used to design an analysis operator φ (a n × m m atrix with m ≥ n ) in which the analysis vector φ ⊤ x of signal x is sparse and the signal with the sparsest analysis coef ficients i s reco vered. The paper by El ad et al. [23] provided deep insight i n to the use of analysis and syn- thesis models as priors for Bayesian methods. They poi n ted out that the two m odels are equiv alent as long as ψ is square and i nv ertible. They also p rovided examples showing the dichotom y between the two models in the case of over -complete dictionaries, where m > n . Finally , they presented theoretical resul ts i ndicating that any sparse analysis- based probl em po s es an equivalent sparse synthes i s-based problem, but the re verse does not hold. They als o demonstrated that th e reformulati o n o f an analysis problem to an identical syn t hesis problem can lead t o an e xponentially large diction ary . Nam et al. [24] demo nstrated that, for over -comp l ete ψ and φ , there is generally a cons iderable diffe rence in the uni o n of subspaces provided by sy nthesis and analysis models . If ψ and φ both hav e sizes of n × m , then the number of atoms of ψ th at can syn thesize a signal is freely from 0 to n − 1 while that of φ that obtains zero coef ficients of φ ⊤ x (co-sparsity) freely and from 0 to n − 1 . The number of zero coefficients of φ ⊤ x is in versely proportio n al to the subsp ace of ψ in which si gnal x = ψ u lies on th e number of n o n-zero coefficients in u . Thus , the algorithm of t h e co-sparsit y analysis model is designed to deri ve zero coefficients whereas the algo ri thm of the synthesis-based s p arse model is d esi gned to deriv e no n -zero coefficients. The above studies demonstrate that analysis and sy nthesis models can be vie wed as complem entary . This , i n turn, sugg ests that analysi s and synthesis operators and the corresponding recovery algorithms could perhaps be designed in pairs. I I I . S P A R S E S Y N T H E S I S C O E FFI C I E N T S V I A D U A L F R A M E The dual frame of ψ that min i mizes the l 2 -norm frame coefficient of any si gnal i s its canonical du al frame. W e sought t o d et erm i ne whether th ere exists a uni versal dual DRAFT 6 frame of ψ that min imizes the l 1 -norm frame coefficient for any signal. The existence of l 1 -norm frame coef ficients for any ψ i s provided by the following proposition. Pr oposition 1 [7]. Let ψ = [ ψ i ∈ C n ] i be a f rame for finit e-dimension vector sp a ce V . Given x ∈ V , ther e exists coef ficients d i ∈ C m such that x = ψ d = m X i =1 d i ψ i , (8) and k d k 1 = inf { d i | x = m X i =1 d i ψ i } , (9) wher e d = [ d 1 , · · · , d m ] ⊤ is the vector of sparse synt hesis coef ficients. If ψ is under-c omplete or s q uare with n ≥ m , then the answer to our question is af firmativ e and th e canonical dual frame is the solution [23]. Unfortunately , Theorem 3 shows that this is not the case if ψ is over -complete. First, we highlight an in teresting result o f Nam et al. [25], [26]: if φ is a frame with the form of an n × m m atrix and with colum n s φ i in general positio n s (any s et of at most n columns are linearly independent), then the number k of non-zero coeffi cients of φ ⊤ x , where x is a sign al in R n , is at least m − n + 1 . Since m − k rows i n φ ⊤ x are zeros and these ro ws take up no more than n − 1 dimensions, k ≥ m − n + 1 . In the case of m > n , th e result excludes the case i n whi ch φ ⊤ x is sparse for any k . Theorem 3 extends this result by showing that th ere does not exist a dual frame φ of over -comp lete ψ for an y pair of x ∈ R n and u ∈ R m with k ≥ m + n − 1 that s atisfies    x = ψ u ; u = φ ⊤ x. (10) Lemma 2 . Let the n × m matrix φ = [ φ 1 · · · φ m ] be a frame in R n with columns φ i in g eneral posi tions, an d φ is a dual frame of ψ (th u s , ψ φ ⊤ = I n × n ). Suppose that k ≥ m − n + 1 . Then, for m > n , ther e exists no pair of φ and ψ that meets (10) for any pair of k -sparse vector u and signal x . Pr oof. S uppose t hat there is a pair of frame ψ and dual frame φ s atisfying conditi ons (10) for all pairs o f signal x and correspond i ng k -sparse vector u . W i t h no los s of generality , we supp o se that the indices of no n -zero coefficients of u are 1 , · · · , k . Because x = ψ u , x is in the subsp ace spanned by ψ 1 , · · · , ψ k (denoted as V ψ a ( k ) ). Meanwhile, because u = φ ⊤ x , x is in the subspace spanned by φ 1 , · · · , φ k (denoted as V φ a ( k ) ), as well as DRAFT 7 in the subspace perpendicular to that spanned by the remaining m − k vectors (denoted as V φ b ( m − k ) ⊥ ). Hence, x ∈ V ψ a ( k ) ∩ V φ a ( k ) ∩ V φ b ( m − k ) ⊥ . (11) Nam et al. [25], [26] reported t hat to ha ve a non -empty intersection of V φ a ( k ) and V φ b ( m − k ) ⊥ , k must be at least m − n + 1 . Because t he set of signals satis fying Equation (11) is a vector sp ace of dimension k ( x = ψ u and the bottom m − k elements in u are zeros.), the dimension of V φ b ( m − k ) is n − k . Let φ = [ φ a φ b ] and ψ = [ ψ a ψ b ] , where φ a = [ φ 1 · · · φ k ] , φ b = [ φ k +1 · · · φ m ] , ψ a = [ ψ 1 · · · ψ k ] , and ψ b = [ ψ k +1 · · · ψ m ] . Equation (10) stipulates that we ha ve u = φ ⊤ ψ u. (12) Therefore, φ ⊤ b ψ a u 1 = 0 m − k and φ ⊤ a ψ a u 1 = u 1 , where u 1 is the first k elements in u . Since u 1 ∈ R k is arbitrary , we obtain the follo wing: φ ⊤ b ψ a = 0 m − k × k ; (13) φ ⊤ a ψ a = I k × k ; (14) which corresponds to V φ b ( m − k ) ⊥ V ψ a ( k ) and V φ a ( k ) = V ψ a ( k ) , respectiv ely . Combining the resul t s of Equations (13) and (14) whi l e consid erin g th at the dimensions of V φ b ( m − k ) and V ψ a ( k ) are respectiv ely n − k and k , we obtain the following: V ψ a ( k ) = V φ a ( k ) = V φ b ( m − k ) ⊥ . ( 15) Based on the fa ct that the dimensi on spanned by m − k columns in φ b is n − k , we know th at m − k = n − k + ( m − n ) . Because m − n > 0 , we can deduce th at there exists a set of n − k + 1 vectors in φ b spanning a dimensi o n n − k s ubspace. This violates the assumption that φ i are in general positions. End of Pr oof. Theor em 3 . Let ψ be an over -complete f rame and let u ∗ 1 ( x ) be the k -sparse minimi zer of l 1 -synthesis pr oblem for sign a l x with    min u k u k 1 x = ψ u ; (16) for signal x . Ther e does not e xist a dual frame φ 1 in general positions that yields u ∗ 1 ( x ) = φ ⊤ 1 x for any x , r e gar dless of the value of k . DRAFT 8 Pr oof. Th results reported by Nam et al. [25], [26] forfeit the existence of φ for k = k u ∗ 1 k 0 ≤ m − n . Meanwhile, Lemma 2 forfeits the existence of φ for k ≥ m − n + 1 . End of Pr oof. Proposition 1 claim s the existence of l 1 -norm mini mizer for any x . Theorem 3 shows that the minimizer cannot be deri ved t h rough the decomposition o f signals with a univ ersal linear operator , which depends exclusi vely on ψ . I V . O P T I M A L P R OX Y T O S P A R S E S Y N T H E S I S C O E FFI C I E N T S Theorem 3 provides a negative answer concerning th e existence of a du al frame φ 1 of over -comp lete frames as th e analysis operator to obtain the l 1 -minimizers of any sig n al. Ne vertheless, this section presents an af firmati ve answer concerning the e xistence of a univ ersal dual frame for t he follo wing approximation problem:    min φ k φ ⊤ x − u ∗ 1 ( x ) k 2 φ is a dual frame of ψ. (17) The foll owing prop o sition shows that the canonical dual frame, φ 2 , is the solution to (17) resulting in frame coeffic ients closest to the l 1 -minimizer of any signal and can be deriv ed by solving th e foll owing ℓ 2 -analytical problem:    min φ k φ ⊤ x k 2 φ is a dual frame of ψ. (18) Theor em 4 . (i) The canonical du al frame yields the op t imal proxy of u ∗ 1 ( x ) , the solution of l 1 -synthesis pr oblem (16), for any signal x . (ii) φ 2 is also the mini mizer of the l 2 -analysis pr oblem (18). Pr oof. (i) If ψ is an und er-complete or squ are fra me, then φ ⊤ 2 x = u ∗ 1 ( x ) is applicalbe t o any x [23]. The canonical dual frame, φ 2 , of th e over -complete frame ψ is ( ψ ψ ⊤ ) − 1 ψ . Since φ ⊤ 2 ψ = ψ ⊤ ( ψ ψ ⊤ ) − 1 ψ = ψ ⊤ φ 2 , the kernel φ ⊤ 2 ψ is th e rank n ort h ogonal projection 1 from R m to R m , due to the f act that    ( φ ⊤ 2 ψ ) 2 = φ ⊤ 2 ( ψ φ ⊤ 2 ) ψ = φ ⊤ 2 ψ ; ( φ ⊤ 2 ψ ) ⊤ = ψ ⊤ φ 2 = φ ⊤ 2 ψ . (19) 1 P i s an orthogonal projection i f and only if P 2 = P and P ⊤ = P . DRAFT 9 Consequently , the optimal v alue of problem (17) for any signal x is k φ ⊤ 2 x − u ∗ 1 ( x ) k 2 = k φ ⊤ 2 ψ u ∗ 1 − u ∗ 1 ( x ) k 2 = k ( I m × m − φ ⊤ 2 ψ ) u ∗ 1 ( x ) k 2 . (20) If φ differs from φ 2 , then kernel φ ⊤ ψ is a projection operator but is not necessarily an orthogonal proj ecti on as ( φ ⊤ ψ ) 2 = φ ⊤ ( ψ φ ⊤ ) ψ = ψ ⊤ ψ . Thus, the canonical dual frame yields the optimal proxy of u ∗ 1 ( x ) . (ii) The case where ψ is a square or under-complete frame: Equation x = ψ u ∗ 2 , wi th u ∗ 2 being the minimi zer of ℓ 2 -synthesis problem (7), implies that ( ψ ψ ⊤ ) − 1 ψ ⊤ x = u ∗ 2 . Since the canonical dual frame φ ⊤ 2 is ( ψ ψ ⊤ ) − 1 ψ ⊤ , we ha ve φ ⊤ 2 x = u ∗ 2 . In the case where ψ is an ove r-complete frame: u ∗ 2 is also the solution of    min u 1 2 k u k 2 2 x = ψ u . (21) T aking partial deriv ativ es with respect to λ and u on the Lagrangian function g ives L ( u, λ ) = 1 2 k u k 2 2 + λ ⊤ ( ψ u − x ) (22) and then setting the results to zero, we obtain s olution ( u ∗ 2 , λ ∗ ) that satisfy , ψ u ∗ 2 = x and u ∗ 2 = ψ ⊤ λ ∗ ; (23) respectiv ely . It foll ows that x = ψ ψ ⊤ λ ∗ . ( ψ ψ ⊤ ) − 1 x = λ ∗ and φ ⊤ 2 x = u ∗ 2 can be deduced from the fac t that ψ is an ov er -complete frame and φ ⊤ 2 = ψ ⊤ ( ψ ψ ⊤ ) − 1 . End of Pr oof. Thus, φ ⊤ 2 x is referred to as the optimal proxy of u ∗ 1 ( x ) . Th e following corollary is summarized from Theorems 3 and 4. Corollary 5 . (i) Let ψ be an over -complete frame; φ be i ts dual frame; and let x be a sig nal. Let S ( x ) = { u | x = ψ u } be the synthesis coefficients of x and let A ( x ) = { φ ⊤ x | φ is a dual frame of ψ } b e the frame coef ficients of the sig n al. Then, A ( x ) ⊆ S ( x ) , u ∗ 1 ( x ) ∈ S ( x ) . The orthogonal pr ojection of u ∗ 1 ( x ) to A ( x ) obt a ins the frame coef ficient u ∗ 2 ( x ) , which is the solution of ℓ 2 -synthesis pr oblem (7) (see F igur e 1). (ii) u ∗ 2 can be obtained linearly by applying φ ⊤ 2 to x (see F igure 2). The fact that u ∗ 1 cannot be ob t ained wi th a li near operator impl i es that s ome non-linear method must be adopted in order to obtain u ∗ 1 analytically . This conclusio n provides theoretical support for the recent trend of adopting non-linear operators to deri ve the solutions to sparse representation and compressed sensing problems. DRAFT 10 Fig. 1. Signal x and frame ψ are fixed. S ( x ) is a linear sub-space. A ( x ) is a con vex set due to the fact that if φ a and φ b are dual frames of ψ , so does αφ a + (1 − α ) φ b for any α ∈ [0 , 1] . I f u ∗ 1 does not belong to A ( x ) , then u ∗ 2 is it s best approximation in A ( x ) ; otherwise, u ∗ 1 = u ∗ 2 , obtainable by a linear operation. Fig. 2. Frame ψ is used as the reco nstruction (synthesis) operator . The decomposition (analysis) operator is the transpose of canonical dual f r ame φ 2 . T he optimal solution of ℓ 2 -sythesis poblem, u ∗ 2 , can be obtained analytically with li near operator φ 2 . V . L E A R N I N G P A R S E V A L F R A M E S F O R S PA R S E R E P R E S E N TA T I O N Aharon et al. [27] address ed the problem o f learning a s y nthesis di ctionary to sparsely represent a set of obs ervations. The K-SVD algorit hm is perhaps the most popul ar algorithm of that kind. Let Y = [ y 1 , · · · , y N ] be a n × N matrix, where each colu m n is an observ ation. The K-SVD algorithm optimizes          min X,ψ k Y − ψ X k 2 F X = [ x 1 , · · · , x N ] k x i k 0 ≤ k for i = 1 , · · · N . (24) Rubinstein et al. [28] addressed a parallel problem, wherein they concentrated on learning an analysis dictionary to produce a sparse outcome from a set of observations. The synthesis and analys i s dictionaries are learned using a s i milar approach; howe ver , they are mu tually exclusi ve. Dictionary learning problems can be brought wi thin th e context of frame theory by controlling th e frame bounds of the learned dictionary . Frame bounds A and B i n (5) are related t o the correlation between basis elem ent s in frames. Higher B A values increase the likelihood of correlation between basis elements in frames. If B A is far from 1 , then un iform distortio n i n a s i gnal yields non-uni form, direction-dependent DRAFT 11 distortion i n i ts frame coefficients. The K-SVD algorit hm is able to learn a fra me ψ k svd ; howe ver , it has no control over the bounds, due to the fact th at it i m poses no constraints o n frame bound s (the l argest and smallest si n gular values of ψ k svd ψ ⊤ k svd ). Parse val tight frames with frame bounds of 1 hav e been widely used in signal processing [29]–[31]. T h us, we re-formulated the K -SVD optimizatio n problem in our deve lopment of a learning algo ri t hm to learn the opt imal Parse v al tight frame as w ell as t he sparse coef ficients from a set of ob serv ations. Note that ψ is a Parse val tight frame if and only if |h ψ ⊤ x, ψ ⊤ x i| 2 = k x k 2 . Thus, ψ ψ ⊤ = I n × n . A. Design of P arse val F rames By imposing a Parse v al frame ψ using a K-SVD-like approach, we obtain the fol- lowing:                min X,ψ k Y − ψ X k 2 F X = [ x 1 , · · · , x N ] k x i k 0 ≤ k for i = 1 , · · · , N ψ is a Parse val tig ht fra me and ψ ψ ⊤ = I n × n . (25) This probl em is difficult to solve because applying the augmented Lagrangian ap- proach would introduce fourth-order polynom ial terms on ψ , due t o the tight frame constraint ψ ψ ⊤ = I n × n . Thus, we take advantage o f the f act that the canonical dual frame of a Parse v al tigh t frame is itself ( φ 2 = [( ψ ψ ⊤ ) − 1 ] ⊤ ψ = ψ ) and form ulate (25) as follows:                      min X,φ,ψ k φ ⊤ Y − φ ⊤ ψ X k 2 F X = [ x 1 , · · · , x N ] k x i k 0 ≤ k for i = 1 , · · · , N φ = ψ ψ φ ⊤ = I n × n . (26) In the following, we present the Parsev al-frame learning algorithm, also referred to as the Parse val K-SVD algorithm, which trains dictionaries as well as sparse coef ficients based on the formulation of (26). B. P arseval K-SVD The l earnin g algorith m is designed through the alternati n g optim ization of ( φ, ψ ) and X . Par ameter ρ 1 is introduced and the objective function of the problem is altered DRAFT 12 to make it the weighted sum of (25) and (26), with the aim of enhancing num erical stability in the optimization process, as follo ws: k φ ⊤ Y − φ ⊤ ψ X k 2 F + ρ 1 k Y − ψ X k 2 F . (27) The alternation does not change the m inimizers, due to the fact that problems (25) and (26) are equivalent 2 . Howe ver , as we s h own below , ρ 1 makes the updated X numerically stable, due to the fa ct that the nominator of (38) is a non-zero v alue. The augmented Lagrangian function of problem (26) is L ρ 2 ,ρ 3 ( X , φ, ψ ; λ 2 , λ 3 ) = ρ 1 k Y − ψ X k 2 F +   φ ⊤ Y − φ ⊤ ψ X   2 F + T race ( λ ⊤ 2 ( ψ φ ⊤ − I n × n )) + ρ 2 2   ψ φ ⊤ − I n × n   2 F (28) + T race ( λ ⊤ 3 ( ψ − φ )) + ρ 3 2 k ψ − φ k 2 F , where ρ 2 and ρ 3 are parameters, and λ 2 and λ 3 are matrices of Lagrangian multipliers of s i zes n × n and n × m , respectiv ely . W e solve the following optim ization problem by minimizing the primal v ariables and maximizing the dual v ariables as follo ws:    max λ 2 ,λ 3 min X,ψ ,φ L ρ 2 ,ρ 3 ( X , φ, ψ ; λ 2 , λ 3 ) , k x i k 0 ≤ k for i = 1 , · · · N . (29) The learning algorithm is deriv ed based on the primal-dual approach, in which the optimal Parse val tight frame is l earned from im age blocks by the alternating direction method of multipli er (ADMM) method [32], [33]. W e adopted the alternating direction method of multiplier (ADMM) due to the robustness of the updates of d u al variables and that f act that this m et h od suppo rt s decompositio n of prim al va riable updates. The updates of φ and ψ are diffe rent from that of X , due to the f act that the latter has sparse constraints on its columns. The updates of φ and ψ are φ k +1 = arg min φ L ρ 2 ,ρ 3 ( X k , φ, ψ k ; λ k 2 , λ k 3 ); (30) ψ k +1 = arg min ψ L ρ 2 ,ρ 3 ( X k , φ k +1 , ψ ; λ k 2 , λ k 3 ); (31) and the update of X is    X k +1 = arg min X L ρ 2 ,ρ 3 ( X , φ k +1 , ψ k +1 ; λ k 2 , λ k 3 ) k x k +1 i k 0 ≤ k for i = 1 , · · · N . (32) 2 φ ⊤ Y = φ ⊤ ψ X i f and only if Y = ψ X because of ψ φ ⊤ = I . DRAFT 13 The dual v ariables λ 2 and λ 3 are updated usi ng the standard gradient ascent m ethod of ADMM ( λ 0 2 and λ 0 3 are initialized using zero matrices): λ k +1 2 ← λ k 2 + ρ 2 ( ψ k +1 ( φ k +1 ) ⊤ − I n × n ); λ k +1 3 ← λ k 3 + ρ 3 ( ψ k +1 − φ k +1 ) . (33) The complexity of updat i ng λ 2 and λ 3 are O ( n 2 ) and O ( nm ) , respectively . In th e following, we detail the procedures in volve d in minimizing the primal variables in (30), (31), and (32), wherein we omit the superscript indices on all variables in order to simply the notation. 1) Up d ate of φ an d ψ : T aking the partial deriv ative of the augmented Lagrangian L ρ 2 ,ρ 3 ( X , φ, ψ ; λ 2 , λ 3 ) with respect to φ and setting the resul t to zero, we obtain the following: (2 Y Y ⊤ − 2 Y X ⊤ ψ ⊤ + 2 ψ X X ⊤ ψ ⊤ − 2 ψ X Y ⊤ ) φ + φ ( ρ 2 ψ ⊤ ψ + ρ 3 I ) = − λ ⊤ 2 ψ + ρ 2 ψ + λ 3 + ρ 3 ψ . (34) Similarly , taking the partial deriv ativ e of the augm ented Lagrangian with respect to ψ and setting the result to zero, we obtain 2 φφ ⊤ ψ + ψ (2 ρ 1 X X ⊤ + ρ 2 φ ⊤ φ + ρ 3 I )( X X ⊤ ) − 1 = (2 ρ 1 Y X ⊤ − λ 2 φ + ρ 2 φ − λ 3 + ρ 3 φ + 2 φφ ⊤ Y X ⊤ )( X X ⊤ ) − 1 . (35) Both (34) and (35) are deri ved from Sylvester matrix equations: A 1 β + β B 1 = C 1 , (36) the solutions of which can be deriv ed by solving linear systems using the least square method obtained by taking the v ec operator on both sides of (36): ( I m × m ⊗ A 1 + B ⊤ 1 ⊗ I n × n ) v ec ( β ) = v ec ( C 1 ) , (37) where ⊗ is the Kronecker produ ct, A 1 , B 1 and C 1 are known m at ri ces, and β is the unknown term. The A 1 , B 1 , C 1 and β corresponding to (34) are 2 Y Y ⊤ − 2 Y X ⊤ ψ ⊤ + 2 ψ X X ⊤ ψ ⊤ − 2 ψ X Y ⊤ , ρ 2 ψ ⊤ ψ + ρ 3 I , − λ ⊤ 2 ψ + ρ 2 ψ + λ 3 + ρ 3 ψ , and φ , respective ly . Like wise, A 1 , B 1 , C 1 , and β of (35) are respectively 2 φφ ⊤ , (2 ρ 1 X X ⊤ + ρ 2 φ ⊤ φ + ρ 3 I )( X X ⊤ ) − 1 , (2 ρ 1 Y X ⊤ − λ 2 φ + ρ 2 φ − λ 3 + ρ 3 φ + 2 φφ ⊤ Y X ⊤ )( X X ⊤ ) − 1 , and ψ . The numbers of constraints and variables for both sy s tems are equal to mn . The complexity of u p dating φ and ψ by solving systems of linear equation s is thus no more than O ( m 3 n 3 ) . DRAFT 14 2) Up d ate of X : Th e objectiv e related t o th e update of X is (27), t he same at that to update φ and ψ . Following the approaches in [27], [34], we iterati vely and exclusi vely update the non-zero coefficients in on e row of X . This ensures that the number of zero coef ficients is increased as the number of iteration s i s increased, because once a coef ficient becomes zero, it remains at zero thereafter . First, for each row of X , we generate a row vector that contains on ly the non-zero entries i n the row . Let ψ k , ( φ ⊤ ψ ) k denote the k -th columns in ψ and φ ⊤ ψ , respectiv ely; and let the 1 × N vector r ⊤ k denote the k -th ro w in X . Furthermore, let p ( i ) be th e i -th non-zero index in r ⊤ k and let   r ⊤ k   0 be t h e number of non-zero coeffi cients. If we let G k be an N ×   r ⊤ k   0 matrix in which ( p ( i ) , i ) is set at one and the other entries are set at zero, t h en r ⊤ k G k is a 1 ×   r ⊤ k   0 vector wi th non-zero entries in row r ⊤ k of X . For example, if r ⊤ k = [0 0 0 2 0 0 1] , then r ⊤ k G k is [2 1] . Next, we consider the update of r ⊤ k G k . Let B = φ ⊤ Y ; let ˆ r ⊤ k = r ⊤ k G k ; and let E k = Y − P j 6 = k ψ j r ⊤ j , F k = B − P j 6 = k ψ j r ⊤ j , ˜ E k = E k G k , and ˜ F k = F k G . E k , F k , and ˜ E k are all known matrices. Then, we obtain the following ρ 1 k Y G k − ψ X G k k 2 F +   φ ⊤ Y G k − φ ⊤ ψ X G k   2 F = ρ 1 k Y G k − ψ X G k k 2 F +   B G k − φ ⊤ ψ X G k   2 F = ρ 1      ( Y − X j 6 = k ψ j r ⊤ j ) G k − ψ k r ⊤ k G k      2 F +      ( B − X j 6 = k ψ j r ⊤ j ) G k − φ ⊤ ψ k r ⊤ k G k      2 F = ρ 1   E k G k − ψ k r ⊤ k G k   2 F +   F k G k − φ ⊤ ψ k r ⊤ k G k   2 F = ρ 1    ˜ E k − ψ k ˆ r ⊤ k    2 F +    ˜ F k − φ ⊤ ψ k ˆ r ⊤ k    2 F . The non-zero coef ficients in the k -th row of X can be updated by solving min ˆ r ⊤ k ρ 1    ˜ E k − ψ k ˆ r ⊤ k    2 F +    ˜ F k − φ ⊤ ψ k ˆ r ⊤ k    2 F . The above has the following clos ed-form solut ion: ˆ r ⊤ k = ρ 1 ψ ⊤ k ˜ E k + ψ ⊤ k φ ˜ F k ρ 1 k ψ k k 2 + k φ ⊤ ψ k k 2 . (38) Thus, the complexity of one update of non-zero entries in a row of X takes O ( n 2 ) . This update of non-zero entries in a ro w is processed from row 1 to row m and repeats j tim es in one update of X . A full description of t h e proposed learning algorithm is giv en in T able I. DRAFT 15 T A B LE I P A R S E V A L K - S V D Parse v al K-SVD Algorithm Deriv ing initial frame ψ and coef ficients X from observation Y using the K-SVD with DCT as the initial dictionary of K- SVD. Input: (i) Set initial ψ (of size n × m ) and X (of size m × N ). (ii) Set the values for ρ 1 , ρ 2 , ρ 3 , and the maximum number of iterations. (iii) The i nitial matrices of Lagrangian multi pliers λ 2 and λ 3 are set t o zero. (iv) T he default iteration number j = 20 for one update of X . For i = 1 to max number of iterations: (v) Update φ by solving (34) followed by updating ψ using (35). (vi) Update Lagrangian matrices λ 2 and λ 3 using (33). (vii) Update X on e row at a time f r om the first to t he last row . This process is repeated j times. end i Output: ψ , φ , and X The complexity associated with one update of φ and ψ of Step (v) in volves solving systems o f mn linear equations, each of whi ch takes at m ost O ( m 3 n 3 ) . M eanwhile, one update o f non-zero entries i n X in Step (vii) takes O ( j n 2 m ) , where j (default is s et at 20 ) is th e nu m ber of it erations associated wi t h one update of X . The complexity of updating λ 2 and λ 3 in Step (vi ) are O ( n 2 ) and O ( n 2 m ) , respectively . T h u s, one update o f Parse val K-SVD (S teps (v), (vi) and (vii )) is O ( m 3 n 3 ) . Let M (typi call y , M = 200 ) be the n u mber of iteratio n s required for th e learning algorit hm to con ver ge. The complexity of l earning the P arse val dictionary i s thus O ( m 3 n 3 M )+ the complexity of learning K - SVD di ct i onary . Howe ver , the Parse v al K-SVD algorithm needs to be compl eted only once, after which th e learned dictionary can b e used for a host of s parse recover y problems. Finally , we address the q u estion of whether th e Parse v al K-SVD algorithm con ver ges. First, let us assum e that X is fixed. Each up date of either φ or ψ would decrease the objective value of (27) or resul t in no change. Then, with fixe d φ and ψ , updati ng non-zero entries in each row of X would reduce this objecti ve value as well as the number of zero elem ents in the row , or result in no change. Executing a series of such st eps would ensure a monotonically non -i n creasing of the objective value, thereby DRAFT 16 guaranteeing con vergence to a lo cal minim um (27). V I . E X P E R I M E N T R E S U L T S W e first applied the P arsev al K-SVD algori t hm to bl ocks of natural images in order to determine wh et h er the learned dictionary is capable of recov ering the desi red results. W e then cond u cted several experiments using natural i mage d at a, in an attempt to demonstrate the applicabilit y of the synt hesis and analysis views of image representatio n when used in conjunctio n with the learned dictionaries (K-SVD and Parse val K-S VD) for the restoration of images. A. Learning P arseval Dictionari es and Reconstruction of Original Imag es Our goal is to derive a Parse val d ictionary capable o f recovering original images using linear operators, based on the tenets of frame t heory . W e consi d er a num ber of implementati o n is s ues. The matrix of observa tions Y is of size n × N , set at 64 × 4096 . The i nitial diction ary ψ and coefficients X in T able I(i) are derived using the K-SVD method. The dictionary ψ is of size n × m , set at 64 × 256 ; and the coeffic ient X of si ze m × N is set at 256 × 4096 . The traini n g data consists of 4 , 096 blocks of size 8 × 8 , which was o btained from generic 512 × 512 gray scale i mages. The 8 × 8 blocks are mapped to 1-D vectors, each of size 64 . The K-SVD removes t he mean values of all blocks and reserves one dictio n ary element exclusi vely for the mean value. Alth o ugh it is preferable that all element s except one have a zero m ean, it is no t a necessary for the design of dictionary . Thu s , we preserv ed the mean v alues of all bl ocks. The initial dictionary for the K-SVD algorithm in T able I i s the d iscrete cosi ne transform (DCT) of size 64 × 256 with the sparsity le vel of K-SVD set at 64 . As previously discussed, the parameter ρ 1 is ro bust because it s value does not change the minimi zer o f (27) and is set to m ake X update num erically stable by not dividing a zero in (38). Thus, ρ 1 is set at 0 . 1 . W e then applied the proposed algorithm to image blocks obtained from the image Barbara, using various s ets of parameters ρ 2 and ρ 3 in conjuncti on with the matrices of Lagrangian multip l iers λ 2 and λ 3 , which were initially set at 0 64 × 64 and 0 64 × 256 , respectiv ely . For each set of parameters, we conducted l iterations and observed the con vergence of functions corresponding to the equality constraints of (26) i n terms of l og 10 ( k ψ l − φ l k 2 F ) , l og 10 (   ψ l φ ⊤ l − I   2 F ) and | l og 10 ( tr ( ψ l φ ⊤ l )) − l og 10 (64) | . The t erm log 10 ( k ψ l − φ l k 2 F ) m easures the closeness of DRAFT 17 the constraint ψ = φ i s reached and the last two t erm s m easure that of t h e const raint ψ l φ ⊤ l = I 64 × 64 . Sparsity le vel k in (26) for each column of X is set at 6 4 . Figure 3 illustrates the distribution of the num ber o f non-zero coefficients i n columns of X . The maximum num ber of iteratio ns i n T able I was set at 200 . As shown i n Figures 4, 5, and 6, this is a sufficient numb er of iterations to achieve numerical con ver g ence for som e parameter sets. As shown in the figures, large ρ 2 and ρ 3 values result in smaller o b jectiv e values. Parameters ρ 2 and ρ 3 determine the s everity of the penalty to the approxim ations of ψ φ ⊤ − I and ψ − φ , respectively , and the sizes of ψ φ ⊤ is 64 × 64 and that of ψ is 64 × 256 . T h e fact that ρ 2 and ρ 3 take higher values, means that the approx i mation become increasingly accurate. Figure 7 dis p lays t he array of reconst ructed images of Barbara, obt ai n ed through decompositio n fol l owed by reconst ruction operators usin g the learned dictionaries, as shown i n Figure 8. The quantity and quality performance depend s on the v alues selected for ρ 2 and ρ 3 . Experimental results indicate t hat respectively setting their values at 10 11 and 10 11 yields optim al PSNR performance in Figure 7. This set of values is then fixed. Figure 9 shows reconstructed images of Lena and Boat with the dictionaries trained using data fr om the im age Ba rbara. Figure 10 p resent s a reconstruction of the image Lena, w h ich was derived using dictionaries learned using 8 , 192 image blocks obtained from Barbara and Boat. This demonstrates that increasing the num ber of observations does not necessarily im prove or quality of reconstructed images. The trained dicti o naries are displ ayed in Figure 1 1. Low v ariance elements (bott o m rows) comp rise separable horizontal and vertical cosi n e wa ves of d i f ferent frequencies, whereas high variance elements (top rows) contain rotation-depend ent v ariations and details. Similar patterns can be found in report s such as [27], [28], [35]. Althoug h they are perc eptually similar , the two dictionaries are dissimi lar in terms of quantity . Figure 12 sho ws the distribution of the distances of the i - th element d i between dictionary D and dictionary E , which is defined as follo ws: 1 − max j | d ⊤ i e j | , (39) where e j is the j -th elements of dictionary E . DRAFT 18 Fig. 3. Bell-shaped distribution of number of non-zero coef ficients in columns of X with a peak at 32 . Numerically , it is impossible to determine whether a number is zero; therefore, we set any number with an absolute value smaller than 10 − 6 as zero. Fig. 4. Comparison of the con vergence rate of l og 10 (   ψ φ ⊤ − I   2 F ) vs. number of iterations using various sets of ρ 1 , ρ 2 , and ρ 3 v alues. Curves corresponding to ρ 2 , ρ 3 ≥ 10 9 sho wed small numerical fl uctuations of approximately 10 − 20 where number of iterations exceeded 100 . B. Application to Ima ge Pr ocessin g: Pre liminar y r esults The optim ization formulations usi ng in th e following experiments were designed to take adv antage of the analysis and synthesis perspecti ves of representation coefficients (deri ved from frame theory) to prob lems in volv ed in image recovery . Our objective in the performance comparison using the learned Parse v al dicti o nary (the ratio o f frame bounds is 1 ) and the learned K-SVD dictionary (the rati o of frame bounds i s 3 . 74 ) DRAFT 19 Fig. 5. Comparison of con ver gence rate of | l og 10 ( tr ( ψ φ ⊤ )) − log 10 (64) | vs. number of iterations using various sets of ρ 1 , ρ 2 , and ρ 3 v alues. Curves corresponding to ρ 2 , ρ 3 ≥ 10 7 sho wed small numerical fl uctuations of approximately 0 where number of iterations e xceeded 100 . First term | log 10 ( tr ( ψ φ ⊤ )) − log 10 (64) | measures the sum of the diagonal elements of ψ φ ⊤ , t he desired solution of which 64 . T hus, the value is subtracted from log 10 (64) . Fig. 6. Comparison of con ver gence rate of log 10 ( k ψ − φ k 2 F ) vs. number of iterations using v arious sets of ρ 1 , ρ 2 , and ρ 3 v alues. Curves corresponding to ρ 2 , ρ 3 ≥ 10 7 sho wed small numerical fluctuations of approximately 10 − 10 where number of it erations exceeded 100 . was not t o identity th e better d i ctionary , but rather to elucidate how dictio naries wi th diffe rent frame bounds would af fect performance. In denoisi n g , tests were designed to le vera ge the sparse synthesi s and redundant DRAFT 20 Fig. 7. Reconstruction of i mage Barbara using flow diagram in Figure 8 wi th analysis and synthesis op erations deri ved using values of ρ 2 and ρ 3 . T op: Original i mage. S econd row f rom l eft to right: Parameters are ρ 2 , ρ 3 = 10 1 , ρ 2 , ρ 3 = 10 3 , ρ 2 , ρ 3 = 10 5 and ρ 2 , ρ 3 = 10 7 . Third row from left to right: Parameters are ρ 2 , ρ 3 = 10 9 , ρ 2 , ρ 3 = 10 11 , ρ 2 , ρ 3 = 10 13 and ρ 2 , ρ 3 = 10 15 . Comparison of peak-signal-to-noise-ratio (P SNR) of reconstructed images: PS NRs of images in the second row are 16.09dB, 27.17dB, 88.27dB and 172.27dB, respective ly; PSNRs of the images i n the third row are 274.32dB, 288.01dB, 271.0dB and 266.51dB, respecti vely . The highest PS NR was from the image derive d using ρ 2 = ρ 3 = 10 11 with ρ 1 fixed at 0 . 1 . Fig. 8. F low diagram of frame-based image reconstruction: Decomposition operation is follo wed by a reconstruction operation. An image is partiti oned into non-ov erlapping blocks. Each block, a column in Y , is decomposed using analysis dictionary φ ⊤ to obtain the f r ame coefficients, a row of Z . The coefficients are then used to reconstruct t he block using synthesis dictionary ψ . representations of an image. In im age compression and filling-the-missing-pi xels, we demonstrate that th e ef fectiv eness of a metho d relies on the frame bound s , which are used to o b tain the redundancy measure of the frame [22] in the learned frame dictionaries. DRAFT 21 Fig. 9. Recon struction of images Boat and Lena. The dicti onaries were t rained, using parameters ρ 1 = 0 . 1 , ρ 2 = 10 11 , and ρ 3 = 10 11 , on image blocks taken from image Barbara. PSNRs of t he reconstructed images Boat and Lena are 287 . 75 dB and 289 . 20 dB, respectiv ely . Fig. 10. Image L ena after reconstruction using dictionaries learned on image blocks of Barbara and Boat, using parameters ρ 1 = 0 . 1 , ρ 2 = 10 11 , and ρ 3 = 10 11 . T he number of training blocks was 8 , 192 , which is double the number of training blocks used to deriv e the dictionaries for Figure 9. PSNR of the reconstructed Lena is 288 . 19 dB. By comparing the P S NR of Lena in Figure 9, we can see t hat increasing the number of tr ai ning blocks does not improv e the P S NR of reconstructed images. Fig. 11. Learned Parse v al frame: (left) trained on Barbara, (right) trained on Boat and Barbara. Columns of the frames are sorted in descendin g order of v ariance and st r etched to their maximal range for the purpose of ill ustration. DRAFT 22 Fig. 12. Distribution of distances between elements in the dictionary learned on Barbara and those in the dictionary learned on Barbara and Boat. (Most of the distances exceed 0 . 5 and the range of distances is [0 , 1] ). 1) Image Denoising: W e address the following image denoi sing problem: Y = X + N , (40) where the ideal im age X is corrupted i n t he presence of additive white and homogeneous Gaussian noise, N , wit h mean zero and standard deviation of σ . W e form u l ate this problem using the frame and the Bayesian approach as a constrained in verse problem. The ob j ectiv e of the o ptimization process is to use the learned dicti onaries to recover t he sparse sy nthesis coeffic ients ( W ) of image with cons trained redun dant frame coefficients ( φ ⊤ ψ W ). min W k W k 1 sub j ect to   Z − φ ⊤ ψ W   2 ≤ ε, (41) where ψ and φ are synthesis and analysis dictionaries, respectively; Z is t he frame coef ficient of Y ; and the v alue of ε is related to σ . In implementation, ε ’ s va lue is manually selected from a set of candidate value rangi n g from 2 to 24 . The ones that yield the best performance are retained for us e in deriving the denoised results from which the a verage performance is then measured. This problem can be sol ved using the basis pursuit denoising (BPDN) method [14], [36]. T able II and T able III compare the av erage PSNR and structural s i milarities (SS IM) between the images recover ed using the proposed meth od with K-SVD and learned Par - DRAFT 23 sev al dictionaries in va rious noi sy en vironments. As shown in th e tables, u sing Parse val dictionaries obtains higher av erage PSNR and SSIM p erformance. In the experiments, the dictionary of K-SVD was u s ed as t he synthesis dictionary and its pseudo-in verse was adopted as the analysis di ctionary . Note that th e performance of the nois y results can be i mproved as long as structure in the frame coef ficients can be explored [37]–[39]. T A B LE II C O M PA R I S O N O F A V E R A G E P S N R F O L L OW I N G N O I S E R E D U C T I O N T A S K . T H E M E A N S O F T H E O R I G NA L I M A G E S W E R E S U B T R AC T E D F RO M O B S E RV A T I O N S I N E X P E R I M E N T S . T H E M E A N S W E R E T H E N A D D E D T O T H E R E S T O R E D I M A G E S F R O M W H I C H P S N R V A L U E S W E R E O B TA I N E D . T H E V A L U E S A R E T H E A V E R AG E S O B T A I N E D F R O M FI V E E X P E R I M E N T S U S I N G T H E S A M E N O I S E L E V E L ( σ ) . E A C H B L O C K O F A FI X E D σ C O M P R I S E S F O U R R O W S . T H E T O P R O W I S T H E P S N R O F A N O I S Y I M A G E . T H E S E C O N D RO W I S T H E P S N R O F A N I M AG E P R O C E S S E D U S I N G T H E K - S V D D I C T I O N A R I E S . T H E T H I R D R O W I S T H E P S N R O F A N I M A G E P RO C E S S E D U S I N G T H E L E A R N E D P A R S E V A L D I C T I O NA R I E S . T H E P S N R U N I T I S dB . σ Barbara Man Lena Hill Boat A verag e noisy image 34.16 3 4.16 34.16 34.16 34.16 34.16 5 K-SVD dictionaries 34.77 3 4.73 35.11 34.73 34.76 34.82 Parse v al dictionaries 35.82 3 5.78 36.53 35.69 35.73 35.91 noisy image 28.14 2 8.14 28.14 28.14 28.15 28.14 10 K-SVD dictionaries 29.82 2 9.98 31.01 30.22 30.14 30.23 Parse v al dictionaries 31.22 3 1.34 32.53 31.44 31.41 31.59 noisy image 24.64 2 4.62 24.62 24.64 24.65 24.63 15 K-SVD dictionaries 27.37 2 7.82 28.98 28.20 27.95 28.07 Parse v al dictionaries 28.78 2 9.10 30.44 29.34 29.21 29.37 noisy image 22.18 2 2.15 22.14 22.17 22.19 22.17 20 K-SVD dictionaries 25.82 2 6.57 27.79 27.04 26.61 26.77 Parse v al dictionaries 27.17 2 7.69 29.06 28.02 27.78 27.94 noisy image 20.31 2 0.25 20.24 20.29 20.29 20.28 25 K-SVD dictionaries 24.73 2 5.66 26.83 26.16 25.62 25.80 Parse v al dictionaries 26.00 2 6.70 28.08 27.11 26.75 26.93 noisy image 18.80 1 8.73 18.71 18.76 18.75 18.75 30 K-SVD dictionaries 23.92 2 5.07 26.25 25.59 24.92 25.15 Parse v al dictionaries 25.11 25.93 27 .24 26.39 25.94 26.12 2) Image Compr ession: W e conducted a comparison of im age compression between the synth esi s dictionary learned using the K-SVD method and its pseudo-in verse as the analys is d ictionary and the proposed Parse v al di ctionaries. The entropy of t h e bit distribution at each bit-plane of frame coefficients is encoded under th e assumption that DRAFT 24 the bi t di stribution i s an independently identically distributed (i.i.d.) random v ariable. The resulting rate-dis tortion graph i s presented in Figure 13. An image was partitioned into 8 × 8 disjoint blocks. All blocks were encoded independently from t he other blocks. The mean of the im age was assumed t o be k n own and subtracted from the im age in the experiments. The mean was t h en added to th e compressed i m age to measure th e PSNR. As shown in Figure 13, using the Parse va l dictionary obtains a better performance than using t he K-SVD dicti o nary . For a fixed PSNR, using t he Parsev al dictionary achieved a sa vings in bits per pixels (bpp) up t o 0 . 7 dB when its value w as abov e 0 . 8 . The frame b o unds can be used as an indicatio n of the correlation between frame coef- T A B LE III C O M PA R I S O N O F A V E R AG E S S I M A S S O C I A T E D WI T H N O I S E R E D U C T I O N TA S K . T H E M E A N S O F T H E O R I G N A L I M AG E S W E R E S U B T R A C T E D F RO M O B S E RV A T I O N S I N E X P E R I M E N T S . T H E M E A N S W E R E T H E N A D D E D T O T H E R E S T O R E D I M A G E S F R O M W H I C H S S I M V A L U E S W E R E O B TA I N E D . T H E V A L U E S A R E T H E A V E R AG E O B TA I N E D F R O M FI V E E X P E R I M E N T S U S I N G T H E S A M E N O I S E L E V E L ( σ ) . E AC H B L O C K O F A FI X E D σ C O N S I S T S O F F O U R R O W S . T H E T O P R O W I S T H E S S I M O F A N O I S Y I M A G E . T H E S E C O N D R O W I S T H E S S I M O F A N I M A G E P R O C E S S E D U S I N G T H E K - S V D D I C T I O N A R I E S . T H E T H I R D R O W I S T H E S S I M O F A N I M A G E P RO C E S S E D U S I N G L E A R N E D P A R S E V A L D I C T I O NA R I E S . σ Barbara Man Lena Hill Bo at A verage noisy image 0.89 0.88 0.85 0.89 0.89 0.88 5 K-SVD dictionaries 0.91 0.90 0.87 0.90 0.90 0.90 Parse v al dictionaries 0.94 0.92 0.91 0.92 0.92 0.92 noisy image 0.72 0.68 0.61 0.69 0.69 0.68 10 K-SVD dictionaries 0.81 0.77 0.76 0.78 0.78 0.78 Parse v al dictionaries 0.86 0.83 0.83 0.82 0.83 0.83 noisy image 0.58 0.52 0.45 0.53 0.54 0.52 15 K-SVD dictionaries 0.74 0.69 0.69 0.69 0.70 0.70 Parse v al dictionaries 0.80 0.75 0.77 0.74 0.76 0.76 noisy image 0.48 0.41 0.34 0.41 0.43 0.42 20 K-SVD dictionaries 0.68 0.63 0.66 0.63 0.65 0.65 Parse v al dictionaries 0.74 0.69 0.72 0.68 0.70 0.71 noisy image 0.40 0.33 0.27 0.33 0.35 0.34 25 K-SVD dictionaries 0.63 0.59 0.62 0.58 0.60 0.60 Parse v al dictionaries 0.70 0.65 0.69 0.63 0.66 0.67 noisy image 0.35 0.27 0.22 0.27 0.29 0.28 30 K-SVD dictionaries 0.59 0.57 0.61 0.55 0.57 0.58 Parse v al dictionaries 0.66 0.61 0.66 0.59 0.62 0.63 DRAFT 25 Fig. 13. Rate-distortion graphs of images: Barbara (T op) and Boat (Bottom). Proposed method uses the learned Parse v al analysis dictionary t o obtain frame coefficients, whereas the K-SVD method is based on the pseudo-in verse of the K-S VD dictionary . ficients, wherein a lower value B A indicates lower correlation between frame coef ficients [40]. The abov e two experiments d em onstrate the advantage of using a Parse v al tight frame to remove the correlation between frame coeffi cients. T h e frame bo u nds A and B of the learned Parse val frame are 1 wherea s those of K-SVD are respectively 1 and 3 . 74 . The smaller B A value of a Par sev al frame is the reason why t he method based on Parserv al diction aries obtains su p erior performance in im age compression. Sim ilarly , a high degree of correlation between frame coeffi cients facilitates the recov ery of m issing pixels, due to the fa ct that m issing coef ficients can b e compensated for using other coef ficients through high frame redundancy . DRAFT 26 3) Filling Missing Pixels: Information p ertain i ng t o an i mage is give n by a s equ ence of irregularly s am pled pixels. Th i s means that the p rob lem lies in determining the means by which to reco ver the image from t he input. W e partitioned an im age into 8 × 8 bl o cks . W e then random l y removed a fraction of the pixels of an image block (between 0 . 1 and 0 . 8 ) by setting their values at ze ro. Let X denote the ideal im age and let Q be a matrix of either 0 or 1 , respectiv ely indi cati ng the positions of t he corrupted and non-corrupted pixels in X . Thus, thei r produ ct, QX , consists exclusively of non-corrupted pixels in X . For example, if X k =[ 2 3 5 0 0 1 ] ⊤ is the k -th column of X , t hen the k -th column of QX is [ 2 3 5 1 ] ⊤ . W e form ulate the prob l em using the frame as a constrained in verse p rob lem, with the aim of recov ering the sparse synthesis coefficients ( W ) of image X : min W k W k 1 sub j ect to k QX − Qψ W k 2 ≤ ε, (42) where ε is a parameter set at 0 . 0 1 for all cases; and ψ is the synthesis dictionary . This problem can be solved usi ng the BPDN algorithm . The reconstructed image is obtained by applying synthesis frame ψ to solution W . The performance of the dictionaries is presented i n T ables IV and V. W e assum ed that the mean of an image i s known. Th e m ean was subtracted in the experiments and th en added to the image to determine the performance. The images in Figure 14 are the results . The perceptual q u ality of the restored images remains goo d as long as no more than 60% of the pixels are missing. This is cons istent with the PS NR and SSIM performance shown in T ables IV and V. The PSNR and SSIM with mis sing fractions below 60% exceed 30 dB and 0 . 88 , respectiv ely . The performance of t he K- SVD dictionary is uniform l y better than that of the P arse val diction ary , as the frame redundancy o f the K-SVD dictionary is higher . DRAFT 27 T A B LE IV C O M PA R I S O N O F A V E R AG E P S N R I N E X P E R I M E N T S I N VO L V I N G M I S S I N G P I X E L S . T H E V A L U E S A R E T H E A V E R AG E O F FI V E E X P E R I M E N T S U S I N G T H E S A M E P E R C E N T AG E O F M I S S I N G P I X E L S . T H E T O P R O W O F E A C H B L O C K I S T H E P S N R O F A N I M AG E D E R I V E D FR O M M I S S I N G P I X E L S . T H E S E C O N D R O W P R E S E N T S T H E R E S U LT S D E R I V E D U S I N G T H E K - S V D D I C T I O N A RY , A N D T H E T H I R D R O W P R E S E N T S R E S U LT S D E R I V E D U S I N G T H E L E A R N E D P A R S E V A L D I C T I O NA RY . T H E P S N R U N I T I S dB . Missing L e vel Barbara Man Lena Hill Boat A verage Corrupted image 31.17 32.87 33.68 34.00 31.95 32.73 10% K-SVD dictionaries 40.37 39.69 42.85 40.48 40.20 40.72 Parse v al dictionaries 38.36 38.11 40.88 39.19 38.60 39.03 Corrupted image 28.16 29.81 30.68 30.96 28.99 29.72 20% K-SVD dictionaries 36.02 35.89 38.95 36.73 36.28 36.78 Parse v al dictionaries 34.30 34.45 37.07 35.58 34.82 35.24 Corrupted image 26.38 28.08 28.92 29.20 27.25 27.97 30% K-SVD dictionaries 33.07 33.42 36.36 34.32 33.62 34.16 Parse v al dictionaries 31.52 32.07 34.54 33.29 32.27 32.74 Corrupted image 25.13 26.82 27.66 27.95 26.01 26.71 40% K-SVD dictionaries 30.68 31.43 34.17 32.40 31.46 32.03 Parse v al dictionaries 29.31 30.10 32.42 31.47 30.19 30.70 Corrupted image 24.16 25.85 26.69 26.98 25.05 25.75 50% K-SVD dictionaries 28.61 29.62 32.19 30.71 29.51 30.13 Parse v al dictionaries 27.37 28.42 30.53 29.86 28.42 28.92 Corrupted image 23.36 25.06 25.90 26.18 24.27 24.95 60% k-SVD dictionaries 26.67 27.92 30.22 29.12 27.67 28.32 Parse v al dictionaries 25.56 26.89 28.73 28.39 26.77 27.27 Corrupted image 22.69 24.39 25.22 25.50 23.60 24.28 70% K-SVD dictionaries 24.82 26.24 28.17 27.53 25.83 26.52 Parse v al dictionaries 23.90 25.44 26.98 26.99 25.17 25.69 Corrupted image 22.11 23.80 24.64 24.92 23.02 23.70 80% K-SVD dictionaries 22.98 24.51 26.02 25.79 23.94 24.65 Parse v al dictionaries 22.36 24.05 25.30 25.51 23.59 24.16 DRAFT 28 T A B LE V C O M PA R I S O N O F A V E R A G E S S I M I N E X P E R I M E N T S I N VO LVI N G M I S S I N G P I X E L S . T H E V A L U E S A R E T H E A V E R AG E O F FI V E E X P E R I M E N T S U S I N G T H E S A M E P E R C E N T AG E O F M I S S I N G P I X E L S . T H E T O P R O W O F E A C H B L O C K I S T H E S S I M O F A N I M A G E D E R I V E D F RO M M I S S I N G P I X E L S . T H E S E C O N D RO W P R E S E N T S T H E R E S U LT S D E R I V E D U S I N G T H E K - S V D D I C T I O N A RY A N D T H E T H I R D R O W P R E S E N T S R E S U LT S D E R I V E D U S I N G T H E L E A R N E D P A R S E V A L D I C T I O NA RY . Missing L e vel Barbara Man Lena Hill Boa t A verage Corrupted image 0.94 0.94 0.95 0.95 0.94 0.95 10% K-SVD dictionaries 0.99 0.98 0.99 0.98 0.99 0.99 Parse v al dictionaries 0.99 0.98 0.98 0.98 0.98 0.98 Corrupted image 0.89 0.89 0.91 0.89 0.89 0.90 20% K-SVD dictionaries 0.98 0.97 0.97 0 .96 0.96 0.97 Parse v al dictionaries 0.97 0.95 0.97 0.95 0.96 0.96 Corrupted image 0.84 0.85 0.87 0.85 0.85 0.85 30% K-SVD dictionaries 0.96 0.94 0.96 0.93 0.94 0.95 Parse v al dictionarie 0.94 0.92 0.94 0.92 0.93 0.93 Corrupted image 0.80 0.80 0.83 0.80 0.80 0.81 40% K-SVD dictionaries 0.93 0.91 0.94 0.90 0.91 0.92 Parse v al dictionaries 0.91 0.89 0.92 0.88 0.89 0.90 Corrupted image 0.75 0.76 0.80 0.75 0.75 0.76 50% K-SVD dictionaries 0.90 0.87 0.91 0.86 0.87 0.88 Parse v al dictionaries 0.87 0.84 0.88 0.84 0.84 0.85 Corrupted image 0.70 0.71 0.77 0.70 0.71 0.72 60% K-SVD dictionaries 0.85 0.82 0.87 0.81 0.81 0.83 Parse v al dictionarie 0.81 0.78 0.84 0.78 0.79 0.80 Corrupted image 0.65 0.67 0.73 0.65 0.66 0.67 70% K-SVD dictionaries 0.77 0.75 0.82 0.74 0.74 0.76 Parse v al dictionaries 0.73 0.71 0.79 0.71 0.72 0.73 Corrupted image 0.59 0.62 0.70 0.60 0.61 0.62 80% K-SVD dictionaries 0.67 0.65 0.75 0.64 0.64 0.67 Parse v al dictionaries 0.63 0.63 0.72 0.63 0.63 0.65 DRAFT 29 Fig. 14. Filling in pixels missing from Barbara. T op: Original image. Second row from left to right respectiv ely presents images with missing pixels percentage 20% and 80% missing pixels. The images i n t he third row were deri ved using the K-SV D dictionary . The images in the f ourth row w ere deri ved using the l earned Parsev al dictionary . DRAFT 30 V I I . C O N C L U S I O N S Frames theory provides the foundation for the desi gn of linear operators used in the decompositio n and reconstructi on of signals. In this paper , we s o ught to formulate a dual frame design in which th e sparse vector obtained t h rough t he decompos ition of a signal is also the sparse so lution representing the signals th at use a reconstruction frame. Our findings demonstrate that such a dual frame does not e xist for over -complete frames. Nonetheless, t h e best approximation t o th e sparse synthesi s s olution can be derived from the analysis coefficient using the canon i cal dual frame. W e took advantage of the analysis and s ynthesis views of signal representation from the frame perspective and proposed optimizatio n form ulations for problems pertaining t o image recovery . W e then compared the performance of the solution s derived using dictionaries with d iff erent frame bounds . Our resul ts revea led a correlation b etween recovered images and the frame bounds o f dictionaries, thereby demonstrating the importance of using different dicti onaries for diffe rent applications. R E F E R E N C E S [1] A. Oppenheim nd A. Willsky , “ Signals and analysis ”, Prentice Hall, 1997. [2] G. Strang and T . Nguyen, “ W avelets and filter banks ”, W ellesley-Cambridg e Press, 1996. [3] S. Mallat, “ The w avelet tour of signal pr ocessing ”, Academic Press, 2008. [4] J. Duffin and A. S chaef fer , “ A class of nonharmonic Fourier series”, Tr ans. Amer . Math. Soc. , V ol. 72, No. 2, pp.341-366 , 1952. [5] R. Y oung, “ An Intr oduction t o nonharmonic F ourier series ”, Academic Press, 1980. [6] I. Duabechies, “ T en lectures on wavelets ”, SIAM, 1992. [7] O. Christensen, “ F rame s and bases: an intro ductory course ”, Birkhaus er, 2008. [8] I. Daubechies, A. Grossmann, and Y . Meyer , “Painless nonorthogonal expansions”, J. Math. Phys. , V ol. 27, pp. 1271-1283, 1986. [9] I. Daubechies, “The wav elet transform, ti me-frequenc y localization and signal analysis”, IEEE Tr ansactions on Information T heory , V ol. 36, No. 5, pp.961-1005, 1990. [10] E. Candes, an d D. Donoho, “Curve lets: a surprising effecti ve nona daptiv e representation for objects with edges”, Curves and Surface Fitting (Edited by Cohen, Rab ut, S chumak er), V anderbilt Uni ve rsity Press, pp. 105-120, 2000. [11] A. Danielyan, V . Katko vnik an d K. Egiazarian, “BM3D F rames and V ariational Image Deblurring”, IEEE T ra nsactions on Image P rocessing , V ol. 21 , No.4, pp. 1715-1728, Nov ember 2011. [12] Y . Qua n and H. Ji an d Z. Shen, “Data-driv en multi-scale non-loca l w ave let frame construction and image recov ery”, J ournal of Scientific Computing , V ol. 63, No. 2, pp. 307-329, 2015. [13] E. Candes and J. Romberg and T . T ao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequenc y information”, IEE E T ransaction s on Information Theory , V ol. 52, No. 2, pp.489-50 9, Feb. 2006. DRAFT 31 [14] S. Chen and D. Donoho and M. Saunders, “ At omic Decomposition by Basic P ursuit”, SIAM Review , V ol. 43, No. 1, pp. 129-159, Apr . 2006. [15] D. Donoho, “Compressed sensing”, IEEE T ransa ctions on Information T heory , V ol. 52, No.4, pp.1289-13 06, Apr . 2001. [16] J. T ropp, “Greedy is good: algorithmic results for sparse representation”, IEE E Tr ansactions on Information Theory , V ol. 50, No. 10, pp.2231-2242, Oct. 2004. [17] A. Cohen, W . Dahmen, and R. Dev ore, “Compressed sensign and best k-term approximation”, J. Amer . Math. Soc. , V ol. 22, pp. 211-231, 2009. [18] E. C andes and T . T ao, “Decoding by Linear Pr ogramming”, IEEE T ransa ctions on Information T heory , V ol. 51, No.12, pp.4203-4215, Dec. 2005. [19] D. Donoho and M. Elad “Optimally sparse representation in general ( nonorthogonal) dictionaries via l 1 minimization”, PNAS , V ol. 100, No. 5, pp.2197-2202, March 2003. [20] A. Mousavi and R. Baraniuk, “Learning t o In vert: Signal Recovery via Deep Conv olutional Networks”, arXiv:1701.038 91 [stat.ML] , 2017. [21] A. Bora and A. Jalal and E. Price and A. G. Dimakis, “Compressed Sensing usin g Generativ e Models”, arXiv:1703.032 08 [stat.ML] , 2017. [22] S. C. Pei and M. H. Y eh, “ An Introduction to Discrete Finite Frames”, IEEE Tr ans. on Signal Pro cessing , V ol. 14, No. 6, pp. 84-96, Nov . 1997. [23] M. Elad, P . Milanfar , and R. Rubinstein, “ Analysis V ersus Synthesis in Signal Priors”, In verse Prob lems , V ol. 23, No.38, pp. 947-968, June 2007. [24] S. Nam, M.E. Davies, M. Elad, and R . Gr i bon val, “Co-sparse Analysis Modeling - Uniqueness and Al gorithms”, ICASSP , May , 2011. [25] S. Nam, M. Davies, M. Elad, and R.Gribon v al, “Cosparse analysis modeling - uniqueness and algorithms”, IEEE IC A SSP , 2001. [26] S. Nam, M. Davies, M. Elad, and R.G r ibon val, “The coparse analysis and model and algorithms”, Appl. Comput. Harmon. Anal. , V ol. 34, No. 1, pp.30-56, Feb . 2013. [27] M. Aharon, M. Elad, and A. Bruckstein, “K-S V D: An Algorithm for Designing Overcomp lete Dictionaries for Sparse Representation”, IEEE T ran s. on Signal P r ocessing , V ol. 54, No. 11, pp. 4311-4322, Nov . 2006. [28] R. Rubinstein, T . Peleg, and M. Elad, “ Analysis K-S V D: A Dictionary-Learning Algorithm for t he Analysis Sparse Model”, IEE E Signal Pr ocessing Magazine , V ol. 61, No. 3, pp. 661-677, Feb . 2013. [29] A Ron and Z Shen, “Frames and stable bases for shift-i n v ariant subspaces of L 2(Rd)”, Canad ian Jou rnal of Mathematics , V ol. 47, No. 5, pp. 1051-1094, 1995. [30] L. Shen and M. Papadak is and I. A. Kakadiaris and I. K onstantinidis and D. Ko uri and D. Hoffman, “Image Denoising using a Tigh t Frame”, IEE E T rans. on Image Pr ocessing , V ol. 15, No. 5, pp. 1254-1263, May 2006. [31] B. Dong and Q. Jiang and C . L iu and Z. Shen, “Mu ltiscale representation of surfaces by tight wav elet frames with applications to denoising”, Applied and C omputational Harmonic Analysis , V ol. 41, No. 2, pp. 561-589, 2016. [32] D. Gabay and B. Mercier, “ A dual al gorithm for the solution of nonlinear v ariational p roblems via finite element approximations” , Computers and Mathematics with Applications , V ol. 2, pp. 17-40, 1976. [33] S. Boy d, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Dist r ibuted Optimization and Statistical Learning via the Alternating Direction Method of Multipliers”, F ounda tions and T r ends in Machine Learning , V ol. 3, No.1, pp.1-122, 2011. DRAFT 32 [34] G.-J. Peng and W .-L . H wang, “ A Proximal Method for Dictionary Updating in S parse Representation”, IE EE T ra ns. on Signal Processin g , V ol. 63, No. 15, pp. 3946-3358, Aug. 2015. [35] M. Lewicki and T . Sejnowski, “Learning ov ercomplete representations”, Neural Comp. , V ol. 12, pp. 337-365, 2000. [36] P . Gill, A. W ang, and A. Molnar , “The I n-Cro wd Algorithm for Fast B asis Pursuit Denoising”, IEEE T ransaction on Signal Processin g , V ol. 59, No. 10, pp. 4595 - 4605, Oct. 2011. [37] M. Elad and M. Aharon, “Image Denoising V ia Sparse and Redundant Representations over Learned Dictionaries”, IEEE T rans. on Signal Processin g , V ol. 15, No. 12, pp. 3738-3745, Dec. 2006. [38] E. B alster , Y . Zheng, and R. Ewing, “Feature-based wav elet shrinkage algorithm for image denoising”, IEEE T ra nsactions on Image P rocessing , V ol. 14, N o. 11, pp. 2024 - 2039, Nov . 2005. [39] J. Ho and W . L. Hwang, “W avelet Bayesian Network Image Denoising”, IEEE T ransa ctions on Image Pr ocessing , V ol. 22, No. 4, pp. 1277-1290, Apr . 2013. [40] K. Gr ¨ ochenig, “ Acceleration of the F rame A l gorithm”, IEEE Tr ans. on Signal Proc essing , V ol. 41, No. 12, pp. 3331-3340, Dec. 1993. DRAFT

Original Paper

Loading high-quality paper...

Comments & Academic Discussion

Loading comments...

Leave a Comment