Table Of Content2011(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.
