Fourth

IEEETRANSACTIONSONSIGNALPROCESSING,VOL.55,NO.6,JUNE20072965

Fourth-OrderCumulant-BadBlindIdentification

ofUnderdeterminedMixtures

LievenDeLathauwer,SeniorMember,IEEE,JoséphineCastaing,andJean-FrançoisCardoso

Abstract—Inthispaperwestudytwofourth-ordercumulant-

badtechniquesfortheestimationofthemixingmatrixinun-

derdeterminedindependentcomponentanalysis.Thefirstmethod

isbadonasimultaneousmatrixdiagonalization.Thecond

isbadonasimultaneousoff-diagonalization.Thenumberof

sourcesthatcanbeallowedisroughlyquadraticinthenumber

ofobrvations.Forbothmethods,explicitexpressionsforthe

maximumnumberofsourcesaregiven.Simulationsillustratethe

performanceofthetechniques.

IndexTerms—Cumulant,higherorderstatistics,higherorder

tensor,independentcomponentanalysis(ICA),parallelfactor

analysis,simultaneousdiagonalization,underdeterminedmixture.

I.I

NTRODUCTION

C

ONSIDERthefollowingbasiclinearmixturemodel:

(1)

reprentsmultichannelob-Thestochasticvector

rvations,thecomponentsofthestochasticvector

correspondtounobrvedsourcesignals,and

denotesadditivenoi.Theaprioriunknownmixingmatrix

characterizesthewaythesourcesarecom-

binedintheobrvations.Thegoalofindependentcomponent

analysis(ICA)[13],[31],orblindsourceparation(BSS),con-

sistsoftheestimationofthesourcesignalsand/orthemixing

matrixfromobrvationsof,assumingthatthesourcesare

statisticallyindependent.TheliteratureonICAaddressfor

.themostparttheso-calledoverdeterminedca,where

Here,weconsidertheunderdeterminedorovercompleteca,

.where

AlargeclassofalgorithmsforunderdeterminedICAstarts

fromtheassumptionthatthesourcesare(quite)spar[5],[26],

[29],[33],[39].Inthisca,thescatterplottypicallyshows

ManuscriptreceivedMarch9,2006;revidSeptember22,2006.This

workwassupportedinpartbytheRearchCouncilK.U.Leuvenunder

GrantGOA-AMBioRICS,CoEEF/05/006OptimizationinEngineering,in

partbytheFlemishGovernmentunderF.W.O.ProjectG.0321.06,Tournesol

2005—ProjectT20013,andF.W.O.rearchcommunitiesICCoS,ANMMM,

andinpartbytheBelgianFederalSciencePolicyOfficeunderIUAPP5/22.

Theassociateeditorcoordinatingthereviewofthismanuscriptandapproving

itforpublicationwasProf.TulayAdali.

L.DeLathauweriswiththeRearchGroupETIS,UMR8051,F95014

Cergy-PontoiCedex,FranceandalsowiththeRearchGroupESAT-SCD,

KatholiekeUniversiteitLeuven,Leuven,Belgium(e-mail:delathau@).

J.CastaingiswiththeRearchGroupETIS,UMR8051,F95014Cergy-

PontoiCedex,France(e-mail:castaing@).

J.-F.CardosoiswiththeTSIDepartment,ÉcoleNationaleSupérieuredes

Télécommunications,75634ParisCedex13,France(e-mail:cardoso@.

fr).

DigitalObjectIdentifier10.1109/TSP.2007.893943

highsignalvaluesinthedirectionsofthemixingvectors.The

extremamaybelocalizedbymaximizationofsomeclustering

measure[5],[26],[39].Someoftheclustering-badtechniques

performanexhaustivearchinthemixingvectorspace,andare

thereforeveryexpensivewhentherearemorethantwoobr-

vationchannels.Inapreprocessingstepalineartransformmay

beappliedsuchthatthenewreprentationofthedataissparr

(e.g.,short-timeFouriertransforminthecaofaudiosignals)

[5].Themethodin[1]onlyrequiresthatforeachsourceonearea

inthetime-frequencyplanecanbefoundwhereonlythatpar-

ticularsourceisactive;thesignalsmayoverlapanywhereel.

In[24]thedifferencebetweenlong-timestationarysourcesand

sourcesthatareonlyshort-timestationaryisexploitedtopa-

ratethelatter.

TherearetwoaspectstoICA:estimationofthemixingma-

trixandparationofthesources.Intheoverdeterminedca,

sourcesareusuallyparatedbymultiplyingtheobrvations

withthepudo-inverofthemixingmatrixestimate.Thisis

nolongerpossibleinthecaofunderdeterminedmixtures:for

eachsample

,thecorrespondingsourcesamplethatsatis-

fiesisonlyknowntobelongtoanaffinevarietyofdi-

—hencetheterm“underdetermined.”However,mension

themixingmatrixandthesourcedensitiesarestilluniqueunder

mildlyrestrictiveconditions[27].Uniquenessofthesourcedis-

tributionsallowsforthejointestimationofsourcesandmixing

matrixinaprobabilisticframework[34].However,eveninthe

caofunderdeterminedmixtures,theestimationofthemixing

matrixisanoverdeterminedproblem(e.g.,eSectionsIIand

III),suchthatitmakesntoestimatethemixingmatrixfirst,

andthenestimatethesources.Thesourcevalues

