ebook img

Iterative Gaussianization: from ICA to Random Rotations PDF

2.3 MB·
Save to my drive
Quick download
Download
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Iterative Gaussianization: from ICA to Random Rotations

2011(cid:13)cIEEETRANSACTIONSONNEURALNETWORKS 1 Iterative Gaussianization: from ICA to Random Rotations Valero Laparra, Gustavo Camps-Valls, Senior Member, IEEE and Jesús Malo Abstract—Most signal processing problems involve the chal- The conventional approach is to transform data into a do- lenging task of multidimensional probability density function main where interesting features can be easily (i.e. marginally) (PDF) estimation. In this work, we propose a solution to this characterized.Inthatcase,onecanapplywell-knownmarginal problembyusingafamilyofRotation-basedIterativeGaussian- techniques to each feature independently and then obtain a ization(RBIG)transforms.Thegeneralframeworkconsistsofthe 6 sequential application of a univariate marginal Gaussianization description of the multidimensional PDF. The most popular 1 transform followed by an orthonormal transform. The proposed approachesrelyonlinearmodelsandstatisticalindependence. 0 procedure looks for differentiable transforms to a known PDF However, they are usually too restrictive to describe general 2 so that the unknown PDF can be estimated at any point of data distributions. For instance, principal component analysis n the original domain. In particular, we aim at a zero mean unit (PCA) [7], that reduces to DCT in many natural signals such covariance Gaussian for convenience. a J RBIG is formally similar to classical iterative Projection as speech, images and video, assumes a Gaussian source [3], Pursuit (PP) algorithms. However, we show that, unlike in [7]. More recently, linear ICA, that reduces to wavelets in 1 PP methods, the particular class of rotations used has no naturalsignals,assumesthatobservationscomefromthelinear 3 special qualitative relevance in this context, since looking for combination of independent non-Gaussian sources [8]. In interestingness is not a critical issue for PDF estimation. The ] general,theseassumptionsmaynotbecompletelycorrect,and L key difference is that our approach focuses on the univariate part (marginal Gaussianization) of the problem rather than on residual dependencies still remain after the linear transform M the multivariate part (rotation). This difference implies that one thatlooksforindependence.Asaresult,anumberofproblem- . may select the most convenient rotation suited to each practical oriented approaches have been developed in the last decade t application. a to either describe or remove the relations remaining in these The differentiability, invertibility and convergence of RBIG t linear domains. For example, parametric models based on s are theoretically and experimentally analyzed. Relation to other [ joint statistics of wavelet coefficients have been successfully methods,suchasRadialGaussianization(RG),one-classsupport vector domain description (SVDD), and deep neural networks proposed for texture analysis and synthesis [5], image coding 1 (DNN)isalsopointedout.ThepracticalperformanceofRBIGis [9] or image denoising [10]. Non-linear methods using non- v 9 successfullyillustratedinanumberofmultidimensionalproblems explicit statistical models have been also proposed to this such as image synthesis, classification, denoising, and multi- 2 end in the denoising context [11], [12] and in the coding information estimation. 2 context[13],[14].Infunctionapproximationandclassification 0 IndexTerms—Gaussianization,IndependentComponentAnal- problems, a common approach is to first linearly transform 0 ysis (ICA), Principal Component Analysis (PCA), Negentropy, the data, e.g. with the most relevant eigenvectors from PCA, . Multi-information, Probability Density Estimation, Projection 2 Pursuit. and then applying nonlinear methods such as artificial neural 0 networks or support vector machines in the reduced dimen- 6 sionality space [3], [4], [7]. 1 : I. INTRODUCTION Identifying the meaningful transform for an easier PDF v MANY signal processing problems such as cod- descriptioninthetransformeddomainstronglydependsonthe i X ing, restoration, classification, regression or synthesis problem at hand. In this work we circumvent this constraint r greatlydependonanappropriatedescriptionoftheunderlying by looking for a transform such that the transformed PDF is a probability density function (PDF) [1]–[5]. However, density known. Even in the case that this transform is qualitatively estimation is a challenging problem when dealing with high- meaningless, being differentiable, allows us to estimate the dimensionalsignalsbecausedirectsamplingoftheinputspace PDF in the original domain. Accordingly, in the proposed is not an easy task due to the curse of dimensionality [6]. As context, the role (meaningfulness) of the transform is not a result, specific problem-oriented PDF models are typically that relevant. Actually, as we will see, an infinite family of developed to be used in the Bayesian framework. transformsmaybesuitabletothisend,soonehasthefreedom to choose the most convenient one. Copyright (c) 2011 IEEE. Personal use of this material is permitted. In this work, we propose to use a unit covariance Gaussian Permission from IEEE must be obtained for all other users, including as target PDF in the transformed domain and iterative trans- reprinting/republishingthismaterialforadvertisingorpromotionalpurposes, creating new collective works for resale or redistribution to servers or lists, formsbasedonarbitraryrotations.Wedosobecausethematch or reuse of any copyrighted components of this work in other works. DOI: between spherical symmetry and rotations makes it possible 10.1109/TNN.2011.2106511. to define a cost function (negentropy) with nice theoretical This work was partially supported by projects CICYT-FEDER TEC2009- 13696,AYA2008-05965-C04-03andCSD2007-00018,andgrantBES-2007- properties.Thepropertiesofnegentropyallowustoshowthat 16125. one Gaussianization transform is always found no matter the Authors are with Image Processing Laboratory (IPL) at the Universitat selected class of rotations. de València, Spain. Address: C/ Catedrático Escardino, 46980 - Paterna (València)Spain.E-mail:{lapeva,gcamps,jmalo}@uv.es,http://ipl.uv.es The remainder of the paper is organized as follows. In 2 2011(cid:13)cIEEETRANSACTIONSONNEURALNETWORKS Section II we present the underlying idea that motivates decomposed as the sum of two non-negative quantities, the the proposed approach to Gaussianization. In Section III, multi-information and the marginal negentropy: we give the formal definition of the Rotation-based Itera- J(x)=I(x)+J (x). (2) tive Gaussianization (RBIG), and show that the scheme is m invertible, differentiable and it converges for a wide class This can be readily derived from Eq. (5) in [26], by con- of orthonormal transforms, even including random rotations. sidering as contrast PDF (cid:81) q (x ) = N(0,I). The multi- i i i Section IV discusses the similarities and differences of the information is [27]: proposedmethodandProjectionPursuit(PP)[15]–[18].Links (cid:81) I(x)=D (p(x)| p (x )) (3) toothertechniques(suchassingle-stepGaussianizationtrans- KL i i i forms[19],[21],one-classsupportvectordomaindescriptions Multi-information measures statistical dependence, and it is (SVDD) [22], and deep neural network architectures [23]) are zero if and only if the different components of x are indepen- alsoexplored.SectionVshowstheexperimentalresults.First, dent. The marginal negentropy is defined as: we experimentally show that the proposed scheme converges d to an appropriate Gaussianization transform for a wide class (cid:88) J (x)= D (p (x )|N(0,1)) (4) of rotations. Then, we illustrate the usefulness of the method m KL i i in a number of high-dimensional problems involving PDF es- i=1 timation: image synthesis, classification, denoising and multi- GivenadatadistributionfromtheunknownPDF,ingeneral information estimation. In all cases, RBIG is compared to both I and Jm will be non-zero. The decomposition in (2) relatedmethodsineachparticularapplication.Finally,Section suggests two alternative approaches to reduce J: VI draws the conclusions of the work. 1) Reducing I: This implies looking for interesting (inde- pendent) components. If one is able to obtain I = 0, II. MOTIVATION thenJ =J ≥0,andthisreducestosolvingamarginal m This section considers a solution to the PDF estimation problem. Marginal negentropy can be set to zero with problem by using a differentiable transform to a domain the appropriate set of dimension-wise Gaussianization with known PDF. In this setting, different approaches can be transforms,Ψ.Thisiseasyaswillbeshowninthenext adopted which will motivate the proposed method. section. Let x be a d-dimensional random variable with (unknown) However, this is an ambitious approach since looking PDF, p (x). Given some bijective, differentiable transform of forindependentcomponentsisanon-trivial(intrinsically x x into y, G :Rd →Rd, such that y=G(x), the PDFs in the multivariate and nonlinear) problem. According to this, original and the transformed domains are related by [24]: linear ICA techniques will not succeed in completely removing the multi-information, and thus a nonlinear (cid:12)(cid:12)dG(x)(cid:12)(cid:12) post-processing is required. px(x)=py(G(x))(cid:12)(cid:12) dx (cid:12)(cid:12)=py(G(x))|∇xG(x)|, (1) 2) Reducing Jm: As stated above, this univariate problem is easy to solve by using the appropriate Ψ. Note where |∇ G| is the determinant of the Jacobian matrix. x that I will remain constant since it is invariant under Therefore, the unknown PDF in the original domain can be dimension-wisetransforms[27].Inthisway,oneensures estimated from a transform of known Jacobian leading to an that the cost function is reduced by J . Then, a further m appropriate(knownorstraightforwardtocompute)targetPDF, processing has to be taken in order to come back to p (y). y a situation in which one may have the opportunity to One could certainly try to figure out direct (or even closed removeJ again.Thisadditionaltransformmayconsist m form) procedures to transform particular PDF classes into a of applying a rotation R to the data, as will be shown target PDF [19], [21]. However, in order to deal with any in the next section. possible PDF, iterative methods seem to be a more reasonable The relevant difference between the approaches is that, in approach. In this case, the initial data distribution should be the first one, the important part is looking for the interesting iteratively transformed in such a way that the target PDF is representation (multivariate problem), while in the second progressively approached in each iteration. approach the important part is the univariate Gaussianization. Theappropriatetransformineachiterationwouldbetheone In this second case, the class of rotations has no special thatmaximizesasimilaritymeasurebetweenPDFs.Asensible qualitative relevance: in fact, marginal Gaussianization is the cost function here is the Kullback-Leibler divergence (KLD) only part reducing the cost function. betweenPDFs.Inordertoapplywell-knownpropertiesofthis The first approach is the underlying idea in Projection measure[25],[26],itisconvenienttochooseaunitcovariance Pursuitmethodsfocusedonlookingforinterestingprojections Gaussian as target PDF, p (y) = N(0,I). With this choice, y [16], [17]. Since the core of these methods is looking for thecostfunctiondescribingthedivergencebetweenthecurrent meaningful projections (usually ICA algorithms), they suffer data,x,andtheunitcovarianceGaussianisthehereaftercalled fromabigcomputationalcomplexity:forexample,robustICA negentropy1, J(x) = D (p(x)|N(0,I)). Negentropy can be KL algorithms such as RADICAL [28] would lead to extremely 1Thisusageofthetermnegentropyslightlydiffersfromtheusualdefinition slow Guassanization algorithms whereas relatively more con- [25]wherenegentropyistakentobeKLDbetweenpx(x)andamultivariate venient alternatives such as FastICA [29] may not converge Gaussianofthesamemeanandcovariance.However,notethatthisdifference in all cases. This may explain why, so far, Gaussianization hasnoconsequenceassumingtheappropriateinputdatastandardization(zero meanandunitcovariance),whichcanbedonewithoutlossofgenerality. techniques have been applied just to low-dimensional (audio) LAPARRAETAL.:ITERATIVEGAUSSIANIZATION:FROMICATORANDOMROTATIONS 3 Fig.1. ExampleofmarginalGaussianizationinsomeparticulardimensioni.Fromlefttoright:marginalPDFofxi,uniformizationtransformu=Ui(xi), PDFoftheuniformizedvariablep(u),GaussianizationtransformG(u),andPDFoftheGaussianizedvariablepi(Ψi(xi)). signals in either simple contexts based on point-wise non- the PDF in the original domain from the Jacobian of the linearities [30], [31], or after ad hoc speech-oriented feature transform in each point, cf. Eq. (1). Invertibility guarantees extraction steps [32]. In this work, we propose following that the transform is bijective which is a necessary condition the simpler second approach using the most computationally to apply Eq. (1). Additionally, it is convenient for generating convenient rotation. Intentionally, we do not pay attention to samples in the original domain by sampling the Gaussianized the meaningfulness of the rotations. domain. Before getting into the details, we take a closer look at the III. ROTATION-BASEDITERATIVEGAUSSIANIZATION basic tool of marginal Gaussianization. Marginal Gaussian- (RBIG) ization in each dimension i and each iteration k, Ψi , can (k) bedecomposedintotwoequalizationtransforms:(1)marginal This section first introduces the basic formulation of the uniformization,Ui ,basedonthecumulativedensityfunction proposed method, and then analyzes the properties of differ- (k) of the marginal PDF, and (2) Gaussianization of a uniform entiability, invertibility, and convergence. Finally, we discuss variable,G(u),basedontheinverseofthecumulativedensity on the role of the rotation matrix used in the scheme. function of a univariate Gaussian: Ψi =G(cid:12)Ui , where: (k) (k) A. Iterative Gaussianization based on arbitrary rotations (cid:90) x(k) Accordingtotheabovereasoning,weproposethefollowing u=U(ik)(x(ik)) = i pi(x(cid:48)i(k))dx(cid:48)i(k) (6) −∞ classofRotation-basedIterativeGuassianization(RBIG)algo- (cid:90) xi rithms:givenad-dimensionalrandomvariablex(0),following G−1(x ) = g(x(cid:48))dx(cid:48) (7) i i i an unknown PDF, p(x(0)), in each iteration k, a two-step −∞ processing is performed: and g(x ) is just a univariate Gaussian. Figure 1 shows an i G :x(k+1) =R ·Ψ (x(k)) (5) exampleofthemarginalGaussianizationofaone-dimensional (k) (k) variable x . i where Ψ(k) is the marginal Gaussianization of each dimen- Onedimensionaldensityestimationisanissuebyitself,and sion of x(k) for the corresponding iteration, and R is it has been widely studied [4], [33]. The selection of the most (k) a generic rotation matrix for the marginally Gaussianized convenient density estimation procedure depends on the par- variable Ψ (x(k)). ticular problem and, of course, the univariate Gaussianization (k) Thefreedominchoosingtherotationsisconsistentwiththe stepintheproposedalgorithmcouldbenefitfromtheextensive intuitive fact that there is an infinite number of ways to twist literatureontheissue.Inourcase,wetakeapracticalapproach a PDF in order to turn it into a unit covariance Gaussian. and no particular model is assumed for the marginal variables In principle, any of these choices is equally useful for our to keep the method as general as possible. Accordingly, purpose, i.e. estimating the PDF in the original domain using the univariate Gaussianization transforms are computed from Eq.(1).Notethatwhenusingdifferentrotations,thequalitative the cumulative histograms. Of course, alternative analytical meaningofthesameregionofthecorrespondingGaussianized approximations could be introduced at the cost of making the domain will be different. As a result, in order to work in the model more rigid. On the positive side, parametric models Gaussianized domain, one has to take into account the value may imply better data regularization and avoid overfitting. of the point-dependent Jacobian. Incidentally, this is also the However,exploringtheeffectofalternativedensityestimators caseinthePPapproach,andmoregenerally,inanynon-linear will not be analyzed here. approach. However, the interpretation of the Gaussianized Let us consider now the issue of invertibility. By simple domain is not an issue when working in the original domain. manipulationof(5),itcanbeshownthattheinversetransform Finally, it is important to note that the method just depends is given by: on univariate (marginal) PDF estimations. Therefore, it does not suffer from the curse of dimensionality. G−1 :x(k) =Ψ−(k1)(R(cid:62)(k)·x(k+1)). (8) The rotation R is not a problem for invertibility since the (k) B. Invertibility and differentiation inverse is just the transpose, R−1 = R(cid:62) . However, the key (k) (k) TheconsideredclassofGaussianizationtransformsisdiffer- to ensure transform inversion is the invertibility of Ψ . This (k) entiable and invertible. Differentiability, allows us to estimate is trivially ensured when the support of each marginal PDF is 4 2011(cid:13)cIEEETRANSACTIONSONNEURALNETWORKS connected,thatis,therearenoholes(zeroprobabilityregions) Property 3.2 (Redundancy reduction): Given a marginally in the support. In this way all the marginal CDFs are strictly Gaussianized variable, Ψ(x), any rotation reduces the redun- monotonic and hence invertible. Note that the existence of dancy among coefficients, holes in the support of the joint PDF is not a problem as long ∆I =I(Ψ(x))−I(RΨ(x))≥0, ∀R (13) as it gives rise to marginal PDFs with a connected support. Problems in inversion will appear only when the joint PDF Note that this property also implies that the combination of gives rise to clusters that are so distant that their projections marginalGaussianizationandrotationgivesrisetoredundancy onto the axes do not overlap. However, in such a situation, reduction since I(Ψ(x))=I(x). it may make more qualitative sense to consider that distinct Proof: Using Eq. (2) on both I(Ψ(x)) and I(RΨ(x)), clusters come from different sources and learn each one with the redundancy reduction is: a different Gaussianization transform. ∆I =J(Ψ(x))−J (Ψ(x))−J(RΨ(x))+J (RΨ(x)). TheJacobianoftheseriesofK iterationsisjusttheproduct m m of the corresponding Jacobian in each iteration: Sincenegentropyisrotationinvariantandthemarginalnegen- ∇ G =(cid:81)K R ·∇ Ψ (9) tropy of a marginally Gaussianized variable is zero, x k=1 (k) x(k) (k) ∆I =J (RΨ(x))≥0, ∀R Marginal Gaussianization, Ψ , is a dimension-wise trans- m (k) form, whose Jacobian is the diagonal matrix, The above properties suggest the convergence of the pro- ∂Ψ1(k)  posed Gaussianization method. Property 3.1 (Eq. (12)) en- ··· 0 ∂x(k)  sures that the distance between the PDF of the transformed ∇x(k)Ψ(k) = ...1 ... ...  (10) visarrieadbulecetdo ainzeeraochmietaenratuinoint.cPorvoapriearntyce3m.2u(ltEivqa.ri(a1t3e))Geanussusiraens    ∂Ψd  that redundancy among coefficients is also reduced after each  0 ··· (k) ∂x(k) iteration. According to this the distance to a Gaussian will d decay to zero for a wide class of rotations. According to the two equalization steps in each marginal Gaussianization, Eq. (7), each element in ∇ Ψ can be x(k) (k) D. On the rotation matrices easily computed by applying the chain rule on u defined in Admissiblerotationsarethosethatchangethesituationafter Eq. (6): marginal Gaussianization in such a way that J is increased. m ∂Ψi ∂G ∂u (cid:18)∂G−1(cid:19)−1 Using different rotation matrices gives rise to different prop- (k) = = p (x(k)) ∂x(k) ∂u∂x(k) ∂xi i i erties of the algorithm. i i The above Properties 3.1 and 3.2 provide some intuition on = g(Ψi(k)(x(ik)))−1pi(x(ik)) (11) the suitable class of rotations. By using (12) and (13) in the sequence (5), one readily obtains the relations: Again, the differentiable nature of the considered Gaussian- ization is independent from the selected rotations R(k). ∆J(k) =Jm(x(k))=∆I(k−1), (14) and thus, interestingly, the amount of negentropy reduction C. Convergence properties (the convergence rate) at some iteration k will be determined Here we prove two general properties of random variables, by the amount of redundancy reduction obtained in the pre- which are useful in the contexts of PDF description and vious iteration, k−1. Since dependence can be analyzed in redundancy reduction. terms of correlation and non-Gaussianity [26], the intuitive Property 3.1 (Negentropy reduction): Marginal Gaussian- candidates for R include orthonormal ICA, hereafter simply ization reduces the negentropy and this is not modified by referred to as ICA, which maximizes the redundancy reduc- any posterior rotation: tion; and PCA, which removes correlation. Random rotations (RND) will be considered here as an extreme case to point ∆J =J(x)−J(RΨ(x))≥0, ∀R (12) out that looking for interesting projections is not critical to Proof: Using Eq. (2), the negentropy reduction due to achieve convergence. Note that other rotations are possible, marginal Gaussianization followed by a rotation is: for instance, a quite sensible choice would be randomly selecting projections that uniformly recover the surface of ∆J =J(x)−J(RΨ(x))=J(x)−J(Ψ(x)) an hypersphere [34]. Other possibilities include extension to since N(0,I) is rotation invariant. Therefore, complex variables [35]. As an illustration, Table I summarizes the main character- ∆J = I(x)+J (x)−I(Ψ(x))−J (Ψ(x)) m m istics of the method when using ICA, PCA and RND. The Since the multi-information is invariant under dimension-wise table analyzes the closed-form nature of each rotation, the transforms [27] (such as Ψ), and the marginal negentropy of theoretical convergence of the method, the convergence rate a marginally Gaussianized variable is zero, (negentropy reduction per iteration), and the computational costofeachrotation.SectionV-Aisdevotedtotheexperimen- ∆J =J (x)≥0, ∀R m tal confirmation of the reported characteristics of convergence presented here. LAPARRAETAL.:ITERATIVEGAUSSIANIZATION:FROMICATORANDOMROTATIONS 5 Using ICA guarantees the theoretical convergence of the computation time of the rotation is negligible compared to Gaussianization process since it seeks for the maximally non- that of the marginal Gaussianization. Gaussian marginal PDFs. Therefore, the negentropy reduction ∆J (Eq. (12)) is always strictly positive, except for the case IV. RELATIONTOOTHERMETHODS that the Gaussian PDF has been finally achieved. This is In this section we discuss the relation of RBIG to pre- consistentwithpreviouslyreportedresults[17].Moreover,the viously reported Gaussianization methods, including iterative convergence rate is optimal for ICA since it gives rise to PP-like techniques [15]–[17] and direct approaches suited for the maximum J (x) (indicated in Table I with ‘Max ∆J’). m particular PDFs [19], [21], [40]. Additionally, relations to However, the main problem of using ICA as the rotation other machine learning tools, such as Support Vector Domain matrixisthatithasnoclosed-formsolution,soICAalgorithms Description [22] and deep neural networks [23] are also typically resort to iterative procedures with either difficulties considered. in convergence or high computational load. Using PCA leads to sub-optimal convergence rate because itremovessecond-orderredundancy(indicatedinTableIwith A. Iterative Projection Pursuit Gaussianization ‘∆J =2ndorder’),butitdoesnotmaximizethemarginalnon- As stated above, the aim of Projection Pursuit (PP) tech- Gaussianity Jm(x). Using PCA guarantees the convergence niques [15], [16] is looking for interesting linear projections for every input PDF except for one singular case: consider a accordingtosomeprojectionindexmeasuringinterestingness, variable x(k) which is not Gaussian but all its marginal PDFs and after, this interestingness is captured by removing it are univariate Gaussian and with a unit covariance matrix. In through the appropriate marginal equalization, thus making a this case, ∆J(k+1) = Jm(x(k)) = 0, i.e. no approximation step from structure to disorder. When interestingness or struc- to the Gaussian in negentropy terms is obtained in the next ture is defined by departure from disorder, non-Gaussianity iteration. Besides, since Ψ (x(k)) = x(k), the next PCA, or negentropy, PP naturally leads to iterative application of (k+1) R , will be the identity matrix, thus x(k+1) = x(k): non-orthogonal ICA transforms followed by marginal Gaus- (k+1) as a result, the algorithm may get stuck into a negentropy sianization, as in [17]: local minimum. In our experience, this undesired effect never G :x(k+1) =Ψ (RICA·x(k)) (15) happened in real datasets. On the other hand, advantages of (k) using PCA is that the solution is closed-form, very fast, and As stated in Section II, this is Approach 1 to the Gaussian even though the convergence rate is lower than for ICA, the goal. Unlike PP, RBIG aims at the Gaussian goal following solution is achieved in a fraction of the time. Approach 2. The differences between (15) and (5) (reverse Using RND transforms guarantees the theoretical conver- order between the multivariate and the univariate transforms) gence of the method since random rotations ensure that, even suggest the different qualitative weight given to each coun- intheaboveconsideredsingularcase,thealgorithmwillnotbe terpart. While PP gives rise to an ordered transition from stuck into this particular non-Gaussian solution. On the con- structure to disorder2, RBIG follows a disordered transition trary,iftheachievedmarginalnon-Gaussianityiszeroafteran to disorder. infinite number of random rotations, it is because the desired Gaussian solution has been finally achieved (Cramer-Wold B. Direct (single-iteration) Gaussianization algorithms Theorem[39]).Inpractice,theabovepropertyofRNDcanbe usedasawaytocheckconvergencewhenusingotherrotations Direct (non-iterative) Gaussianization approaches are pos- (e.g. PCA): when the zero marginal non-Gaussianity situation sible if the method has to be applied to restricted classes is achieved, a useful safety check consists of including RND- of PDFs, for example: (1) PDFs that can be marginally based iterations. In the RND case, the convergence rate is Gaussianized in the appropriate axes [19], or (2) elliptically clearly sub-optimal, yet non-negative (∆J ≥ 0): the amount symmetric PDFs so that the final Gaussian can be achieved ofnegentropyreductionmaytakeanyvaluebetweenzeroand byequalizingthelength(norm)ofthewhitenedsamples[21], themaximumachievedbyICA.However,themethodismuch [40]. faster in practice: even though it may take more iterations Themethodproposedin[19]isusefulwhencombinedwith to converge, the cost of each transform does not depend on toolsthatcanidentifymarginallyGaussianizablecomponents, the number of samples. The rotation matrix can be computed somewhat related to ICA transforms. Nevertheless, the use of by fast orthonormalization techniques [38]. In this case, the alternative transformations is still an open issue. Erdogmus et al. proposed PCA, vector quantization or clustering as alternatives to ICA in order to find the most potentially TABLEI ‘Gaussianizable’ components. In this sense, the method could PROPERTIESOFTHEGAUSSIANIZATIONMETHODFORDIFFERENT be seen as a particular case of PP in that it only uses one ROTATIONS(SEECOMMENTSINTHETEXT). iteration: first finding the most appropriate representation and Closed Theoretical Convergence CPUcost† then using marginal Gaussianization. Elliptically symmetric Rotation -form convergence rate [36]–[38] √ PDFs constitute a relevant class of PDFs in image processing ICA × Max∆J O(2md(d+1)n) √ √ PCA ∆J =2ndorder O(d2(d+1)n) √ √ RND ∆J ≥0 O(d3) 2InPPthestructureoftheunknownPDFintheinputdomainisprogres- sivelyremovedineachiterationstartingfromthemostrelevantprojectionand † Computationalcostconsidersnsamplesofdimensiond.Thecostforthe continuingbythesecondone,andsoon,untiltotaldisorder(Gaussianity)is ICAtransformisthatofFastICArunningmiterations. achieved. 6 2011(cid:13)cIEEETRANSACTIONSONNEURALNETWORKS Image Image PDF RG RBIG C. Relation to Support Vector Domain Description (0.34) (0.04) (0.0006) The Support Vector Domain Description (SVDD) is a one- class classification method that finds a minimum volume sphereinakernelfeaturespacethatcontains1−ν fractionof the target training samples [22]. The method tries to find the (0.59) (0.034) (0.0002) transformation (implicit in the kernel function) that maps the target data into a hypersphere. The proposed RBIG method and the SVDD method are conceptually similar due to their apparent geometrical similarity. However, RBIG and SVDD represent two different approaches to the one-class classifi- (0.066) (0.052) (0.0001) cation problem: PDF estimation versus separation boundary estimation. RBIG for one-class problems may be naively seen as if test samples were transformed and classified as target if lying inside the sphere containing 1−ν fraction of the learned Gaussian distribution. According to this interpretation, both methods reduce to the computation Fig.2. Gaussianizationofpairsofneighborpixelsfromdifferentimageswith of spherical boundaries in different feature spaces. However, RGandRBIG:naturalimage(toprow),remotesensingLandsatchannelinthe optical range(middle row)and intensityof aERS2 synthetic apertureradar this is not true in the RBIG case: note that the value of (SAR)image(bottomrow).ContourplotsshowthePDFsinthecorresponding the RBIG Jacobian is not the same at every location in domains.Theestimatedmutualinformation(inbits)isgiveninparenthesis. the Gaussianized domain. Therefore, the optimal boundary to reject a ν fraction of the training data is not necessarily a sphere in the Gaussianized domain. In the case of the SVDD, though,byusinganisotropicRBFkernel,alldirectionsinthe applications since this kind of functions is an accurate model kernel feature spaces are treated in the same way. of natural images (e.g. Gaussian Scale Mixtures [10] and related models [20] share this symmetry). Radial Gaussian- ization (RG) was specifically developed to deal with these D. Relation to Deep Neural Networks particular kind of models [21]. This transform consists of a nonlinear function that acts radially, equalizing the histogram RBIG is essentially an iterated sequence of two operations: of the magnitude (energy) of the data to obtain the histogram non-linear dimension-wise squashing functions and linear ofthemagnitudeofaGaussian.Othermethodshaveexploited transforms. Intuitively, these are the same processing blocks this kind of transformation to generalize it to Lp symmetric used in a feedforward neural network (linear transform plus distributions [40]. Obviously, elliptical symmetry is a fair sigmoid-shapedfunctionineachhiddenlayer).Therefore,one assumption for natural images, but it may not be appropriate could see each iteration as one hidden layer processing of for other problems. Even in the image context, particular the data, and thus argue that complex (highly non-Gaussian) images may not strictly follow distributions with elliptical tasksshouldrequiremorehiddenlayers(iterations).Thisview symmetry,thereforeifRG-liketransformsareappliedtothese is in line with the field of deep learning in neural networks images, they will give rise to non-Gaussianized data. [23],whichconsistsoflearningamodelwithseverallayersof nonlinearmappings.Thefieldisveryactivenowadaysbecause Figure2showsthiseffectinthreetypesofacquiredimages: some tasks are highly nonlinear and require accurate design (1) a standard grayscale image, i.e. a typical example of a of processing steps of different complexity. Note, that it may natural photographic image, (2) a band (in the visible range) appear counterintuitive the fact that full Gaussianization of a of a remote sensing multispectral image acquired by the dataset is eventually achieved with a large enough number of Landsatsensor,and(3)aERS2syntheticapertureradar(SAR) iterations, thus leading to overfitting in the case of a neural intensityimageforthesamescene(ofcourseoutofthevisible network with such number of layers. Nevertheless, note that range). In these illustrative examples, RG and RBIG were capacity control also applies in RBIG: we have observed trainedwiththedatadistributionofpairsofneighborpixelsfor that early-stopping criteria must be applied to allow good each image, and RBIG was implemented using PCA rotations generalization properties. In this setting, one can see early according to the results in Section V-A. Both RG and RBIG stopping in the Gaussianization method as a form of model strongly reduce the mutual-information of pairs of neighbor regularization. This is certainly an interesting research line to pixels (see the mutual information values, in bits), but it is be pursued in the future. noticeablethatRGismoreeffective,higherI reduction,inthe naturalimagecases(photographicandvisiblechannelimages), inwhichtheassumptionofellipticallysymmetricPDFismore Finally, we would like to note that it does not escape our reliable. However, it obviously fails when considering non- notice that the exploitation of the RBIG framework in the natural (radar) images, far from the visible range (I is not previous contexts might eventually be helpful in designing significantlyreduced).Theproposedmethodismorerobustto new algorithms or helping understanding them from different these changes in the underlying PDF because no assumption theoretical perspectives. This is of course out of the scope of is made. this paper. LAPARRAETAL.:ITERATIVEGAUSSIANIZATION:FROMICATORANDOMROTATIONS 7 k 1 5 10 20 35 RND PCA Fig. 4. Cumulative negentropy reduction (top) and multivariate Gaussian ICA significance test (bottom) for each iteration in 2D (left) and 4D (right) dimensional synthetic problem. Average and standard deviation results from 5realizationsisshown. Fig.3. Scatterplots ofa2Ddataindifferentiterationsfortheconsidered rotationmatrices:RND(top),PCA(middle)andICA(bottom). TABLEII AVERAGE(±STD.DEV.)CONVERGENCERESULTS. RND PCA ICA V. EXPERIMENTALRESULTS Dim. iterations time[s] iterations time[s] iterations time[s] 2 14±3 0.01±0.01 7±3 0.005±0.002 3±1 6±5 This section shows the capabilities of the proposed RBIG 4 44±6 0.06±0.01 33±6 0.05±0.01 11±1 564±223 methods in some illustrative examples. We start by experi- 6 68±7 0.17±0.01 43±12 0.1±0 11±2 966±373 8 92±4 0.3±0.1 54±23 0.2±0 16±1 1905±534 mentally analyzing the convergence of the method depend- 10 106±10 0.4±0 58±25 0.3±0.1 19±1 2774±775 ing on the rotation matrix in a controlled toy dataset, and 12 118±10 0.5±0.2 44±5 0.2±0.1 21±2 3619±323 14 130±8 0.7±0.1 52±21 0.4±0.1 19±1 4296±328 give useful criteria for early-stopping. Then, method’s perfor- 16 139±10 0.7±0 73±36 0.4±0.2 22±1 4603±932 mance is illustrated for mutual information estimation, image synthesis, classification and denoising. In each application, results are compared to standard methods in the particular field. A documented Matlab implementation is available at method converges to a multivariate Gaussian independently http://www.uv.es/vista/vistavalencia/RBIG.htm. of the rotation matrix; (2) ICA requires a less number of iterations to converge, but it is closely followed by PCA; (3) random rotations take a higher number of iterations to A. Method convergence and early-stopping converge and show high-variance in the earlier iterations; and TheRBIGmethodisanalyzedhereintermsofconvergence (4) convergence in cumulative negentropy is consistent with rateandcomputationalcostfordifferentrotations:orthonormal the parametric estimator in [41] which, in turn, confirms the ICA, PCA and RND. Synthetic data of varying dimensions analysis in Table I. (d = 2,...,16) was generated by first sampling from a Despitethepreviousconclusions,andaspointedoutbefore, uniform distribution hypercube and then applying a rotation in practical applications, it is not the length of the path to the transform. This way we can compute the ground-truth negen- Gaussiangoalwhatmatters,butthetimerequiredtocomplete tropiesoftheinitialdistributions,andestimatethereductionin this path. Table II compares the number of iterations for negentropies in every iteration by estimating the difference in appropriateconvergenceandtheCPUtimeof5realizationsof marginal negentropies, cf. Eq. (13). A total of 10000 samples RBIG with different matrix rotations (RND, PCA and ICA) was used for the methods, and we show average and standard in several dimensions. While, in general, CPU time results deviation results for 5 independent random realizations. are obviously implementation dependent, note that results in Two-dimensionalscatterplotsinFigure3qualitativelyshow Table II are fairly consistent with the computational burden that different rotation matrices give rise to different solutions per iteration shown in Table I since each ICA computation is ineachiterationbut,afterasufficientnumberofiterations,all an iterative procedure itself which needs m iterations. of them transform the data into a Gaussian independently of The use of ICA rotations critically increases the conver- the rotation matrix. gence time. This effect is more noticeable as the dimension RBIG convergence rates are illustrated in Fig. 4. Top plots increases, thus making the use of ICA computationally unfea- show the negentropy reduction for the different rotations as a siblewhenthenumberofdimensionsismoderateorhigh.The function of the number of iterations and data dimension. We useofPCAinRBIGisconsequentlyagoodtrade-offbetween alsogivetheactualnegentropyestimatedfromthesamples,is Gaussianizationerrorandcomputationalcostifthenumberof an univariate population estimate since Eq. (12) can be used. iterationsisproperlychosen.Anearly-stoppingcriterioncould Successful convergence is obtained when the accumulated bebasedontheevolutionofthecumulativenegentropyreduc- reduction in negentropy tends to the actual negentropy value tion, or of a multivariate test of Gaussianity such as the one (cyan line). Discrepancies are due to the accumulation of usedhere[41].Botharesensiblestrategiesforearly-stopping. computational errors in the negentropy reduction estimation According to the observed performance, we restrict ourselves in each iteration. to the use of PCA as the rotation matrix in the experiments Bottom plots in Fig. 4 give the result of the multivariate hereafter. Note that by using PCA, the algorithm might not Gaussianity test in [41]: when the outcome of the test is converge in a singular situation (see Section III-D). However, 1, it means accepting the hypothesis of multidimensional we checked that such singular situation never happened by Gaussianity. Several conclusions can be extracted: (1) the jointly using both criteria in each iteration. 8 2011(cid:13)cIEEETRANSACTIONSONNEURALNETWORKS TABLEIII Original data Gaussianized data Synthesized data AVERAGE(±STD.DEV.)MULTI-INFORMATION(INBITS)FORTHE DIFFERENTESTIMATORSIN2DPROBLEMS. DIST EG GG UU RBIG 0.49 ± 0.01 1.38 ± 0.004 0.36 ± 0.03 NE 0.35 ± 0.02 1.35 ± 0.006 0.39 ± 0.002 RE 0.32 ± 0.01 1.29 ± 0.004 0.30 ± 0.002 Actual 0.51 1.38 0.45 TABLEIV MULTI-INFORMATION(INBITS)WITHRBIGINDIFFERENT d-DIMENSIONALPROBLEMS. Dim. EG GG UU d RBIG Actual RBIG Actual RBIG Actual 3 1.12±0.03 1.07 1.91±0.01 1.9 1.6±0.1 1.6 4 5±0.1 5.04 1.88±0.02 1.86 2.2±0.1 2.2 5 4.7±0.1 4.82 1.77±0.02 1.75 2.7±0.1 2.73 6 7.8±0.1 7.9 2.11±0.01 2.08 3.5±0.1 3.72 7 6.2±0.1 6.33 2.68±0.03 2.65 3.6±0.1 3.92 8 8.1±0.1 8.19 2.72±0.02 2.68 4.1±0.1 4.29 9 9.5±0.1 9.6 3.22±0.02 3.18 5.3±0.1 5.69 10 12.7±0.1 13.3 3.45±0.03 3.4 5.8±0.2 6.24 B. Multi-information Estimation As previously shown, RBIG can be used to estimate the Fig.5. ToydataexamplessynthesizedusingRBIG. negentropy, and therefore could be used to compute multi- information(I)ofhighdimensionaldata(Eq.(2)).Essentially, random Gaussian samples in the transformed domain inverted one learns the sequence of transforms to Gaussianize a given back to the original domain using G−1. Two examples are dataset, and the I estimate reduces to compute the cumulative given here to illustrate the capabilities of the method. ∆I since, at convergence, full independence is supposedly 1) Toydata: Figure5showsexamplesof2Dnon-Gaussian achieved. We illustrate the ability of RBIG in this context by distributions(leftcolumn)transformedintoaGaussian(center estimating multi-information in three different synthetic dis- column). The right column was obtained sampling data from tributions with known I: uniform distribution (UU), Gaussian a zero mean unit covariance Gaussian and inverting back the distribution(GG),andamarginallycomposedexponentialand transform. This example visually illustrates that synthesized Gaussian distribution (EG). An arbitrary rotation was applied data approximately follow the original PDF. ineachcasetoobtainnon-zeromulti-information.Inallcases, 2) Face synthesis: In this experiment, 2,500 face images we used 10,000 samples and repeated the experiments for 10 were extracted from [43], eye-centered, cropped to have the realizations. Two kinds of experiments were performed: same dimensions, mean and variance adjusted, and resized to • A 2D experiment, where RBIG results can be compared 17×15pixels.Imageswerethenreshapedto255-dimensional to the results of naive (histogram-based) mutual infor- vectors,andGaussianizedwithRGandRBIG.Figure6shows mation estimates (NE), and to previously reported 2D illustrativeexamplesoforiginalandsynthesizedfaceswithRG estimatessuchastheRudyestimate(RE)[42](seeTable and RBIG. III). Note that both methods achieve good visual qualitative • A set of d-dimensional experiments, where RBIG results performance.Inordertoassessperformancequantitatively,we are compared to actual values (see Table IV). compared 200 actual and synthesized images using the inner TableIIIshowstheresults(inbits)forthemutualinformation productasameasureoflocalsimilarity.Weaveragedthissim- estimation in the 2D experiment to standard approaches. The ilarity measure over 300 realizations and show the histograms ground-truth result is also given for comparison purposes. for RG and RBIG. Results suggest that the distribution of the For Gaussian and exponential-Gaussian data distributions, samples generated with RBIG is more realistic (similar to the RBIG outperforms the rest of methods, but when data are original dataset) than the obtained with RG. marginally uniform, NE yields better estimates. Table IV extends the previous results to multidimensional cases, and D. One-class Classification compares RBIG to the actual I. Good results are obtained in all cases. Absolute errors slightly increase with data dimen- In this experiment, we assess the performance of the RBIG sionality. method as one-class classifier. Performace is illustrated in the challenging problem of detecting urban areas from mul- tispectral and SAR images. The ground-truth data for the C. Data Synthesis images used in this section were collected in the Urban RBIG obtains an invertible Gaussianization transform that Expansion Monitoring (UrbEx) ESA-ESRIN DUP project3 canbeusedtogenerate(orsynthesize)samples.Theapproach [44]. The considered test sites were the cities of Rome and is simple: the transform G is learned from the available training data, and then synthesized samples are obtained from 3http://dup.esrin.esa.int/ionia/projects/summaryp30.asp LAPARRAETAL.:ITERATIVEGAUSSIANIZATION:FROMICATORANDOMROTATIONS 9 Original RG RBIG Fig.6. Exampleofreal(top)andsynthesizedfaceswithRG(middle)andRBIG(bottom). Fig.7. Histogramofthesimilarity(innerproduct)betweenthedistributionof originalandsynthesizedfaceimagesfor300realizations.Forreference,the averageimageenergy(averageinnerproduct)intheoriginalsetis1.81·104. Naples, Italy, for two acquisitions dates (1995 and 1999). The available features were the seven Landsat bands, two SAR backscattering intensities (0–35 days), and the SAR interferometric coherence. We also used a spatial version of the coherence specially designed to increase the urban areas discrimination[44].Afterthispreprocessing,allfeatureswere stacked at a pixel level, and each feature was standardized. Fig.8. Overallaccuracy(left)andkappastatistic,κ(right)forRBIG(solid We compared the RBIG classifier based on the estimated line)andSVDD(dashedline)indifferentscenes:Naples1995(top),Naples PDF for urban areas with the SVDD classifier [22]. We used 1999(center)andRome1995(bottom). the RBF kernel for the SVDD whose width was varied in the range σ ∈ [10−2,...,102]. The fraction rejection parameter was varied in ν ∈ [10−2,0.5] for both methods. The optimal signatures. Results show that SVDD behavior is similar to the parameters were selected through 3-fold cross-validation in proposed method for small size training sets. This is because the training set optimizing the κ statistic [45]. Training sets more target samples are needed by the RBIG for an accurate of different size for the target class were used in the range PDFestimation.However,formoderateandlargetrainingsets [500,2500].Weassumedascarceknowledgeofthenon-target the proposed method substantially outperforms SVDD. Note class: 10 outlier examples were used in all cases. The test set thattrainingsizerequirementsofRBIGarenottoodemanding: was constituted by 10,000 pixels of each considered image. using 750 samples in a 10-dimensional problem is enough for Trainingandtestsampleswererandomlytakenfromthewhole RBIG to outperform SVDD when very little is known about spatialextentofeachimage.Theexperimentwasrepeatedfor the non-target class. 10 different random realizations in the three considered test Figure9showstheclassificationmapsfortherepresentative sites. Naples95 scene for SVDD and RBIG. Note that RBIG better Figure 8 shows the estimated κ statistic and the overall rejects the ‘non-urban’ areas (in black). This may be because accuracy (OA) in the test set achieved by SVDD and RBIG SVDD training with few non-target data gives rise to a too in the three images. The κ scores are relatively small because broad boundary. As a result, too many pixels are identified samples were taken from a large spatial area thus giving rise as belonging to the target class (in white). Another relevant to a challenging problem due to the variance of the spectral observation is the noise in neighboring pixels, which may 10 2011(cid:13)cIEEETRANSACTIONSONNEURALNETWORKS GT SVDD(0.67,0.32) RBIG(0.79,0.42) Fig. 9. Ground-truth (GT) and classification maps obtained with SVDD and RBIG for the Naples 1995 scene. The white points represent urban area and theblackpointsrepresentnon-urbanarea.Thecorrespondingoverallaccuracyandκ-statisticaregiveninparenthesis. come from the fact that no spatial information was used. This points from the neighborhood of each wavelet coefficient problem could be easily alleviated by imposing some post- by generating samples with the PDF of the noise model classificationsmoothnessconstraintorbyincorporatingspatial (p(x |x)), and evaluated the probability for each sample with n texture features. the PDF obtained in the training step p(x). The estimated coefficient xˆ is obtained as the expected value over the 8000 samples of the posterior PDF. Obtaining the expected value is E. Image Denoising equivalent to using the L norm [51]. Note that the classical Image denoising tackles the problem of estimating the 2 hard-thresholding(HT)andsoft-thresholding(ST)results[46] underlying image, x, from a noisy observation, x , assuming n areausefulreferencesincetheycanbeinterpretedassolutions an additive degradation model: x = x+n. Many methods n to the same problem with a marginal Laplacian image model haveexploitedtheBayesianframeworktothisend[10],[46]– and L and L norms respectively [47]. [48]: 1 2 (cid:26)(cid:90) (cid:27) Figure 10 shows the denoising results for the ‘Barbara’ xˆ =argmin L(x,x∗)p(x|xn)dx , (16) image corrupted with Gaussian noise of σn2 = 100 using x∗ marginal models (HT and ST), and using a RBIG as the PDF wherex∗ isthecandidateimage,L(x,x∗)isthecostfunction, estimator. Accuracy of the results is measured in Euclidean andp(x|xn)istheposteriorprobabilityoftheoriginalsample terms (RMSE), and using a perceptually meaningful image xgiventhenoisysamplexn.Thislasttermplaysanimportant quality metric such as the Structural Similarity Index (SSIM) role since it can be decomposed (using the Bayes rule) as [52]. Note that RBIG method obtains better results (numeri- cally and visually) than the classical methods due to the more p(x|x )=Z−1p(x |x)p(x), (17) n n accurate PDF estimation. whereZ−1isanormalizationterm,p(x |x)isthenoisemodel n (probability of the noisy sample given the original one), and VI. CONCLUSIONS p(x) is the prior (marginal) sample model. In this work, we proposed an alternative solution to the Note that, in this framework, the inclusion of a feasible PDF estimation problem by using a family of Rotation-based image model, p(x), is critical in order to obtain a good Iterative Gaussianization (RBIG) transforms. The proposed estimationoftheoriginalimage.Imagesaremultidimensional procedure looks for differentiable transforms to a Gaussian signals whose PDF p(x) is hard to estimate with traditional sothattheunknownPDFcanbecomputedatanypointofthe methods. The conventional approach consists of using para- original domain using the Jacobian of the transform. metric models to be plugged into Eq. (17) in such a way that The RBIG transform consists of the iterative application of theproblemcanbesolvedanalytically.However,mathematical univariatemarginalGaussianizationfollowedbyarotation.We convenience leads to the use of too rigid image models. Here show that a wide class of orthonormal transforms (including we use RBIG in order to estimate the probability model of trivial random rotations) is well suited to Gaussianization natural images p(x). purposes.Thefreedomtochoosethemostconvenientrotation In this illustrative example, we use the L -norm as cost 2 function, L(x,x∗) = ||x−x∗|| , and an additive Gaussian isthedifferencewithformallysimilartechniques,suchasPro- 2 noise model, p(x |x)=N(0,σ2I). We estimated p(x) using jection Pursuit, focused on looking for interesting projections n n (which is an intrinsically more difficult problem). In this way, 100 achromatic images of size 256 × 256 extracted from the McGill Calibrated Colour Image Database4. To do this, here we propose to shift the focus from ICA to a wider class of rotations since interesting projections as found by ICA images were transformed using orthonormal QMF wavelet are not critical to solve the PDF estimation problem in the domainwithfourfrequencyscales[49],antheneachsubband original domain. The suitability of multiple rotations to solve was converted to patches in order to obtain different PDF the PDF estimation problem may help to revive the interest models for each subband according to well-known properties of classical iterative Gaussianization in practical applications. of natural images in wavelet domains [12], [50]. In order to As an illustration, we showed promising results in a number evaluate Eq. (16), we sampled the posterior PDF at 8,000 of multidimensional problems such as image synthesis, clas- 4http://tabby.vision.mcgill.ca/ sification, denoising, and multi-information estimation.

See more

The list of books you might like

Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.