Table Of ContentA PROBABILISTIC LEAST-MEAN-SQUARES FILTER
Jesu´sFerna´ndez-Bes,V´ıctorElvira StevenVan Vaerenbergh
DepartmentofSignalTheoryandCommunications DepartmentofCommunicationsEngineering
UniversidadCarlosIIIdeMadrid UniversityofCantabria
5 Legane´s(Spain) Santander(Spain)
1
0
2
ABSTRACT Gaussianstate-spacemodel(i.e. Kalmanfilter)andtheRLS
n filter,bySayedandKailath[7]andthenbyHaykinetal. [8].
a We introduce a probabilistic approachto the LMS filter.
TheRLSadaptivefilteringalgorithmemergesnaturallywhen
J
By means of an efficient approximation, this approach pro-
one defines a particular state-space model (SSM) and then
7 vides an adaptable step-size LMS algorithm together with a
performs exact inference in that model. This approach was
2
measureofuncertaintyabouttheestimation. Inaddition,the
laterexploitedin[9]todesignakernelRLSalgorithmbased
proposed approximation preserves the linear complexity of
] onGaussianprocesses.
L thestandardLMS.Numericalresultsshowtheimprovedper-
A first attempt to approximate the LMS filter from a
M formanceofthealgorithmwithrespecttostandardLMSand
probabilistic perspective was presented in [10], focusing on
state-of-the-artalgorithmswithsimilarcomplexity. Thegoal
t. ofthiswork,therefore,istoopenthedoortobringsomemore akernel-basedimplementation. Thealgorithmof[10]makes
a
useofaMaximumaPosteriori(MAP)estimateasanapprox-
t Bayesianmachinelearningtechniquestoadaptivefiltering.
s imationforthepredictivestep. However,thisapproximation
[
IndexTerms— probabilisticmodels,least-mean-squares, doesnotpreservetheestimateoftheuncertaintyineachstep,
1 adaptivefiltering,state-spacemodels thereforedegradingtheperformanceofthealgorithm.
v
In this work, we provide a similar connection between
9
state-space models and least-mean-squares (LMS). Our ap-
2 1. INTRODUCTION
9 proach is based on approximating the posterior distribution
6 Probabilistic models have proven to be very useful in a lot with an isotropic Gaussian distribution. We show how the
0 of applications in signal processing where signal estimation computationofthisapproximatedposteriorleadstoalinear-
.
1 isneeded[1,2,3]. Someoftheiradvantagesarethat1)they complexityalgorithm,comparabletothestandardLMS.Sim-
0 forcethedesignertospecifyalltheassumptionsofthemodel, ilarapproacheshavealreadybeendevelopedfora varietyof
5
2)theyprovideaclearseparationbetweenthemodelandthe problems such as channel equalization using recurrent RBF
1
: algorithmusedto solveit, and3)theyusuallyprovidesome neuralnetworks[11],orBayesianforecasting[12]. Here,we
v measureofuncertaintyabouttheestimation. show the usefulness of this probabilistic approach for adap-
i
X On the other hand, adaptive filtering is a standard ap- tivefiltering.
r proach in estimation problems when the input is received The probabilistic perspective we adopt throughout this
a
as a stream of data that is potentially non-stationary. This work presents two main advantages. Firstly, a novel LMS
approach is widely understood and applied to several prob- algorithm with adaptable step size emerges naturally with
lemssuchasechocancellation[4],noisecancellation[5],and thisapproach,makingitsuitableforbothstationaryandnon-
channelequalization[6]. stationary environments. The proposed algorithm has less
Althoughthesetwoapproachessharesomeunderlyingre- freeparametersthanpreviousLMSalgorithmswithvariable
lations, thereare veryfew connectionsin the literature. The stepsize[13,14,15],anditsparametersareeasiertobetuned
firstimportantattemptinthesignalprocessingcommunityto w.r.t. thesealgorithmsandstandardLMS.Secondly,the use
relate these two fields was the connection between a linear of a probabilistic model providesus with an estimate of the
errorvariance,whichisusefulinmanyapplications.
TheworkofJesu´sFerna´ndez-BeswassupportedbytheSpanishMECD
Experiments with simulated and real data show the ad-
FPU program and MICINN project PRI-PIBIN-2011-1266. The work of
V´ıctorElvirawassupportedbytheSpanishGovernment(MINECO)under vantages of the presented approach with respect to previous
project OTOSiS (TEC2013-41718-R), the Community of Madrid Govern- works.However,weremarkthatthemaincontributionofthis
mentunderproject AGES(S2010/BMD-2422), andtheEuropean Union’s paper is that it opens the door to introduce more Bayesian
7thFPthroughtheMarieCurieITNMLPM2012(GrantNo.316861).
machine learning techniques, such as variational inference
TheworkofStevenVanVaerenberghwassupportedbytheSpanishGov-
ernment(MINECO)underprojectRACHEL(TEC2013-47141-C4-3-R). andMonteCarlosamplingmethods[16],toadaptivefiltering.
scalar summarizingthe varianceof the estimate couldprove
to be sufficiently useful. In the next section, we show how
suchascalarisobtainednaturallywhenp(w |y )isapprox-
2. PROBABILISTICMODEL k 1:k
imatedwithanisotropicGaussiandistribution.Wealsoshow
Throughoutthis work, we assume the observation model to thatthisapproximationleadstoanLMS-likeestimation.
belinear-Gaussianwiththefollowingdistribution,
4. APPROXIMATINGTHEPOSTERIOR
p(yk|wk)=N(yk;xTkwk,σn2), (1) DISTRIBUTION:LMSFILTER
where σ2 is the variance of the observationnoise, x is the
n k Theproposedapproachconsistsin approximatingtheposte-
regression vector and w is the parameter vector to be se-
k rior distribution p(w |y ), in general a multivariate Gaus-
quentiallyestimated,bothM-dimensionalcolumnvectors. k 1:k
siandistributionwithafullcovariancematrix,byanisotropic
In a non-stationaryscenario, w followsa dynamicpro-
k sphericalGaussiandistribution
cess. Inparticular,weconsideradiffusionprocess(random-
walkmodel)withvarianceσ2forthisparametervector:
d pˆ(w |y )=N(w ;µˆ ,σˆ2I). (4)
k 1:k k k k
p(w |w )=N(w ;w ,σ2I), (2) In order to estimate the mean and covariance of the ap-
k k−1 k k−1 d proximatedistributionpˆ(w |y ),weproposetoselectthose
where I denotes the identity matrix. In order to initiate the k 1:k
that minimize the Kullback-Leibler divergence with respect
recursion,weassumethefollowingpriordistributiononw
k totheoriginaldistribution,i.e.,
p(w )=N(w ;0,σ2I).
0 0 d
{µˆ ,σˆ }=arg min{D (p(w |y ))kpˆ(w |y ))}.
k k KL k 1:k k 1:k
µˆk,σˆk
3. EXACTINFERENCEINTHISMODEL:
The derivation of the corresponding minimization prob-
REVISITINGTHERLSFILTER
lem can be found in Appendix A. In particular, the optimal
meanandthecovariancearefoundas
GiventhedescribedprobabilisticSSM,wewouldliketoinfer
the posterior probability distribution p(w |y ). Since all Tr{Σ }
k 1:k µˆ =µ ; σˆ2 = k . (5)
involved distributions are Gaussian, one can perform exact k k k M
inference,leveragingtheprobabilityrulesinastraightforward
We now show that by using (4) in the recursive pre-
manner.Theresultingprobabilitydistributionis
dictive and filtering expressions we obtain an LMS-like
adaptive rule. First, let us assume that we have an approx-
p(w |y )=N(w ;µ ,Σ ),
k 1:k k k k imate posterior distribution at k − 1, pˆ(w |y ) =
k−1 1:k−1
inwhichthemeanvectorµk isgivenby N(wk−1;µˆk−1,σˆk2−1I). Sinceallinvolveddistributionsare
Gaussian,thepredictivedistributionisobtainedas
µ =µ +K (y −xTµ )x ,
k k−1 k k k k−1 k
pˆ(w |y ) = p(w |w )pˆ(w |y )dw
k 1:k−1 k k−1 k−1 1:k−1 k−1
Z
wherewehaveintroducedtheauxiliaryvariable
= N(w ;µ ,Σ ), (6)
k k|k−1 k|k−1
Σ +σ2I
K = k−1 d , wherethemeanvectorandcovariancematrixaregivenby
k xT (Σ(cid:0) +σ2I)x(cid:1) +σ2
k k−1 d k n
µˆ = µˆ
andthecovariancematrixΣ isobtainedas k|k−1 k−1
k Σˆ = (σˆ2 +σ2)I.
k|k−1 k−1 d
Σ = I−K x xT (Σ +σ2),
k k k k k−1 d From(6),theposteriordistributionattimekcanbecom-
(cid:0) (cid:1)
putedusingBayes’TheoremandstandardGaussianmanipu-
Note that the mode of p(w |y ), i.e. the maximum-a-
k 1:k lations(seeforinstance[17,Ch. 4]). Then,weapproximate
posterioriestimate (MAP), coincideswith the RLS adaptive
theposteriorp(w |y )withanisotropicGaussian,
rule k 1:k
pˆ(w |y )=N(w ;µˆ ,σˆ2I),
w(RLS) =w(RLS)+K (y −xTw(RLS))x . (3) k 1:k k k k
k k−1 k k k k−1 k
where
Thisruleissimilartotheoneintroducedin[8].
(σˆ2 +σ2)
ofthFeinuanllcye,rntaoitnetythoafttthheeecsotvimaraitaencwekmcaotnridxitΣiokneidsoanmtheaesoubre- µˆk = µˆk−1+ (σˆk2−1+kσ−d21)kxkdk2+σn2(yk−xTkµˆk−1)xk
serveddatay . Nevertheless,formanyapplicationsasingle = µˆ +η (y −xTµˆ )x . (7)
1:k k−1 k k k k−1 k
Notethat,insteadofagainmatrixK asinEq.(3),wenow In stationary environments, the proposed algorithm has
k
haveascalargainη thatoperatesasavariablestepsize. onlyoneparameter,σ2. Wesimulateboththescenariowhere
k n
Finally,toobtaintheposteriorvariance,whichisourmea- we have perfectly knowledge of the amount of noise (prob-
sureof uncertainty,we apply(5) and the trickTr{x xT} = LMS1)andthecasewherethevalueσ2 is100timessmaller
k k n
xTx =kx k2, thantheactualvalue(probLMS2). TheMean-SquareDevia-
k k k
tion(MSD=Ekw −w k2),averagedoutover50indepen-
0 k
Tr(Σ ) dentsimulations,ispresentedinFig. 1.
σˆ2 = k (8)
k M
1 −20
= Tr I−η x xT (σˆ2 +σ2) (9) lms
= (cid:18)M1−(cid:8)η(cid:0)kkMxkkk2(cid:19)k(σˆkk2(cid:1)−1+k−σ1d2). d (cid:9) (10) 2w − w ||)0k−−−332505 nvpprlsslrrmoosbbnsLLlmMMsSS12
IfMAPestimationisperformed,weobtainanadaptablestep- E || −40
sizeLMSestimation B) ( −45
d
D ( −50
S
M
w(LMS) =w(LMS)+η (y −xTw(LMS))x , (11) −55
k k−1 k k k k−1 k
−60
0 1000 2000 3000 4000 5000 6000
with iterations
(σˆ2 +σ2)
η = k−1 d .
k (σˆ2 +σ2)kx k2+σ2 Fig. 1. Performancein terms of MSD of probabilistic LMS
k−1 d k n
with bothoptimal(probLMS1)and suboptimal(probLMS2)
Atthispoint,severalinterestingremarkscanbemade:
comparedtoLMS,NLMS,VS-LMS,andRLS.
• The adaptive rule (11) has linear complexity since it
doesnotrequireustocomputethefullmatrixΣ . The performance of probabilistic LMS is close to RLS
k
(obviously at a much lower computationalcost) and largely
• Forastationarymodel,wehaveσ2 =0in(7)and(10).
d outperformspreviousvariablestep-sizeLMSalgorithmspro-
In this case, the algorithm remains valid and both the
posedintheliterature. Notethat,whenthemodelisstation-
stepsizeandtheerrorvariance,σˆ ,vanishovertimek.
k ary,i.e. σ2 = 0in(2),boththeuncertaintyσˆ2,andtheadap-
d k
• Finally,theproposedadaptablestep-sizeLMShasonly tivestepsizeηk,vanishovertime. Thisimpliesthattheerror
two parameters, σ2 and σ2, (andonlyone, σ2, in sta- tendstozerowhenk goestoinfinity. Fig. 1alsoshowsthat
d n n
tionaryscenarios)incontrasttoothervariablestep-size the proposed approach is not very sensitive to a bad choice
algorithms [13, 14, 15]. More interestingly, both σ2 ofitsonlyparameter,asdemonstratedbythegoodresultsof
d
andσn2 haveaclearunderlyingphysicalmeaning,and probLMS2,whichusesaσn2 thatis100timessmallerthanthe
theycanbeestimatedinmanycases. Wewillcomment optimalvalue.
moreaboutthisinthenextsection.
MISO channel (real)
0.3
5. EXPERIMENTS Estimated channel (real)
0.2
0.1
Wboethesvtaaltuioanteartyheanpdertfroacrmkianngceexopfertihmeepnrtso.pIonsetdheafilgrostriethxmperiin- e (V) 0
ment,weestimate afixedvectorwo ofdimensionM = 50. oltag−0.1
V−0.2
The entries of the vector are independently and uniformly
−0.3
chosen in the range [−1,1]. Then, the vector is normalized
−0.4
so that kwok = 1. Regressors x are zero-mean Gaussian
k −0.5
vectors with identity covariance matrix. The additive noise
0 200 400 600 800 1000 1200 1400 1600 1800 2000
varianceis such that the SNR is 20 dB. We compareoural- Samples
gorithmwithstandardRLSandthreeotherLMS-basedalgo-
rithms: LMS,NLMS[18],VSS-LMS[15].1 Theprobabilis- Fig.2. Realpartofonecoefficientofthemeasuredandesti-
matedchannelinexperimenttwo.Theshadedarearepresents
tic LMS algorithm in [10] is not simulated because it is not
twostandarddeviationsfromtheprediction(themeanofthe
suitableforstationaryenvironments.
posteriordistribution).
1Theusedparametersforeachalgorithmare: forRLSλ = 1,ǫ−1 =
0.01;forLMSµ=0.01;forNLMSµ=0.5;andforVSS-LMSµmax =
1,α=0.95,C=1e−4.
Method LMS NLMS LMS-2013 VSSNLMS probLMS RLS • A similar approximation technique could be applied
MSD(dB) -28.45 -21.07 -14.36 -26.90 -28.36 -25.97 to more complex dynamical models, i.e. switching
dynamical models [20]. The derivation of efficient
adaptive algorithmsthat explicitly take into accounta
Table1. Steady-stateMSDofthedifferentalgorithmsforthe
switchinthedynamicsoftheparametersofinterestis
trackingofarealMISOchannel.
a non-trivial and open problem, though the proposed
approachcouldbeuseful.
In a second experiment, we test the tracking capabilities
• Finally, like standardLMS, this algorithmcan be ker-
oftheproposedalgorithmwithrealdataofawirelessMISO
nelized for its application in estimation under non-
channelacquiredin a realisticindoorscenario. Moredetails
linearscenarios.
onthesetupcanbefoundin[19]. Fig. 2showstherealpart
ofoneofthechannels,andtheestimateoftheproposedalgo-
rithm. The shaded area represents the estimated uncertainty A. KLDIVERGENCEBETWEENAGENERAL
for each prediction, i.e. µˆ ±2σˆ . Since the experimental GAUSSIANDISTRIBUTIONANDANISOTROPIC
k k
setupdoesnotallowustoobtaintheoptimalvaluesforthepa- GAUSSIAN
rameters,wefixtheseparameterstotheirvaluesthatoptimize
thesteady-statemeansquaredeviation(MSD).Table1shows Wewanttoapproximatepx1(x)=N(x;µ1,Σ1)bypx2(x)=
this steady-state MSD of the estimate of the MISO channel N(x;µ ,σ2I). In order to do so, we have to compute the
2 2
with different methods. As can be seen, the best tracking parametersofpx2(x),µ2andσ22,thatminimizethefollowing
performanceis obtainedby standard LMS and the proposed Kullback-Leiblerdivergence,
method.
6. CONCLUSIONSANDOPENEDEXTENSIONS DKL(px1kpx2) = Z−∞∞px1(x)lnppxx21((xx))dx
1
We havepresenteda probabilisticinterpretationof theleast- = 2{−M +Tr(σ2−2I·Σ−11)
mean-square filter. The resulting algorithm is an adaptable
+(µ −µ )Tσ−2I(µ −µ )
step-sizeLMSthatperformswellbothinstationaryandtrack- 2 1 2 2 1
ing scenarios. Moreover, it has fewer free parameters than σ2M
+ln 2 }. (12)
previousapproachesandtheseparametershaveaclearphysi- detΣ
1
calmeaning. Finally,asstatedintheintroduction,oneofthe
Usingsymmetryarguments,weobtain
advantagesofhavingaprobabilisticmodelisthatitiseasily
extensible: µ∗2 =argmin{DKL(px1kpx2)}=µ1. (13)
µ2
• If,insteadofusinganisotropicGaussiandistributionin Then,(12)getssimplifiedinto
the approximation, we used a Gaussian with diagonal
covariancematrix,wewouldobtainasimilaralgorithm
1 Σ σ2M
wfoirtheadcihffceoremnptosnteepntsoizfewska.nAdlmtheoausguhrethsisofmuondceelrctaainnbtye, DKL(px1kpx2)= 2{−M +Tr(σ221)+lndet2Σ1}. (14)
moredescriptive,itneedsmoreparameterstobetuned, The variance σ2 is computed in order to minimize this
2
andtheparallelismwithLMSvanishes. Kullback-Leiblerdivergenceas
• Similarly,ifwesubstitutethetransitionmodelof(2)by
anOrnstein-Uhlenbeckprocess, σ2∗ = argminD (P kP )
2 σ22 KL x1 x2
p(wk|wk−1)=N(wk;λwk−1,σd2), = argmσi22n{σ2−2Tr{Σ1}+Mlnσ22}. (15)
asimilaralgorithmisobtainedbutwithaforgettingfac-
Derivingandmakingitequalzeroleadsto
torλmultiplyingw(LMS) in(11). Thisalgorithmmay
k−1
have improved performance under such a kind of au-
toregresivedynamicsofwk,though,again,theconnec- ∂ Tr{Σ1} +Mlnσ2 = M − Tr{Σ1} =0.
tionwithstandardLMSbecomesdimmer. ∂σ22 (cid:20) σ22 2(cid:21) σ22 (σ22)2 (cid:12)(cid:12)σ22=σ22∗
(cid:12)
• Asin[10],themeasurementmodel(1)canbechanged Finally,sincethedivergencehasasingleextrem(cid:12)uminR+,
toobtainsimilaradaptivealgorithmsforclassification, Tr{Σ }
ordinalregression,andDirichletregressionforcompo- σ2∗ = 1 . (16)
2 M
sitionaldata.
7. REFERENCES [13] R.H.KwongandE.W.Johnston, “Avariablestepsize
LMSalgorithm,” IEEETrans.SignalProcess.,vol.40,
[1] L. Rabiner, “A tutorial on hidden Markov models no.7,pp.1633–1642,1992.
andselectedapplicationsinspeechrecognition,” Proc.
[14] T. Aboulnasrand K. Mayyas, “A robustvariable step-
IEEE,vol.77,no.2,pp.257–286,1989.
size LMS-type algorithm: analysis and simulations,”
[2] M. S. Arulampalam, S. Maskell, N. Gordon, and IEEE Trans. Signal Process., vol. 45, no. 3, pp. 631–
T. Clapp, “A tutorial on particle filters for on- 639,1997.
line nonlinear/non-gaussianBayesian tracking,” IEEE
[15] H.-C. Shin, A. H. Sayed, and W.-J. Song, “Variable
Trans. Signal Process., vol. 50, no. 2, pp. 174–188,
step-sizeNLMSandaffineprojectionalgorithms,”IEEE
2002.
Signalprocessing Letters, vol. 11, no. 2, pp. 132–135,
2004.
[3] S.Ji,Y.Xue,andL.Carin,“Bayesiancompressivesens-
ing,” IEEE Trans. Signal Process., vol. 56, no. 6, pp.
[16] D. Barber, Bayesian reasoning and machine learning,
2346–2356,2008.
CambridgeUniversityPress,2012.
[4] A. Gilloire and M. Vetterli, “Adaptivefiltering in sub- [17] K. P. Murphy, Machine learning: a probabilistic per-
bandswithcriticalsampling:analysis,experiments,and spective, MITpress,2012.
applicationtoacousticechocancellation,” IEEETrans.
SignalProcess.,vol.40,no.8,pp.1862–1875,1992. [18] A.H.Sayed,Adaptivefilters, JohnWiley&Sons,2008.
[19] J. Gutie´rrez, O´. Gonza´lez, J. Pe´rez, D. Ram´ırez,
[5] P. A. Nelson and S. J. Elliott, Active controlof sound,
L. Vielva, J. Iba´n˜ez, and I. Santamar´ıa, “Frequency-
Academicpress,1991.
domain methodology for measuring MIMO channels
using a generic test bed,” IEEE Trans. on Instrum.
[6] D.Falconer,S.L.Ariyavisitakul,A.Benyamin-Seeyar,
Meas.,vol.60,no.3,pp.827–838,Mar.2011.
and B. Eidson, “Frequency domain equalization for
single-carrierbroadbandwirelesssystems,” IEEECom-
[20] D. Barber and A. T. Cemgil, “Graphical models for
mun.Mag.,vol.40,no.4,pp.58–66,2002.
time-series,” IEEE Signal Process. Mag., vol. 27, no.
6,pp.18–28,2010.
[7] A. H. Sayed and T. Kailath, “A state-space approach
toadaptiveRLSfiltering,” IEEESignalProcess.Mag.,
vol.11,no.3,pp.18–60,1994.
[8] S.Haykin,A.H.Sayed,J.R.Zeidler,P.Yee,andP.C.
Wei, “Adaptivetrackingoflineartime-variantsystems
byextendedRLSalgorithms,” IEEETrans.SignalPro-
cess.,vol.45,no.5,pp.1118–1128,1997.
[9] S. Van Vaerenbergh,M. La´zaro-Gredilla, and I. Santa-
mar´ıa, “Kernelrecursiveleast-squarestrackerfortime-
varying regression,” IEEE Trans. Neural Netw. and
Learn.Syst.,vol.23,no.8,pp.1313–1326,2012.
[10] I. M. Park, S. Seth, and S. Van Vaerenbergh, “Proba-
bilistic kernelleast mean squaresalgorithms,” in 2014
IEEE International Conference on Acoustics, Speech
andSignalProcessing(ICASSP),Florence(Italy),May
2014.
[11] J.Cid-Sueiro,A.Arte´s-Rodr´ıguez,andA.R.Figueiras-
Vidal, “Recurrentradialbasisfunctionnetworksforop-
timalsymbol-by-symbolequalization,”SignalProcess.,
vol.40,no.1,pp.53–63,1994.
[12] J. Harrison, BayesianForecasting&DynamicModels,
Springer,1999.