maysub-

quentlybeestimatedbymaximizingthelogposteriorlikelihood

[34].Inthecaofsparsources,characterizedbyLaplacian

densities,thiscanbeformulatedintermsofalinearprogram-

sourcescanbemingproblem[5],[10],[33].Ifatmost

activeatthesametime,thenforeachsampletheactivemixing

vectorsmaybedeterminedandthecorrespondingmixturein-

verted[29].Inthecaoffinitealphabetsignalsintelecommu-

nication,onemayperformanexhaustivearchoverallpossible

combinations.Inthispaperwefocusontheestimationofthe

mixingmatrix.Theestimateofthemixingmatrixmaysub-

quentlybeudtoparatethesourcesbymeansofthetech-

niquesmentionedearlier.

Thispaperprentsnewcontributionstotheclassofalge-

braicalgorithmsforunderdeterminedICA.In[15],[18],and

[19]algorithmsarederivedforthespecificcaoftwomix-

turesandthreesources.Anarbitrarynumberofmixingvectors

canbeestimatedfromtwoobrvationchannelsbysampling

derivativesofsufficientlyhighorderofthecondcharacter-

isticfunction[38].Amorestableversionof[38]isprented

1053-587X/$25.00©2007IEEE

2966IEEETRANSACTIONSONSIGNALPROCESSING,VOL.55,NO.6,JUNE2007

in[17].AlgebraicunderdeterminedICAisbadonthede-Hermiteanmatrices.Finally,werecallthedefinitionofthe

compositionofahigherordertensorinasumofrank-1terms.KroneckerproductandtheKhatri–Raoproduct[35]

Somelinkswiththeliteratureonhomogeneouspolynomialsare

discusdin[14].Asimultaneousmatrixdiagonalizationtech-

niquethatmaystillbeudwhenthenumberofsourcesex-

ceedsthenumberofnsorsisprentedin[42].In[3]theal-

gebraicstructureofthesixth-ordercumulanttensorispartially

exploited.Asimilarideacanbeappliedtoatoffourth-order

cumulanttensors,correspondingtodifferenttimelags,when

theindividualsourcesignalsaredependentoversometimein-

terval[28].Inthispaperwemerelyassumethatthesources

havenonzerokurtosis.Forconvenience,wealsoassumethat

thenoiisGaussian.Non-Gaussiannoileadstoaperturba-

tionoftheequations.Thisisadmissibleaslongastheperturba-

tionstaysrelativelysmall,i.e.,thesignal-to-noiratio(SNR)

hastobesufficientlyhigh.

Thepaperisorganizedasfollows.Afirstfourth-ordercu-

mulant-badapproachisdiscusdinSectionII.Theresulting

algorithmisbadonasimultaneousmatrixdiagonalization.A

condapproach,leadingtoasimultaneousoff-diagonalization,

isdiscusdinSectionIII.Simulationresultsareprentedin

SectionIV.SectionVistheconclusion.Theproofsofthethe-

oremsaregivenintheAppendix.Theprentationisinterms

ofcomplexsignals.Whenevertheresultscannotbedirectlyap-

pliedtorealdata,thisisexplicitlymentioned.

ThefoundationsofSectionIIwerelaidin[6].Somemath-

ematicalaspectsaredevelopedinmoredetailin[21].In[22],

[23]avariantofthetechniqueisprentedthatgeneralizessi-

multaneousmatrixdiagonalization-badmethods(involvinga

tofcorrelationmatrices,forinstance)totheunderdetermined

ca.

Notation:Scalarsaredenotedbylowercaitalicletters

,vectorsbylowercaboldfaceletters,

matricesbyboldfacecapitals,andtensorsbycal-

ligraphicletters.Italiccapitalsareudtodenote

indexupperbounds.Theentrywithrowindex

andcolumnindexinamatrix,i.e.,,issymbol-

izedby.Likewi,wehave.The

columnsofaredenotedbyWewillfrequentlyu

matrixreprentationsoftensors.

Tothisend,wedefine

Analogously,matriceswilloftenbestackedin-di-

mensionalvectors

Theinverofthelatteroperationisdenotedby

.Vectorizationofantensorisdone

asfollows:

ThesymbolstandsfortheKroneckerdelta,i.e.,Thelatterequationisesntiallyamatrixeigenvaluedecompo-

ifand0otherwi.denotesthespaceofsition(EVD),whichmayeasilybecomputed.Theeigenvectors

..

..

..

II.FOOBIA

LGORITHM

Considerthequadricovariance

.Duetothemultilinearitypropertyofcumulantten-

sors,wehave

(2)

inwhichisthekurtosisofthethsource.Thisisadecom-

positionofasymmetricfourth-ordertensorinasumofsym-

metricrank-1terms,cf.[9],[14],[16],[20],[21],[30],and

thereferencestherein.Theminimalnumberofrank-1tensors

inwhichahigherordertensorcanbedecompodiscalled

itsrank.Intermsof

and

,(2)canbewrittenas

(3)

Notethateachtermin(2)consistsofthecontributionofone

distinctsourceto.Hence,intermsofthequadricovariance,

“mixtureidentification”amountstothecomputationofdecom-

position(2)–(3).Wewillworkviaaconddecomposition,

whichisintroducedinthefollowingtheorem.

Theorem1:Atensor,satisfyingthesym-

