Applying alternating direction method of multipliers for constrained dictionary learning

Applying Alternating Direction Method of Multipliers for ConstrainedDictionary LearningA.Rakotomamonjy1,∗LITIS EA4108,University of RouenAbstractThis paper revisits the problem of dictionary learning on several points although we still consider an alter-nating optimization scheme.Ourfirst contribution consists in providing a simple proof of convergence for this numerical scheme for a large class of constraints and regularizers on the dictionary atoms.We also investigate the use of a well-known optimization method named alternating direction method of multipliers for solving each of the alternate step of the dictionary learning problem.We show that such an algorithm yields to several benefits.Indeed,it can be more efficient than other competing algorithms such as Iterative Shrinkage Thresholding approach and besides,it allows one to easily deal with mixed constraints or regular-izers over the dictionary atoms or approximation coefficients.For instance,we have induced joint sparsity, positivity and smoothness on the dictionary atoms by means of some total variation and sparsity-inducing regularizers.Our experimental results prove that using these mixed constraints helps in achieving better learned dictionary elements especially when learning from noisy signals.1.IntroductionIn the recent years,a lot of works have been devoted to the problem of sparse representation of signals and images.Several research communities[24,36,15]have focused on this problem in order to develop some new tools for analyzing signals or images,to select features for discrimination tasks or to study the theoretical properties of sparse representation.This large success of“sparsity”is essentially due to the fact that many real world signals or natural images can be represented as a linear combination of few representative elements, denoted as atoms or dictionary elements.One of the key problem related to sparse representation is the choice of the dictionary on which the signal of∗Corresponding authorEmail address:alain.rakotomamonjy@insa-rouen.fr(A.Rakotomamonjy)1Complete postal address:LITIS EA4108,University of Rouen,Avenue de l’universit´e,76800Saint Etienne du Rouvray, France.Preprint submitted to Elsevier August6,2013interest is decomposed.One simple approach is to consider an off-the-shelf dictionary such as wavelet basis, wavelet packet basis,Gabor atoms or Discrete Cosine Basis.A recent trend which achieved state-of-the art results on many low-level signal and image processing tasks is to learn the dictionary itself from the data [13,20,25,19,12].Several algorithms for dictionary learning have been proposed and most of them are based on an alternating optimization scheme which involves a signal sparse coding step and a dictionary optimization step[2,20,13,24,18,33].Dictionary learning algorithms have been successfully applied to a wide range of applications including face recognition[35],image reconstruction[21,34]and image denoising [12].Recently,Yaghoobi et al.[33]have proposed a novel algorithm for dictionary learning.Their approach consists in majorizing the datafidelity term and then in performing alternating optimization of this novel majorizing objective function.Owing to a clever majorization,their algorithm is simply based on several resolutions of constrained or regularized matrix-value approximations.When looked at closer details,at each alternating optimization stage,the algorithm they propose is actually equivalent to a Landweber iteration [17]or a Iterative Shrinkhage-Thresholding algorithm(ISTA)[16].The very nice point of their dictionary learning method is that it can easily handle constraints over dictionary atoms.In this paper,we focus on the problem of constrained dictionary learning that we also address through an alternating optimization scheme.One of our contribution is to provide a simple and graceful proof of convergence of such an alternate scheme for a large class of constraints and regularizers,independently of the algorithms used for solving each alternating stage.Then,we propose to tackle the dictionary learning problem through efficient algorithms that can also easily handle constraints or regularizers on the dictionary elements.More specifically,we investigate novel algo-rithms for solving the dictionary learning problem in which constraints on the dictionary atoms actually involve several constraints or regularizers.The approach we promote is based on alternate direction method of multipliers(ADMM)[11,1].Unlike Yaghoobi’s et al which algorithm takes into account onlyfirst-order gradient-based information,by doing so,we involve at each alternating optimization step,some second-order information(related to the Hessian’s inverse)on the problem which,as experimentally shown,helps us in solving the dictionary learning problem in an efficient way,for medium-scale problems and when several complex constraints are considered.As made clear in the sequel,the gain in efficiency brought by the ADMM approach can be substantial.Another of our contribution is to provide evidence on the benefits of using mixed-constraints for dictionary learning.While several works on dictionary learning acknowledge the need of learning atoms[13],under some norm constraints,few works have addressed the situation where several mixed constraints or regularizers are in play[33,24].In other application domains[1,28],algorithms dealing with compound constraints or2regularizers have recently been proposed,we provide in this paper a simple instantiation of the ADMM approach for coping with several constraints that can be useful for a dictionary learning problem.Then, owing to this algorithm,we show that alternating optimization methods for dictionary learning can handle constraints and regularization that induce the learned dictionary atoms to be for instance jointly sparse, smooth,positive and unit-norm.Experimental results show the benefit of these mixed constraints especially when learning dictionary from noisy signals.The paper is organized as follows.Section2introduces the dictionary learning framework we are dealing with and reviews works related to the one we present.The ADMM algorithm we promote for solving the problem is presented and discussed in Section3.Numerical experiments on toy signal denoising and image denoising problem are presented in Section4.They clearly show that our ADMM approach is faster thanfirst-order algorithms and that using appropriate constraints on the dictionary atoms leads to better denoising.The code that has been used for producing results andfigures is available for reviewing purpose at http://asi.insa-rouen.fr/enseignant/~arakoto/code/dicolearn.html.2.Dictionary learningIn this section,we formally state the dictionary learning problem we are interested in and present the current state of the art methods for solving this problem.2.1.Simultaneous Sparse approximation and dictionary learningThe problem of simultaneous sparse approximation(SSA)[31,29]consists in looking for a sparse decom-position over afixed dictionary of a set of signals under the constraint that,to some extents,they all share the same sparsity profile.Formally,this translates as follows.Suppose we have a set of L signals X=[x1,x2,···,x L],with X∈R N×L,a dictionary D∈R N×M,then the SSA problem is:min A∈R M×L 12X−DA 2F+λΩ(A)(1)whereΩ(A)is a sparsity-inducing regularizer on A andλa trade-offparameter that balances the data fitting term(the square loss)and the regularization term.Typically,for inducing a shared sparsity profile on the signal approximation,which means that all x i should preferably be approximated by the same set of dictionary elements,one can consider a mixed-norm regularizer over the rows of A of the form:Ω(A)= Mi=1A i,· p q1/p(2)where typically p=1and q∈{2,∞}[10,31].Note that for all matrices,we will denote as A i,·and A·,j the i-th row and the j-th column.When the sparsity profile is known to differ for all signals to be approximated,3one instead may consider p=1and q=1[6]which allowsΩ(A)to decouple and thus untie coefficients of A.Dictionary learning problems go beyond simultaneous sparse approximation by jointly optimizing over the coefficient matrix A and the dictionary D leading to the optimization problemmin A,D∈D 12X−DA 2F+λΩ(A)(3)where D is a set that imposes some constraints on the dictionary elements chosen according to prior knowledge on the problem or chosen so as to help in resolving the scale invariance of the problem.The most frequent constraint imposed to D,and already used in other works,is that each column of D has a unit 2norm or that the Frobenius norm on D is also unit[13,20,33].In a more general form,we consider the following problem for dictionary learningmin A,D Φ(A,D)=12X−DA 2F+λAΩA(A)+λDΩD(D)(4)whereΩA(·)andΩD(·)are general regularizers on A and D that can involve mixed terms(e.g non-negativity and sparsity)or indicator functions that are related to projections on convex set defining some constraints (e.g unit-norm dictionary elements).We will provide more details on these regularizers in the sequel.2.2.Related worksSeveral methods for solving the dictionary learning problem given in Equation(4)have been investigated in the last decade.These algorithms usually solve problem(4)by means of an alternating optimization procedure which is summarized in Algorithm1.Basically,this algorithm consists in alternatively learning the sparse approximation coefficient matrix A when the dictionary is consideredfixed and then in updating the dictionary withfixed matrix approximation A.For instance,Engan et al.[13]consider the problem with no constraints on A and impose unit-norm dictionary elements.Accordingly,the solution at each alternate step can be straightforwardly obtained by solving the related least-square problem over A and D.For the D update,assuming that(A k A Tk)−1exists(for instance when L≥M),this gives:D k+1=XA T k(A k A T k)−1(5) and unit-norm atoms are obtained by normalizing each column of D k+1.In a similar way,Kreutz-Delgado et al.have proposed an alternating optimization algorithm that solves a Bayesian interpretation model of the dictionary learning problem[20].The recent work of Yaghoobi et al.investigates problem(4)when several simple types of constraints and regularizations on D(typically unit-norm constraints or Frobenius norm boundedness)as well as sparsity-inducing regularizers on A are in play.For instance they consider both unit-norm Frobenius or unit 2norm4Algorithm 1:Alternate optimization for dictionary learning1:set k=1,initialize A (1),D (1)2:repeat3:A k +1=min A Φ(A ,D k )4:D k +1=min D Φ(A k +1,D )5:k ←k +16:until stopping criterion is metconstraints,jointly to sparsity-inducing regularizers of the same form of the one given in Equation (2)on the dictionary atoms.The main idea of their algorithm is based on the use of a surrogate function that majorizes the data fidelity term X −DA 2F in problem (4).For instance,when updating the dictionary D ,they replace the problem of minimizing Φ(A k ,D )with respects to D by the following surrogate problem :min D 1 X −DA k 2F +λD ΩD (D )+C D 1 D −D k −1 2F −1 DA k −D k −1A k 2F (6)where C D > A k A 2k 2is a constant ensuring that the majorization holds.In order to get a better insight on how this surrogate approach works,we provide below the solution of this minimization problem .Let us denote as J the objective function of Equation (6).By expanding all the square norms,we can show that2·J = X 2F −2tr(X DA k )+ DA k 2F +λD ΩD (D )+C D D 2F −2C D tr(D D k −1)+C D D k −1 2F − DA k 2F +2tr(A k D D k −1A k )+ D k −1A k 2F=C D D 2F −2C D tr(D T (D k −1+1C D (X −D k −1A k )A k ))+λD ΩD (D )+G where G is a term that englobes all components of the equation that do not depend on D .Note that the second equality is obtained because we have tr(A k D D k −1A k )=tr(D D k −1A k A k )and tr(X DA k )=tr((DA k ) X )=tr(D XA k ).According to this equality,we can deduce that problem given in Equation (6)is equivalent to the following one :min D 12 D −(D k −1+1C D (X −D k −1A k )A k ) 2F +λC D ΩD (D )(7)This problem can be easily solved by introducing the notion of proximal operator ,a key element of several algorithms such as the forward-backward splitting [8,7]or the Iterative Shrinkage/Thresholding algorithm (ISTA)[16]known.The proximal operator [26]of a lower semi-continuous convex function Ωis defined asprox Ω(ˆD )=arg min D 12D −ˆD 2F +Ω(D )and is unique and well-defined.For interested readers,a comprehensive review on basic properties of proximal operators and optimization methods based on proximal operators have recently been given in [9].Note that for several functions Ω(D ),this proximal operator has a closed-form solution.For instance,when Ω(D )=λ i,j |d i,j |,which corresponds to a soft thresholding operation for matrices,then prox λ i,j |d i ,j |(ˆD)5is the matrix whose(i,j)entry is sign(ˆd i,j)max(|ˆd i,j|−λ,0).WhenΩ(D)takes the form of an indicator over a convex set,then proxΩ(ˆD)is the projection ofˆD on that set.From this definition of the proximal operator,the solution of problem(7),as proposed by Yaghoobi et al., can now be expressed asD k=proxλC D ΩD(D k−1+1C D(X−D k−1A k)A k))(8)Hence,Yaghoobi et al.solve problem(6)byfirst updating D k−1according to a gradient descent step,since(X−D k−1A k)Ak is the negative gradient of12X−D k−1A k) 2F,and then in“projecting”this novel updateon the convex set imposed by the regularizers or the constraintsΩD.We can highlight that this surrogate approach can help in dealing with several types of constraints over the dictionary atoms,as long as the proximal operator ofΩD can be(easily)computed.We also note that this update equation is also the update we would have obtained if an iterative shrinkage thresholding algorithm(ISTA)[3]has been used for solving problem min DΦ(A k,D).Regarding ISTA,while this algorithm is mainly devoted to sparse coding approach based on a 1regularizer,it is generic enough to handle a large class of regularizers and thus can also be used for the dictionary update.Still in the framework of alternating optimization on A and D,several interesting research outcomes have also recently been proposed by the machine learning community[22,24].Efforts have essentially been spent for designing appropriate regularizers or constraints on the matrix code A so as to take into account some knowledge over the learning problems.For instance,Jenatton et al.have investigated methods that are able to exploit a hierarchical relationship between the dictionary elements through the definition of appropriate regularizersΩA(A)[18].They have also proposed methods that are able to learn some structures on the dictionary elements[19].As an example,they imposed atoms to have contiguous non-zeros patterns.These works are relevant to ours since the optimization strategy they consider is based,as in Yaghoobi et al.[33]on ISTA and its fast version FISTA.Indeed,they use an alternating optimization scheme in which the building block of the sparse coding and dictionary update involves afirst-order gradient step and the application of a proximal operator.In addition,they introduce some novel regularizers and constraints that can be mixed with other ones through the framework we promote.The K-SVD of Aharon et al.[2]presents a different perspective on the dictionary learning problem.Indeed, their dictionary update step also involves some updates of the approximation coefficients.This update is performed through a SVD decomposition on a representation error matrix of the set of signal X.The resulting decomposition defines a novel dictionary element and the approximation coefficients related to that atom. This SVD decomposition is performed M times(M being the dictionary size)for each dictionary update in the alternating optimization algorithm(1).6It is also interesting to mention that two works have lately investigated the problem of online dictionary learning[30,24].At the contrary of the batch approaches presented above,in the online framework it is supposed that signal examples x i are available on thefly.Skretting et al.proposes an extension of the method of optimal direction based on recursive least-squares for dealing with this online framework while Mairal et al.propose a stochastic approximation algorithm that is able to deal with constraints on A and D. These methods are those tailored to very large-scale(in the number of signals)problems.The approach we describe in this paper is more suitable for medium-scale problems since it needs the computation of an inverse Hessian matrix which dimension depends on the number of dictionary elements.However,our experimental results show that our approach is still competitive compared to other batch algorithms such as the one of Yaghoobi et al.for these medium-scale problems.In this paper,we consider the problem of dictionary learning in a batch framework as given by Equation(4). In particular,we are interested in problems where several constraints on the dictionary atoms are imposed. Our objective is to efficiently learn dictionary atoms while taking into account some prior knowledge on the dictionary.We note that the MOD algorithm of Engan et al.[13]indeed uses second-order information )−1in the dictionary update but it can not deal with specific constraints on the dictionary atoms.(A k AkContrarily,the approach of Yaghoobi et al.[33]can easily deal with constraints but it only usesfirst-order information in their update stages.We achieve these two objectives(efficiency and mixed-constraints)by means of an alternating direction of multipliers method(ADMM)approach.Our efficient dictionary learning algorithm is still based on the alternating optimization approach as described in Algorithm 1.The novelty we bring is that we solve each of the alternating problem using an ADMM method,which involves some second-order information in the optimization scheme,making it efficient for medium-scale problems.Since the alternating optimization approach nicely decouples the dictionary atom learning stage and the sparse coding step,This allows us to deal with several constraints imposed on the dictionary atoms.Since these mixed-constraints actually define a novel(global)constraint/regularizerΩ(D),we deal with this newΩ(D)by numerically computing its proximal operator.For this purpose,we again use an algorithm which is based on an ADMM.Interestingly,owing to a clever implementation,computing the proximal operator ofΩ(D)only needs the proximal operator of each single constraint composingΩ(D).From our numerical results,we show that our algorithm is indeed more efficient than otherfirst-order methods and that substantial improvements of the learned dictionary can be achieved by using mixed constraints,in particular in a noisy signal setting.73.Alternating Direction Method of Multiplier and Application to Dictionary LearningIn this section,wefirst review augmented Lagrangian method and the alternating direction method of multipliers.Then,we detail how such a framework can be instantiate for the dictionary learning problem and how it can be used to handle constraints over the dictionary elements.3.1.Augmented Lagrangian and Alternating Direction Method of MultipliersWe are interested in an optimization problem which general form is:minUf1(U)+f2(GU)(9) where f1(·)and f2(·)are convex lower semi-continuous functions over the set of matrices,U∈R n1×n2 and G∈R n ×posite optimization problems based on the sum of two functions are emcompassed into this framework by choosing G as the identity matrix.By introducing an auxiliary variable,the above unconstrained problem is equivalent to the following constrained oneminU,Vf1(U)+f2(V)s.t.GU−V=0(10) One advantage of introducing this auxiliary variable V is that it allows the decoupling of the two functions f1and f2since each of them now applies to one specific optimization variable U or V.Because of this property,in some situations,it thus may be easier to solve the equality-constrained problem(10)than the unconstrained one(9).This occurs for instance when f1or f2has a closed-form proximal operator.One typical way for solving such an equality-constrained problem(10)is to deploy an augmented Lagrangian approach[27].Basically,this consists in defining the problem’s Lagrangian augmented with a quadratic term that penalizes U and V not satisfying the equality contraint:L(U,V,Λ)=f1(U)+f2(V)−tr(ΛT(GU−V))+µ2GU−V 2FwithΛbeing a dense matrix of Lagrangian multipliers related to the equality constraints andµa positive parameter that balances the quadratic penalization.As described by Nocedal and Wright[27],the idea of augmented Lagrangian method is tofind a saddle point of the Lagrangian L(U,V,Λ)which coincides with the solution of the original problem(10).This is done byfinding a minimizer of the augmented Lagrangian withΛandµfixed,and then in updating these parameters and repeating these two steps until convergence. ForfixedΛandµ,after some simple algebras,minimizing the Lagrangian can be shown to be equivalentmin U,V f1(U)+f2(V)+µ2GU−V−Λµ2F(11)leading thus to the iterative algorithm given in Algorithm2with Z being1µΛ.8Algorithm2:Augmented Lagrangian method1:set k=1,initialize U1,V1,Z12:repeat3:(U k+1,V k+1)=argmin U,V f1(U)+f2(V)+µ2 GU−V−Z k 24:Z k+1=Z k+V k+1−GU k+15:k←k+16:until stopping criterion is metAlgorithm3:Alternating Direction Method of Multipliers1:set k=1,initialize V(1),Z(1)2:repeat3:U k+1=argmin U f1(U)+µ2 GU−V k−Z k 2F4:V k+1=argmin V f2(V)+µ2 GU k+1−Z k−V 2F5:Z k+1=Z k+V k+1−GU k+16:k←k+17:until stopping criterion is met8:output:U k+1This minimization problem is an unconstrained optimization problemjointly on U and V,involving in the objective function both f1(·),f2(·)and a quadratic term coupling U and V.Hence,it seems that little has been gained by introducing the auxiliary variable V and by considering Augmented Lagrangian approach. However,one should note that problem(11)provides us with an useful feature since f1and f2are now separable in U and V.This separability suggests thus the use of a block-coordinate descent(BCD)for solving problem(11).Alternating direction method of multipliers exploits this idea of block-coordinate optimization instead of a joint optimization on the two variables.The idea of ADMM is to optimize alternatively over U and V. Interestingly,ADMM considers only a single iteration of a block-coordinate descent over U and V and then updates the Lagrangian multipliers.This is different from an exact implementation of a BCD algorithm which would have alternated over U and V until convergence before updating the multipliersΛ.This results in the ADMM approach described in Algorithm3.More details on ADMM can be found in[11,14,4].Note that thefirst minimization problem in the ADMM algorithm involves a regularized least-square problem, where the design matrix is G and the regularizer is f1(U),its resolution thus can be expensive.Instead, the minimization of the second problem makes use of the proximal operator of1µf2(V),which can be cheap to compute if it has a closed-form solution.To make clear this connection,we can simply rewrite Line4of Algorithm3)asargminV 12V−(GU k+1−Z k) 2F+1µf2(V)which by definition is prox1µf2(V)(GU k+1−Z k).One of the very nice feature of this alternating direction method of multipliers is that it converges towards9a minimizer of problem(9)or(10)under the simple hypotheses that G is full-column rank and f1and f2 are proper convex functions[11]and this for an arbitrary choice ofµ.In addition,the convergence proof of ADMM methods given by Eckstein and Bertsekas states that,even if the two minimization steps(Lines 3and4)in the algorithm are not solved exactly,then as long as the two sequences of optimization errors are absolutely summable[11],the sequence generated by the ADMM algorithm still converges towards a solution of the problem(10).This implies that in practice,if the proximal operator of f2has to be computed numerically,the ADMM algorithm is still guaranteed to convergence as long as the numerical optimization error sufficiently decreases,which means that the stopping criterion of the proximal operator computation has to be tighter as ADMM iterations go.3.2.Application of ADMM to dictionary learning:global frameworkThe approach we propose for solving the dictionary learning problem is still based on the alternating approach presented in Algorithm1.Each of these minimization problems is here solved by means of an ADMM algorithm.According to the general problem of dictionary learning given in Equation(4),we define the functions f1(·), f2(·)asf1(A,D)=12X−DA 2F and f2(A,D)=λAΩA(A)+λDΩD(D)andΦ(A,D)=f1(A,D)+f2(A,D).NowΦ(A,D k)andΦ(A k,D)are clearly defined for each minimization step of Algorithm1.Before delving into the details of how ADMM have been implemented for solving these problems,we provide a proof of convergence of the alternating minimization scheme given in Algorithm1for dictionary learning.Proposition1.Suppose that the constraints and regularizers on the dictionaryΩD(D)and on the coefficient matrixΩA(A)are convex,possibly non-smooth and lower semicontinuous functions,leading to the dictionary learning problem given in Equation(4),then Algorithm1converges towards a stationary point of Equation (4).Proof.Note that as a quadratic function f1(A,D)is continuous and differentiable on its domain and Φ(A,D k)andΦ(A k,D)are convex respectively in A and D.Furthermore,it can be easily shown that f1(A,D)is coercive and thus the level sets of f1(A,D)+f2(A,D)are bounded.Owing to all these prop-erties,we invoke Theorem5.1of[32]and we can conclude that the sequence of{A k,D k}generated by Algorithm1,which is a block-coordinate descent with a cyclic rule,converges towards a stationary point of problem(4).10We want to emphasize that this proof of convergence is independent to how each alternate step is solved (by means of ADMM or any other methods).Besides being simpler,it is thus more general than the one provided by Yaghoobi et al.[33]which is specific to their algorithm.Interestingly,convergence holds for a large class of contraints and regularizers on the dictionary and the coefficient matrices.In accordance with the alternating optimization scheme given in Algorithm1,we now want to solve either minΦ(A,D k)or minΦ(A k+1,D)using the ADMM scheme described in Algorithm3.In order to make the parallel between the ADMM notations in Equation(9)and the dictionary learning context,we rewrite the dictionary update minimization problem as(with D=U))min U 12X−UA k+1 2Ff1(A k+1,D)+λAΩ(A k+1)+λDΩD(U)f2(A k+1,D)(12)and the weight update problem asmin U 12X−D k U 2Ff1(A,D k)+λAΩ(U)+λDΩD(D k)f2(A,D k)(13)The next subsections describe how we solve each of these alternating step through the ADMM scheme given in Algorithm3.3.3.Solving thefirst ADMM stepFor both the dictionary and weight update problems given in Equations(12)and(13),we naturally set G=I so that the auxiliary variable V is a direct copy of the variable to optimize A or D denoted as U in the sequel,regardless of the problem solved.Hence,the minimization problem which considers f1(·,·)in the ADMM approach(line3of Algorithm3), boils down to be a least-square problem.For the D update problem,i.e U=D,the problem ismin U 12X−UA k+1 2F+µ2U−V k−Z k 2FA closed-form solution can be straightforwardly obtained by setting the gradient of this problem to zero leading toU k+1=(XA T k+1+V k+Z k)(A k+1A T k+1+µI)−1(14) Similarly,when A is updated while D being keptfixed,we have:U k+1=(D T k D k+µI)−1(D T k X+V k+Z k)(15)In these two update equations,the terms A k+1A Tk+1+µI and D TkD k+µI can be interpreted as the regularizedHessian of respectively f1(A k+1,D)and f1(A,D k+1).Compared to the dictionary learning algorithm of11。