Table Of ContentMon.Not.R.Astron.Soc.000,000–000 (0000) Printed26January2013 (MNLATEXstylefilev2.2)
Constraints on f from Wilkinson Anisotropy Probe
NL
7-year data using a neural network classifier
B. Casaponsa,1 2⋆ M. Bridges3, A. Curto1, R.B. Barreiro1, M.P. Hobson3,
E.Mart´ınez-Gonz´alez1
2
1 1 Institutode F´ısica de Cantabria, CSIC-Universidad de Cantabria, Avda. de los Castros s/n, 39005 Santander, Spain.
0 2 Dpto. de F´ısica Moderna, Universidad de Cantabria, Avda. los Castros s/n, 39005 Santander, Spain.
2 3 Astrophysics Group, Cavendish Laboratory, Madingley Road, Cambridge, CB3 0HE, U.K.
n
a
J Accepted —.Received—;inoriginalform26January2013
1
3
ABSTRACT
] We present a multi-class neural network (NN) classifier as a method to measure non-
O Gaussianity, characterised by the local non-linear coupling parameter fNL, in maps
C of the cosmic microwave background (CMB) radiation. The classifier is trained on
. simulatednon-GaussianCMBmaps witharangeofknownfNL valuesby providingit
h with wavelet coefficients of the maps; we consider both the HealPix (HW) wavelet
p
andthesphericalMexicanhatwavelet(SMHW).Whenappliedtosimulatedtestmaps,
-
o the NN classfier produces results in very good agreement with those obtained using
r standardχ2 minimization. The standarddeviations of the fNL estimates for WMAP-
t
s like simulations were σ = 22 and σ = 33 for the SMHW and the HW, respectively,
a whichareextremelyclosetothoseobtainedusingclassicalstatisticalmethodsinCurto
[
et al. and Casaponsa et al. Moreover,the NN classifier does not require the inversion
2 of a large covariance matrix, thus avoiding any need to regularise the matrix when it
v is not directly invertible, and is considerably faster.
6
Key words: methods: data analysis — cosmic microwavebackground
1
1
6
.
5
0 1 INTRODUCTION mator(Bucheret al.2010).Othermethodswhichhavealso
1 beenappliedtonon-GaussianityanalysisincludeMinkowski
1 Artificialintelligencealgorithmsarebeingusedincreasingly functionals(Hikageet al.2008;Natoli et al.2010),wavelet-
: to improve the efficiency of computationally intensive data
v based methods (Cayo´n et al. 2001; Mukherjee & Wang
analysis.Inparticular,neuralnetworks(NN)havebeensuc-
i 2004;Curto et al.2009a,b;Pietrobon2010;Casaponsa et al.
X cessfully appliedtopatternrecognition, classification ofob-
2010), a Bayesian approach (Elsner & Wandelt 2010) and
jectsandparameterestimationinanumberoffields,includ-
r the analysis of the N-dimensional probability density func-
a ing cosmology (see e.g. Auld et al. 2007).
tion (Vielva& Sanz 2010).
Cosmological analysistypicallyinvolvestheuseoflarge
datasets and high precision numerical tools. In particular,
This paper introduces an approach based on a neural
the study of deviations from Gaussianity in the distribu-
network classifier which, after training on third order mo-
tion of temperature anisotropies in the cosmic microwave
ments of wavelet coefficients derived from simulated Gaus-
background (CMB) require very demanding computational
sianandnon-GaussianCMBrealisations,canbeusedtoes-
methods.Thesimplest waytocharacterisesuchadeviation
timate the presence and degree of non-Gaussianity for any
isthroughthirdordermoments,asthesevanishintheGaus-
given data map. We havechosen to estimate the local non-
siancase.ItisnowcommonplaceinCMBanalysistoworkin
linearcouplingparameterf ,whichparameterisesthelocal
NL
spherical harmonic space, where computing the three point
non-Gaussianity as a quadratic term in the primordial cur-
correlation functionorbispectrumcanprovedifficult,orin-
vatureperturbation.Moreprecisely,fNListheamplitudeof
deedimpossible,duetonumericalinstability.Somerecentef-
the corrections at second order of the primordial curvature
fortshavebeenappliedtolessenthecomputationaldemand
perturbations (Salopek & Bond 1990; Gangui et al. 1994;
without reducing efficiency; see for example the KSW bis-
Verdeet al. 2000; Komatsu & Spergel 2001). This type of
pectrumestimator(Komatsu et al.2005),orthebinnedesti-
non-Gaussianity is predicted even in the simplest slow-roll
inflationary scenario, albeit at a very low level f < 1,
NL
whereas a wide range of non-standard inflationary models
⋆ e-mail:[email protected] predict much larger typical fNL values (for a more com-
(cid:13)c 0000RAS
2 B. Casaponsa et al.
plete review see Bartolo et al. (2004),Babich et al. (2004) the wavelet transform, one produces both a high- and low-
andYadav& Wandelt(2010)).Estimating thevalueoff resolution map. The low-resolution map for the HW is ob-
NL
fromagivendatamapusingexistingmethodstypicallyhas tained simply by averaging over 4-pixel blocks, and the
a high computational cost and usually numerical problems high-resolution map isjust theoriginal map minusthelow-
arise (e.g. matrix inversion). As we will show, the use of resolution map. One begins with the original map at reso-
neural networks bypassesthese difficulties. lution J =9 (N =512) and performs successive wavelet
side
In principle, one could use the pixel temperatures in decompositions until resolution J =2 (Nside = 2), thereby
theCMBmapdirectly,oritssphericalharmoniccoefficients, constructing 7 sets of high- and low-resolution maps. Al-
as the inputs to the neural network classifier. Nonetheless, though the original map is fully represented by the 7 high-
we perform a pre-processing step in which we decompose resolution maps plus the low-resolution map at J = 2, in
the temperature maps into their wavelet coefficients, which our analysis we have used all the high- and low-resolution
haveshown themselves to be sensitive to non-Gaussian sig- maps, plus the original map, since this has been shown to
nals(e.g.Curto et al.2009b,2011a;Casaponsa et al.2010). improveresults (see Casaponsa et al. 2010, for details).
In particular, we consider the HealPix wavelet (HW) and Using all these maps, the third order moments of the
a spherical Mexican hat wavelet (SMHW), and compute wavelet coefficients are computed as follows:
tvhailrude-oorfdwerhimchomisenptrsopoofrtthioensealwtaoveflNeLt.coTehffiecnieentwtso,rtkheismtheaenn Sjkl = N1 PNi=l1σwiσ,jwσi,kwi,l, (1)
trained so that when presented with these cubic statistics l j k l
for a new (data) map, it can estimate the fNL value and where wi,j is the ith wavelet coefficient of the map at res-
its error bar. We apply this method to estimate the degree olution j, σj is the dispersion of wi,j, and Nl is the num-
of non-Gaussianity in the Wilkinson microwave anisotropy ber of pixels in the map at resolution l (since one requires
probe (WMAP) 7-year data release. j 6k 6l). Some of these statistics are redundant (linearly
This paper is organized as follows. In Section 2, we dependencyexists between them), so we restrict our analy-
give a brief introduction to the wavelet analysis used. An sis to theset of non-redudant statistics, which gives a total
overview of the type of neural network employed and our of232quantities;thesearethencomputedfornon-Gaussian
training procedure follows in Section 3. In Section 4 we simulations with a range of known valuesof f .
NL
explain the generation of our Gaussian and non-Gaussian Theexpectedvaluesofthesestatistics areproportional
simulations, and thespecificcharacteristics of ourfNL clas- to the non-linear coupling parameter, and they have pre-
sification network. We present the results of applying our viously been used to estimate the best fit f value for
NL
classifier to simulations and to WMAP 7-year data in Sec- the data by weighted least squares parameter estimation
tion 5. Ourconclusions are summarised in Section 6. (Casaponsa et al. 2010). In this case, the dispersion in the
estimated f value for Gaussian simulations and is found
NL
tobe σ(f )=34, which is slightly larger that theoptimal
NL
value. The main advantage of the HW is the computing ef-
2 WAVELETS ficiency;for example, thethird-orderstatistics construction
is 103 times faster than for the KSW bispectrum estimator
Wavelet methods have seen increasing usage in cosmol-
(Komatsu et al.2005)and102 timesfasterthantheSMHW
ogy. This has been particularly marked in CMB non-
(seebelow).Thisprocedure(forboththeHWandSMHW)
Gaussianityanalyses,inwhichcompetitiveresultshavebeen
does,however,includetheestimationandinversionofacor-
obtained using wavelets such as the SMHW (Cayo´n et al.
relation matrix, which can be computationally demanding
2003; Vielva et al. 2004; Cruz et al. 2005; Curto et al.
and,insomecases,closetosingular.Aswewillshowbelow,
2011a), directional spherical wavelets (McEwen et al.
thisstep is avoided with theuse of a NN classifier.
2008), spherical Haar wavelet (SHW) (Tenorio et al. 1999;
Barreiro et al. 2000), and recently the HealPix wavelet
(HW)(Casaponsa et al.2010).Forareviewof waveletsap- 2.2 Spherical Mexican Hat Wavelet
pliedonthesphere,see,forexample,McEwen et al.(2007).
The spherical Mexican hat wavelet (SMHW)
Inessence,decomposingaCMBmapintoitswaveletcoeffi-
(Antoine& Vandergheynst 1998; Mart´ınez-Gonz´alez et al.
cientsallowsonetoseparatethesearchfornon-Gaussianity
2002) has produced competitive results in constraining
ondifferentlength-scales,whileretainingpositionalinforma-
primordial non-Gaussianity (Mukherjee& Wang 2004;
tion.Inthissectionwewillbrieflydiscussthecharacteristics
Curto et al.2009a,b,2011a).Itisacontinuouswaveletthat
ofboththeHWandSMHWanddescribehowweconstruct
has better frequency localization than the HW, although
thestatistics which are used in our analysis.
the computing efficiency is lower. Curto et al. (2011a) use
the SMHW to constrain f with an accuracy equivalent
NL
2.1 HealPix wavelet to the bispectrum estimators (see for example Smith et al.
2009; Fergusson & Shellard 2009; Fergusson et al. 2010;
The HealPix wavelet is very similar to that presented by Komatsu et al. 2011; Bucher et al. 2010). The definition of
Shahram et al. (2007). Casaponsa et al. (2010) haveused a thethird-ordermomentsis thesame asfor theHW.Inthis
revised version of this wavelet and perform the first cos- case, however, there are more inter-scale combinations be-
mological application. In both papers, the central idea is causethescalesinvolvedarenotrestrictedbytheHealPix
the construction of a fast wavelet method adapted to the pixelization. The total number of non-redundant statistics
HealPixpixelization scheme(G´orski et al.2005).TheHW for the SMHW wavelet coefficients is 680. Using the mean
is similar to the SHW in the sense that, at each level of values and covariances of these statistics computed from
(cid:13)c 0000RAS,MNRAS000,000–000
WMAP 7yr local f constraints using neural networks 3
NL
indexirunsoverallinputnodes.Thefunctionsg(1)andg(2)
x y arecalledactivationfunctions.Thenon-linearnatureofg(1)
i k
is a key ingredient in constructing a viable and practically
hj useful network. This non-linear function must be bounded,
smooth and monotonic; we use g(1)(x) = tanhx. For g(2)
wesimplyuseg(2)(x)=x.Thelayoutandnumberofnodes
are collectively termed the architecture of the network. For
a basic introduction to artificial neural networksthe reader
is directed toMacKay (2003).
b b Foragivenarchitecture,theweightswandbiasesbde-
b
b bb b fine the operation of the network and are the quantities we
b b
wish to determine by some training algorithm. We denote
w and bcollectively bya.Astheseparameters varyduring
w ij w
jk training,averywiderange ofnon-linearmappingsbetween
inputs and outputs is possible. In fact, according to a ‘uni-
Figure 1.Schematic ofa3-layerfeed-forwardneuralnetwork. versalapproximation theorem’Leshnoet al. (1993),astan-
dardthree-layerfeed-forward network canapproximateany
continuous function to any degree of accuracy with appro-
non-Gaussian simulations, Curto et al. (2011a) applied a priatelychosenactivationfunctionsandasufficientnumber
χ2 minimisation method to obtain optimal uncertainties of hidden nodes.
on the fNL estimates of σ ≈ 21. However, this method Inourapplication,wewillconstructaclassificationnet-
requires a principal component analysis to deal with the work.Theaimofanyclassificationmethodistoplacemem-
degenerancies present in the covariance matrix. As we will bersofasetintogroupsbasedoninherentpropertiesorfea-
see, this problem is avoided with the use of the multi-class tures of the individuals, given some pre-classified training
neural network classifier. data. Formally, classification can be summarised as finding
aclassifierC :x→C whichmapsanobjectfromsome(typ-
icallymulti-dimensional)featurespacextoitsclassification
label C, which is typically taken as one of {1,...,N} where
3 NEURAL NETWORKS
N is the number of distinct classes. Thus the problem of
Artificial neural networks are a methodology for comput- classification is to partition feature space into regions (not
ing, based on massive parallelism and redundancy, which necessarily contiguous), assigning each region a label corre-
are features also found in animal brains. They consist of sponding to the appropriate classification. In our context,
a number of interconnected processors each of which pro- theaim istoclassify setsof third-orderstatistics ofwavelet
cesses information and passes it to other processors in the coefficients of (possibly) non-Gaussian CMB maps (assem-
network. Well-designed networks are able to ‘learn’ from bled into an input feature vector x) into classes defined by
a set of training data and to make predictions when pre- ranges of f ; this is discussed in more detail below.
NL
sented with new, possibly incomplete, data. These algo- Inbuildingaclassifierusinganeuralnetwork,itiscon-
rithmshavebeensuccessfullyappliedinseveralareas,inpar- venienttoviewtheproblemprobabilistically.Tothisendwe
ticular, we note the following applications in astrophysics: consider a 3-layer MLP (multi-layer percepton) consisting
Storrie-Lombardi et al. (1992); Baccigalupi et al. (2000); of an input layer (xi), a hidden layer (hj), and an output
Vanzella et al.(2004);Auldet al.(2007)andCarballo et al. layer (yi). In classification networks, however, the outputs
(2008). are transformed according to thesoftmax procedure
The basic building block of an ANN is the neuron. In-
eyk
formationispassedasinputstotheneuron,whichprocesses p = , (4)
k eym
themandproducesanoutput.Theoutputistypicallyasim- Pm
ple mathematical function of the inputs. The power of the such that they are all non-negative and sum to unity. In
ANNcomes from assembling many neuronsintoanetwork. this way p can be interpreted as the probability that the
k
The network is able to model very complex behaviour from input feature vector x belongs to the kth class. A suitable
input to output.Weuse a three-layer feed-forward network objective function for the classification problem is then
consisting of a layer of input neurons, a layer of ‘hidden’
neurons and a layer of output neurons. In such an arrange- L(a)=XXt(kl)lnpk(x(l),a), (5)
menteachneuronisreferred toasanode.Figure1showsa l k
schematic design of such a network. where the index l runs over the training dataset D =
The outputs of the hidden layer and the output layer {x(l),t(l)}, in which the target vector t(l) for the network
are related to their inputsas follows:
outputs has unity in the element corresponding to the true
hidden layer: hj =g(1)(fj(1)); fj(1) =Xwj(1i)xi+b(j1),(2) tchlaesnsowfisthheesltthofcehaotuorseevneecttwoorrxk(lp)aarnadmzeetreorsesaelssoewashetroe.mOanxe-
i
imisethisobjectivefunctionasthetrainingprogresses. The
output layer: yk =g(2)(fk(2)); fk(2) =Xwk(2j)hj+b(k2),(3) advantageofthisprobabilistic approach isthat wegain the
j abilitytomakestatistical decisionsontheappropriateclas-
wheretheoutputofthehiddenlayerhandoutputlayeryare sification in very large feature spaces where a direct linear
given for each hidden node j and each output node k. The partition would not be feasible.
(cid:13)c 0000RAS,MNRAS000,000–000
4 B. Casaponsa et al.
One wishes to choose network parameters a so as to (0,0,1). In this example, the output given by the network
maximise the objective function L(a) as the training pro- forsometestinputvectorxwouldbea3-dimensionalvector
gresses. This is, however, a highly non-linear, multi-modal p=(p ,p ,p ),where p =1andp canbeinterpreted
1 2 3 Pk k k
function in many dimensions whose optimisation poses a as the probability that the input vector belongs to class k.
non-trivial problem. We perform this optimisation using Thediscrepancybetweenthetargetsandtheoutputsduring
theMemSyspackage(Gull & Skilling1999).Thisalgorithm trainingcanbemeasuredbythetruepositiverate,whichis
considers the parameters a to have prior probabilities pro- simply the fraction of the training input vectors for which
portional to eαS(a), where S(a) is the positive-negative en- thenetworkassignsthemaximumprobabilitytothecorrect
tropyfunctional(Hobson & Lasenby1998).αistreatedasa class.
hyper-parameter of the prior, and sets the scale over which From the output values p obtained for each map, we
k
variations in a are expected. α is chosen to maximise its definetheestimatorofthelocalnon-Gaussianityparameter
marginal posterior probability whose value is inversely pro- to besimply
portional to the standard deviation of the prior. Thus for
aL(gai)ve+nααS,(ath).eFloorg-epaocshtecrhioorsepnroαbathbeilrietyisisapsrooluptoirotnioaˆnatlhtaot fˆNL =nXclasshfNLikpk (6)
k=1
maximises the posterior. As α varies, the set of solutions aˆ
is called the maximum-entropy trajectory. We wish to find wherehf i isthemean valueof f in thekth class. The
NL k NL
the solution for which L is maximised which occurs at the statisticalpropertiesofthisestimator,namelyitsmeanand
endofthetrajectorywhereα=0.Forpracticalpurposeswe dispersion, determinethe accuracy of themethod.
start at a large value of α and iterate downwards untilα is
sufficiently small so that the posterior is dominated by the
L term. MemSys performs this algorithm using conjugate
4.1 Training data
gradient descent at each step to convergeto themaximum-
entropy trajectory. The required matrix of second deriva- The training input vectors x(l) were generated as fol-
tives of L is approximated using vector routines only, thus lows. We began with a set of 1000 non-Gaussian CMB
circumventing the need for O(N3) operations required for realisations from which aNG and aG were generated by
lm lm
exactcalculations.TheapplicationofMemSystotheprob- Elsner & Wandelt (2009) and normalized to the WMAP
lem of network training allows for the fast efficient train- 7-year concordance model power spectrum generated by
ing of relatively large network structures on large data sets CAMB.Thesea arepubliclyavailable1.Theultimateac-
lm
that would otherwise be difficult to perform in a reason- curacyofthenetworkclassifierisimproved,however,bythe
able time. Moreover the MemSys package also computes inclusion of further training data. Given the finite number
the Bayesian evidence for the model (i.e. network) under of available simulations, we thus created a further set by
consideration, (see for example Jaynes 2003, for a review), rotation of the original maps by 90◦ perpendicular to the
whichprovidesapowerfulmodelselectiontool.Inprinciple, galacticplane.Thisrotationcreatesroughly20percentex-
values of the evidence computed for each possible architec- tra map area based on the original mask; we verified that
tureofthenetwork(andtrainingdata)provideamechanism its inclusion improves the results. Using this procedure we
toselect themost appropriatearchitecture, whichissimply generate a further 1000 non-Gaussian simulations. Of the
theonethatmaximisestheevidence(althoughwewillusea 2000 non-Gaussian maps, 1800 were used for training and
more prosaic method below for decidingon thenetwork ar- theremainder were set aside for testing of thenetworks.
chitecture). The MemSys algorithm is described in greater Foreachnon-Gaussiansimulationusedfortraining,sets
detail in (Gull & Skilling 1999). of alm were then generated with varying fNL using thefol-
lowing prescription
a =aG +f aNG, (7)
lm lm NL lm
4 THE F CLASSIFICATION NETWORK
NL
with 20differentf random valuesbetween−120 and120
Totrainourf classificationnetworkweprovideitwithan NL
NL fortheHWdecomposition andbetween −76and 76for the
ensembleoftrainingdataD={x(l),t(l)}.Thelthinputvec-
SMHWanalysis.Thus,foreachnon-Gaussiansimulation,20
torx(l) containsthethird-orderstatisticsofthewaveletco-
setsofa weregenerated.Hencethetotalnumberofavail-
efficientsofthelth simulatedCMBmap.Theoutputclasses lm
abletrainingdatasetsis36000.Noise-weightedV+W-band
of our network correspond to contiguous ranges of f val-
NL WMAP realizations were then constructed as explained in
ues.Thus,thetargetvectort(l) forthenetworkoutputshas
Curto et al. (2009a) and Casaponsa et al. (2010), and the
zeroeseverywhereexceptforaunitentryintheelementcor-
KQ75 mask was then applied, which covers roughly 29%
responding to the class in which the true f value of the
NL of the sky. A wavelet decomposition for both the HW and
lth simulated CMB map falls.
SMHW was performed to determine the wavelet coefficents
The N output classes of the network were defined by
for each a set, and their third-order moments computed.
dividing some initial (anticipated) range of f values into lm
NL These statistics were provided as inputs to the neural net-
N equal-width subranges. For example, for a total range of
work. Each input vector contained 232 values for the HW
−30 6 f < 30 and a network with just 3 output classes,
NL and 680 for the SMHW.
inputvectorsconstructedfrommapswith−306f <−10
NL
were ascribed to class=1 with an associated target vector
t=(1,0,0),mapswith−106f <10toclass=2witht=
NL
(0,1,0),andthosewith 106fNL <30toclass=3witht= 1 http://planck.mpa-garching.mpg.de/cmb/fnl-simulations/
(cid:13)c 0000RAS,MNRAS000,000–000
WMAP 7yr local f constraints using neural networks 5
NL
4.2 Network architecture
45
HW
The architecture of our 3-layer neural networks are defined
SMHW
by two free parameters: the number of hidden nodes nhid NL40
f
andthenumberofoutputclasses,n ,intowhichthef
class NL n
rangeisdivided.Afurtherparameter,whichdeterminesthe o 35
n
accuracy of the network classifier, is the quantity of train- o
i
s
ing data ndata. Variation of theseparameters can affect the er 30
training efficiency so it is desirable to explore this training sp
i
spaceadequately in ordertofindan optimalset ofparame- d 25
ters.
Although the MemSys algorithm provides routines to 20
0 0.5 1 n 1.5 2 2.5
determinetheoptimal valueofthenumberofhiddennodes data x 104
using the Bayesian evidence Gull & Skilling (1999), in this
applicationn isdeterminedsimplybymeasuringtraining Figure 2. Results of the dispersion of fˆNL for 1000 Gaussian
hid simulationsfordifferentvaluesofndata.
times and the accuracy of the trained networks on an in-
dependent testing set. In this example, we have found that
in fact the optimal architecture contains no hidden nodes, ingset.ThisisplottedinthebottomtwopanelsofFigure3
resultingin whatiseffectively alinearclassifier. Thisis not fortheHWandSMHW,respectively.Wenotethat,ineach
surprising, since we are effectively ‘asking’ the network to case, this dispersion increases noticeable beyond the point
learn the mean value and dispersion of the third-order mo- where theTPRs on thetraining and testing sets diverge.
mentsofthewavelet coefficients foreach fNL;since theex-
pectation value is linearly dependent on the fNL, this net-
work architecture trivially satisfies thisrequirment.Indeed,
5 RESULTS
networks of this sort provide a simple way of obtaining the
(pseudo)inverseof any matrix. 5.1 Application to WMAP simulations
The number of output classes, n , of the network
class
We first applied our classifiers to 1000 WMAP-7yr simu-
is clearly related to the total range of f considered and
NL
lations made from Gaussian CMB maps (f = 0). For
size of the subranges into which this range is divided. Here NL
the HW classifier, we obtained hfˆ i = −1, which indi-
we consider networks with n = 9 (an odd number en- NL
class
cates the estimator is essentially unbiassed. Moreover, the
sures that f =0 does not lie on the boundary of a class)
NL dispersionoftheestimatorσ(fˆ )=33isextremelysimilar
The range of f was chosen a priori to correspond to NL
NL
to that obtained with the weighted least-squares method
approximately ±3σ, where σ is the dispersion in the f
estimates obtained previously using the standard χ2 miNnL- (σ(fˆNL) = 34). The full distribution of the estimator is
showninthetoppanelofFig.4.FortheSMHWclassifier,we
imisation method. Thus, the full range was taken to be
again found hfˆ i=−1, with a dispersion of σ(fˆ )=22,
−120 6 f < 120 for the HW and −76 6 f < 76 for NL NL
NL NL whichisveryclosetotheoptimalvalueofσ(fˆ )=21.The
the SMHW, resulting in subranges per class of width 27 NL
distributionoftheestimatorfortheSMHWisshowninthe
and17 units,respectively.This combination fulfilled all the
bottom panelof Fig. 4.
requirementsofclassificationovertherangeofoursimulated
The histogram bins in Fig. 4 have the same size and
data.
central values as those used to define the network classes.
Thequantityoftrainingdata,n ,determinestheac-
data We see that the classes at extremal f values are empty,
curacyoftheresultingclassification network.Naturally,the NL
indicating that the network placed no maps in these f
network accuracy increases with n , but it typically sat- NL
data ranges. Thus for estimating f from Gaussian or nearly
urates after a given number. We found that the quantity NL
Gaussian maps therange in f used is sufficiently wide.
of data required saturated at roughly n ∼ 10000 (see NL
data
We next applied our estimator to sets of non-Gaussian
Fig. 2).
simulations, each with a different non-zero f value. For
NL
eachtruef value,weanalysedthecorrespondingWMAP
NL
simulations and calculated the mean and dispersion of our
4.3 Training evolution
estimator fˆ for boththeHWand SMHWclassifiers. The
NL
Figure 3 illustrates the training evolution for the classifica- resultsare shown in fig.5,in which we plot themean value
tion network with n = 0 and n = 9. In the top two offˆ againstthetruef value.Weseethattheclassifiers
hid class NL NL
panelsweplotthetruepositiverates(TPR)ofthenetwork areunbiassedfor |f |.σ withan almost constantdisper-
NL
onthetrainingsetandthetestset,fortheHWandSHMW sion.Forlarger|f |values,however,theestimatorbecomes
NL
respectively; in each plot, the TPR on the training set has progressively more biassed and its dispersion decreases.
been mutliplied by a factor less than unity to highlight the The latter behaviour is simply understood as an edge
divergence with the TPR for the test set. We see that this effect dueto thefinitetotal range of f considered by the
NL
divergence occurs after ∼ 100 and ∼ 500 iterations of the networks. This point is illustrated in Fig. 6 in which we
MemSys optimiser for the HW and SMHW, respectively. plot the full distributions of fˆ obtained for a number of
NL
Thusthetrainingwasterminatedatthispointtoconstruct representativevaluesofthetruef .Weseethatfor|f |.
NL NL
our finalclassification networks. σ,we obtain close tosymmetricdistribution centredon the
Akeycriterion indeterminingthequalityofourclassi- true f value, with no maps being placed in the extreme
NL
fiersisthedispersionofthef valuesobtainedinthetest- classes. As |f | > σ, however, we see that the classifier
NL NL
(cid:13)c 0000RAS,MNRAS000,000–000
6 B. Casaponsa et al.
) )
(%40 (%35
e HW Training*3/4 e SMHW Training*2/3
at Testing at Testing
r35 r
s s
e e
v v30
i i
t t
i30 i
s s
o o
p p
e e
tru250 100 i2t0e0rati3o0n0s 400 500 tru250 500 iter1a0t00ions1500 2000
35 23
NL HW NL SMHW
n f 34 n f 22.5
o o
33
n n
o o 22
i i
s 32 s
r r
spe 31 spe 21.5
i i
d d
30 21
0 100 200 300 400 500 0 500 1000 1500 2000
iterations iterations
Figure 3. Evolution of the true positive rate for each iteration of the training process with a neural network with nhid = 0 and
ndata=10000.NotethattheTPRofthetrainingsethavebeenmultipliedbyafactorlessthanunityinordertohighlightthedivergence
ofthebehaviours.Thebottom panelsshowthevariationofthedispersionontheestimatefˆNL duringthetraining.LeftpanelsforHW
andrightpanelsforSMHW.
does begin to place maps in the extreme classes, resulting
0.015
in thedistribution offˆNL becoming severely skewed andno HW
longer centred on the true value. Of course, if one were to
encounter this behaviour in the analysis of a real data set, 0.01
one could simply alter the range of f considered by the
NL
network and retrain.
0.005
In any case, the above results show that both the HW
andSMHWnetworkclassifiersproduceunbiassedestimates
fˆNL provided −σ < fNL < σ. Moreover, the dispersions 0−100 −50 0 50 100
on these estimators are very similar to those obtained with
the classical weighted least squares (WLS) method, indi-
cating that neural networks can produce very accurate re- 0.02
sults within the limitations described above. In the case of SMHW
theSMHW,thisisaparticularly importantresultsincethe 0.015
complexity of the covariance matrix inversion required in
0.01
thestandardapproach isbypassedviatheuseoftheneural
networkclassifier.Curto et al.(2011a)usedaprincipalcom-
0.005
ponent analysis to reduce the covariance matrix dimension
to allow inversion.
0
−50 0 50
Figure 4.Distribution of fˆNL obtained from1000 Gaussian re-
alizationsforHW(top) andSMHW(bottom).
5.2 Application to WMAP 7-year data
Applying the neural network classifiers to real data (the
V+W WMAP 7-year data map), we obtain fˆ = −12
NL
for the HW and fˆ = 19 for the SMHW. Both these val-
NL
ues lie well within the corresponding dispersion of the esti-
mator. From the corresponding fˆ distributions obtained
NL
on simulated data, we find that 95% confidence limits are
−78 < f < 51 for the HW and −24 < f < 61 for the
NL NL
(cid:13)c 0000RAS,MNRAS000,000–000
WMAP 7yr local f constraints using neural networks 7
NL
fˆNL,data σ(fˆNL) hfˆNL,gaussi P2.5 P97.5
SMHW(NN) 19 22 −1 −43 42
SMHW(WLS)
37 21 0 −42 46
Curtoetal.2011b
HW(NN) −12 33 −1 −66 63
HW(WLS)
6 34 1 −68 67
Casaponsaetal.2011
Table 1. Results obtained with neural networks (NN) and weighted least squares (WLS). fˆNL,data is the best fitting value for V+W
WMAP data, hfˆNL,gaussi and σ(fˆNL) are the expected value and the standard deviation for Gaussian simulations. P2.5 and P97.5
representthepercentilevalues at95%confidence leveloffˆNL forGaussianrealizations.
f =−5 f =−20 f =−50 f =−120 f =−2 f =−10 f =−30 f =−70
NL NL NL NL NL NL NL NL
0.02 0.02 0.02 0.02 0.02 0.02 0.02
0.03
0.015 0.015 0.015 0.015 0.015 0.015 0.015
0.02
0.01 0.01 0.01 0.01 0.01 0.01 0.01
0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.01
0 0 0 0 0 0 0 0
−1000100 −1000100 −1000100 −1000100 −100 0 100 −100 0 100 −100 0 100 −100 0 100
f =5 f =20 f =50 f =120 f =2 f =10 f =30 f =70
NL NL NL NL NL NL NL NL
0.02 0.02 0.02 0.025 0.02 0.02 0.02
0.03
0.015 0.015 0.015 0.02 0.015 0.015 0.015
0.015 0.02
0.01 0.01 0.01 0.01 0.01 0.01
0.01
0.005 0.005 0.005 0.005 0.005 0.005 0.01
0.005
0 0 0 0 0 0 0 0
−1000100 −1000100 −1000100 −1000100 −100 0 100 −100 0 100 −100 0 100 −100 0 100
Figure 6. Distribution of fˆNL obtained from 200 non-Gaussian realizations with representative true fNL values, for HW (left) and
SMHW(right).
SMHW.2Wethereforeconcludethatthedataareconsistent 6 CONCLUSIONS
withtheGaussianhypothesis.WenotethattheSMHWcon-
fidence limits are very similar to those obtained with the Wehavetrainedamulti-classneuralnetworkclassifier with
optimal f estimator (Komatsu et al. 2011; Smith et al. third-order moments of the HW and SMHW coefficients
NL
2009). of non-Gaussian realizations in order to set constraints on
These results are summarised in Table 1, along with thelocalnon-linearcouplingtermfNL usingWMAP7-year
found via the weighted least squares (WLS) method. The data. We found that with a very simple network and with
latterresultsarealsoconsistentwithGaussianity.Itisworth fewiterations(requiringjustafewsecsCPUtime)itispos-
mentioning,however,thedifferentvaluesoffˆ obtainedby sible to produce the same results as those obtained with
NL
theneuralnetworkandtheWLSmethods,forbothHWand the weighted least squares method. This is an interesting
SMHW. Although all four values lie well within their cor- achievement, as it bypasses any covariance matrix related
respondingdispersions,eachmethodreturnsadifferentfˆ computations and their associated regularisation problems.
NL
value when applied to the same WMAP-7yr dataset. This Theestimationofthecovariancematrixwithbothwavelets
behaviour is to be expected, however, since these are four requirestheanalysis of at least 10000 Gaussian simulations
different estimators of f . Therefore, in general, they will which involves a huge demand in CPU time, in particular
NL
not be equal, even when applied to the same input data. withtheSMHWstatisitcs.Theerrorbarsontheestimation
Only the statistical properties (e.g. mean, dispersion) of offNLcomputedwithGaussiansimulationsareσ(fˆNL =33)
their sampling distributions are important. for HW and σ(fˆNL) =22 for SMHW, which are extremely
similar totheonesobtained in Casaponsa et al. (2010) and
Curto et al. (2011a) using the same statisitcs but a differ-
2 Note that the constraints are not corrected forthe unresolved ent estimator based on the weighted least squares method
pointsourcescontribution. (σ = 34, σ = 21 for HW and SMHW respectively). The
(cid:13)c 0000RAS,MNRAS000,000–000
8 B. Casaponsa et al.
AuldT.,BridgesM., HobsonM.P.,GullS.F.,2007, MN-
100
RAS,376, L11
HW
50 Babich D.,Creminelli P.,Zaldarriaga M., 2004, Journal of
Cosmology and Astro-Particle Physics, 8, 9
i
L
ˆfN 0 BaccigalupiC.,BediniL.,BuriganaC.,DeZottiG.,Farusi
h A.,MainoD.,MarisM.,PerrottaF.,SalernoE.,Toffolatti
−50 L., Tonazzini A., 2000, MNRAS,318, 769
BarreiroR.B.,HobsonM.P.,LasenbyA.N.,BandayA.J.,
−100 G´orski K. M., Hinshaw G., 2000, MNRAS,318, 475
−100 −50 0 50 100
fNL Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004,
Phys.Rev.D,402, 103
60 Bucher M., van Tent B., Carvalho C. S., 2010, MNRAS,
40 SMHW 407, 2193
Carballo R., Gonz´alez-Serrano J. I., Benn C. R., Jim´enez-
20
Luj´an F., 2008, MNRAS,391, 369
i
ˆfNL 0 CasaponsaB.,BarreiroR.B.,CurtoA.,Mart´ınez-Gonz´alez
h
−20 E., Vielva P., 2010, ArXiv e-prints
Cayo´nL.,Mart´ınez-Gonz´alezE.,Argu¨esoF.,BandayA.J.,
−40
G´orski K. M., 2003, MNRAS,339, 1189
−60
−60 −40 −20 0 20 40 60 Cayo´nL.,SanzJ.L.,Martinez-GonzalezE.,BandayA.J.,
fNL Argueso F., Gallegos J. E., Gorski K. M., Hinshaw G.,
2001, MNRAS,326, 1243
Figure5.ThemeananddispersionoffˆNLobtainedforanumber CruzM.,Mart´ınez-Gonz´alezE.,VielvaP.,Cayo´nL.,2005,
ofrepresentativevaluesofthetruefNLfortheHWnetwork(top) MNRAS,356, 29
andtheSMHWnetwork(bottom). Curto A., Mart´ınez-Gonz´alez E., Barreiro R. B., 2009b,
ApJ,706, 399
Curto A., Martinez-Gonzalez E., Barreiro R. B., 2011a,
constraintsforWMAP7-yeardatawerefoundtobe−78<
MNRAS,412, 1023
f <51 for the HW and −24<f <61 for the SMHW,
NL NL Curto A., Martinez-Gonzalez E., Barreiro R. B., Hobson
which are compatible to a Gaussian distribution as found
M. P., 2011b, submitted to MNRAS
in Smith et al. (2009); Curto et al. (2009b); Komatsu et al.
Curto A., Mart´ınez-Gonz´alez E., Mukherjee P., Barreiro
(2011); Casaponsa et al. (2010) and Curto et al. (2011b).
R. B., Hansen F. K., Liguori M., Matarrese S., 2009a,
The results obtained with the SMHW statistics are similar
MNRAS,393, 615
totheonesfoundinSmith et al.(2009)andKomatsu et al.
Elsner F., Wandelt B. D., 2009, ApJS,184, 264
(2011),whicharethemoststringentonescurrentlyavailable
Elsner F., Wandelt B. D., 2010, ApJ,724, 1262
atthelimitoftheWMAPresolution.Furtheranalysis,asto
Fergusson J. R.,Liguori M.,Shellard E. P.S.,2010, Phys-
thecontribution to f of unresolved point sources or fore-
NL ical Review D, 82, 023502
groundscanbeperformedbyapplyingthelinearclassifierto
Fergusson J.R.,Shellard E. P.S.,2009, Phys.Rev.D,80,
thestatisticsofnewsimulatedmapswiththischaracteristic
043510
signal.
Gangui A., Lucchin F., Matarrese S., Mollerach S., 1994,
ApJ,430, 447
G´orski K. M., Hivon E., Banday A. J., Wandelt B. D.,
ACKNOWLEDGMENTS Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ,
622, 759
We acknowledge partial financial support from the Spanish
GullS.F.,SkillingJ.,1999,Quantifiedmaximumentropy:
MinisteriodeCienciaeInnovaci´onproject AYA2010-21766-
Mem-Sys 5 users manual. Maximum Entropy Data Con-
C03-01 and from the CSIC-The Royal Society joint project
sultantsLtd, Royston
withreference2008GB0012.B.CasaponsathankstheSpan-
Hikage C., Matsubara T., Coles P., Liguori M., Hansen
ishMinisteriodeCienciaeInnovaci´onforapre-doctoralfel-
F. K., Matarrese S., 2008, MNRAS,389, 1439
lowship. The authors acknowledge the computer resources,
Hobson M., Lasenby A., 1998, 298, 905
technical expertise and assistance provided by the Span-
JaynesE.,2003,ProbabilityTheory:TheLogicofScience.
ishSupercomputingNetwork(RES)nodeatUniversidadde
Cambridge University Press
Cantabria. We acknowledge the use of Legacy Archive for
KomatsuE.,SmithK.M.,DunkleyJ.,BennettC.L.,Gold
MicrowaveBackgroundDataAnalysis(LAMBDA).Support
B.,HinshawG.,JarosikN.,LarsonD.,NoltaM.R.,Page
foritisprovidedbytheNASAOfficeofSpaceScience.The
HealPix package was used throughout the data analysis L., Spergel D. N., Halpern M. an Wright E. L., 2011,
ApJS,192, 18
(G´orski et al. 2005).
Komatsu E., Spergel D., 2001, Phys.Rev.D,63, 063002
KomatsuE.,SpergelD.N.,WandeltB.D.,2005,ApJ,634,
14
REFERENCES
LeshnoM., V.Y.,A.P., SchockenS.,1993, Neural Netw.,
Antoine J.-P., Vandergheynst P., 1998, Journal of Mathe- 6, 861
matical Physics, 39, 3987 MacKay D., 2003, Information Theory, Inference and
(cid:13)c 0000RAS,MNRAS000,000–000
WMAP 7yr local f constraints using neural networks 9
NL
Learning Algorithms. Cambridge University Press
Mart´ınez-Gonz´alez E., Gallegos J. E., Argueso F., Cayon
L., SanzJ. L., 2002, MNRAS,336, 22
McEwen J. D., Hobson M. P., Lasenby A. N., Mortlock
D.J., 2008, MNRAS,388, 659
McEwenJ.D.,VielvaP.,WiauxY.,BarreiroR.B.,Cayon
L., Hobson M. P., Lasenby A. N., Martinez-Gonzalez E.,
Sanz J. L., 2007, Journal of Fourier Analysis and Appli-
cations, 13, 495
Mukherjee P., Wang Y., 2004, ApJ, 613, 51
Natoli P., de Troia G., Hikage C., Komatsu E., Migliac-
cio M., Ade P. A. R., Bock J. J., Bond J. R., Borrill J.,
Boscaleri A.,ContaldiC. R.,Crill B.P.,deBernardis P.,
deGasperisG.,deOliveira-CostaA.,diStefanoG.,Hivon
E., 2010, MNRAS,408, 1658
Pietrobon D., 2010, Memorie della Societa Astronomica
Italiana Supplementi,14, 278
Salopek D. S., Bond J. R., 1990, Phys.Rev.D,42, 3936
Shahram M., Donoho D., Starck J., 2007, in Society of
Photo-OpticalInstrumentationEngineers(SPIE)Confer-
ence Series Vol. 6701 of Society of Photo-Optical Instru-
mentationEngineers(SPIE)ConferenceSeries,Multiscale
representation for data on thesphereand applications to
geopotential data
Smith K. M., Senatore L., Zaldarriaga M., 2009, Journal
of Cosmology and Astro-Particle Physics, 9, 6
Storrie-Lombardi M. C., Lahav O., Sodre Jr. L., Storrie-
Lombardi L. J., 1992, MNRAS,259, 8P
TenorioL.,JaffeA.H.,HananyS.,LineweaverC.H.,1999,
MNRAS,310, 823
Vanzella E., Cristiani S., FontanaA., NoninoM., Arnouts
S., Giallongo E., Grazian A., Fasano G., Popesso P.,
Saracco P., Zaggia S., 2004, A&A,423, 761
VerdeL.,WangL.,HeavensA.F.,KamionkowskiM.,2000,
MNRAS,313, 141
VielvaP.,Martinez-GonzalezE.,BarreiroR.B.,SanzJ.L.,
Cayon L., 2004, ApJ, 609, 22
Vielva P., Sanz J. L., 2010, MNRAS,404, 895
YadavA.P.S.,WandeltB.D.,2010, Advancesin Astron-
omy,2010
(cid:13)c 0000RAS,MNRAS000,000–000