metriesand,canbeeigendecompod

as

(4)

inwhichthematricesareHermiteanandmutuallyor-

thonormaltheEuclideaninnerproduct,andinwhich

thevaluesarerealandnonzero.istherankof

.Denote

and.

Then(4)isequivalentto

(5)

inwhichiscolumnwiorthonormal,with

,andinwhichthevaluesarereal

andnonzero.

FromTheorem1wehave

(6)

(7)

DELATHAUWERetal.:FOURTH-ORDERCUMULANT-BASEDBLINDIDENTIFICATIONOFUNDERDETERMINEDMIXTURES2967

havetobenormalizedinordertomakethematricesHer-

mitean;intheAppendixitisexplainedhowthiscanbedone.

Notethat,ifisfullcolumnrankandif(com-

plexca)or(realca),thenumberofsources

isgivenbytherankof.(Wenoticethatinarrayprocessing

applications,propertiesofthearraymaycautobe

columnrankdeficient—werefertoRemark2.)Wewillnow

showhow(2)–(3)and(6)–(7)arelinked.Assumeatthispoint

thatallsourcesaresuper-Gaussian,i.e.,

.

ThemoregeneralsituationwillbeaddresdinRemark1.From

(3)itisclearthatispositive(mi)definite.Wehavethefol-

lowingtheorem.

Theorem2:Let

bepositive(mi)definiteandassume

thatitcanbedecompodasin(3)and(7).Thenwehave

(8)

inwhich

isrealorthogonal.

Afterthecomputationofandfrom(7),thenextstep

isthecomputationof.Equation(8)showsthatmultiplication

ofbyyieldsamatrixofwhichthecolumnsare

Kroneckerproducts.TheKroneckerstructureofcor-

respondstotherank-1structureof.

Thismaybeexploited.Whatweneedisatoolthatallowsusto

distinguishbetweenHermiteanmatricesthatareatmostrank-1

andHermiteanmatricesofwhichtherankisgreaterthanone.

Sucha“rank-1detectingdevice”isintroducedinthefollowing

theorem.

Theorem3:Considerthemapping

definedby

(9)

Thenwehavethatifandonlyifisatmost

rank-1.

Definematrix,Hermiteanmatrices

andfourth-ordertensors.Now,

letbeanydiagonalmatrixandlet.Then,

usingthebilinearityof,itsrank-onedetectingfeature,and(8),

itisreadilyfoundthat.Thissuggeststo

determineamatrixfromthelatterequation,andfindasits

eigenmatrix.Morespecifically,wehavethefollowingtheorem.

Theorem4:Assumethatthetensors

,arelinearlyindependent.Thenthereexistprecily

linearlyindependentrealsymmetricmatrices

thatsatisfy

(10)

Thematriceshaveasacommoneigenmatrix,i.e.,

.

.

.

(11)

inwhicharediagonal.sourcesaresuper-Gaussian.Ifallsourcesaresub-Gaussian,i.e.,

Wecannowproceedasfollows.Given,thenwesimplyprocess.Inca

linearlyindependentmatricesarecomputedfrom(10),notallkurtosisvalueshavethesamesign,isindefinite.The

TABLEI

FOOBIA

LGORITHM

whichisjustahomogeneoustoflinearequations.Thenthe

matrixfollowsfromthesimultaneousEVD(11).

Inpractice,weworkwithnoisycumulantestimates,suchthat

(10)willonlyapproximatelybesatisfied.Thematrices

are

thendeterminedasfollows.Duetothesymmetryof,and

thefactthat,(10)canbewrittenas

(12)

Intheusualformofatofhomogeneouslinearequations,we

have

(13)

inwhichthecoefficientmatrixisgivenby

(14)

Theleast-squaressolutionof(13)consistsoftherightsin-

gularvectorsofthatcorrespondtothesmallestsingular

values.Afterstackingthevectorsinuppertriangularmatrices

,inthemannersuggestedby(13),thematrices

areobtainedas.Thefollowing

theoremguaranteesthatthevectorsarereal,evenintheca

ofnoisycumulantestimates.

Theorem5:Therightsingularvectorsofthematrixin(14)

arereal.

Aftercomputationofthematrices,thecommoneigen-

matrixin(11)canbeobtainedbymeansoftheJacobialgo-

rithmdevelopedin[7]and[8].Multiplicationby,

asin(8),yieldsamatrix

ofwhichthecolumns

aretheoreticallyproportionalto.Inpractice,weesti-

matefromthebestrank-1approximationof.The

overallFourth-Order-OnlyBlindIdentification(FOOBI)algo-

rithmisoutlinedinTableI.

Remark1:Inthederivationabove,wehaveassumedthatall

2968IEEETRANSACTIONSONSIGNALPROCESSING,VOL.55,NO.6,JUNE2007

derivationthenstillapplies,withtheexceptionthatis-or-

thogonal[41]insteadoforthogonal.Thisimpliesthat,forthe

simultaneousdiagonalization(11),avariantofthealgorithmin

[7]and[8]hastobeworkedout,thatinvolves-orthogonalma-

trices.Wecanalsoworkasfollows.Insteadofimposing-or-

thogonality,wesimplystartfrom

(15)

withsomerealnonsingularmatrix.Theprocedureisesn-

tiallythesame,butthesimultaneousdiagonalizationin(11)now

involvesarealnonsingularmatrix.Foralgorithmsforthistype

ofsimultaneousdiagonalizationwereferto[20],[40],[42],and

thereferencestherein.

Theconditionon

inTheorem4yieldsan

upperboundonthenumberofsourcesFOOBIcanhandle.We

callaproperty“generic”whenitholdswithprobabilityone

whentheentriesofthemixingmatrixaresampledfromcon-

tinuousprobabilitydensityfunctions.Wehavethefollowing

theorem.

Theorem6:Inthecomplexca,linearindependenceof

isgenericallyguaranteedif

.Intherealca,isboundedas

follows:

Remark2:In[11],[12],and[25]itisexplainedthatinan-

tennaarrayapplications,thecharacteristicsoftheantennasand

thegeometryofthearraymayinduceastructureintheentries

ofthehigherordercumulantthatlimitsthenumberofsources

thatcaneffectivelybedealtwith.Suchastructureisneglected

inTheorem6.Asaresult,thenumberofsourcesthatcanbeal-

lowedisboundedbytheminimumof:1)thenumberofsources

inTheorem6and2)themaximalnumberofvirtualnsors

(VSs),derivedin[11],[12],and[25].

III.FOOBI-2A

LGORITHM

Likeinthepreviousction,westartfromtheEVD

(7).Generically,aslongas(complexca)or

(realca),thenumberofsourcescorre-

spondstotherankof.Inthisctionweassumethatall

sourcesaresuper-Gaussian.(Ifallsourcesaresub-Gaussian,

thenweprocessinsteadof.)ThismeansthatTheorem2

stillapplies.Wenowintroduceanewrank-1detectingdevice.

Theorem7:Considerthemapping

definedby

(16)

Thenwehavethatifandonlyifisatmost

rank-1.

Letanddefinesymmetric

matricesby

Thefollowingtheoremsuggestsanewalgorithmforthecom-hand,thenumberofdistinctrealparametersinthedecompo-

putationof

.sition,withreal,isequalto

TABLEII

FOOBI-2A

LGORITHM

Theorem8:Thematrixin(8)satisfies

(17)

Thistheoremshowsthatthematrix

canbecomputedby

meansofsimultaneousoff-diagonalizationof(real

ca)or(complexca)realsymmetricmatrices.Thesimul-

taneousoff-diagonalizationcanberealizedbymeansofasimple

variantoftheJacobialgorithmderivedin[7],[8].Itsufficesto

choineachsteptheJacobirotationthatminimizes(instead

ofmaximizes)thesumofthesquareddiagonalentries.Simul-

taneousoff-diagonalizationalsoappearedin[4].TheFOOBI-2

algorithmissummarizedinTableII.

TheiterationthatformsthecoreofFOOBI-2(step4inAl-

gorithmII)iscomputationallymoreexpensivethanFOOBI’s

coreiteration(step5inAlgorithmI),becauthesimultaneous

off-diagonalizationinvolvesmorematrices.Ontheotherhand,

FOOBIrequiresthecomputationofpartoftheSVDofthe

matrix(step4inAlgorithmI).Also,

FOOBI-2islessrestrictiveintermsofthenumberofsources

thatcanbeallowed.Itonlyrequiresthat(complex

ca)or(realca),provideddecomposition

(2)isunique.Conquently,weinvestigateunderwhichcondi-

tionsgenericuniquenessholds.(Wenoticethatinnongeneric

casRemark2stillholds.)

In[16]itisstatedthatadecompositioninrank-1termsis

genericallyuniquewhenthenumberofparametersinthede-

compositionisstrictlysmallerthanthenumberofdistincttensor

entries.Whenbothnumbersareequal,thengenericallyonly

afinitenumberofdecompositionsarepossible.Inthecom-

plexca,thetotalnumberofdistinctrealanddistinctimagi-

narypartsoftheentriesofageneric

tensor

satisfyingthesymmetriesand

isgivenby

(18)

whereweassumethat

when.Ontheother

DELATHAUWERetal.:FOURTH-ORDERCUMULANT-BASEDBLINDIDENTIFICATIONOFUNDERDETERMINEDMIXTURES2969

.Themaximalrankforwhichthedecompositionisunique

isthengivenbyinthefollowingtable:

Notethatcanbegreaterthan.Intherealca,the

numberofdistinctentriesofagenericsymmetric

tensorisequalto,whilethenumber

ofparametersinthedecompositionequals.Themaximal

rankforwhichthedecompositionisuniqueisthengivenby

inthetableabove.Notethatcanbegreaterthan

.

IV.S

IMULATIONS

Inthefirstsimulation,narrow-bandsourcesarereceived

byauniformcirculararray(UCA)ofidenticalnsorsof

radius.Weassumefree-spacepropagation.Thismeansthat

theentriesofthemixingmatrixbeforenormalizationaregiven

by

where

and.Wehave

.Themixingmatrixisobtainedbydividing

thecolumnsofbytheirFrobeniusnorm.Weconsidertwo

cas:and.Thevaluesofarenotgreater

thanthenumberoffourth-orderVSsoftheUCA[11],[12],

[25].Thedirections-of-arrival(DOAs)ofthesourcesaregiven

by

and

.Intheca,

weconsiderthefirstfiveDOAs.Thesourcesareunit-variance

QAM4inbaband,whichmeansthattheytaketheirvalues

equallylikelyinthet.Additivezero-mean

complexGaussiannoiisaddedtothedata.Themixing

matrixisestimatedbymeansof:1)theFOOBIalgorithm;

2)theFOOBI-2algorithm;and3)theBIRTHalgorithm[2]

(or6-BIOME1algorithm,intheterminologyof[3]),which

usthesixth-ordercumulantoftheobrvations.(Wenote

thatthe6-BIOME3algorithmissomewhatmoreaccurate

than6-BIOME1,attheexpenofahighercomputationcost

[3].)Theprecisionismeasuredintermsofthemeanrelative

error

,inwhichthenormistheFrobenius

normandinwhichreprentstheoptimallyorderedand

scaledestimateof.WeconductMonteCarloexperiments

consistingof100runs.

Fig.1showstheaccuracyasafunctionoftheSNR,when

5000samplesareud.TheFOOBIandFOOBI-2curvesprac-

ticallycoincide.BIRTHislessaccurate.Wehavealsocompared

totheAC-DCalgorithm[42],appliedtothedominant

Her-

miteaneigenmatricesofthefourth-ordercumulant.Thismeans

thatexactlythesamestatisticsasinFOOBIandFOOBI-2are

ud.However,AC-DCfailedtoreliablyestimatethemixture.

Fig.2showstheaccuracyasafunctionofthenumberofdataInFig.3wecomparethecomputationalcostofthealgo-

samples

,fortheca.TheSNRwastakenequalto16rithms.FOOBIandFOOBI-2areaboutequallyexpensivein

dB.Again,theFOOBIandFOOBI-2curvespracticallycoincidethissimulation.BIRTHisaboutafactor40moreexpensivethan

andBIRTHislessaccurate.FOOBIandFOOBI-2.ThereasonisthatBIRTHrequiresthe

Fig.1.AccuracyasafunctionofSNRinthefirstexperiment(

J=4;R=

5;6;5000

samples).

Fig.2.Accuracyasafunctionofdatalengthinthefirstexperiment(

J=

4;R=5;16

dB).

Fig.3.Computationtimeasafunctionofdatalengthinthefirstexperiment

(

J=4;R=5;16

dB).

2970IEEETRANSACTIONSONSIGNALPROCESSING,VOL.55,NO.6,JUNE2007

Fig.4.Accuracyasafunctionofangleoffirstmixingvector(

J=4;R=

5;16

dB).

estimationofcumulantentries.Theestimationofthe

sixth-ordercumulantaccountsformorethan90%ofthetotal

computationalcost.Theestimationofthefourth-ordercumu-

lantaccountsforabout10%ofthetotalFOOBIorFOOBI-2

costwhen200samplesareud;respectively,70%when5000

samplesareud.Thecomputationtimevarieslittleasafunc-

tionoftheSNR.

Fig.4showstheaccuracyasafunctionoftheconditionofthe

problem,fortheca

.TheSNRwastakenequalto16dB,

and5000sampleswereud.Thefigureshowswhathappens

whenisvariedfromto(thelattervalue

isverycloto).Again,theFOOBIandFOOBI-2curves

practicallycoincide.FOOBIandFOOBI-2aremoreaccurate

thanBIRTHwhentheproblemiswellconditioned.Ontheother

hand,whenthetwomixingvectorsareveryclo,BIRTHis

morereliablethanFOOBIandFOOBI-2.Thereasonisthatthe

vectors

andarelessclothanthe

vectorsand,asexplainedin[12].

Inthecondsimulation,narrow-bandsourcesare

receivedbyaUCAofidenticalnsors.Thenumber

ofsourcesislessthanthenumberoffourth-orderVSsforthis

array,whichisequalto6[12].However,thenumberofsources

isabovetheFOOBIbound(Theorem6).Conquently,the

mixingmatrixisonlyestimatedbymeansof:1)theFOOBI-2

algorithmand2)theBIRTHalgorithm.TheDOAsofthe

sourcesareequaltothefirstfiveDOAsinthefirstexperiment.

Allotherparametersareasinthefirstexperiment.

Fig.5showstheaccuracyasafunctionoftheSNR.FOOBI-2

turnsouttobemoreaccuratethanBIRTH.Similarcurvesasin

thefirstexperimenthavebeenobtainedfortheaccuracyandthe

computationalcostasafunctionofthenumberofsamples.

Finally,weshowtheresultsofasimulationwithentirelysyn-

theticdata.Inthissimulation,thereare18obrvationchan-

nelsand25sources.Thesourceshaveunitkurtosis.Theentries

ofthemixingmatrixaredrawnfromazero-meanunit-vari-

ancecomplexGaussiandistribution.Thecolumnsaresub-

quentlyscaledtounitlength.Thenoi-freecumulantiscom-

puteddirectlyfrom(3).Whenevertheconditionnumberof

isgreaterthan100,anewmixingmatrixisgenerated,sothat

Fig.5.AccuracyasafunctionofSNRinthecondexperiment(

J=

3;R=

5;5000

samples).

Fig.6.AccuracyasafunctionofSNRinthethirdexperiment

(J=18;R=

25)

.

wedonotconsiderverelyill-conditioneddata.Additivezero-

meanGaussiannoiisaddeddirectlyonthecumulant.Withthe

noitermreprentedby,theSNRisdefined

as,inwhichthenormsareFrobeniusnorms.The

mixingmatrixisestimatedbymeansoftheFOOBIalgorithm.

Fig.6showstheresultsofaMonteCarloexperimentconsisting

of100runs.

V.C

ONCLUSION

Inthispaperwehavestudiedtheestimationofthemixing

matrixfromtheobrvedfourth-ordercumulanttensorinunder-

determinedICA.Aslongasthenumberofsourcesislessthan

thenumberofVSsoftheantennaarray(iftheresultsof[11],

[12],and[25]apply),itcanbefoundastherankofamatrixrep-

rentationofthecumulant.Underaspecificconditiononthe

mixingvectors,allowingforanumberofsourcesthatincreas

quadraticallywiththenumberofobrvations,thenoi-free

solutionmaybefoundfromanEVD.Fornoisydatawepro-

podtheFOOBIalgorithm,whichcomputesthesolutionby

meansofasimultaneousmatrixdiagonalization.Acondal-

gorithm,calledFOOBI-2,wasbadonasimultaneousoff-di-

agonalization.FOOBI-2isevenlessrestrictiveinthenumberof

DELATHAUWERetal.:FOURTH-ORDERCUMULANT-BASEDBLINDIDENTIFICATIONOFUNDERDETERMINEDMIXTURES2971

sourcesthanFOOBI.Thealgorithmsarebadonnewresultseven

thedecompositionofafourth-ordersymmetrictensorinayieldsanti-Hermiteaneigenmatrices.Wealsomentionthat,in

sumofsymmetricrank-1terms.Wehavedeterminedthemax-thecomplexca,therankcanbeaslargeas,themax-

imumnumberoftermssuchthatthisdecompositionisgeneri-

callyunique.Throughoutthepaper,boththerealandthecom-

plexcahavebeenaddresd.Theperformanceofthealgo-

rithmshasbeenillustratedbymeansofsimulations.

A

PPENDIX

I

P

ROOFS

Theorem1:

Proof:Duetothesymmetry,thematrix

isHermitean.Hence,itsEVDtakestheformof(5),with

columnwiorthonormalandreal.Thetensorcanthus

bedecompodasin(4),withmutuallyorthonormal

thescalarproductofmatrices,andreal.Considerthetensor

,definedby,anditsmatrixreprentation

.Wehave

(19)

inwhich

isdefinedby

.Ontheotherhand,

wehave

(20)

Becauofthesymmetry

,wehaveand

.Inthecawherealleigenvaluesaredifferent,projec-

torsandcorrespondingtothesameeigenvalueare

equal.Hence

(21)

If

isamultipleeigenvalue,thenthecorrespondingrank-1pro-

jectorscanbechonequal.Weconcludethattheprojectorssat-

isfythesamesymmetriesasitlf.

Nowweshowthatthisimpliesthatthematricescanal-

waysbenormalizedtoHermiteanmatrices.Notethatthepro-

jectordoesnotchangewhenismultipliedbyaunit-

modulusscalar.Let.Ifsomediagonalentryof

,say,isnonzero,thenwechoosuchthat

isreal.Sincewehaveforall

isHermitean.Ifallthediagonalentriesofarezero,

thenweproceedasfollows.Firstnoticethat(21)impliesthat

all.If,say,,thenwemultiply

bychonsuchthat.Sincewehave

isHermitean.

Thecomputationoftensordecomposition(4)amountstothe

computationoftheclassicalmatrixEVDin(5),inwhichthe

eigenvectorsarenormalizedinordertomakethematrices

Hermitean,asexplainedintheproof.Weemphasizethatthe

eigenmatricesarenotHermiteanbydefault,astheymaybemul-

tipliedbyanyunit-modulusscalar.Multiplicationby

imalrankofmatrices.TheequivalentofTheorem1

forreal-valuedtensorsissimplyobtainedbydroppingcomplex

conjugations.Theproofistrivial.Sincetheeigenmatricesare

realsymmetrichere,

isboundedby,thedimen-

sionofthevectorspaceofrealsymmetricmatrices.

Theorem2:

Proof:Bothandaresquare

rootsofthepositive(mi)definitematrix.Hence,theyare

relatedasin(8),withunitary.Wewillnowshowthat

isinfactreal.Considerthepermutationmatrix,

definedby

elwhere.

Fromthesymmetrypropertiesofthecolumnsof

and

,wehave

(22)

Combinationof(8)and(22)showsthat

isreal.

Theorem3:

Proof:The“if”partisobvious.Forthe“onlyif”part,we

startfrom

whichimplies

Thelatterequationcanbewritteninmatrixtermsas

(23)

Theca

canbediscarded,sinceitimpliesthat

andhence,sinceisHermitean.Dividing

(23)byshowsthattheunittraceHermiteanmatrix

satisfies.Hence,isanorthogonal

projector.Moreover,sincethedimensionoftheimagespaceof

anorthogonalprojectorisequaltoitstrace,therankofis

equaltoone.Weconcludethatisrank-1.Thetheoremcan

alsobeprovedinanalogywith[21,Th.2.1].

Theorem4:

Proof:Wefirstshowthateveryrealsymmetricmatrix

thatsatisfies(10),hasaseigenmatrix.Duetothebilinearity

of,wehavefrom(8)

(24)

Substitutionof(24)in(10)yields

2972IEEETRANSACTIONSONSIGNALPROCESSING,VOL.55,NO.6,JUNE2007

AccordingtoTheorem3,wehaveasfollows:

.Additionallytakingintoaccountthesymmetryof

andthefactthatissymmetricinitsarguments,weobtain

(25)

Ifthetensors,arelinearly

independent,thenthecoefficientsin(25)havetobezero

(26)

Thiscanbewritteninamatrixformatas

(27)

inwhichisdiagonal.Sincelinearindependenceofmatrices

amountstolinearindependenceofdiagonalmatrices,at

mostrealsymmetricmatricescansatisfy(10).Ontheother

hand,itiseasytoverifythatanyrealdiagonalmatrixgen-

eratesarealsymmetricmatrixthatdoessatisfy(10).This

provesthetheorem.

Theorem5:

Proof:ItsufficestoprovethatisHermitean.This

canbedonebycomputingitsentriesandtakingintoaccount

thatthematricesareHermitean.

Theorem6:

Proof:Thecomplexcaisatechnicalvariantoftheproof

of[21,Th.2.5].Therealcaisanalyzedin[37].Analgorithm

isdescribedthatallowstocomputeforanygiven.Itis

conjecturedthatintherealcatheboundisoftheform

where

if

if

Theorem7:

Proof:Itiseasytoverifythatifisrank-1.

Forthe“onlyif”part,lettheEVDofbegivenby.

Wehaveifandonlyif

Hence,atmostoneeigenvaluecanbedifferentfromzero.

Theorem8:

Proof:Duetothebilinearityof,wehave

(28)

Thisequationcanbewrittenintermsof

Thisequationisequivalentwith(17)becauofthelinkbe-

tweenand.

R

EFERENCES

[1]F.AbrardandY.Deville,“Atime-frequencyblindsignalparation

methodapplicabletounderdeterminedmixturesofdependentsources,”

SignalProcess.,vol.85,no.7,pp.1389–1403,Jul.2005.

[2]L.Albera,A.Ferréol,P.Comon,andP.Chevalier,“Sixthorderblind

identificationofunderdeterminedmixtures(BIRTH)ofsources,”in

Proc.4thInt.Symp.IndependentComponentAnalysisandBlindSignal

Separation(ICA’03),Nara,Japan,Apr.2003.

[3]L.Albera,A.Ferréol,P.Comon,andP.Chevalier,“Blindidentification

ofovercompletemixturesofsources(BIOME),”Lin.Alg.Appl.,vol.

391,pp.1–30,Nov.2004.

[4]A.Belouchrani,K.Abed-Meraim,M.G.Amin,andA.M.Zoubir,

“Jointantidiagonalizationforblindsourceparation,”inProc.IEEE

Int.Conf.onAcoustics,Speech,andSignalProcessing(ICASSP-01),

SaltLakeCity,UT,May2001,vol.5,pp.2789–2792.

[5]P.BofillandM.Zibulevsky,“Underdeterminedblindsourcepara-

tionusingsparreprentations,”SignalProcess.,vol.81,no.11,pp.

2353–2362,Nov.2001.

[6]J.-F.Cardoso,“Super-symmetricdecompositionofthefourth-ordercu-

mulanttensor.Blindidentificationofmoresourcesthannsors,”in

Proc.IEEEInt.Conf.onAcoustics,Speech,andSignalProcessing

(ICASSP-91),Toronto,ON,Canada,1991,pp.3109–3112.

[7]J.-F.CardosoandA.Souloumiac,“Blindbeamformingfor

non-Gaussiansignals,”Proc.Inst.Elect.Eng.,F,vol.140,no.6,

pp.362–370,Dec.1993.

[8]J.-F.CardosoandA.Souloumiac,“Jacobianglesforsimultaneousdi-

agonalization,”SIAMJ.MatrixAnal.Appl.,vol.17,no.1,pp.161–164,

Jan.1996.

[9]J.CarrollandJ.Chang,“Analysisofindividualdifferencesinmulti-

dimensionalscalingviaan-waygeneralizationof‘Eckart-Young’

N

decomposition,”Psychometrika,vol.35,pp.283–319,1970.

[10]S.S.Chen,D.L.Donoho,andM.A.Saunders,“Atomicdecomposition

bybasispursuit,”Dept.Statist.,StanfordUniv.,Stanford,CA,Tech.

Rep.30401,1996.

[11]P.ChevalierandA.Ferréol,“Onthevirtualarrayconceptforthe

fourth-orderdirectionfindingproblem,”IEEETrans.SignalProcess.,

vol.47,no.9,pp.2592–2595,Sep.1999.

[12]P.Chevalier,L.Albera,A.Ferréol,andP.Comon,“Onthevirtual

arrayconceptforhigherorderarrayprocessing,”IEEETrans.Signal

Process.,vol.53,no.4,pp.1254–1271,Apr.2005.

[13]P.Comon,“Independentcomponentanalysis,Anewconcept?,”Signal

Process.,vol.36,no.3,pp.287–314,Apr.1994.

[14]P.ComonandB.Mourrain,“Decompositionofquanticsinsumsof

powersoflinearforms,”SignalProcess.,vol.53,no.2,pp.93–107,

Sep.1996.

[15]P.Comon,“Blindidentificationandsourceparationin2

2

3under-

determinedmixtures,”IEEETrans.SignalProcess.,vol.52,no.1,pp.

11–22,Jan.2004.

[16]P.Comon,“Canonicaltensordecompositions,”Lab.I3S,Sophia-An-

tipolis,France,Tech.Rep.RR-2004-17,Jun.2004.

[17]P.ComonandM.Rajih,“Blindidentificationofunder-determined

mixturesbadonthecharacteristicfunction,”inProc.IEEEInt.

Conf.onAcoustics,Speech,andSignalProcessing(ICASSP’05),

Philadelphia,PA,Mar.18–23,2005,vol.4,pp.1005–1008.

[18]L.DeLathauwer,P.Comon,andB.DeMoor,“ICAalgorithmsfor

3sourcesand2nsors,”inProc.6thSignalProcessingWorkshop

onHigherOrderStatistics,Caesarea,Israel,Jun.14–16,1999,pp.

116–120.

[19]L.DeLathauwer,B.DeMoor,andJ.Vandewalle,“AnalgebraicICA

algorithmfor3sourcesand2nsors,”inProc.10thEur.SignalPro-

cessingConf.(EUSIPCO2000),Tampere,Finland,Sep.5–8,2000.

[20]L.DeLathauwer,B.DeMoor,andJ.Vandewalle,“Computationof

thecanonicaldecompositionbymeansofasimultaneousgeneralized

Schurdecomposition,”SIAMJ.MatrixAnal.Appl.,vol.26,pp.

295–327,2004.

[21]L.DeLathauwer,“Alinkbetweenthecanonicaldecompositioninmul-

tilinearalgebraandsimultaneousmatrixdiagonalization,”SIAMJ.Ma-

trixAnal.Appl.,vol.28,no.3,pp.642–666,2006.

DELATHAUWERetal.:FOURTH-ORDERCUMULANT-BASEDBLINDIDENTIFICATIONOFUNDERDETERMINEDMIXTURES2973

[22]L.DeLathauwerandJ.Castaing,“Second-orderblindidentificationof

underdeterminedmixtures,”inIndependentComponentAnalysisand

BlindSourceSeparation,Proc.6thInt.Conf.IndependentComponent

Analysis(ICA2006).Berlin,Germany:Springer-Verlag,Mar.5–8,

2006,vol.3889,LectureNotesinComputerScience,pp.40–47.

[23]L.DeLathauwerandJ.Castaing,“Independentcomponentanalysis

badonsimultaneousmatrixdiagonalization:Theunderdetermined

ca,”Lab.ETIS,Cergy-Pontoi,France,Tech.Rep.,submittedfor[40]A.-J.vanderVeenandA.Paulraj,“Ananalyticalconstantmodulus

publicationtoIEEETrans.SignalProcess.

[24]Y.Deville,M.Benali,andF.Abrard,“Differentialsourceparation

forunderdeterminedinstantaneousorconvolutivemixtures:Concept

andalgorithms,”SignalProcess.,vol.84,no.10,pp.1759–1776,Oct.

2004.

[25]M.C.DoganandJ.M.Mendel,“Applicationsofcumulantstoarraynwithapplicationinblindsourceparation,”IEEETrans.Signal

processing—PartI:Apertureextensionandarraycalibration,”IEEEProcess.,vol.50,no.7,pp.1545–1553,Jul.2002.

Trans.SignalProcess.,vol.43,no.5,pp.1200–1216,May1995.

[26]D.Erdogmus,L.Vielva,andJ.C.Principe,“Nonparametricestimation

andtrackingofthemixingmatrixforunderdeterminedblindsource

paration,”inProc.ICA’01,SanDiego,CA,Dec.2001,pp.189–194.

[27]J.ErikssonandV.Koivunen,“Identifiability,parability,andunique-

nessoflinearICAmodels,”IEEESignalProcess.Lett.,vol.11,no.7,

pp.601–604,Jul.2004.

[28]A.Ferréol,L.Albera,andP.Chevalier,“Fourthorderblindidentifica-

tionofunderdeterminedmixturesofsources(FOBIUM),”IEEETrans.

SignalProcess.,vol.53,no.5,pp.1640–1653,May2005.

[29]P.Georgiev,F.Theis,andA.Cichocki,“Sparcomponentanalysis

andblindsourceparationofunderdeterminedmixtures,”IEEETrans.

NeuralNetw.,vol.16,no.5,pp.992–996,Jul.2005.

[30]R.A.Harshman,“FoundationsofthePARAFACprocedure:Model

andconditionsforan‘explanatory’multi-modefactoranalysis,”UCLA

WorkingPapersinPhonetics,vol.16,pp.1–84,1970.

[31]A.Hyvärinen,J.Karhunen,andE.Oja,IndependentComponentAnal-

ysis.NewYork:Wiley,2001.

[32]J.B.Kruskal,“Rank,decomposition,anduniquenessfor3-wayand

N

Eds.Amsterdam,TheNetherlands:North-Holland:ElvierScience,

-wayarrays,”inMultiwayDataAnalysis,R.CoppiandS.Bolasco,

1989.

[33]Y.Li,A.Cichocki,andS.-I.Amari,“Analysisofsparreprenta-

tionandblindsourceparation,”NeuralComput.,vol.16,no.6,pp.

1193–1234,Jun.2004.

[34]M.LewickiandT.J.Sejnowski,“Learningovercompletereprenta-